-
(10)
I developed transient 1-D transport-reaction models to analyze the hydrological systems of the Costa Rica subduction zone. The aim was to quantify vertical and horizontal fluid flow as well as the respective time scales of the flow systems in the oceanic crust and the prism wedge. In addition, the AMO reaction front across the boundary between prism and underthrust sediments was simulated to investigate the dynamics in the sedimentary carbon cycle during subduction.
The most obvious geochemical characteristic of sediments on the incoming Cocos plate is the observed reversal of pore water constituents to their respective bottom seawater concentrations at ~400 m sediment depth, just above the first sill complex (Fig. F3). This has been interpreted as the influence of horizontal fluid flow of cold modern seawater in the oceanic crust (Kastner et al., 2000; Kastner and Rudnicki, 2004; Kimura, Silver, Blum, et al., 1997; Silver et al., 2000).
Model runs were performed to estimate vertical flux of solutes from the basement aquifer into the overlying oceanic sediments as well as the timescale of this event. Whereas seawater acts as a source or a sink for solutes at the surface as well as at the bottom of the sediment column, sulfate reduction of organic material is the dominating diagenetic reaction within sediments. Hence, SO42– is consumed below the sediment surface, whereas NH4+, PO43–, and alkalinity are produced (Fig. F3). Therefore, the biogeochemical reaction term in equation 1 consists of only one term, the Monod-type rate law for OM degradation utilizing sulfate (equation 5). At Site 1039, the sediments contain >2 wt% organic matter at the surface and concentrations decrease exponentially to an average of 0.2 wt% below 200 mbsf (Fig. F4). This profile shape was least-squares fitted to the analytical solution of steady-state OM degradation (equation 6) and prescribed in the model runs to reduce computational costs. Finally, upper and lower boundary conditions for sulfate are set to bottom-water concentrations to adequately mimic the two seawater sources.
Two scenarios can explain the observed solute profiles (Fig. F3): (1) upward-directed advection of pore water and (2) purely diffusional transport of solutes from the lateral basement flow into the overlying sediment column. The model results of these two possibilities are discussed below.
First, a sensitivity analysis on the sulfate initial condition was performed because the sulfate distribution prior to establishing the lateral flow system in the crust is unknown but critical to the modeling. The tested initial profiles covered the range from a sediment column fully at seawater concentration to completely sulfate depleted. These tests clearly prove (1) upward advection will only lead to the observed sulfate profile if initial concentrations were lower than today, whereas (2) purely diffusional transport requires the initial concentrations to have been close to seawater values in the entire sediment column. In addition, analysis also revealed that the curved shape of the observed sulfate profile indicates a transient situation. These observations cannot be explained by a steady-state situation. Both scenarios predict steady-state profiles with a linear increase with depth below a near-surface sulfate minimum.
The initial sulfate condition can vary between an almost completely sulfate-free sediment column and an exponentially decaying profile with 27.5 mM at the surface and 14 mM at the bottom. In addition, the rate constant for OM degradation, kG, was chosen such as to fit the observed near-surface gradient when pore water advection is applied. A kG of 1.6 x 10–8/yr sufficiently fits this surface gradient.
Starting from such an initial sulfate distribution and applying the above kG of 1.6 x 10–8/yr, an upward advection rate, u0, of 0.15 cm/yr (i.e., 1500 m/m.y.) for the past 200–240 k.y. will lead to a good fit of the observed sulfate data (solid line; Fig. F4). Thus, the total flux of sulfate into the sediment column from below is
(10)
which is equal to 2.79 µmol/(cm2 yr). This is equivalent to a total flux of 5586–6703 mol/m2 within a time span of 200–240 k.y.
This scenario requires a sulfate-rich sediment column as initial condition (i.e., a constant concentration of 27.5 mM). In addition, the model suggests a reactive surface fraction of OM in the upper 30 mbsf degrading with a rate constant, kG, of 2.0 x 10–3/yr, whereas nonreactive OM is buried below this depth. The combination of anoxic OM degradation and upward diffusing SO42– results in a good fit of the observed sulfate data within a predicted time of 1 m.y. (dotted line; Fig. F4). About 0.025 µmol/(cm2 yr) sulfate is diffusing from the basement aquifer into the sediment column, corresponding to a total flux of 249 mol/m2 within the modeled time span of 1 m.y.
The diffusion scenario suggests a drastic change in organic matter reactivity that somewhat correlates with the observed change in the sedimentation regime from 30 m/m.y. to 105 m/m.y. about 500 k.y. ago (Kimura, Silver, Blum, et al., 1997), and this model result also indicates that the current situation has lasted for the past 1 m.y. Upward advection from the basement into the overlying sediment column, on the other hand, requires sufficient permeability. The total of ~150 m of gabbroic sill found at both drilling sites, Sites 1039 and 1253, shows an increasing number of fractures with depth (Morris, Villinger, Klaus, et al., 2003). However, these fractures are small, and it is not clear if they can maintain the necessary pressure difference to drive the upward advection.
Unfortunately, the other pore water constituents produced during anoxic organic matter degradation, such as NH4+, PO43–, and alkalinity, do not clarify either scenario. The results of these variables are highly dependent on the chosen initial profile, and as a consequence I have omitted their presentation here.
Because the lateral flow of seawater in the oceanic basement replenishes the sedimentary pore water with solutes, the resulting vertical sulfate flux must be balanced by the flow of seawater in the crust. Thus, it is possible to provide a minimum estimate of the lateral seawater flow based on the sulfate flux calculations.
In the advection scenario, the sulfate flux is equivalent to a seawater flux of 5905 µmol/(cm2 yr) or a volume flux of 1 L/(m2 yr) at a minimum lateral flow rate of 0.15 cm/yr. This results in a total water volume flux of 203,117–243,740 L/m2 within 200–240 k.y.
The diffusive sulfate flux of the second scenario requires a lateral basement flow of only 52.5 µmol/(cm2 yr). This is equivalent to a volume flux of 0.009 L/(m2 yr) at a lateral flow rate of 0.001 cm/yr or a total volume flux of 9037 L/m2 in 1 m.y.
Silver et al. (2000) estimated that a lateral basement flow of cold seawater at discharge rates of 1–5 m/yr would sufficiently cool the oceanic plate to explain the observed low heat flow data of 12 mW/m2. At these high lateral flow rates, therefore, the basement aquifer can easily maintain the required constant sulfate concentration at the interface to the sediment column, also suggesting that most seawater circulating through the oceanic crust is subducted.
Fluid flow in the prism sediments is primarily driven by sediment compaction during the subduction process and chemical dewatering of minerals at elevated temperatures. The latter reactions occur at ~15–20 km sediment depth (Silver et al., 2000), and fluids are channelled through prism sediments parallel to the décollement by faults and conduits. Hence, these pore fluids show a characteristic geochemical composition: distinct peaks of higher alkanes (Kimura, Silver, Blum, et al., 1997) and lithium (Chan and Kastner, 2000) are observed as well as enrichments in Ca2+ and depletions in K+ (Fig. F5) (Kimura, Silver, Blum, et al., 1997). During Leg 205, I collected piercing tool gas samples at Site 1254 at higher spatial resolution (i.e., three samples per core) (Morris, Villinger, Klaus, et al., 2003). The measured propane data show two major peaks: (1) directly above the décollement at ~360 mbsf and (2) at ~220 m sediment depth (Fig. F7). In addition, the asymmetric shape of these peaks also indicates upward advection of pore water due to compaction of the underthrust sediments during subduction.
Both flow systems were quantified during the numerical model runs. Whereas the vertical advection rate can be parameterized directly, the horizontal fluid flow reduces to a source/sink term in equation 1. Both peaks were treated independently using a constant source rate over a distinct interval representing the flow conduit. Less intense fluid flow between these two peaks led to elevated propane levels without any distinct peaks (see data from Leg 170 in Fig. F5). During Leg 205 this interval was not cored, and it is therefore not considered in the model. Because propane is not produced in the ambient sediment, zero initial and boundary conditions were chosen. Figure F7 depicts the model results with the best fit to the measured propane data.
A constant propane source, Rprop, of 22 ppmv/yr reproduced the peak maxima best, whereas the width of the propane peak is indicative of the transient character of both flow systems. The model results suggest that the upper conduit has been active for ~2000 yr and the décollement conduit for ~4000 yr. A thickness of 10 cm was chosen for the conduits as suggested by the stratigraphic results from Leg 205 (Morris, Villinger, Klaus, et al., 2003). However, a sensitivity analysis revealed that fluid flow through, for example, a 10-m-thick conduit combined with a propane rate of 0.22 ppmv/yr leads to an identical profile. Thus, the peak height and width depends on the depth-integrated propane rate, that is,
My model result suggests fairly continuous lateral fluid flow in the prism sediments over a couple of thousand years. Previous studies have suggested that episodic events rather than a stationary process dominate fluid flow in subduction zones. For the Nankai accretionary complex, for example, Saffer and Bekins (1998) suggested transient dewatering on timescales of 80–160 k.y. resulting from 2.5- to 5-m.y. cycles of pore pressure buildup and drainage. From the scatter of the chloride and ammonium profiles of Site 1040, Hensen and Wallmann (2005) inferred these lateral flow pulses need to be shorter than 10 k.y. Even short fluid pulses caused by, for example, seismic activity (Brown et al., 2005) cannot be excluded either, but they likely do not significantly affect the pore water signals. These short pulses will be levelled out by diffusion and advection.
Based on the 1-D modeling approach, it is difficult to determine lateral flow rate in the conduits because the source/sink term also includes the rate of propane production at depth.
In addition, an upward advection rate of 0.4 cm/yr was necessary to properly reproduce the asymmetry of both peaks. This value is in good agreement with the pore water flow resulting from the compaction of the underthrust sediments (i.e., the 25% reduction in overall thickness) (Kimura, Silver, Blum, et al., 1997). This is equivalent to an advection rate of 0.55 cm/yr. Screaton and Saffer (2005) discussed that the hemipelagic sediments dewater vertically, whereas the pelagic carbonate-rich section dewaters parallel to the décollement. However, the 29% reduction of the sediment Units U1 and U2 (Fig. F2) would still contribute to an advection rate of 0.26 cm/yr, a value close to my result. Hensen and Wallmann (2005) modeled the gas hydrate distribution in the prism sediments and determined a much lower, but steady-state, advection rate (0.03 cm/yr) when simulating the steep sulfate gradient near the sediment/water interface.
On subduction of the oceanic sediments, these sulfate-rich pore waters come into contact with the completely sulfate-depleted, methane-rich sediments of the prism wedge. As a consequence, a reaction front of AMO develops across the décollement. Whereas close to the deformation front (Sites 1043/1255) the sulfate/methane boundary is still very close to the décollement, the reaction front has already moved ~20 m into the underthrust sediments at Sites 1040/1254 (1.5 km further arcward) (Kimura, Silver, Blum, et al., 1997; Morris, Villinger, Klaus, et al., 2003). Figure F6 depicts the downward progressing reaction front combining data from ODP sites in a graph. A horizontal black line indicates the respective depths of the décollement.
A combination of the following three factors explains the downward progression of the AMO reaction front: (1) porosity values are 10%–15% higher in the underthrust sediments than in the overlying prism sediments (Fig. F6), thus facilitating downward transport of solutes, (2) methane diffusivity is ~60% higher than that of sulfate, and (3) the methane flux from the gas hydrate reservoir and microbial production is potentially larger than the upward-directed replenishment of sulfate from the basement aquifer.
The third numerical transport-reaction model was developed to simulate the downward-progressing AMO reaction front. Because the subduction process cannot be modeled in a 1-D model, the sediment depth is treated relative to the lithologic boundary (represented by the horizontal black lines in Fig. F6). The upper boundary of methane was set to 60 mM, the theoretical equilibrium concentration with the hydrate phase at 300 mbsf, and the measured CH4 concentrations were normalized to this value. Initially, concentrations above the décollement were set to 60 mM, whereas they were set to zero in the underthrust sediments. In contrast, sulfate concentrations were set to zero in the prism wedge, whereas the observed profile at Site 1039 was prescribed below the décollement. Additionally, porosity distribution at Site 1040, including the step-like change across the décollement, was used in the model calculations.
An AMO rate constant kAMO of 1 x 103/(mM yr) is sufficient to reproduce the measured sulfate data ~1500 m arcward of the deformation front at Sites 1040/1254 (solid lines and open circles, Fig. F6). At a subduction rate of 88 mm/yr (Kimura, Silver, Blum, et al., 1997), this takes 17,000 yr. Thus, 7.4 µmol/(cm2 yr) of methane are oxidized by AMO, or a total amount of 1262 mol/m2. This is equivalent to ~15 kg of carbon per square meter in 17 k.y. These are minimum values because the methane gradient decreases with time as the AMO front continues to progress deeper. For Sites 1043/1255, closer to the deformation front, the respective rates are 23.9 µmol/(cm2 yr) of methane oxidized or 12.9 kg of carbon per square meter after 4.5 k.y. AMO (equation 7) also produces hydrogen sulfide and bicarbonate, as indicated by the alkalinity peak below the décollement (Fig. F6). The model profile fits the observed alkalinity data without any further need to adjust, suggesting an appropriate model formulation. In addition, as SO42– becomes depleted in the underthrust sediments, barite starts to dissolve, thereby releasing barium and sulfate into the pore water (Fig. F6). The numerical model determines a barite dissolution rate of 1.2 x 10–5 mM/yr, which is equivalent to a Ba2+ flux of 23.6 µmol/(cm2 yr). Hence, a total amount of 130 kg of sulfur per square meter was dissolved within 17 k.y.
The model parameters (i.e., the AMO rate constant and barite dissolution rate) were adjusted solely to match the observed profiles at Sites 1040/1254. This parameterization was then used to successfully reproduce pore water distributions after 4500 yr of subduction at the drilling sites closer to the deformation front (i.e., Sites 1043/1255) without any further model adjustments. However, geochemical fluxes during these first 4500 yr of subduction are slightly underestimated because I used the porosity distribution of Sites 1040/1254 in this model run.
