Large variations in iron input to an oligotrophic Baltic Sea estuary: impact on sedimentary phosphorus burial

Estuarine sediments are key sites for removal of phosphorus (P) from rivers and the open sea. Vivianite, an Fe(II)-P mineral, can act as a major sink for P in Fe-rich coastal sediments. In this study, we investigate the burial of P in the Öre Estuary in the northern Baltic Sea. We find much higher rates of P burial at our five study sites (up to ∼0.145 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) when compared to more southern coastal areas in the Baltic Sea with similar rates of sedimentation. Detailed study of the sediment P forms at our site with the highest rate of sedimentation reveals a major role for P associated with Fe and the presence of vivianite crystals below the sulfate methane transition zone. By applying a reactive transport model to sediment and porewater profiles for this site, we show that vivianite may account for up to ∼40 % of total P burial. With the model, we demonstrate that vivianite formation is promoted in sediments with a low bottom water salinity and high rates of sedimentation and Fe oxide input. While high rates of organic matter input are also required, there is an optimum rate above which vivianite formation declines. Distinct enrichments in sediment Fe and sulfur at depth in the sediment are attributed to short periods of enhanced input of riverine Fe and organic matter. These periods of enhanced input are linked to variations in rainfall on land and follow dry periods. Most of the P associated with the Fe in the sediment is likely imported from the adjacent eutrophic Baltic Proper. Our work demonstrates that variations in land-to-sea transfer of Fe may act as a key control on burial of P in coastal sediments. Ongoing climate change is expected to lead to a decrease in bottom water salinity and contribute to continued high inputs of Fe oxides from land, further promoting P burial as vivianite in the coastal zone of the northern Baltic Sea. This may enhance the role of this oligotrophic area as a sink for P imported from eutrophic parts of the Baltic Sea.

1 Introduction

Phosphorus (P) is an important nutrient for primary producers. Burial of reactive P (i.e., bioavailable P) in coastal sediments can permanently remove P from the water column (Delaney, 1998; Froelich et al., 1982; Ruttenberg, 2003). This removal allows coastal systems to act as filters for P, reducing the flux of land-derived P to the open sea (Asmala et al., 2017; Bouwman et al., 2013; Froelich, 1988), and as sinks of P imported from the open sea (Asmala et al., 2017). The major sedimentary phases contributing to the removal of reactive P in coastal sediments are organic P, iron-bound P and authigenic calcium (Ca)-bound P (Ruttenberg and Berner, 1993; Slomp et al., 1996a).

In most coastal systems, organic matter is a major carrier of P to the sediment (Slomp, 2011). The subsequent microbial degradation of part of the organic matter can lead to elevated concentrations of phosphate (${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$) in the porewater and precipitation of authigenic P minerals (Froelich et al., 1988; Ruttenberg, 2003). While, until recently, carbonate fluorapatite was thought to act as the major sink of P in coastal systems (Ruttenberg and Berner, 1993; Slomp et al., 1996b; Ruttenberg, 2003), increasing evidence points towards a significant role for the Fe(II)-phosphate mineral vivianite (Dijkstra et al., 2016; Egger et al., 2015a; Li et al., 2015; Slomp et al., 2013).

Vivianite can form in sediments where the supply of sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) is low relative to the input of Fe(oxyhydr)oxides (henceforth termed Fe oxides). As a consequence, not all Fe oxides are converted to Fe sulfides (FeS and FeS₂; Ruttenberg, 2003) and Fe oxides can continue to dissolve at depth through Fe reduction coupled to organic matter degradation or anaerobic oxidation of methane (Beal et al., 2009; Egger et al., 2015b; Ettwig et al., 2016). In coastal sediments, conditions for the formation of vivianite are especially favorable below the sulfate methane transition zone (SMTZ; Egger et al., 2015a; Hsu et al., 2014; März et al., 2008; Roden and Edmonds, 1997). Here, both porewater Fe²⁺ and ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$ can accumulate in the absence of sulfide and precipitate as vivianite (Dijkstra et al., 2016; Egger et al., 2015a; Ruttenberg, 2003). Therefore, sediments with high rates of sedimentation and inputs of organic matter and Fe oxides and with an SMTZ that is close to the sediment–water interface are thought to be particularly conducive to vivianite formation (Egger et al., 2015a; Slomp et al., 2013).

Recent studies of diagenesis in sediments of the northern Baltic Sea indicate that variations in the input of organic matter and Fe oxides may contribute to temporal variations in vivianite formation in coastal sediments (Egger et al., 2015a; Rooze et al., 2016). While Egger et al. (2015a) attributed this to eutrophication, the large magnitude of the change in input of organic matter and of Fe required to describe these data with a reactive transport model (factors of 15 and 2, respectively; Rooze et al., 2016) suggests that additional processes are active. We hypothesize that variations in river input of both organic matter and Fe, linked to variations in river discharge, play a role. This hypothesis is based on observed changes in Fe dynamics in soils and rivers in Finland (Sarkkola et al., 2013). More specifically, dry periods may lower groundwater levels exposing Fe minerals in previously reducing soil layers to oxidation (Laudon et al., 2011). After dry periods, the freshly formed Fe oxides can be reduced again, after which Fe can be mobilized as Fe²⁺ and transported to rivers and the coastal zone (Sarkkola et al., 2013). While some of the Fe²⁺ will be oxidized to form amorphous and crystalline Fe oxides (Pokrovsky and Schott, 2002; Schwertmann and Taylor, 1989) and will be transported downstream in particulate form, some of the Fe²⁺ can also form complexes with dissolved organic compounds and remain in solution (Kuma et al., 1996; Rue and Bruland, 1995; Wu and Luther, 1995). Both the particulate and dissolved Fe and the associated organic matter may be transported to the coastal zone through river flow where they may impact P burial.

Future climate change is expected to alter precipitation patterns in northern Europe and change both the magnitude and seasonality of river discharge (Harley et al., 2006). In the northern Baltic Sea, an increase in the input of freshwater is expected over the coming century (Meier et al., 2006). Such a change may further promote vivianite formation by reducing the availability of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (Hartzell et al., 2017; Jordan et al., 2008). Environmental change may also alter riverine input of Fe. In the northern Baltic Sea, the input of dissolved Fe from rivers has increased over the past decades (Björnerås et al., 2017; Kritzberg and Ekström, 2012; Sarkkola et al., 2013). The associated increased input of Fe oxides to coastal sediments may also enhance vivianite formation.

