Modeling Benthic–pelagic Nutrient Exchange Processes and Porewater Distributions in a Seasonally–hypoxic Sediment: Evidence for Massive Phosphate Release by Beggiatoa?

Supplementary Material Parameters used in the biogeochemical model not listed in Table 3 are provided in Table A1. Boundary conditions are given in Table A2. Constitutive equations which define the transport terms, temperature corrections and conversion factors are given in Table A3. General details of these functions are presented in Dale et al. (2011). The parameters for the porosity function were determined from the measured porosity values for all geochemical cores analyzed. Porosity was used, amongst others, to calculate in situ diffusion coefficients, D, from the molecular diffusion coefficients in porewater, D W , for the bottom water temperature and mean salinity of 20 using data provided by Boudreau (1997) and Schulz (2006). D was therefore time dependent through the temperature.


Introduction
The Baltic Sea is a marginal brackish sea in northern Europe consisting of deep and shallow basins interconnected by sills and channels.The only exchange with the ocean occurs via the Kattegat where saline water from the North Sea flows into the Baltic Sea underneath a fresher outgoing water mass.This results in a stratified water column in the central Baltic Sea with hypoxic (O 2 < 63 µM) or even anoxic bottom waters existing in the deep water areas (Conley et al., 2002).Widespread increases in the surface area affected by seasonal and permanent hypoxia and eutrophication are of Published by Copernicus Publications on behalf of the European Geosciences Union.
prime concern for policy making and environmental management strategies by Baltic Sea states (HELCOM, 2009).
Sediments underlying seasonally hypoxic waters are important sites of nutrient regeneration to the water column (Conley and Johnstone, 1995;Cowan and Boynton, 1996;Rozan et al., 2002;Reed et al., 2011).These fluxes are generally enhanced by bioturbation and bioirrigation; processes that describe the mixing of sediment and porewater by animals (Aller and Aller, 1998).Hypoxia generally decreases biodiversity and macrofaunal abundance (Levin et al., 2009), thus altering the processing of organic material in time and space.Yet, the change in depth and intensity of sediment mixing by animals at the community level is extremely complex to predict due to the huge variability in ecology and tolerance levels of animals to low oxygen (see reviews by Levin et al., 2009 andMiddelburg andLevin, 2009).Consequently, a key challenge in understanding the impacts of hypoxia in shallow or semi-enclosed aquatic ecosystems, such as the Baltic Sea, is to predict the response of the benthic compartment to decreasing levels of oxygen.Nonetheless, knowledge on the sensitivity of sediments to periodic and episodic fluctuations in bottom water O 2 still remains very fragmentary.This is partly because there are very few time series field observations in such settings (e.g.Rozan et al., 2002).Better incorporation of time series benthic field data into models should greatly improve the predictive capacity of feedbacks operating between the sediment and water column in areas prone to hypoxia (e.g.Eilola et al., 2011).
Reduction of ferric oxyhydroxide phases and release of iron-bound P to the porewater during hypoxia has been proposed as a major source of phosphate to the overlying water column by some authors (e.g.Sundby et al., 1992;Mc-Manus et al., 1997).Other authors highlight the potential role of microorganisms (Ingall and Jahnke, 1997;Hupfer et al., 2004;Sannigrahi and Ingall, 2005).Under oxic conditions, microorganisms can perform "luxury" uptake of phosphate and store P as intracellular polyphosphates.Under hypoxic or anoxic conditions, the polyphosphate pool is hydrolyzed, releasing phosphate back to the porewater.The role of vacuolated nitrate-storing sulfide oxidizing bacteria, such as Thiomargarita and Beggiatoa, in benthic phosphorus cycling is even more striking (Schulz and Schulz, 2005;Brock and Schulz-Vogt, 2010;Goldhammer et al., 2010).These microorganisms perform dissimilatory nitrate reduction to ammonium (DNRA) but also appear to uptake and release phosphorus.Working with Thiomargarita, Schulz and Schulz (2005) showed that subsurface total dissolved phosphate (TPO 4 ) peaks of > 200 µM could be sustained by exceptional bacterial TPO 4 release rates of 500-1700 nmol cm 3 d −1 .Moreover, Goldhammer et al. (2010) provided direct evidence for uptake and internal storage of TPO 4 by Thiomargarita and Beggiatoa and also showed that these bacteria mediate authigenic apatite precipitation under anoxic conditions.However, the phylogenetically and metabolically related bacteria Thioploca that inhabits the seasonally anoxic sediments off the Chilean coast appears not to perform luxury P uptake and release, or at least not at the time that the sediments were sampled (Holmkvist et al., 2010).
In this paper we present a suite of biogeochemical data, including measured rates of bioirrigation, collected over a year at roughly monthly intervals at the Boknis Eck time series station located at the entrance of Eckernförde Bay (Germany) in the SW Baltic Sea.The highly transient features of the data are analyzed using a previous steady-state model that has been modified to account for the temporal changes in bottom water chemistry, organic matter fluxes and sediment mixing rates by animals (Dale et al., 2011).The principle motive for the modeling work is to interpret the major features of the data and quantify the fluxes of bioavailable chemical compounds to and from the sediment.We provide circumstantial evidence for TPO 4 release by Beggiatoa at rates that are in good agreement with the few literature data available.The abundance and distribution of these bacteria in the sediments has been documented previously by Preisler et al. (2007), and the cycling of nitrogen by Beggiatoa at Boknis Eck was investigated by Dale et al. (2011).
The present work exemplifies how the combination of time series data and modeling can be a powerful tool for quantifying benthic solute exchange fluxes and identifying the response times of biogeochemical processes to O 2 depletion.

Study site
Boknis Eck is a narrow channel located at the northern entrance of Eckernförde Bay (54 • 31 N, 10 • 20 E) and has a maximum water depth of about 28 m.The deep waters are frequently ventilated in winter due to saline inflows from the Great Belt (Hansen et al., 1999).From mid-March until mid-September, vertical mixing is restricted by density stratification of the water column, which leads to pronounced periods of hypoxia during late summer caused by microbial respiration of organic material in the bottom waters and sediment (Graf et al., 1983;Meyer-Reil, 1983;Hansen et al., 1999).Phytoplankton blooms dominated by diatoms and dinoflagellates generally occur in spring (March to April) and autumn (September to November).This causes an increase in benthic activity (Graf et al., 1983) and elevated fluxes of phosphate and ammonium from the sediment (Balzer et al., 1983).Autumn storms and a decrease in surface water temperature lead to an overturn of the water column and ventilation of the deeper water layers (Hansen et al., 1999).The oxygen penetration depth into the sediments is limited to the upper millimeters when bottom waters are oxic (Preisler et al., 2007), yet may be as high as 10 mm in the shallower sandier sediments at 18 m depth (Graf et al., 1983).
Table 1.Sampling dates where multiple cores (MUC) were retrieved for geochemical analysis and bioirrigation experiments.Temperature, salinity and dissolved oxygen data in the bottom water (28 m) are also listed, along with the temperature-correction factor for reaction rates, f (T ).