Figure 1 (a) The Öre Estuary on the Swedish coast in the Bothnian Sea; (b) Locations of the five sampling sites. Site NB8, in the deepest part of the estuary, is indicated with a blue star. The mouth of the Öre River is indicated with a red star. Figure drawn using Ocean Data View (Schlitzer, 2015).

Figure 2 Öre River discharge (km³ yr⁻¹) from 1970 to 2012 (Swedish Meteorological and Hydrological Institute, http://SMHI.se, last access: 9 May 2016). The red dots indicate two years with an unusually low river runoff, 1976 and 1996, respectively.

In this study, we investigate the burial of P in the Öre Estuary in the northern Baltic Sea, with a specific focus on the factors contributing to vivianite formation. Porewater and sediment geochemical depth profiles for five sites characterized by varying rates of sedimentation are presented. We apply a reactive transport model to the data for one of the sites and investigate how vivianite formation and P burial in the sediment responds to changes in (1) bottom water salinity, (2) sedimentation rate, (3) the input of Fe oxides and (4) the input of organic matter. We also discuss the expected changes in the burial of P in coastal sediments upon ongoing climate change. Our model results suggest that vivianite will become a more important sink for P in the coastal zone of the Baltic Sea in future because of the expected decline in salinity and the continued high input of Fe from land.

2 Methods

2.1 Study area and sampling

The Öre Estuary is located along the Swedish coast in the Bothnian Sea, in the northern Baltic Sea (Fig. 1). The estuary is oligotrophic and has a surface area of approximately 70 km², a mean depth of 10 m and a bottom water salinity of ca. 5 (http://SMHI.se, last access: 8 September 2016). The estuary is fed by the Öre River, which has a strongly varying rate of annual discharge (Fig. 2). While surface water salinity varies with time and distance from the river, there is no salinity gradient in the bottom water of this estuary. The spatial extent of the particle plume from the river varies with time but is known to reach all of our sample sites (Malmgren and Brydsten, 1992). In this region, the spring flood (April, May) is the major annual hydrological event. The spring flood results in a brief period of enhanced water flow in streams and rivers (Björkvald et al., 2008; Hölemann et al., 2005) and export of Fe and terrestrial organic matter to the coastal zone (Algesten et al., 2006; Rember and Trefry, 2004). During high flow, total and dissolved Fe and Mn in streams can increase by a factor of up to 10, with concentrations being highly variable between years (Björkvald et al., 2008).

This study focuses on sediments from five sites in the Öre Estuary (Fig. 1b), with water depths varying from 10 to 33 m (Table 1). Sediments at all sites are fine-grained and rich in organic matter (Hellemann et al., 2017; Table 1). The major macrofauna species in the area are the mollusk Limecola balthica (formerly named Macoma balthica), the amphipod Monoporeia affenis, the spionid Marenzelleria spp. and the isopod Saduria entomon (Fig. S1). Marenzelleria first appeared in 1995, but only became abundant after 2002/2003, with densities at our study sites ranging from 187 to 780 ind m⁻² (Kauppi et al., 2015). The density of Marenzelleria spp. in the sediment is water depth dependent with maximum densities at water depths between 30 and 50 m and lower densities both above and below (Kauppi et al., 2015).

Sediment was collected during two field campaigns with R/V Lotty in April and August 2015 using a Gemini gravity corer (8 cm inner diameter). In April, three cores were taken per site. The first core was used for porewater and solid-phase analyses, the second core was used for methane (CH₄) sampling and the third core for sediment dating using ²¹⁰Pb. In August, two cores were collected per site for porewater analyses and for CH₄. At site NB8, a third and a fourth core were collected in August for mineral analyses and additional ²¹⁰Pb measurements.

Table 1 Water depth, temperature, coordinates and organic carbon at five sites sampled in the Öre Estuary in April and August 2015. Depth unit mbss is meters below sea surface. Bottom water temperatures were derived from CTD data. Total organic carbon (C_org) is given as the average for the top 2 cm (n=2), with standard deviation between parentheses. Sedimentation rates for sites NB1, N6, N10 and N7 were determined on cores collected in April 2015. The sedimentation rate at site NB8 was determined on a core collected in August 2015.

Methane was sampled directly after core recovery using a core liner with pre-drilled holes with a 2.5 cm depth spacing. Samples of 10 mL were taken with cutoff syringes from each hole and immediately transferred to a 65 mL glass bottle filled with saturated salt solution. The bottles were stoppered, capped and stored upside down until analysis.

After core retrieval, two bottom water samples were taken and the cores were subsequently sliced into intervals of 1 to 4 cm under a nitrogen atmosphere at bottom water temperature (Table 1). Each sediment sample was divided over a pre-weighed glass vial for solid-phase analysis and determination of water content, and a 50 mL centrifuge tube. The glass vials were stored in nitrogen flushed gas-tight aluminum bags at −20 ^∘C until analysis. The 50 mL centrifuge tubes were centrifuged at 3500 rpm for 20 min to extract porewater.

2.2 Bottom water and porewater analyses

Bottom and porewater samples were filtered through 0.45 µm pore size filters and subsampled under a nitrogen atmosphere. Subsamples were taken for analysis of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, hydrogen sulfide (where H₂S represents the sum of H₂S, HS⁻ and S²⁻), Fe, manganese (Mn), ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$, ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) and dissolved inorganic carbon (DIC). Total Fe and Mn are assumed to represent Fe²⁺ and Mn²⁺; although in the latter case some Mn³⁺ may also be included (Madison et al., 2011). Subsamples for ${\mathrm{NH}}_{\mathrm{4}}^{+}$ were stored frozen at −20 ^∘C. All other subsamples were stored at 4 ^∘C until analysis. Samples for ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ were analyzed with ion chromatography (detection limit of <75 µmol L⁻¹; average analytical uncertainty based on duplicates and triplicates <2 %). For H₂S, 0.5 mL of porewater was immediately transferred into a 4 mL glass vial containing 2 mL of 2 % zinc acetate solution to trap the H₂S as ZnS. Sulfide was determined spectrophotometrically by complexing of the ZnS precipitate in an acidified solution of phenylenediamine and ferric chloride (Cline, 1969). The sulfide standard was verified by titration with thiosulfate. Subsamples taken for porewater Fe²⁺, Mn²⁺ and ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$ were acidified with 10 µL 30 % Suprapur^® HCl per milliliter of sample and were analyzed by inductively coupled plasma optical emission spectrometry (ICP-OES; Spectro Arcos).

Concentrations of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ were determined with the modified indophenol blue method (Solorzano, 1969). Samples for DIC analyses were collected in 5 mL vials without a headspace and poisoned with 10 µL of saturated HgCl₂. Analysis for DIC was performed using an AS-C3 analyzer (Apollo SciTech), consisting of an acidification and purging unit in combination with a LI-COR 7000 CO₂∕H₂O gas analyzer.

Samples for CH₄ were prepared for measurement by injecting a 10 mL nitrogen headspace into the bottle. Subsequently, the CH₄ concentrations in the headspace were determined by injection of a subsample (50–200 µL) into a Thermo Finnigan Trace GC gas chromatograph (flame ionization detector), after which calculated CH₄ concentrations were corrected for sediment porosity. The average analytical uncertainty based on duplicates and triplicates was <3 % for DIC and <5 % for CH₄.

2.3 Solid-phase analyses

Sediments were freeze-dried and the porosity was determined from the weight loss upon freeze-drying. Freeze-dried sediments were ground in an agate mortar under a nitrogen atmosphere and were split into subsamples that were stored either under oxic or anoxic conditions. For site NB8, the speciation of solid-phase Fe, S and P was determined on the anoxic subsamples to avoid oxidation artifacts (Kraal and Slomp, 2014; Kraal et al., 2009). All other analyses were performed on the oxic subsamples.

Samples of ca. 300 mg of sediment was decalcified with two wash steps of 1 M HCl (Van Santvoort et al., 2002) and subsequently dried, powdered and analyzed for carbon using an elemental analyzer (Fison Instruments model NA 1500 NCS). Organic carbon (C_org) contents were determined after correction for the weight loss during decalcification. The average analytical uncertainty based on duplicates and triplicates was <5 %. A second subsample of ca. 125 mg was digested in 2.5 mL of HClO₄ and HNO₃ (ratio 3 : 2) and 2.5 mL 40 % HF in a Teflon vessel at 90 ^∘C overnight. The acid was evaporated at 140 ^∘C until a gel was formed, which was subsequently dissolved in 25 mL of 4.5 % HNO₃ at 90 ^∘C overnight. Total elemental concentrations in the HNO₃ solution of Al, Fe, Mn, P and S were determined by ICP-OES, from which total concentrations in the sediment samples were calculated. The average analytical uncertainty based on duplicates and triplicates was <3 % for all reported elements.

Sediment samples of ca. 50 mg were subjected to the four-step Fe speciation procedure of Poulton and Canfield (2005) that targets the following Fe phases: (i) carbonate-associated Fe (including siderite and ankerite), extracted for 24 h with 1 M Na acetate (pH brought to 4.5 with acetic acid); (ii) easily reducible Fe oxides and vivianite (including ferrihydrite and lepidocrocite), extracted for 24 h with 1 M hydroxylamine-HCl in 25 % v∕v acetic acid; (iii) reducible (crystalline) Fe oxides and vivianite (including goethite, hematite and akagenéite), extracted for 2 h with Na-dithionite (pH 4.8); (iv) recalcitrant Fe oxides (mainly magnetite), extracted for 2 h with 0.2 M ammonium oxalate/0.17 M oxalic acid. Samples were measured colorimetrically using the 1,10-phenanthroline method (APHA, 2005). For simplicity, fractions (ii) and (iii) were summed and henceforth referred to as total Fe oxides. Average analytical uncertainty, based on duplicates, was <5 % for all fractions.

Sediment samples of ca. 500 mg were subjected to the three-step S speciation procedure of (Burton et al., 2006, 2008) that targets the following S phases: (i) acid-volatile sulfur (AVS, i.e., mostly FeS) by addition of 2 mL of ascorbic acid and 10 mL of 6 M HCl to the sediment and trapping of the released H₂S into a tube filled with 7 mL Zn-acetate; (ii) elemental sulfur (S⁰) by extracting overnight with 25 mL methanol; (iii) chromium-reducible sulfur (CRS, i.e., mostly FeS₂) by addition of acidic chromium chloride solution and trapping of the released H₂S into a tube filled with 7 mL Zn-acetate. Samples for AVS and CRS were analyzed by iodometric titration of the alkaline Zn-acetate trap. Elemental sulfur was measured colorimetrically according to Bartlett and Skoog (1954). Average analytical uncertainty, based on duplicates, was <3 % for all fractions.

Sediment samples of ca. 100 mg were subjected to the SEDEX method described by Ruttenberg (1992) as modified by Slomp et al. (1996a), but including the exchangeable P step. Five P phases were distinguished: (i) exchangeable-P, extracted for 30 min with 1 M MgCl₂ (pH 8); (ii) Fe-bound P fraction (including Fe-oxide-bound P and vivianite; Dijkstra et al., 2014; Nembrini et al., 1983) extracted for 8 h with citrate-dithionite-bicarbonate (CDB; pH 7.5) followed by extraction for 30 min with 1 M MgCl₂ (pH 8); (iii) authigenic Ca-P (including carbonate fluorapatite, biogenic hydroxyapatite and carbonate-bound P), extracted for 6 h with Na acetate buffer (pH 4) followed by extraction for 30 min with 1 M MgCl₂ (pH 8); (iv) detrital Ca-P, extracted for 24 h with 1 M HCl; (v) organic P, extracted for 24 h with 1 M HCl after ashing at 550 ^∘C for 2 h. Steps (i)–(iv) were performed under an argon atmosphere to avoid oxidation artifacts (Kraal and Slomp, 2014; Kraal et al., 2009). CDB solutions were analyzed for P with ICP-OES. All other solutions were measured colorimetrically according to Strickland and Parsons (1972). Average analytical uncertainty, based on duplicates, was <5 % for all fractions.

2.4 Sedimentation rate and P burial

Sediment accumulation rates at all 5 sites were determined from depth profiles of ²¹⁰Pb for April (NB1, N6, N10 and N7) or August (NB8). ²¹⁰Pb was measured on freeze dried sediment by direct gamma counting at 46.5 keV using a high-purity germanium detector (Ortec GEM-FX8530P4-RB). Self-absorption was measured directly and the detector efficiency was determined by counting a National Institute of Standards and Technology sediment standard. Excess ²¹⁰Pb was calculated as the difference between the measured total ²¹⁰Pb and the estimate of the supported ²¹⁰Pb activity as given by ²¹⁴Pb (²¹⁰Pb_exc = ²¹⁰Pb_total − ²¹⁴Pb). Sediment accumulation rates at each site were estimated by fitting a reactive transport model (Soetaert and Herman, 2008) to the ²¹⁰Pb depth profiles assuming depth dependent changes in porosity and bioturbation (Fig. S2).

Total P burial ($\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) for all sites was calculated as follows:

$$\begin{array}{}\text{(1)}& \text{P burial}={\mathrm{P}}_{\mathrm{total}}\times \text{sed. rate}\times \mathit{\rho }\times (\mathrm{1}-\mathit{\varphi })\times {\mathrm{10}}^{\mathrm{4}},\end{array}$$

where P_total is the averaged concentration of total P (mol g⁻¹) over the deepest 10 cm of the sediment sampled, ϕ is the porosity in the same interval (vol vol⁻¹), sed. rate is the sedimentation rate (cm yr⁻¹) and ρ is the density of dry sediment, 2.65 g cm⁻³ (Burdige, 2006).

2.5 Scanning electron microscopy (SEM)

Sediment samples from site NB8 were analyzed by SEM to determine whether large vivianite crystals (≥38 µm size) were present. Wet sediments from five depths (sediment depth intervals: 3–4, 14–16, 34–36, 49–52 and 55–58 cm) were sieved through a 38 µm mesh size sieve under an argon atmosphere with deoxygenated ultraclean water. The sieved material was washed seven times (5 min) with deoxygenated ultraclean water in a sonic bath. After washing, the samples were dried in an argon-filled glove box at ambient temperature. A subsample of the sieved and dried material was mounted on an aluminum sample holder using double sided carbon tape and subsequently coated with 0.8 mm platinum. Samples were analyzed using SEM energy dispersive X-ray spectroscopy (EDS; JCM 6000PLUS NeoScope Benchtop SEM) with 15 kV accelerating voltage using a Si/Li detector, in scanning electron mode. Measurements with EDS were performed in the 0–20 keV energy range for elemental quantification (probe current: 1 nA; acquisition time: 50 s, live time). SEM-EDS software was used to estimate the relative abundances in mol % for the major elements (oxygen, O; sodium, Na; magnesium, Mg; Al; silica, Si; P; calcium, Ca; Mn; Fe). Samples were detected and photographed in secondary electron imaging (SEI) mode. Measurements were performed with a 1 µm beam in backscattered electron imaging (BEI).

Figure 3 Transient fluxes at the sediment–water interface from 1965 to 2015 as used in the reactive transport model; (a) sedimentation rate; (b) organic matter input; (c) Fe oxide input; (d) total Mn input (Mn oxide and Mn carbonate). Partitioning of different phases of organic matter, Fe oxides and manganese are shown in Fig. S4a–c.

3 Reactive transport modeling

3.1 General model description

To investigate the mechanisms that control the formation of vivianite in the Öre Estuary a reactive transport model was applied to key porewater and solid-phase depth profiles for site NB8. Site NB8 was characterized by a SMTZ close to the sediment–water interface and a relatively high sedimentation rate (Table 1). The geochemistry at site NB8 is strikingly similar to that at another site with a high sedimentation rate in the Bothnian Sea (Egger et al., 2015a, b). The model describes the mass balance of 11 dissolved and 18 particulate species (Table S1) and is a modified version of that of Rooze et al. (2016), extended here to include the sedimentary Mn cycle. This is a standard multicomponent reactive transport model based on principles outlined by, for example, Wang and Van Cappellen (1996). The large number of model components, the parameters specific for the field site, and the reaction kinetics and parameters from the literature greatly limit the degrees of freedom while fitting the model to the data (Van Cappellen and Wang, 1996). Three forms of Mn are distinguished, namely Mn oxides, Mn carbonate and Mn²⁺. The transformations of Mn that are included are Mn carbonate and Mn oxide formation and Mn oxide reduction coupled to either H₂S, Fe²⁺ or CH₄ oxidation (Table S2; Reactions R25–R32).

The model domain consists of a one-dimensional grid of 900 evenly distributed cells that captures the interval from the sediment–water interface to a depth of 90 cm. All chemical species are subject to biogeochemical reactions (Table S2). Solids and solutes are transported by sediment accumulation and bioturbation. Solutes are additionally transported by molecular diffusion and bioirrigation (Boudreau, 1997; Soetaert et al., 1996; Wang and Van Cappellen, 1996). Bioirrigation is modeled as a nonlocal exchange process (Boudreau, 1984; Emerson et al., 1984).

3.2 Model parameterization

The model was parameterized using data for site NB8, information from the literature, and by visually fitting modeled porewater and solid-phase depth profiles to the measured data (model constrained; Table S3). Environmental parameters such as porosity, temperature, salinity and bottom water solute concentrations were measured (Table S4; Fig. S3a). The density of dry sediment and the C : N ratio for organic matter were taken from literature (Burdige, 2006; Redfield, 1958). Bioturbation and bioirrigation coefficients were used as fitting parameters (Fig. S3b and c), taking into account when the most important bioirrigator in these sediments, Marenzelleria spp., established (Fig. S1) and the typical range in these coefficients at similar macrofaunal densities (Renz and Forster, 2013). Fluxes of solids at the sediment–water interface were model constrained. A detailed description of the reactive transport model is given in the Supplement.

3.3 Iron, manganese and phosphorus dynamics

In the model, Fe oxides and Mn oxides are assumed to consist of fractions with different crystallinities, which affect their reactivity towards organic matter and H₂S (Table S1). In both cases, a highly reactive fraction (α) and a less reactive fraction (β) are assumed, whereas for Fe oxides also a refractory fraction (γ) is included. Only the α fraction of the Fe oxides and Mn oxides is susceptible to reductive dissolution linked to degradation of organic matter. This allows the β fraction to be buried below the zone of organic matter degradation (Rooze et al., 2016).

Reactive P is deposited at the sediment–water interface in the form of organic P and Fe-oxide-bound P. For organic matter three fractions are assumed: a highly reactive (α), less reactive (β) and refractory (γ) fraction, following the multi-G approach (Jørgensen, 1978; Middelburg, 1989; Westrich and Berner, 1984). The C : P ratio is 300 : 1 for all three fractions of organic matter and is model constrained. This relatively high and constant value for the C : P ratio is based on the refractory nature of the organic matter in the area (Algesten et al., 2006; Leipe et al., 2011; Stockenberg and Johnstone, 1997). The three fractions of Fe oxides are assumed to have different P : Fe ratios, with highly reactive Fe oxides having a higher P : Fe ratio then less reactive crystalline Fe oxides assuming the former can bind more P (Table S3; Gunnars et al., 2002). In the bioirrigated zone, porewater ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$ is assumed to bind to the Fe oxide^β fraction present (Fig. S3c and d) to a maximum P : Fe ratio of 0.28 (mol mol⁻¹; Table S3). This allows for enhanced formation of Fe-oxide-bound P due to bioirrigation (Norkko et al., 2012). Nonreactive P is deposited at the sediment–water interface as detrital P, authigenic Ca–P, P bound to nonreactive Fe oxides (γ fraction) and P bound to nonreactive organic matter (γ fraction).

In the sediment, the rate of vivianite formation (R_viv) is modeled by means of Michaelis–Menten kinetics for dissolved Fe²⁺ and ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$, which implies that the rate of vivianite formation depends on both porewater species and that there is a maximum rate of formation (Reed et al., 2016):

$$\begin{array}{}\text{(2)}& R_{\text{viv}}=V_{max}\left({\displaystyle \frac{\left[{\mathrm{Fe}}^{\mathrm{2}+}\right]}{\left[{\mathrm{Fe}}^{\mathrm{2}+}\right]+K_{{\mathrm{Fe}}^{\mathrm{2}+}}}}\right)\left({\displaystyle \frac{\left[{\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}{\left[{\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}\right]+K_{{\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}}}}\right),\end{array}$$

where V_max is the maximum rate ($\mathrm{mol}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$), $K_{{\mathrm{Fe}}^{\mathrm{2}+}}$ and $K_{{\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}}$ are the half saturation constants, [Fe²⁺] and [${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$] are the porewater concentrations of Fe²⁺ and ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$.

3.4 Transient modeling scenario

The model was run to steady state for 200 years. Subsequently, temporal changes in the sedimentation rate and the input of organic matter, Fe oxides, Mn oxides and Mn carbonates were implemented to fit key porewater and solid-phase depth profiles for site NB8 (Figs. 3 and S4a–c). The modeled sedimentation rate in the first 200 years was set to 1 cm yr⁻¹ (0.34 $\mathrm{g}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), based on ²¹⁰Pb dating (Fig. S2; Table 1). From 2003 onwards sedimentation was assumed to decrease to 0.6 cm yr⁻¹ (0.21 $\mathrm{g}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$; Fig. 3a). The modeled organic matter deposition in the first 200 years was 7.65 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. From 1987 to 2003, the deposition increased, whereas in the last 12 years the deposition decreased together with the sedimentation rate (Fig. 3b). Our modeled organic matter fluxes are at the high end of the range estimated for the region (1.1–8.2 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$; Algesten et al., 2006). The qualitative trend in organic matter loading assumed in the model is in accordance with results from previous studies on eutrophication in the Bothnian Sea (Fleming-Lehtinen and Laamanen, 2012; Fleming-Lehtinen et al., 2008; Rooze et al., 2016).

The assumption is made that during strongly enhanced riverine Fe oxide input the P : Fe ratio of Fe oxides can be lower because on land the availability of ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$ is not high enough to maintain a high P : Fe ratio. This assumption is based on the dominance of natural forest and wetlands in the region (Björkvald et al., 2008). We note that the model results were very sensitive to the P : Fe ratio and that we could not model the key profiles when assuming a higher P : Fe ratio. To be able to model a transient P : Fe ratio in the Fe oxide^β fraction a fifth Fe oxide fraction, Fe oxide^β(pulse), was added to the model. This fraction has the same reactivity as Fe oxide^β towards organic matter and H₂S but has a P : Fe ratio of 0. By varying the ratio between Fe oxide^β and Fe oxide^β(pulse) the P : Fe ratio of the Fe oxide^β can be adjusted over time.

Figure 4 Porewater depth profiles of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, CH₄, H₂S, Fe²⁺, Mn²⁺, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, DIC and ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$for sites NB1, N6, N10, N7 and NB8 in April (blue) and August (red) 2015. Porewater profiles of station NB8 include model fits to the April 2015 data for the baseline scenario (black lines). The SMTZ at site NB8 is indicated by the green shaded area.

Figure 5 Solid-phase depth profiles of organic carbon, total S, total P, total Mn, total Fe and Fe/Al for sites NB1, N6, N10, N7 and NB8 in April 2015. Solid-phase profiles of station NB8 include model fits to all profiles (black lines) except total Fe and Fe/Al for the baseline scenario.

In the first 200 years, the Fe oxide input was set to 0.87 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. The Fe oxide input was then assumed to increase from 1988 onwards until 1997 when a maximum in incoming Fe oxides was assumed, followed by a decrease in Fe oxide input (Fig. 3c). The strong pulse-type increase in Fe oxide input at this maximum was assumed to consist of mainly the β(pulse) fraction of Fe oxides (Fig. S4b). Despite this lower P : Fe ratio, the total flux of P bound to Fe oxides was highest during this period of enhanced Fe oxide input (Fig. S4d and e). The transient scenario of Fe oxide loading compares well to qualitative trends used in a previous model study for the same region (Rooze et al., 2016). The initial Mn oxide flux was set to 0.085 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Fig. 3d). In the last 17 years of the scenario, the Mn oxide flux was set to a much lower value. The input flux of Mn carbonate was always low compared to that of Mn oxides (4–29 % of the total incoming Mn flux; Fig. S4c). Processes coupled to the presence of bioirrigating macrofauna (i.e., bioirrigation and the binding of P onto Fe oxides in the bioirrigated zone (Reaction R32 in Table S2), were only implemented in the last 12 years of the run, from 2003 onwards.

A model sensitivity analysis was performed to investigate the impact of changes in bottom water salinity, the rate of sedimentation and input of organic matter and Fe oxides on P burial rates and forms. During these sensitivity analyses the transient baseline scenario was used (Fig. 3) and one initial factor was varied per run. Salinity was varied over a range of 0 to 20. The sedimentation rate was varied from 0.25 to 2 cm yr⁻¹. The input of organic matter was varied between −70 % and +25 %; and the input of Fe oxides was varied between −25 % and +25 % from the baseline scenario. Subsequently, a run was performed in which the input of organic matter, Fe oxides and the sedimentation rate were all changed by the same factor as the sedimentation rate to account for the role of rivers as the main source of material in the region (Algesten et al., 2006; Björkvald et al., 2008).

4 Results

4.1 Porewater profiles

At all sites, depth profiles of porewater constituents showed relatively little difference between April and August 2015 (Fig. 4). Porewater ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ decreased with depth, with a distinct SMTZ only being present in the sampled depth interval at sites N7 and NB8 (Fig. 4). Concentrations of CH₄ were low at all other sites. Appreciable H₂S was only present at site NB8 in the SMTZ. Dissolved Fe²⁺ and Mn²⁺ profiles generally showed a maximum near the sediment–water interface and increased concentrations below ca. 5 to 25 cm depth in the sediment. Porewater profiles of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and DIC increased with depth at all sites. Porewater ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations were low close to the sediment–water interface and increased with depth at all sites. At sites N7 and NB8, ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations remained constant or decreased again below the SMTZ. Concentrations of Fe²⁺, Mn²⁺, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and DIC were highest at depth at site NB8 when compared to the other sites.

Figure 6 Solid-phase depth profiles of Fe oxides, Fe-bound P (Fe–P), AVS, organic P (Org. P), exchangeable P (Exch. P), detrital Ca-P (Detr. P) and authigenic Ca-P (Ca–P) for site NB8 in April 2015. Solid-phase profiles include model fits of the baseline scenario for all profiles except exchangeable P (black lines).

4.2 Solid-phase profiles and sedimentation rate

Organic carbon contents in the surface sediment ranged between 2 and 4 wt % (Table 1) and decreased with depth at all sites (Fig. 5). Total S concentrations increased with depth in the upper 10 to 20 cm of the sediment and subsequently remained constant or varied with depth. Total P, total Mn, total Fe and Fe/Al were generally highest close to the sediment–water interface and decreased with depth at all sites except site NB8. Here, several subsurface enrichments in P and Mn were observed directly below maxima in total S, Fe and Fe/Al at depths of 21, 42 and 60 cm, respectively.

The Fe and P speciation for site NB8 reveals that the upper 10 cm of the sediment was strongly enriched in Fe oxides and Fe-bound P. Minima in both Fe oxides and Fe-bound P were found between depths of 10 and 20 cm followed by either constant (Fe oxides) or varying (Fe-bound P) concentrations at greater depth (Figs. 6 and S5). Maxima in AVS coincided with maxima in total S, total Fe and Fe/Al (Figs. 5 and 6), and minima in total Mn and Fe-bound P. Concentrations of exchangeable P were highest close to the sediment–water interface but account for only ∼2 % of total P. Concentrations of detrital P and authigenic Ca-P showed little change with depth.

Sedimentation rates at our study sites varied between 0.225 and 1 cm yr⁻¹ (Table 1). Rates of P burial ranged from 0.026 to 0.145 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Table 1). The rates of sedimentation and P burial were highest at site NB8.

4.3 SEM-EDS

Examination of sieved sediment fractions for site NB8 with a light microscope revealed the presence of transparent to light blue crystals in sediment intervals below the SMTZ (55–58 cm). Scanning electron micrographs of the sieved fractions revealed crystals with a similar shape to that of vivianite crystals found in Lake Ørn, Denmark, (O'Connell et al., 2015) and the deepest part of the Bothnian Sea (site US5B; Egger et al., 2015a; Fig. 7a) at all sieved intervals below the SMTZ (34–36, 49–52 and 55–58 cm). Quantitative analysis of the crystals with electron microprobe EDS revealed an average Fe : P ratio of 3.3 mol mol⁻¹ and the presence of significant amounts of Mn, Mg, Si and Al (Fig. 7b; Table S5). In samples analyzed above and within the SMTZ such crystals were not observed.

Figure 7 (a) Scanning electron micrograph: example of a vivianite crystal from the 34–36 cm depth interval at site NB8 sampled in August 2015. (b) Electron microprobe EDS spectrum of a spot measurement on the surface of the vivianite crystal.

4.4 Model results

Modeled porewater and solid-phase depth profiles for site NB8 capture the main trends in the measured profiles (Figs. 4, 5 and 6). We note that the model slightly underestimates the amount of Fe-bound P in the upper 15 cm of the sediment and overestimates organic P. The observed variations in total S, total Mn and Fe-P below 30 cm depth in the sediment were not targeted with the model.

Depth integrated rates of P cycling in 2015 (Fig. 8) show that the release of ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$ from organic P accounts for only 2.3 % of the ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$ release to the porewater at site NB8. The main source of ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$ is the release from Fe oxides, with H₂S-driven reductive dissolution accounting for 61 % of the ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$ release and dissolution coupled to organic matter degradation accounting for 28 % (Fig. 8). In the model, most of the Fe oxide reduction involves the less reactive beta fraction and takes place in the upper 20 cm of the sediment (Fig. S6; Reactions R14 and R15).

Approximately 30 % of the ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$ removal from the porewater takes place through release to the overlying water by bioirrigation and diffusion, with the former being the most important process. Binding of ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$ to Fe oxides accounts for 60 % of the removal while 10 % is precipitated as vivianite (Fig. 8). However, Fe oxides that are formed in situ are not a permanent sink for ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$ as the majority is dissolved through reductive dissolution. Burial of P as vivianite accounts for 41 % of the total P burial at site NB8 (Fig. 8). At depth, P bound to vivianite becomes more important compared to P bound to Fe oxides (Fig. S7).

Figure 8 Mass balance of P as calculated with the transient reactive transport model for 2015 (Fig. 3). Results are in non-steady state. Fluxes are in $\mathrm{mol}\mathrm{P}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$.

In the model, total P burial rates increased from 0.12 to 0.34 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ upon a decrease in salinity from 20 to 0 (Fig. 9). Changes in salinity did not affect the burial rates of Ca-P, organic P and P associated with Fe oxides. In the latter case, this was because only the unreactive fraction of Fe oxides remained. The percentage of P buried as vivianite increased strongly with decreasing salinity.

In the model scenario in which changes in the rates of sedimentation and organic matter and Fe oxide input were coupled, the total P burial rate increased from 0.025 to 0.45 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ when assuming sedimentation rates between 0.25 and 2 cm yr⁻¹ (Fig. 9). At a sedimentation rate of 0.5 cm yr⁻¹, little P was buried in the form of vivianite. At higher sedimentation rates, formation of vivianite in the sediment increased and the total burial of P was enhanced. Model runs in which changes in the same factors were evaluated separately revealed that increasing the rate of sedimentation or Fe oxide input both enhance vivianite formation and burial of P (Fig. 10a and b). The response to an increase in the rate of sedimentation is nonlinear, with the greatest increase in P burial between 0.5 and 0.85 cm yr⁻¹. For organic matter, the model runs reveal an optimum input rate beyond which P and vivianite burial rates decrease (Fig. 10c).

Figure 9 Sensitivity of (a) total phosphorus burial rate and (b) the burial forms of P at site NB8 to changes in salinity, as well as in sedimentation rate and inputs of organic matter and Fe oxides as calculated in the model. Four burial forms of P are distinguished: authigenic and detrital apatite (Ca-P), organic P (Org-P), P bound to Fe oxides (FeOx-P) and vivianite. The dashed line indicates the total P burial as a function of sedimentation rate in the coastal zone of the Baltic Sea as suggested in Asmala et al. (2017). Red stars are the burial rates of P as measured at our five sites (Table 1). Sedimentation rates are initial sedimentation rates at t=0.

5 Discussion

5.1 Phosphorus burial in the Öre Estuary

Sediments in the coastal zone of the Bothnian Sea have been suggested to act as an efficient sink for P from rivers and the open sea based on budget calculations and an assumed average burial rate of P of 0.007 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Asmala et al., 2017). We find substantially higher rates of P burial in the Öre Estuary ranging from 0.026 to 0.145 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Table 1). If our data are representative for the wider area, sediments in the coastal zone of the northern Baltic Sea may be an even more efficient sink for P than previously thought.

In their study, Asmala et al. (2017) identified a linear relationship between P burial and rates of sedimentation for a range of coastal systems around the Baltic Sea, but excluding the northern areas because of a lack of data (Fig. 9). Strikingly, our study sites in the Öre Estuary are characterized by a higher burial of P than predicted by this relationship. This particularly holds for site NB8. In the following, we will assess whether vivianite formation can play a role in explaining this enhanced P retention. As discussed in detail by Egger et al. (2015a), vivianite formation is generally most pronounced in sediments below the SMTZ, where both Fe²⁺ and ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$ accumulate in the porewater due to dissolution of Fe oxides and release of associated P. While we find both solutes in the porewater at depth at all sites, the highest concentrations are observed at sites N7 and NB8 where a distinct SMTZ was present between ca. 10 and 20 cm depth in the sediment. The position of this SMTZ shows little seasonality, likely because of the refractory nature of the sediment organic matter, which is mostly of terrestrial origin, and the relatively limited seasonal change in bottom water temperature. This implies that any vivianite formed below the SMTZ will likely be preserved and buried, because there will be no further exposure to sulfide.

Detailed analyses of the sediments at site NB8 with SEM-EDS confirm the presence of vivianite crystals below the SMTZ and a lack thereof in and above the SMTZ (Fig. 7). This trend with depth points towards an authigenic origin. The average Fe : P ratio of the crystals of 3.3 mol mol⁻¹ (Table S5) is higher than the 2 : 1 stoichiometric ratio of vivianite. Similar high ratios were also observed for vivianite crystals in sediments of the Landsort Deep in the Baltic Proper, where they were explained by surface enrichments of Fe (Dijkstra et al., 2016). The vivianite was also enriched in Mn and Mg, which are both elements that are known to be included in the structure of the mineral (Dijkstra et al., 2016, 2018; Egger et al., 2015a). The presence of Al and Si (Table S5) likely reflects the presence of clay particles that were not removed prior to the SEM-EDS analysis (Dijkstra et al., 2018; Egger et al., 2015a).

The Fe-P below the SMTZ (Fig. 6) is likely a mixture of both Fe-oxide-bound P and vivianite, as shown previously for sediments in the deepest part of the Bothnian Sea (Egger et al., 2015a). Due to dissolution of Fe oxides and vivianite formation in the sediment, the fraction of P bound to vivianite increases with depth (Fig. S7). The results of the reactive transport model suggest that Fe-oxide-bound P and vivianite account for ca. 15 % and 40 % of the total burial of P, respectively. Authigenic Ca-P and organic P concentrations are relatively low and each account for ca. 15 % of total P burial (Fig. S5). The remainder of the P is buried as nonreactive P.

Figure 10 Sensitivity analysis of P burial as a function of (a) changing sedimentation rate between 0.25 and 2 cm yr⁻¹ and burial forms of P; (b) changing input of Fe oxides ranging between −25 % and 25 % and burial forms of P; (c) organic matter (OM) input ranging between −70 % and 25 % and burial forms of P.

5.2 Vivianite formation in coastal sediments in the northern Baltic Sea

Vivianite formation in sediments strongly depends on the balance between the formation of H₂S and the input of Fe oxides (Ruttenberg, 2003). When there is an excess of Fe oxides over H₂S, not all Fe oxides will be converted to Fe sulfides, and more Fe²⁺ will be available to precipitate with ${\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}$ as vivianite (Gächter and Müller, 2003; Hartzell et al., 2017; Rozan et al., 2002). The results of the sensitivity analysis (Figs. 9 and 10) show that the bottom water salinity and input of Fe oxides play a critical role in controlling this balance, with vivianite formation and total P burial increasing strongly at lower salinities and with high Fe oxide inputs. A higher sedimentation rate also enhances the formation of vivianite and P burial because of more rapid burial of Fe oxides below the SMTZ. In our model scenario, sedimentation rates greater than 0.5 cm yr⁻¹ were particularly conducive to vivianite formation. Vivianite formation is also highly sensitive to the input of organic matter. When the input of organic matter is lowered to −75 % of that assumed in the baseline scenario, the release of Fe²⁺ to the porewater due to Fe oxide reduction decreases to such an extent that formation of vivianite is limited by the availability of Fe²⁺. Consequently, P burial decreases. An increased input of organic matter relative to the baseline scenario enhances the formation of H₂S and leads to a decline in vivianite formation and P burial. When assuming that changes in the rate of sedimentation and the input of Fe oxides and organic matter are coupled (Fig. 9), the net effect is an increase in vivianite formation and P burial. In summary, conditions for vivianite formation are most favorable in sediments with a low bottom water salinity, a high sedimentation rate and a high input of Fe oxides. While high rates of organic matter input are also required, there is an optimum rate above which vivianite formation declines.

5.3 Variations in riverine Fe and organic matter fluxes

In the Bothnian Sea, boreal rivers are the main source of Fe and organic matter to the coastal zone (Algesten et al., 2006; Björkvald et al., 2008; Palviainen et al., 2015). The input from rivers can be highly variable between years (Hölemann et al., 2005; Rember and Trefry, 2004; Sarkkola et al., 2013). The variations in sediment Fe, S and Fe/Al with depth, with maxima at depths of 21, 42 and 60 cm at site NB8, are in accordance with such a variation in input (Fig. 5). This is confirmed by the model, since a strong increase in input of Fe oxides and organic matter is required to describe the maxima in Fe and S in the SMTZ (Fig. 5) and the maximum in vivianite directly below it (Fig. 6). This increased input of Fe oxides and organic matter leads to enhanced conversion of Fe oxides to Fe sulfides. We suggest that these high inputs of Fe and organic matter are directly linked to variations in river discharge, that may be linked to negative phases in the North Atlantic Oscillation (Gillibrand et al., 2005). In 1996 a dry period affected the entire Bothnian Sea region (Fig. S8; Marttila et al., 2016). In that year, the Öre River discharge was unusually low (Fig. 2). In the following years, i.e., starting in 1997, this then led to an increased input of Fe downstream (Sarkkola et al., 2013). This timing coincides with the maximum in Fe and organic matter input assumed in our model scenario (Fig. 3). In the modeling study for the deep basin site in the Bothnian Sea the timing of the enhanced Fe input also corresponds to the year 1997 (Rooze et al., 2016), implying that a wide area was affected. Although not included in the model, a similar scenario would explain the enrichments of Fe, S and Fe-bound P at greater depth in the sediment, with the 1976 dry period leading to enhanced input of Fe and organic matter starting in 1977 and the Fe enrichment at 42 cm depth. In contrast to Fe, limited reoxidation of Mn is expected during a dry period because of less formation of reduced Mn species (Burdige, 1993). Therefore, no strong increase in riverine Mn input is expected after a dry period. In summary, our results show that variations in rainfall on land play a critical role in the transport of Fe and organic matter in rivers to the coastal zone. This riverine transport ultimately determines how much P is buried in the coastal zone. In the coastal Bothnian Sea, P burial is twice as high as the input of P from land, implying that large amounts of the buried P are imported from the open sea (Asmala et al., 2017). Part of the P in the open sea is known to be imported from the adjacent eutrophic Baltic Proper (Savchuk, 2005). The high P burial in the coastal zone coupled to the input of Fe from rivers likely contributes to the seasonal P limitation observed in the Bothnian Sea (Tamminen and Andersen, 2007).

5.4 Implications

In most marine systems, the balance between H₂S formation and Fe oxide input is such that most Fe oxides are converted to Fe sulfides, thereby allowing only limited vivianite formation (Ruttenberg, 2003). Here we show that in estuaries with high inputs of Fe oxides from rivers, high sedimentation rates and sufficient input of organic matter this balance can become very favorable for the formation of vivianite. Such estuaries can act as a highly effective “coastal filter” for P from rivers and for P imported from the open sea (Asmala et al., 2017; Bouwman et al., 2013).

Our study also highlights the role of bottom water salinity in vivianite formation. While vivianite is a key sink for P in the Öre Estuary at a salinity of 5, we show that its formation is enhanced further at low salinity. Many studies focus only on the relatively high salinity parts of estuaries (e.g., as compiled for the Baltic Sea in Asmala et al., 2017). This may explain why the role of vivianite as a sink for P in these systems has been largely overlooked.

Sea level rise is expected to increase bottom water salinity in coastal areas in future (Scavia et al., 2002). In the Baltic Sea, however, salinity is expected to decrease because of enhanced precipitation on land and associated higher river runoff (Graham, 2004; Meier et al., 2006). Recent studies also indicate that Fe²⁺ input from rivers in Europe and North America increased over the past few decades, possibly due to changes in land use and expanded forestry (Björnerås et al., 2017; Kritzberg and Ekström, 2012). For the coastal zone of the Baltic Sea, both a decline in salinity and a continued enhanced input of Fe from rivers will promote conditions for vivianite formation in the sediment and increase its role as a sink for P. Detailed studies of other coastal areas are required to allow a prediction of the role of vivianite formation in the global coastal ocean, since an increased salinity and increased Fe input have opposite effects on the formation of vivianite.

6 Conclusions

Our study reveals a very high burial of phosphorus (P) in sediments of an oligotrophic estuary in the northern Bothnian Sea (up to 0.145 $\mathrm{mol}\mathrm{P}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). These sediments act as a highly efficient sink for P from rivers and the open sea. We demonstrate that the high P retention in the sediment is related to the formation of vivianite, which is visible in the form of discrete crystals below the sulfate methane transition zone at the site with the highest sedimentation rate. By combining P extractions and reactive transport modeling of solid-phase and porewater profiles, we demonstrate that vivianite accounts for 40 % of total P burial at the same site. With the model, we also assess the sensitivity of vivianite formation to the key drivers of variations in the availability of H₂S versus Fe oxides. We show that vivianite formation is promoted when bottom water salinity is low and rates of sedimentation and inputs of Fe oxides are high. High organic matter inputs are also a requirement, but in our scenario there is a distinct optimum for vivianite formation.

We suggest that enrichments of Fe, S and P in the sediment are linked to periods of enhanced riverine input of Fe and organic matter. Two sets of enrichments can be coupled to the years 1977 and 1997, when riverine fluxes of Fe were likely enhanced in a period following a dry year on land in 1976 and 1996, respectively. The enhanced input of Fe and organic matter likely increased the formation of vivianite and the burial of P in the sediment.

In the future, climate change is expected to enhance the freshwater input to the Baltic Sea and thereby decrease its salinity. Continued elevated Fe²⁺ input from rivers is also expected. Both a decreasing salinity and an increased Fe input will create more favorable conditions for vivianite formation in the sediment. We therefore expect that, in future, vivianite will become more important as a sink for P in the coastal zone of the Baltic Sea.