Sampling no. Date in 2010
Geochemistry Bioirrigation T ( a Formula for f (T ) is given in the Supplement.
b Two cores on 28 April were retrieved using a benthic lander.
The sediments at the present study site (28 m water depth) are classified as fine-grained muds (< 40 µm) with a carbon content of 3 to 5 wt % (Balzer et al., 1986).Calcium carbonate content is < 2 % (Balzer et al., 1983;Orsi et al., 1996).With no significant terrestrial runoff, the bulk of organic matter within Eckernförde Bay sediments originates from marine plankton and macroalgal sources (Balzer et al., 1987;Orsi et al., 1996).Degradation and fermentation of organic matter produces free methane gas at Boknis Eck but fluid seepage from groundwater discharge has not been observed (Schlüter et al., 2000).The dominant fauna in the sediments in winter/spring are the polychaetes Pectinaria koreni and Nephtys ciliate with recorded abundances of 201-476 and 63-122 individuals m −2 , respectively (Graf et al., 1982).Specimens up to 10-15 cm long were present during winter in the cores retrieved in this paper.Mats of Beggiatoa were absent on the surface of Boknis Eck sediments during the period with oxygenated bottom waters, but are present below the sediment surface at the redox interface at the top of the sulfide layer (Preisler et al., 2007).In late summer 2010 when the bottom waters became almost anoxic, we observed a blackening of the surface sediments and colonization by Beggiatoa filaments.

Sampling and analytical methods
Sediment samples were obtained on 12 occasions between February and December 2010, which covers the major events of the year, that is, spring and autumn blooms, stagnation and severe hypoxia in late summer and winter dormancy (Graf et al., 1983).The sampling dates as well as temperature, salinity and dissolved oxygen (O 2 ) concentrations in the bottom water (25-28 m) measured using a CTD (Hydro-Bios, Kiel, Germany) are shown in Table 1.The sediment sampling device was a mini-multiple corer (MUC) equipped with 4 core liners (60 cm long) which retrieved sediment cores with a diameter of 10 cm and a maximal length of ca.35 cm.In general, 2 cores were sampled on each occasion for geochemical analysis, although 2 additional cores were retrieved on 28 April using a benthic lander (Sommer et al., 2008).Upon return to the GEOMAR laboratory 2-3 h later, the cores were immediately transferred to a cool room (4 • C).Porewater from one core was sub-sampled anaerobically using rhizons and from the other(s) by sectioning into 1-2 cm thick slices under ambient air and squeezing the sediment with low pressure argon (1 to 5 bar) through 0.2 µm cellulose acetate nucleopore ® filters.Prior to use, the rhizons were preconditioned in an oxygen-free water bath.To ensure that the samples had not been in contact with air, the first 0.5 cm 3 of extracted porewater was discarded.All filtrates were collected in acid-cleaned recipient vessels and further manipulated for analytical determinations.Wet sediment samples (∼ 5 mL) were freeze-dried for the determination of porosity from the water content using a dry sediment density of 2.5 g cm −3 .
The squeezed porewater samples were analyzed for the following solutes: sulfate (SO 2− 4 ), nitrate (NO − 3 ), nitrite (NO − 2 ), ammonium (NH + 4 ), total hydrogen sulfide (TH 2 S ≈ H 2 S + HS − ), ferrous iron (Fe 2+ ), total phosphate (TPO 4 ≈ H 2 PO − 4 + HPO 2− 4 ), bromide (Br − ), dissolved inorganic carbon (TCO 2 = HCO − 3 + CO 2− 3 + CO 2 ) and total alkalinity (TA).The porewater sampled using rhizons was simultaneously analyzed for NH + 4 , Fe 2+ , TPO 4 and TA.Analytical details for SO 2− 4 , NO − 3 , NO − 2 , NH + 4 , TH 2 S, TCO 2 and TA are described in Dale et al. (2011).Br − was determined in the squeezed porewater by ion chromatography using a Metrohm ion chromatograph with a conventional anion exchange column with IAPSO seawater standard for calibration (relative precision of < 2 % for natural seawater samples).For the determination of Fe 2+ , 10 µL of concentrated HNO 3 was added to 1 mL of porewater, which was then analyzed by ICP-AES with a precision of < 2 %.All dissolved iron in the porewater was assumed to be present as Fe 2+ .TPO 4 was determined using photometry with a precision of 5 %.Systematic differences in the measured concentrations of the samples obtained by rhizons and squeezed sediments could not be discerned.Therefore, although the data are presented according to the extraction method, they are not discussed individually.
In this paper, we also present rates of benthic sulfate reduction (SR) measured in parallel cores from April to December presented by Bertics et al. (2012).SR data from 5 and 6 March 2002 measured at the same location and published by Treude et al. (2005) are used to help interpret the data from 22 March 2010 in this study.Analytical details are detailed in those manuscripts.

Bioirrigation experiments
Throughout the duration of the sampling campaign, additional multi-cores (10 cm internal diameter) were taken to the home laboratory for manipulation experiments to determine bioirrigation rates.The procedure involves adding a large excess of a dissolved conservative tracer, in this case Br − , to the water overlying the sediment and allowing the core to incubate for a known period of time.Posterior analytical determination of the depth distribution of Br − allows the rate of irrigation to be approximated by modeling the transient infiltration of Br − into the sediment using a numerical model (Martin and Banta, 1992).The experimental and simulation procedure is described fully in the Supplement and a brief description is provided below.
The measured Br − concentration profiles at the end of the incubation were simulated using a non-steady state model which accounts for solute transport into the sediments by molecular diffusion and bioirrigation using a nonlocal source-sink function (Boudreau, 1984;Emerson et al., 1984).Assuming that porosity does not change during the incubation period, the 1-D mass conservation equation which describes the change in Br − concentration (M) in the sediment with time is expressed as where z (cm) is depth in the sediment, t (d) is time, α bi (d −1 ) is the depth-dependent bioirrigation coefficient describing solute pumping through animal burrows, and Br olw (M) is the time-dependent Br − concentration in the well-mixed overlying water.The sediment porosity, ϕ, was defined using a depth-dependent function (see Supplement).
The depth-dependence of α bi was described using where α bi1 (d −1 ) is approximately equal to the rate of bioirrigation at the sediment surface and α bi2 (cm) controls the irrigation depth.This type of formulation allows for a constant rate of irrigation in the sediment down to a depth equal to ca. α bi2 (Dale et al., 2011) and better describes the data from Boknis Eck compared to the more commonly employed single exponential function (e.g.Martin and Banta, 1992).

Transport processes and biogeochemical reactions
A steady state 1-D numerical reaction-transport model described previously (Dale et al., 2011)  , TPO 4 , Br − and methane (CH 4 ).TA was not modeled explicitly but calculated from the TH 2 S and TCO 2 profiles and a mean measured porewater pH of 7.7 (not shown).Solid components consist of particulate iron oxyhydroxides (Fe(OH) 3 ) with an associated fraction of iron-bound phosphorus (Fe-P), particulate sulfide as pyrite (FeS 2 ) and particulate organic matter (POM).Previously, 3 POM pools, "Gi", of differing reactivity, i (i = 1, 2, 3, with reactivity 1 > 2 > 3) were considered (Dale et al., 2011).Here, a fourth POM pool (G0) is included to account for fresh phytodetritus deposited during the spring and autumn blooms.
Geochemical species are transported through the sediment by advective and diffusive processes.Bioirrigation is also considered for solutes and transport of NO − 3 bac into the sediment is described analogously to bioirrigation (Dale et al., 2011).In this study, we also consider an additional nonlocal solute exchange due to rising gas bubbles upon which we elaborate further below.Using the traditional mathematical approach (Berner, 1980;Boudreau, 1997), the concentration of aqueous species, C a (mmol cm −3 ), solid species, C s (dry wt %), and intracellular nitrate, C b (mmol cm −3 ), in 1-D along the vertical z-axis with time, t, can be described as where ϕ is porosity, v a (cm d −1 ) and v s (cm d −1 ) are the burial velocities of porewater and solids, respectively, D (cm 2 d −1 ) is the tortuosity-corrected molecular diffusion coefficient, D b (cm 2 d −1 ) is the bioturbation coefficient, α bi (d −1 ) is the solute bioirrigation coefficient determined from the Br − incubations, C a (0) is the solute concentration at the sediment-water interface, α b (d −1 ) is the coefficient for non-local transport of nitrate by filamentous sulfide oxidizing bacteria, C b (0) is the intracellular nitrate concentration at the sediment-water interface, α bu (d −1 ) is the non-local exchange coefficient due to irrigation by rising methane gas bubbles, and γ bi is a dimensionless factor that scales the irrigation coefficients of the different solutes to that of Br − .R is the sum of the rate of change of concentration due to biogeochemical reactions and R DNRA is the rate of DNRA.A detailed description of the model is provided in the Supplement.
The biogeochemical reaction network is shown in Table 2 and parameters are given in Table 3.Most of these reactions have been discussed by Dale et al. (2011).Those which are newly described here relate to iron and phosphorus cycling (Fig. 1).POM is a composite of organic C, N and P fractions defined as CN (N : C)i P (P : C)i , where (N : C) i and (P : C) i are the ratios of organic N and organic P to particulate organic carbon (POC) in pool Gi. (N : C) i is specific to the Gi pool (see Supplement), whereas a constant Redfield (P : C) i value of 1/106 determined from in situ flux experiments was prescribed for all POM fractions (Balzer et al., 1983).POM is degraded by aerobic respiration, denitrification, dissimilatory iron reduction (DIR), sulfate reduction (SR) and methanogenesis.The rate of POM mineralization at each depth is determined by the concentration of available electron acceptors using Michaelis-Menten kinetic limitation terms (Table 2).A nominal oxidation state of zero for POC is assumed.
Concentration profiles of Fe 2+ and Fe(OH) 3 indicate that reduction of Fe(OH) 3 is an important component of the iron cycle in Boknis Eck sediments (Preisler et al., 2007).Therefore, and in line with previous studies (Berg et al., 2003;Dale et al., 2009a), Fe(OH) 3 was partitioned into a fraction termed Fe(OH) 3-A that is reduced during DIR and by reaction with dissolved sulfide and a fraction that reacts with sulfide only (Fe(OH) 3-B ).The iron cycle also includes precipitation of iron sulfides (assumed to be FeS 2 ) and the oxidation of FeS 2 and Fe 2+ by O 2 .Bimolecular rate expressions that are first-order in reactants are prescribed for these reactions, with the exception of the sulfidic reduction of iron oxyhydroxides which is half-order in sulfide (Poulton, 2003).Reported rate constant values for these reactions span many ordersof-magnitude and the values used here are based on studies where iron cycling is well-documented (compiled by Dale et al., 2012).
Fe-P is co-deposited to the sediment with iron oxyhydroxides (Fig. 1).The P : Fe molar ratio in this particulate material (ε) is taken to be 0.1.Slomp et al. (1996) argued that this was a typical ratio for natural and synthetic poorly crys-   1.Schematic of the coupled Fe and P cycles considered in the model.Solid species are in squares and solutes in circles.P is added to the sediment as organic P and with iron oxyhydroxides.Phosphate is released to the porewater upon reductive dissolution of ferric iron and oxidation of organic matter.This dissolved P pool is free to be transported out of the sediment, or removed from the porewater by co-precipitation with ferrous iron oxidation.P uptake and release by Beggiatoa is explored in Sect.4.4.talline iron oxides (see references therein) and obtained good modeling results using this value.Dissolution of iron oxyhydroxides releases both Fe 2+ and TPO 4 to the porewater.Authigenic iron oxyhydroxide precipitation (R Fe2ox , Table 2) removes a fraction of dissolved TPO 4 and increases the Fe-P pool.However, this reaction mainly takes place in the thin oxidized sediment layer where the ratio of TPO 4 to Fe 2+ is often lower than ε.Hence, uptake of TPO 4 into Fe-P cannot simply be equal to the rate of iron oxidation multiplied by ε because this would result in negative TPO 4 concentrations at the surface.To circumvent this problem, co-precipitation of TPO 4 with iron oxyhydroxides is limited by the ambient TPO 4 concentration (Table 2).The mass conservative rate of TPO 4 release during dissolution of iron oxyhydroxides is thus based on the relevant depth-dependent fraction of Fe-P to Fe(OH) 3 (θ in Table 2).Note that ε and (P : C) i are not treated as fitting parameters initially, in order to highlight the potential shortcomings of the P cycle included in the model and to identify a possible source of TPO 4 from Beggiatoa.
The presence of uniformly dispersed authigenic carbonate fluorapatite (CFA) is probably widespread in continental margin sediments (Ruttenburg and Berner, 1993;Slomp et al., 1998).Apatite precipitation is favored thermodynamically by high concentrations of TPO 4 in the porewater.Given that TPO 4 concentrations in the surface sediments are low for much of the year at Boknis Eck (see Results), apatite precipitation is also likely to be low.Reed et al. (2011) concluded that apatite precipitation in the seasonally hypoxic Arkona Basin (S Baltic Sea) was unimportant, and that authigenic P www.biogeosciences.net/10/629/2013/Biogeosciences, 10, 629-651, 2013

POM mineralization reactions
a Rates of POM degradation are dependent on the specific organic matter pool, Gi, under consideration.POM is defined as CN (N : C) i P (P : C) i .b For clarity, the reactions are not proton balanced and H 2 O is omitted.c TPO 4 release during DIR (R Fe ) and reductive iron dissolution (R Fe3red ) is calculated using the fraction of iron-associated TPO 4 (Fe-P) bound to Fe(OH) 3 , θ , where K j is the half-saturation constant for species j .e Factors to convert between particulate and dissolved C, Fe and P (f POC , f Fe and f P ) are given in the Supplement.f Rate units for the other species can be determined by the reaction stoichiometry.A separate equation for TPO 4 is provided for R FeS2ox (see text).In situ fluxes and sediment profiles (Balzer et al., 1983;Balzer, 1984) a Bottom water SO 2− 4 and Br − are scaled to salinity and shown in Fig. 2.
is mainly deposited from the water column and buried conservatively.

Time-dependent functions and boundary conditions
This section details the temporally variable boundary conditions and functions used in the model.Further details as well as information on the time-invariable boundary conditions are given in the Supplement.
A major advance of this model compared to the previous steady state application (Dale et al., 2011) is that the entire time series of porewater data and SR rates was used to constrain the biogeochemical reactions.This provides a unique opportunity to more accurately parameterize several key transport and reaction terms compared to transient model applications constrained with data from a single sampling only.The temporal variability in the data is driven by seasonal changes at the seafloor.To capture this dynamic, the model boundary conditions and forcing functions were varied intra-annually but not inter-annually until a new seasonally cycling steady state model simulation was obtained.
To start with, measured bottom water temperatures were imposed in the model (Fig. 2a).Heat transfer is rapid through muddy sediment and temperature changes at the sediment surface were assumed to immediately penetrate the modeled sediment column (40 cm), although this is only an approximation (Dale et al., 2008).Molecular diffusion coefficients and chemical reaction rates were dynamically updated for the temperature changes (Boudreau, 1997;Schulz, 2006), the latter by multiplying the rate by an Arrhenius term, f (T ) (Table 2).The temperature response was broadly defined using a Q 10 value of 3 for POM mineralization and 2 for secondary redox reactions.However, differential temperature responses for some reactions were included if supported by experimental evidence, for example denitrification and anammox (see Supplement).The change in bottom water salinity (21 ± 2.5; Table 1) was much less than for temperature (ca. 5 ± 4 • C) and the effect of salinity on diffusion coefficients was ignored.Salinity changes were nonetheless considered for major seawater ions (Br − , SO 2− 4 ) whose upper boundary concentrations scaled to standard (IAPSO) seawater composition within a margin of ca. 10 % (Fig. 2b).As an upper boundary condition for O 2 , the measured bottom water concentrations based on the data in Table 1 were used, supplemented with O 2 data gathered from the time series monitoring program at the same location (Fig. 2c; Bange et al., 2011).Although bottom water inorganic N concentrations also fluctuate throughout the year (Dale et al., 2011), the trends are not obvious and ignored in the model.Bottom water TPO 4 concentrations show a clear trend, but are also not imposed as time-variable boundary conditions.This will introduce a small error to the calculated fluxes, yet one that is certainly much smaller than the uncertainty in other model parameters.Thus, all bottom water solute concentrations except O 2 , Br − and SO 2− 4 were fixed.Fixed concentrations were defined at the lower boundary for SO 2− 4 , TCO 2 , NH + 4 , TH 2 S, TPO 4 and CH 4 because their distribution in the modeled domain is influenced by POM mineralization at greater depth (Dale et al., 2011).All other solutes at the lower boundary were defined with a zero  concentration gradient (Neumann condition).Lower boundary conditions for all variables were time-invariable.
The upper boundary concentration for intracellular nitrate, NO 3 − bac , was set to the porewater normalized concentration (150 µM) measured in sediment samples from Boknis Eck that were repeatedly frozen and thawed to release nitrate stored within Beggiatoa (Preisler et al., 2007).Previously, we set this boundary condition to bottom water nitrate concentrations (Dale et al., 2011).We have revised this approach for consistency with modeling studies at other sites where porewater normalized NO 3 − bac concentrations were used (Dale et al., 2009a;Bohlen et al., 2011).
In the previous model (Dale et al., 2011), DNRA showed a maximum rate in the upper 5 cm as shown by Preisler et al. (2007), with an attenuated rate down to the bioirrigated depth (10 cm).When the complete set of data was analyzed, however, we were able to refine this further by restricting the depth of intracellular nitrate transport to 5 cm.This better agrees with measured Beggiatoa biomass concentrations at different times of the year (Preisler et al., 2007).We observed Beggiatoa filaments at the sediment surface during hypoxia and very low sulfide concentrations in the upper 5 cm in Oc-tober and November.Based on previous results in sediments with Thioploca (Holmkvist et al., 2010), it is likely that the higher bacterial biomass at this time maintained these low sulfide concentrations.The rate of nitrate transport was therefore increased in October and November to simulate the sulfide profiles.The sulfide depletion lags behind the minimum in O 2 by around 1 month (see Sect. "Results"), which may reflect the time delay for the Beggiatoa biomass to respond to increasing hypoxia at the sediment surface.Thus, nitrate transport (Fig. 2d) was dependent on the time-lagged O 2 concentrations (O 2 ) shown in Fig. 2c.
The rate constant for reductive dissolution of iron oxide particles, k Fe3red , was made to be dependent on O 2 (Table 2).An explanation is given in the Discussion.
POM fluxes to the sea floor in the deep channel at Boknis Eck have been reported to be 0.4 to 0.59 mmol C cm −2 yr −1 based on in situ total oxygen uptake measurements (Balzer et al., 1986).This flux is composed as a continuum of reactive fractions, parameterized in this study using a multi-G approach consisting of 4 POM pools.The total flux of the less reactive G2 and G3 pools has been quantified as 0.12 mmol C cm −2 yr −1 (Dale et al., 2011).Due to the recalcitrant nature of these fractions, seasonal changes in their flux will have a negligible effect on solute dynamics in the upper 40 cm.G2 and G3 fluxes are therefore time-invariable, which leaves up to 0.47 mmol C cm −2 yr −1 to be split between the reactive G0 and G1 fractions.G0 and G1 fluxes were time-dependent.The temporal variability in G1 flux was calculated by first integrating measured sulfate reduction rates over the upper 15 cm (Bertics et al., 2012).Multiplying these rates by 2 provides the carbon mineralization rate, assuming a carbon oxidation state of zero.Dividing these rates by f (T ) (Table 1) corrects for changes in temperature and provides an estimate of the G1 flux.Yet, these fluxes are minimum values since G1 is also oxidized by other electron acceptors such as O 2 and NO − 3 .Therefore, the flux of G1 was constrained with the model by increasing the G1 flux proportionally over the year until a good correspondence with the measured SR rates was achieved.When integrated over the year, this amounted to 0.42 mmol C cm −2 yr −1 .This leaves 0.05 mmol C cm −2 yr −1 for G0, which was split evenly over the spring (March-April) and autumn (∼ September) blooms.Fluxes of G0 and G1 (in mmol C m −2 d −1 ) used in the model are shown in Fig. 2e.
As a first approximation, the fluxes of particulate iron pools to the seafloor were time-invariable.FeS 2 was set to zero.The total flux of iron oxyhydroxides was estimated from mass accumulation rates and iron concentrations at a nearby site (Balzer, 1982;Lapp and Balzer, 1993).This flux was divided equally into Fe(OH) 3-A and Fe(OH) 3-B .
Bioturbation is limited to the upper few cm in Boknis Eck (Dale et al., 2011).The seasonal variability is unknown, but the rate of mixing is likely diminished in autumn when hypoxia is established and macrofauna are absent.Thus, the biodiffusion constant, D b (0), was scaled to O 2 (Fig. 2f) so that in February D b (0) equaled the value determined for winter, ca.0.08 cm 2 d −1 (Dale et al., 2011).
Bioirrigation coefficients α bi1 and α bi2 derived from the tracer incubations were imposed directly in the model and are described in the Results.The parameter γ bi that defines solute specificity of the irrigation parameters is detailed in the Supplement.Note that because of rapid oxidation of Fe 2+ and precipitation of authigenic iron oxyhydroxides on burrow walls, Fe 2+ is assumed not to be transported by bioirrigation (Meile et al., 2005).
The bubble irrigation coefficient, α bu , was prescribed a constant value over a defined depth: where z bu (30 cm) is the maximum depth where the porewater is strongly affected by bubble irrigation.Irrigation of sediment porewater by gas rising through open bubble tubes or vacant animal burrows was proposed by Martens (1976), and the theory was developed further by Haeckel et al. (2007).These authors described how bubbles enhance diffusive mixing in their wake as they rise through tube structures and act to abiotically irrigate the porewater by homogenizing concentration gradients between the sediment and the overlying seawater.The idea of using constant rates of irrigation down to a specific depth is based on simulations of sediments affected by CH 4 gas venting in the Sea of Okhotsk (Haeckel et al., 2007).Although a mechanistic theory which explains this phenomenon has not been reported, the gas bubbles may manage to form their own small fractures and tubes in the muddy, unconsolidated near-surface sediments (Algar et al., 2011).Once these high permeability conduits are formed, rising bubbles reach a terminal rise velocity and are released rapidly from the sediment to the water column (Haeckel et al., 2007).Bubble irrigation is defined analogously to bioirrigation using the same solute specific factors, γ bi (Eq.3a).
The depth, intensity and duration of bubble irrigation events in Boknis Eck were estimated from the entire suite of data, but mainly from SO 2− 4 , NH + 4 , TH 2 S, TCO 2 and TA.The temporal change in α bu is shown in Fig. 2h.
Following the gas irrigation events, a rapid steepening of the SO 2− 4 , H 2 S, TCO 2 and TA gradients, indicative of enhanced sulfate reduction, was observed in the following months (see Results).We argue that the newly opened bubble conduits provide an avenue for a residual reservoir of smaller methane bubbles to travel slowly upward, dissolve and consume SO 2− 4 by anaerobic oxidation of methane (AOM, Table 2) producing H 2 S and TCO 2 .We simulated the net effect of this process by imposing an upward advection term for CH 4 only (Fig. 2h).The rate of advection was constrained from the SO 2− 4 , H 2 S, TCO 2 and TA concentrations.A time lag of around 1 month for this process following bubble irrigation provided the best correspondence with the data.

Model uncertainties and solution
Although the data in Fig. 2 are resolved at approximately monthly intervals, there are several caveats to consider assuming that these forcing functions are constant from yearto-year.Firstly, random fluctuations in the boundary conditions are very likely occurring on shorter timescales.Secondly, the progressive decline in mean O 2 levels and increase in mean temperature since the 1950s is not accounted for (Bange et al., 2011).The effect of these changes is considered to be small compared to those arising from the intraannual variations and to have no significant bearing on the results presented here.Thirdly, Arntz (1981) argued that the benthic community assemblage at any given time is determined by the length and severity of hypoxia over the previous years.This influences the rate and depth of bioturbation, bioirrigation and the redox potential of the surface layers, potentially leading to a variable oxygen debt from year-to-year via the accumulation of reduced geochemical species (Graf et al., 1983).
In view of these uncertainties, it is not our intention to simulate every detail of the measured geochemical data.This would require an excess of tunable parameters, lead to a www.biogeosciences.net/10/629/2013/Biogeosciences, 10, 629-651, 2013 The model presented consists of a complex set of reactions and many adjustable parameterizations.In actuality, most of the parameters are constrainable from field observations at Boknis Eck or from literature studies that demonstrate either experimentally or theoretically the assignment of a parameter value.The sensitivity of the model nitrogen cycle has been dealt with previously (Dale et al., 2011), and in this work we detail the areas where more research effort is needed to more accurately quantify iron and phosphorus cycling.
To solve the model, the continuous spatial derivatives in Eq. ( 3) were replaced with finite differences (Boudreau, 1997).The resulting set of ODEs (Ordinary Differential Equation) was solved using the NDSolve algorithm in MATHEMATICA 7.0 using the method of lines (Boudreau, 1996).A high resolution of 0.06 cm was used at the surface to minimize numerical errors and determine reliable gradients within these highly reactive layers.At the base of the sediment the grid thickness increased to 0.7 cm.A centered finite difference scheme was used for dissolved species and solid species within the bioturbated zone.An upward scheme was applied for the transport of solids below this depth to avoid numerical instabilities (Boudreau, 1996).Adopting this approach, the model simulations were stable despite the often abrupt change in boundary conditions and forcings.The model was run in a seasonally cycling transient model that required ca.600 yr to reach a dynamic steady state.Mass conservation was checked by annually integrating daily fluxes, resulting in a mass balance generally > 99 % for all species.

Results
Measured and modeled Br − concentrations in the incubated cores indicate a general trend of high bioirrigation rates in winter with well-oxygenated bottom waters and low or zero bioirrigation during the hypoxic summer and autumn period (Fig. 3).In winter (23 February), Br − concentrations in the upper 10 cm at the end of the incubation varied little from the Br − concentration measured at the end of the incubation in the overlying water (16 mM), coincident with the greatest depth where polychaetes were observed, and then decreased steadily down to 20 cm.The corresponding modeled Br − tracer profile using Eq. ( 1) provides a very close correspondence with the measured data with values of α bi1 and α bi2 of 0.28 d −1 and 12.8 cm, respectively, indicating near constant irrigation rates to almost 13 cm.This is very different from the theoretical profile considering transport of Br − into the sediment by molecular diffusion only where Br − diffuses down to ca. 8 cm.Deep irrigation did not reoccur over the sampling period.The same trends are mirrored in the concentrations of all solutes (with the exception of Fe 2+ ), which display near-constant concentrations in the upper 10 cm in February at values equal to those measured in the bottom water (Fig. 4).
Between February and April, α bi1 remained high yet then fell dramatically to 0.03 d −1 in July.The porewater concentration gradients in the upper 10 cm also became more pronounced over this period (again excluding Fe 2+ ).The measured Br − tracer concentration in July showed a similar distribution to the diffusion-only profile, coincident with a drop in bottom water O 2 to 70 µM (Table 1).This value is very close to the reported threshold for hypoxia (∼ 63 µM O 2 ), at which large changes are engendered in community structure by oxygen stress (Diaz and Rosenburg, 2008).Although O 2 continued to fall in August (47 µM), α bi1 and α bi2 inexplicably increased again to 0.08 d −1 and 9.6 cm, respectively.Whilst it is possible that this could highlight increased pumping by the benthic community to compensate for low oxygen levels (Forster al., 2003), we are hesitant about these values since they could reflect artifacts caused by removal of the core under hypoxic bottom waters and incubation under atmospheric conditions.
From September, α bi1 fell to low levels of < 0.01 d −1 , essentially indicating a complete absence of bioirrigation as bottom waters became almost anoxic (2 µM O 2 ).No polychaetes were observed in these cores due to mortality or migration to shallower waters and dead polychaetes were observed on the sediment surface.The Br − tracer concentration at this time was indistinguishable from the diffusiononly profiles, as observed on a previous visit to this site (Forster et al., 1999).Interestingly, the sediments became more compacted between August and October with porosity at 2 cm decreasing from an average of 0.91 in winter to 0.88 (Fig. 3).This 30 % increase in solid volume fraction (1 − ϕ) is presumably a result of the partial collapse of relic burrows.Bioirrigation rates did not recover in November and December even though the bottom waters were again well oxygenated due to mixing of the water column by autumnal storms (Table 1).No polychaetes were recovered from the sediments in November and December.This confirms previous observations that a time lag of at least 2 months is required following hypoxia before the sediments at Boknis Eck become habitable to macrofauna (Meyer-Reil, 1983).
Despite the lack of irrigation from September onwards, the concentration gradients of TCO 2 , NH + 4 and TA did not accumulate in response to temperature-induced increase in organic matter degradation as expected from time series data elsewhere (e.g.Klump and Martens, 1989).On the contrary, TCO 2 concentrations decreased from 12 to 7 mM at 20 cm depth between September and October (Fig. 4a), NH + 4 fell by a factor of 2 from 1500 µ M to 700 µM (Fig. 4f), and TH 2 S decreased by at least a factor of 4 (Fig. 4i).Conversely, SO 2− 4 concentrations at depth began to increase at this time (Fig. 4c).This rate of concentration change is much faster than can be explained by molecular diffusion and resembles deep irrigation.Shortly afterwards, extremely high SR rates were measured in two of the triplicate cores in November (Fig. 4e).This perturbation in solute concentrations was not, however, an isolated incident, and similar trends were observed for NH + 4 , TPO 4 and SO 2− 4 in May.Fe 2+ and TPO 4 behaved quasi-independently of the other solutes (Fig. 4g, h).Fe 2+ exhibited a peak approaching 194 µM on 23 February, which then diminished in March and disappeared completely by July.In contrast to Fe 2+ , TPO 4 concentrations were maintained at low levels in the bioirrigated zone in winter.From August onwards, the Fe 2+ peak slowly reappeared and extended downwards, reaching 180 µM at 1 cm depth in October and 311 µM at 3 cm depth in November.Furthermore, the entire Fe 2+ peak extended over the upper 10 cm in November compared to 2 cm in October.Similarly, extremely high TPO 4 concentrations in excess of 437 µM were measured in the upper 6 cm and 10 cm in October and November, respectively, when O 2 concentrations fell below 10 µM.The increase in vertical extension of the Fe 2+ and TPO 4 peaks between October and November coincides with a depletion of sulfide to levels below detection limit (Fig. 4i).
The modeled concentration and SR rate profiles are shown alongside the measured data in Fig. 4 (solid curves).Imposing the transient forcings and boundary conditions in Fig. 2, the model reproduces the major trends in the data to a high degree with the exception of the TPO 4 concentrations in the upper 10 cm during the hypoxic season.These include (1) the constant concentrations in the upper 10 cm in winter, (2) the weakening of the concentration gradients in May and October/November, and (3) the iron peaks in winter and autumn.The first two of these observations are consistent with nonlocal solute exchange by bioirrigation and bubble irrigation, respectively.

Benthic-pelagic solute exchange processes and rates
The near-constant concentrations of SO 2− 4 , TA, TCO 2 , NH + 4 , TPO 4 and TH 2 S in the upper 10 cm in February at values equal to those measured in the bottom waters is strong evidence for non-local exchange of porewater with seawater by the pumping action of indwelling fauna.This process ameliorates the classic biogeochemical fingerprint of organic matter degradation in the surface sediment layers (Dale et al., 2011).Sediments only became sulfidic below 10 cm in the sulfate reduction zone, which corroborates previous measurements using microelectrodes (Preisler et al., 2007).Between February and April, the depth of irrigation decreased by more than a factor of 2, allowing TCO 2 , NH + 4 , PO 3− 4 and TH 2 S to accumulate and SO 2− 4 to become more extensively depleted.More intense bioirrigation at the start of winter may be explained by tube building activities of the pioneering polychaete population following the previous hypoxic season (Meyer-Reil, 1983).
Attributing the transient nature of these concentration profiles to bioirrigation rather than geochemical processes is consistent with the model results that closely simulate the data using the rates of bioirrigation derived from the Br − incubations.The irrigation coefficients determined for winter and spring are similar to the higher end irrigation coefficients of ca.0.38 d −1 reported for the southern North Sea and Kiel Bight (Forster et al., 2003;Schlüter et al., 2000).The biogeochemical model also simulates the natural Br − abundance in February (Fig. 4d).The data-model disagreement for Br − at 10 cm depth is likely due to spatial differences in biomass abundances in the cores sampled for geochemistry and incubation experiments.In addition, NH + 4 and NO − 3 fluxes measured ex situ during winter by Dale et al. (2011) were manyfold higher than the molecular diffusive fluxes.Enhanced diffusion of this type is often ascribed to bioirrigation (Glud et al., 1994).
When the modeled daily solute fluxes across the sedimentwater interface are compared, it becomes clear that bioirrigation accounts for most of the transport to and from the sediment in winter at Boknis Eck (Fig. 5).Even though the relative importance of bioirrigation tails off as the year progresses, non-local exchange by gas bubbling clearly dominates transport for some species at discrete intervals beginning in May (e.g.TCO 2 , NH + 4 and TH 2 S), with two further ebullition events in September and late October.The bubble irrigation rate constant determined from the solute concentration profiles is an order-of-magnitude lower than that for bioirrigation (Fig. 2h).Yet, because it extends over a much greater depth (∼ 30 cm) where TCO 2 , NH + 4 and TH 2 S accumulate at high concentrations, it has a profound effect on their fluxes across the sediment-water interface.The occurrence of TH 2 S in the bottom waters during the anoxic period (Balzer et al., 1983;Ehrhardt and Wenck, 1984;Bange et  al., 2010) could well be a result of bubble-associated release.Porewater-seawater exchange by rising bubbles has the potential to shift a large part of the oxygen demand, supplied initially as POM, from the sediment to the water column in the form of TH 2 S and CH 4 .How these fluxes will develop in the future and their potential to intensify or prolong hypoxia is a topic for future investigation.
Highly reactive species such as O 2 and NO − 3 are mainly consumed in the surface sediment layers, which creates a high diffusive flux through the diffusive sub-layer.Averaged over the year, 60 % of O 2 and NO − 3 enter the sediment by diffusion, compared to < 10 % for bubble irrigation.Up to 75 % of POC is degraded aerobically in winter, which reflects the presence of reactive carbon within the surface oxidized layer (Fig. 6a).Overall, POC mineralization oscillates between aerobic respiration in winter and spring to sulfate reduction in summer and autumn, with other carbon degradation pathways playing a minor role.
A large spike in O 2 consumption of 26 mmol m −2 d −1 occurs prior to the spring bloom (Fig. 6b).Following the spring bloom, oxygen uptake is lower (15 mmol m −2 d −1 ).Interestingly, total oxygen uptake (TOU) from 1 January to the start of the spring bloom in March is 1100 mmol m −2 , compared to a POC deposition flux of only 800 mmol m −2 .This difference represents the oxygen-debt carried over from the previous hypoxic season when organic matter, pyrite and other reduced substances accumulated.This debt is quickly burned off in the following winter when deepwater renewal and invasion of irrigating organisms once again ventilate the upper sediments.The modeled TOU at this time of year agree with benthic respiration rates of 16 to 24 mmol m −2 d −1 made at the same time using landers (Sommer et al., unpublished data).
When O 2 levels are diminished during hypoxia, diffusive fluxes of nitrate become elevated (Fig. 5).At this time, the most reactive POM fractions are consumed by denitrification, which reaches maximum levels of 0.25 mmol m −2 d −1 in September when bottom water O 2 is 2 µM (Fig. 6c).Denitrification consumes only a few percent of POC (Fig. 6a) compared to 10 % on the continental shelf (Middelburg et al., 1996).This is because bottom water NO − 3 concentrations are only 5-6 µM at Boknis Eck compared to > 30 µM on the shelf.Furthermore, the low O 2 concentrations limit the potential for in situ NO − 3 production by coupled nitrificationdenitrification.
Following the peak in denitrification, DNRA rates increase sharply to 0.28 mmol m −2 d −1 in November (Fig. 6c).This is far less than 3.0 mmol m −2 d −1 for the Peruvian OMZ and 16 mmol m −2 d −1 for the highly reactive sediments of the Namibian shelf inhabited by giant sulfide oxidizers (Bohlen et al., 2011;Dale et al., 2009a).Presumably, more extensive periods of anoxia in these environments permit a higher build-up of biomass, which could explain why higher nitrate transport coefficients of 0.7 and 1.1 d −1 were calculated for

Beggiatoa
Aerobic respiration Denitrification + DIR Sulfate reduction Methanogenesis these environments (Bohlen et al., 2011;Dale et al., 2009a) compared to maximum values of 0.12 d −1 in Boknis Eck (Fig. 2d).The NH + 4 produced by DNRA exits the sediment by diffusion, leading to maximum simulated diffusive fluxes at this time of year which agree well with the measured values (Fig. 5).

Evidence for solute exchange by bubble irrigation
The experiments and modeling show that bioirrigation can be ruled out as the principle cause for the unusual shallowing of the solute concentration gradients in October/November and, by extension, in May.The alternative hypothesis concerning bubble irrigation obviously requires independent experimental verification, but adds to a growing body of evidence which suggests that this process could be widespread in gassy coastal sediments.Irrigation of sediment porewater by rising gas bubbles was identified and quantified in Cape Lookout Bight sediments several decades ago (Martens, 1976;Martens and Klump, 1980;Klump and Martens, 1981).These studies have highlighted bubble transport as the major source of pelagic CH 4 in shallow water columns during summer.Roden and Tuttle (1992) also suggested that bubble irrigation was responsible for rapidly changing sulfide concentration in anoxic Chesapeake Bay sediments.More recently, meter-scale irrigation patterns of SO 2− 4 in sediments overlying gas hydrates have been attributed to ascending gas bubbles (Haeckel et al., 2007;Schwalenberg et al., 2010).It is important to note that hydrodynamically induced flushing of relic burrows by bottom currents (Ray and Aller, 1985) can be ruled out for Boknis Eck since the rapid changes in solute concentrations extend much deeper than the bioirrigation depth (ca. 10 cm).
There are several lines of evidence in support of CH 4 gas escape at Boknis Eck.To begin with, free gas is a common occurrence in the sediments throughout Eckernförde Bay, including Boknis Eck (Abegg and Anderson, 1997;Orsi et al., 1996), and bubble release has been observed in some parts of the bay (Jackson et al., 1998).Free gas occurs when the dissolved concentration exceeds the in situ solubility, and the depth at which gas forms depends on the mud thickness and the amount and reactivity of organic carbon contained within it (Dale et al., 2009b;Jensen and Bennike, 2009).Model simulations suggest that CH 4 saturation is reached at ca. 1 m below the seafloor at Boknis Eck (Dale et al., 2011), although some degree of seasonal variability in this depth is to be expected (Wever and Fiedler, 1995;Mogollón et al., 2011).
In agreement with these studies, Bange et al. (2010) reported episodic CH 4 concentrations in excess of 230 nM in the bottom waters of Boknis Eck (25 m) during 2006-2008.Interestingly, when aggregating the data, a pattern emerged of highest CH 4 concentrations in early autumn (September-October) with another, albeit less obvious, concentration spike in late spring (March-April).This dovetails with the porewater data that also shows evidence for strong bubble release during the same periods, and may further suggest that CH 4 release is triggered by a periodic forcing.Bange et al. (2010) postulated that the deposited spring and autumn blooms was the source of the released CH 4 .The porewater data presented here instead suggest that the CH 4 originates from the gas-bearing layers located well below the surface sediment layers where fresh phytodetritus is degraded.The autumn CH 4 release is more likely to be instigated by the high bottom water temperature, which rose to maximum values of 10 • C at this time of year (Table 1).If the mean annual temperature is 5 • C, then a 5 • C temperature increase in October would decrease CH 4 solubility by roughly 1 mM (Duan et al., 1992).Using the equations given in Martens et al. (1998) and a porosity of 0.8, the evolved CH 4 gas would increase the volume of free gas by ca.0.5 % for the ambient conditions of Boknis Eck.This is around one third of the peak gas volume of 1.5 % predicted for this site (Dale et al., 2011).Such an increase could lead to gas overpressure, degassing and non-local mixing of porewater with bottom water.In May the temperature was relatively low (3.8 • C, Table 1), and we postulate that excess CH 4 gas that accumulated over the winter may be released in spring as a sediment "burp" as temperatures start to rise.Changes in hydrostatic pressure due short-term sea level variations are also known to trigger CH 4 release from the sea floor (Martens and Klump, 1980).
Following the irrigation event in October, the SR rates in two of the triplicate sub-cores in November were extremely elevated below 20 cm (Bertics et al., 2012), being up to 10 times higher than those measured at the same depth in previous months (Fig. 4e).We speculate that a fraction of this SR is associated with anaerobic oxidation of methane (AOM, R AOM in Table 2).A large discharge of gas through sediment conduits in October could have opened a pathway for smaller pockets of CH 4 gas to advect slowly upwards as the system relaxed following the initial eruption.Lateral dissolution of gas through the walls of cracks and fissures would be assisted by CH 4 consumption by AOM which causes the porewater to become undersaturated with CH 4 (Mogollón et al., 2011).Further recent evidence from methane vents suggests that methanotrophic microbial biofilms could exist directly on the walls of sediment bubble tubes (Briggs et al., 2011).Potential AOM rates of 100 nmol cm 3 d −1 have been measured at this site (Treude et al., 2005) indicating that CH 4 is not the only electron donor for SR.A pocket of extremely labile organic material may have been sampled, possibly consisting of dead polychaetes (Bertics et al., 2012).Despite these uncertainties, it is important to note that the steepening of the SO 2− 4 , TA and TCO 2 gradients between November and December occurs far too quickly to be explained by molecular diffusion.This alludes to a modification of solute concentrations by a rapid increase in the activity of sulfate reducing bacteria.

Iron and phosphorus remobilization and fluxes
The transient Fe 2+ peaks in winter and autumn represent a decoupling between Fe 2+ source-sink processes (Fig. 4h).The winter maximum was predicted to be driven by reoxidation of iron sulfides.Renewed deep ventilation of the sediments in winter by faunal activity exposes to oxidation reduced iron that had accumulated during the previous hypoxic season (Fig. 6b).Given that O 2 penetration into the sediment is < 1 cm in winter (Dale et al., 2011), deep ventilation by bioirrigation drives most of the pyrite oxidation in the upper 10 cm.Because pyrite oxidation is not associated with TPO 4 release, TPO 4 and Fe 2+ concentrations below the oxidized layer were decoupled in winter and no TPO 4 peak was observed (Fig. 4g).
An additional source of Fe 2+ that was not considered in the model is from siderite (FeCO 3 ).Siderite can be associated with pyrite in fully marine and brackish sediments in micro-environments that allow both minerals to be thermodynamically stable (Haese et al., 1997;Postma, 1982).The decrease in alkalinity during winter could have induced siderite to dissolve.The quantitative significance of siderite in surface marine sediments is generally believed to be minor, however, since its formation is inhibited by low levels of sulfide (Haese, 2006).
Despite the source of Fe 2+ , the diffusive fluxes calculated from the porewater gradients at the sediment surface remain relatively low in winter at 0.2 ± 0.05 mmol m −2 d −1 (Fig. 5).Fe 2+ that diffuses into the surface oxidized layer is rapidly re-oxidized to iron oxyhydroxides (Pakhomova et al., 2007;Severmann et al., 2010).Sequestration of TPO 4 onto freshly precipitated iron oxyhydroxides is a major control on TPO 4 fluxes in environments with high bottom water O 2 concentrations (Froelich et al., 1979(Froelich et al., , 1988;;Krom and Berner, 1981;McManus et al., 1997).This authigenic P sink explains the low diffusive TPO 4 fluxes in winter (Fig. 5) and leads to extremely high C/P ratios of 710 for the diffusive fluxes across the sediment-water interface (Fig. 7a).The C/P ratio still remains above Redfield (113) in March.Despite this efficient trapping of TPO 4 , the intense bioirrigation occurring at Boknis Eck in winter creates a bypass for TPO 4 to enter the water column with maximum predicted rates of 0.5 mmol m −2 d −1 (Fig. 5).This is likely a maximum value due to the potential adsorption of TPO 4 onto iron hydroxides lining animal burrow walls (see below).
In autumn, when the surface oxidized layer was almost absent, diffusive Fe 2+ fluxes increased by a factor of 4 from winter to 0.8 ± 0.4 mmol m −2 d −1 (Fig. 5).The TPO 4 diffusive fluxes calculated from the porewater gradients in October and November were 0.7 ± 0.2 mmol m −2 d −1 and 0.6 ± 0.1 mmol m −2 d −1 , respectively (Fig. 5).These are identical to fluxes reported by Balzer et al. (1983) for low O 2 conditions (< 1 µM) using a benthic chamber deployed at the same site.Enhanced P release leads to low C/P diffusive flux ratios of ca. 9 in September (Fig. 7a), indicating extreme   2c) and the stoichiometry of the TCO 2 to TPO 4 flux (C/P) by molecular diffusion at the sediment-water interface calculated using the measured data.The horizontal black line is the Redfield C/P ratio (106).Note the axis break in the ordinate.(b) TPO 4 concentrations and dissolved inorganic nitrogen to phosphorus atomic ratios (N/P) in the bottom water over the study period (data from Bange et al., 2011).The period with hypoxic bottom water (O 2 < 63 µM) is indicated with the grey shaded area.
P enrichment of the fluxes relative to TCO 2 .A TPO 4 flux of 0.7 mmol m −2 d −1 over the 40 day hypoxic period (Fig. 5) would raise TPO 4 concentrations in the overlying water column (28 m) by around 1 µM.A potential further increase in TPO 4 flux is provided by the bubble release pathway.In fact, the TPO 4 concentration measured in the bottom water increases by 3 to 4 µM between May and September (Fig. 7b) coincident with the 2 major bubble release events (Fig. 2h) and the high benthic TPO 4 diffusive fluxes.Excess benthic P regeneration decreases the dissolved inorganic N/P ratio in the water column to extremely low levels of 1 to 2 (Fig. 7b), which may affect local productivity patterns as studies in other shallow ecosystems indicate (Rozan et al., 2002).The feedbacks on primary production and the timing of the autumn bloom in Boknis Eck by benthic P release are suspected to be important but have not been investigated quantitatively (Smetacek et al., 1987).
As in winter, the presence of a peak in free Fe 2+ during the hypoxic period suggests that the rate of pyrite formation is slow relative to reductive iron dissolution (Fig. 4h).We were only able to simulate the autumn Fe 2+ peak by scaling the rate for reductive dissolution of iron oxyhydroxides to the time-lagged O 2 concentration (Fig. 2c, Table 2), effectively decoupling this reaction from FeS 2 precipitation.Otherwise, FeS 2 precipitation keeps pace with reductive dissolution of iron oxyhydroxides and maintains Fe 2+ at low concentrations (data not shown).
This dependence of iron reduction on O 2 could reflect a missing or poorly represented process in the model.Although the time-lagged O 2 concentration appears to be related to enhanced nitrate transport by Beggiatoa, we are hesitant about linking these two aspects without further supporting data.The apparent dependence on O 2 may be an artifact reflecting the simple 1-D bioirrigation function that fails to account for geochemical gradients around animal burrows (Meile et al., 2005;Meysman et al., 2006).Inhabited burrows are ventilated by the pumping action of the host organism, leading to a thin veneer of iron oxyhydroxides on the exposed burrow walls (Meile et al., 2005).Dissolution of iron oxyhydroxides following vacation of the burrows in response to hypoxia could increase Fe 2+ concentrations in the burrow water, leading to the situation where horizontal chemical gradients exceed those in the vertical (Lewandoski et al., 2007;Meile et al., 2005).This heterogeneity is destroyed when porewater for geochemical analysis is extracted by squeezing the sediment or by using rhizons inserted laterally through the core.A test of the potential importance of iron cycling within animal burrows requires a higher dimensional model that explicitly considers tube geometry (Meile et al., 2005).
At this stage, only a causal relationship between the rate of iron reduction and O 2 concentrations can be established.Nevertheless, the important point is that the simulated Fe 2+ peak in November represents the net uptake of Fe 2+ into sulfidic minerals and Fe 2+ release from dissolution of iron oxyhydroxides.Reductive dissolution also releases Fe-P.Yet, as the model simulation shows (Fig. 4g, solid curves), the TPO 4 profiles are not reproduced using literature P : Fe ratios of ε = 0.1 (Slomp et al., 1996;Anschutz et al., 1998).The integrated TPO 4 concentration contained within the November peak is 2.5 times as high as the Fe 2+ peak, implying that ε ≈ 2.5.There is obviously a much more significant source of TPO 4 than the fraction which can be reasonably explained as Fe-P.We were unable to simulate the TPO 4 peak assuming enhanced release of Fe-P or preferential mineralization of organic P relative to organic C (not shown).A similarly sized peak of 250 µM measured in Namibian shelf sediments was explained by TPO 4 release rates of 500-1700 nmol cm 3 d −1 from Thiomargarita bacteria (Schulz and Schulz, 2005).As we suggest below, our simulations also point toward an intense yet ephemeral microbial source of TPO 4 in the upper 5 cm when O 2 concentrations dropped below 10 µM.

Massive phosphate release by Beggiatoa?
Luxury uptake and release of TPO 4 under oscillating redox conditions by the giant sulfide oxidizing bacteria Beggiatoa and Thiomargarita has been reported (Schulz and Schulz, 2005;Brock and Schulz-Vogt, 2010;Goldhammer et al., 2010).These microorganisms convert porewater phosphate into intracellular polyphosphate granules under oxic conditions.Under anoxic conditions these compounds are hydrolyzed to release energy as well as TPO 4 to the surrounding porewater (Sannigrahi and Ingall, 2005;Schulz and Schulz 2005;Brock and Schulz-Vogt, 2010).The metabolic stimulus and motive for phosphate release is poorly understood and beyond the scope of this study.Some laboratory studies indicate that polyphosphate breakdown may be part of an auxiliary metabolism arising from increased exposure to sulfide due to higher rates of sulfate reduction (Brock and Schulz-Vogt, 2010).
Considering that nitrate-storing Beggiatoa filaments were observed in great abundance during hypoxia in Boknis Eck, it seems reasonable to consider the possibility that the TPO 4 peak in October/November may be driven by these bacteria in addition to the release of Fe-P and organic P. Ingall and Jahnke (1994) and others (e.g.Sannigrahi and Ingall, 2005;Hupfer and Lewandowski, 2008) noted early on that TPO 4 released from iron oxyhydroxides and organic P under low oxygen conditions could only account for a fraction of the total rate of TPO 4 regeneration in sediments, and ascribed the rest to TPO 4 release by microorganisms.More generally, microorganisms are believed to play a major role in benthic phosphorus cycling (Gächter and Meyer, 1993;Hupfer et al., 2004).
To test this idea, we modified the model to include intracellular P as a new dynamic variable (Fig. 1).The basic premise is that luxury uptake of TPO 4 from the porewater by Beggiatoa takes place during periods of oxygenated bottom waters, whereas net release of intracellular P to the porewater as TPO 4 occurs during hypoxia.Although Beggiatoa could conceivably take up TPO 4 directly from the bottom water by extension of their filaments through the diffusive boundary layer, it seems more plausible that they instead capitalize on the porewater TPO 4 pool where concentrations are at least an order-of-magnitude higher.Thus, if Beggiatoa can be considered as viable candidates for enhanced phosphorus cycling in Boknis Eck, then the TPO 4 added to the porewater in the uppermost sediment layers by mineralization processes over the year has to be sufficient to sustain the TPO 4 peak when taken up by the bacteria and then released again in autumn.
As a minimum requirement, we assumed that microbial phosphate uptake (R P up , M yr −1 ) and release (R P rel , M yr −1 ) only occurs if the Beggiatoa are active.This was achieved by making R P up and R P rel dependent on the rate of DNRA (R DNRA ): where [TPO 4 ] bac is the intracellular phosphate concentration and K TPO 4 and K bac TPO 4 are kinetic constants for phosphate uptake and release.These constants are assigned a value of 10 µM to ensure uptake and release of TPO 4 at the low porewater concentrations observed.Because these processes relate to Beggiatoa, they employ the time-lagged O 2 concentration, O 2 , in the same way as bacterial nitrate transport (Fig. 2d).The oxygen limiting term includes the halfsaturation constant K P O 2 , which determines the concentration at which Beggiatoa become net releasers or assimilators of TPO 4 .The rates are mostly sensitive to this parameter, and a value of 0.8 µM best describes the data.χ up and χ rel are dimensionless parameters equal to 10 and 100, respectively, that simply decouple intracellular nitrate and phosphate dynamics and allow for a faster rate of intracellular TPO 4 release during hypoxia.
Model simulations considering an intracellular P pool provide a much better correspondence with the TPO 4 concentrations in November, but there is no improvement in October (Fig. 4g, dashed lines).However, since the rates of P uptake and release are sensitive to O 2 concentration, the October TPO 4 peak can be simulated by allowing for a faster rate of O 2 decrease between August and September than shown by the linear interpolation used in Fig. 2c (not shown).The simulations would therefore benefit from a higher sampling resolution of bottom water O 2 during this time of year.
The model reproduces the TPO 4 diffusive effluxes that accompany the peaks in October and November, with maximum fluxes of 1.2 mmol m −2 d −1 (Fig. 5).This maintains low C/P ratios of the diffusive fluxes of 20 even though bottom water O 2 returns to high levels (Fig. 7a).TPO 4 release rates of 0.8 and 0.6 mmol m −2 d −1 have been reported for the seasonally hypoxic Arkona Basin and Bay of Concepción (Holmkvist et al., 2010;Mort et al., 2010).Benthic fluxes reported for quasi-permanent hypoxic shelf regions are at least an order of magnitude lower, including California borderland basins (Berelson et al., 1987;Jahnke, 1990;McManus et al., 1997) and the Washington and Oregon margins (Devol and Christensen, 1993;Severmann et al., 2010).At the highly productive Peruvian oxygen minimum zone, TPO 4 fluxes are higher at 0.4 mmol m −2 d −1 .In general, it is the oscillation between oxic and hypoxic/anoxic conditions rather than permanent anoxic conditions that favors periods of short but intense bursts of TPO 4 release to the water column (Balzar et al., 1983;Koop et al., 1990;Cowan and Boynton, 1996;Mort , 2010).In this sense, Beggiatoa and/or iron-associated P may act as phosphorus capacitors in systems with oscillating redox conditions, releasing large amounts of TPO 4 in a short space of time and dramatically increasing the internal loading of TPO 4 to overlying waters.The pulsed release of TPO 4 by Beggiatoa (137 mmol m −2 yr −1 ) dominates the annual benthic P budget (Table 4).P release from organic matter and Fe-P contribute similar but much smaller amounts (42-47 mmol m −2 yr −1 ).Peak TPO 4 release rates from Beggiatoa of 4.5 mmol m −2 d −1 are about 10 times larger than the TPO 4 flux supplied by Fe-P and organic P combined (Fig. 6d).The bacteria need only relatively low TPO 4 assimilation rates of 0.2 mmol m −2 d −1 during winter and spring to accumulate enough TPO 4 to sustain the autumn release (not shown).This could be easily supplied by organic P and Fe-P (Fig. 6d; Table 4).For comparison, depthintegrated TPO 4 uptake rates of 0.3 to 1.2 mmol m −2 d −1 have been reported for anoxic laboratory incubations using Beggitaota and Thiomargarita (Goldhammer et al., 2010).Our model predicted rates are thus very much in agreement with the few data available on P cycling by giant sulfide oxidizers and provide more evidence for a significant role of these microorganisms in moderating P fluxes in seasonally hypoxic settings.

Conclusions
This study presents a time series of porewater data from a seasonally hypoxic coastal setting which is analyzed using a numerical model to quantify benthic processes and fluxes.Clear transient signals in the data are interpreted as changes in both reaction rates and transport processes (mainly bioirrigation and bubble irrigation).If the model predictions are correct, then the patterns of nutrient cycling in Boknis Eck sediments may be more widely applicable.
Firstly, we have shown that solute exchange between the sediments and the water column is dominated by bioirriga-tion in winter and spring and by bubble-induced irrigation at irregular intervals from late spring to autumn.Molecular diffusion remains important for highly reactive species such as O 2 and NO − 3 .Secondly, a very large transient peak in TPO 4 in excess of 400 µM develops in subsurface sediments at the end of the hypoxic period in October.This leads to an extreme enrichment of P relative to C in the fluxes from the sediment and coincides with a dramatic lowering of the N/P ratio in the overlying water column.This is mainly driven by a large pulsed release of TPO 4 to the porewater within a short space of time (∼ weeks).Based on published field evidence and laboratory experiments, we were able to reproduce this peak using new kinetic expressions that allow the giant sulfide-oxidizing bacteria Beggiatoa to take up TPO 4 during periods of oxic bottom waters and release TPO 4 when oxygen drops to low levels (∼ 10 µM).This does not require an external source of TPO 4 and agrees with the known distribution of Beggiatoa in these sediments.Alternative sources of TPO 4 , such as from Fe-P and organic P, do not provide a satisfactory simulation of the data, but cannot be ruled out indefinitely at this stage.
Further work is required to substantiate our findings.Although there are several pieces of circumstantial evidence in support of pseudo-irrigation by gas bubble transport, this has yet to be shown experimentally.Given the apparently episodic nature of gas eruption, continuous monitoring of gas escape and sediment porewater would be needed to capture these events and corroborate the hypothesis.Bubble irrigation may be more common than currently presumed and is likely to be overlooked if organic-rich sediments are sampled at a single time rather than as part of a time series study.Our suggestion that Beggiatoa are releasing TPO 4 would be strengthened with data on the temporal variability of the fraction of TPO 4 bound to iron oxyhydroxides as well as intracellular P concentrations.The mechanism and regulation of P storage in bacteria is poorly understood, and at this stage it is not known whether the strain of Beggiatoa at Boknis Eck is capable of intermediate polyphosphate storage.Sampling at bi-weekly intervals with additional support from benthic chamber deployments during the hypoxic period would also help to refine future modeling work.The limitation of the 1-D model approach to simulate iron precipitation and dissolution on the walls of animal burrows is an additional aspect which will be addressed in the future.
The strength of this study lies in the quantification of transient fluxes and process rates that are supported by a set of temporally resolved benthic data.We demonstrate that the accurate nutrient budgeting in muddy coastal sediments requires a careful consideration of temporal dynamics.Nevertheless, there are very few published modeling studies that explore seasonal benthic cycling that are constrained using data collected at sampling frequencies adequate to capture the transient dynamics (e.g.Aller, 1977;Fossing et al., 2004;Klump and Martens, 1989).This is especially true for phosphorus dynamics in seasonally hypoxic or anoxic systems www.biogeosciences.net/10/629/2013/Biogeosciences, 10, 629-651, 2013 with pronounced redox oscillations.These sediments may behave as phosphate capacitors, whereby particles (e.g.iron oxyhydroxides and/or bacteria) are charged with TPO 4 during oxic conditions and then discharge TPO 4 when hypoxia ensues.More studies in environments such as Boknis Eck with a predictable hypoxic period could prove invaluable to better constrain the kinetics of key reactions and improve the predictive power of basin scale models of systems with seasonal or episodic hypoxia such as the Baltic Sea (e.g.Eilola et al., 2011).

Fig.
Fig.1.Schematic of the coupled Fe and P cycles considered in the model.Solid species are in squares and solutes in circles.P is added to the sediment as organic P and with iron oxyhydroxides.Phosphate is released to the porewater upon reductive dissolution of ferric iron and oxidation of organic matter.This dissolved P pool is free to be transported out of the sediment, or removed from the porewater by co-precipitation with ferrous iron oxidation.P uptake and release by Beggiatoa is explored in Sect.4.4.

Fig. 2 .
Fig. 2. Time-dependent boundary conditions and forcing functions used in the model.(a) Temperature, (b) bottom water sulfate and bromide, (c) bottom water dissolved oxygen and the time-lagged values (O 2 ) used in several rate expressions (P uptake and release by Beggiatoa, Fe(OH) 3 reduction by sulfide and nitrate transport by sulfide oxidizing bacteria), (d) nitrate transport by sulfide oxidizing bacteria defined by the equation shown.Panel (e)shows the G0 and G1 organic matter fluxes, (f) is the bioturbation coefficient defined by the equation given in the panel, (g) shows the bioirrigation coefficients (α bi1 and α bi2 ) determined from the ex situ experiments, and (h) is the bubble irrigation coefficient (α bu ) and enhanced advection of dissolved methane (v a (0)).The symbols in (a), (b), (c) and (g) denote values that have been measured except from those at the beginning and the end of the year which are interpolated.Linear interpolations were used between the data for modeling purposes.

Fig. 3 .
Fig. 3. (a) Bromide tracer concentrations (symbols) measured in cores incubated for bioirrigation experiments sampled on the day in 2010 shown on the top of the figure.The curve represents the modeled tracer concentration at the end of the incubation period using the parameters indicated (Eq.2), and the shaded interval shows the expected bromide profile assuming transport by molecular diffusion only.For each model run, the modeled porosity curves shown in (b) were used, defined by the parameters given in each panel and the equation ϕ = ϕ(L) + (ϕ(0) − ϕ(L)) exp (−z/z por ) (see Supplement).

Fig. 4 .
Fig. 4. Simulated (curves) and measured (symbols) geochemical concentrations and SR rates in Boknis Eck sediments sampled on the dates indicated at the top of the figure.Open and closed circles indicate porewater measurements from adjacent multi-cores sampled by squeezing the sediments.Stars indicate porewater extracted by rhizons.The error bars in (e) indicate the standard deviation of 3 replicate SR rates sampled from the same core.SR rates in November are shown for the individual triplicate cores to highlight the very high rates measured in two of the sub-cores (note the change in scale).The modeled SR is the sum of sulfate consumption by organic carbon and methane oxidation (0.5 × R SO 4 + R AOM ).The numbers 1 to 3 in black circles highlight examples of the salient seasonal features (1 = constant solute concentrations in the upper 10 cm in winter, 2 = weakening of the concentration gradients in May and October/November, and 3 = the Fe 2+ and TPO 4 peaks in winter and autumn).For an explanation of the dashed curves in (g) see Sect.4.4.

Fig. 5 .
Fig. 5. Modeled solute fluxes across the sediment-water interface (mmol m −2 d −1 ) by molecular diffusion, bioirrigation and bubble irrigation.The black circles denote the mean (±1σ ) diffusive fluxes calculated by applying Fick's law to the measured data.The dashed line in the NO − 3 plot denotes intracellular nitrate transport by Beggiatoa.Negative fluxes are into the sediment.Modeled phosphate fluxes correspond to the simulation considering P uptake and release by Beggiatoa (see text).

Fig. 6 .
Fig. 6.(a) The daily fraction (%) of POM which is mineralized by each diagenetic pathway and the total rate of POM mineralization ( R POC in mmol m −2 d −1 of C shown by the white line).(b) Rates of O 2 consumption by different electron donors and TOU (mmol m −2 d −1 of O 2 ).(c) Selected N cycling pathways (mmol m −2 d −1 of N).Denitrification corresponds to the N 2 -producing pathway.(d) Biogeochemical reactions producing TPO 4 in the sediment (mmol m −2 d −1 of TPO 4 ).The period with hypoxic bottom water (O 2 < 63 µM) is indicated with the grey shaded area and the black arrows indicate the start of the spring and autumn blooms.

Fig. 7 .
Fig. 7. (a) Oxygen concentrations (from Fig.2c) and the stoichiometry of the TCO 2 to TPO 4 flux (C/P) by molecular diffusion at the sediment-water interface calculated using the measured data.The horizontal black line is the Redfield C/P ratio (106).Note the axis break in the ordinate.(b) TPO 4 concentrations and dissolved inorganic nitrogen to phosphorus atomic ratios (N/P) in the bottom water over the study period (data fromBange et al., 2011).The period with hypoxic bottom water (O 2 < 63 µM) is indicated with the grey shaded area.

Table 2 .
Reaction network used in the biogeochemical model.

Table 3 .
Important parameters used in the biogeochemical model.Time-dependent parameters denoted "variable" are shown in Fig.2.All other parameters are listed in the Supplement.

Table 4 .
Yearly integrated benthic sources and sinks of phosphate for the upper 40 cm in Boknis Eck (mmol m −2 yr −1 ) listed from largest to smallest.Sources are positive and sinks are negative.Values have been rounded to the nearest integer.