Variations in river input of iron impact sedimentary phosphorus burial in an oligotrophic Baltic Sea estuary

Estuarine sediments are key sites for removal of phosphorus (P) from rivers and the open sea. Vivianite, an iron (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 mol m−2 yr−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 5 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 10 enhanced riverine Fe and organic matter input linked to variations in rainfall on land. 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 15 imported from eutrophic parts of the Baltic Sea.


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 (Froelich et al., 1982;Delaney, 1998;Ruttenberg, 2003). This removal allows coastal systems to act as a filters for P, reducing the flux of land-derived P to the open sea (Froelich, 1988;Bouwman et al., 20 2013; Asmala et al., 2017), and as sinks of P imported from the open sea (Asmala et al., 2017). The major sedimentary phases the northern Baltic Sea, the input of dissolved Fe from rivers has increased over the past decades (Kritzberg and Ekström, 2012;Sarkkola et al., 2013;Björnerås et al., 2017). The associated increased input of Fe oxides to coastal sediments may also enhance vivianite formation.
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 5 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 10 land.

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 2 , a mean depth of 10 m and a bottom water salinity of ca. 5 15 (www.SMHI.se). The estuary is fed by the Öre River, which has a strongly varying rate of annual discharge (Fig. 2). In this region, the spring flood (April, May) is the major annual hydrological event resulting in a brief period of enhanced input of water, dissolved Fe (Hölemann et al., 2005;Björkvald et al., 2008) and terrestrial organic matter (Rember and Trefry, 2004;Algesten et al., 2006). This study focuses on sediments from five sites in the Öre Estuary (Fig. 1B). With water depths varying from 10 to 33 20 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. S.1). 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 −2 (Kauppi et al., 2015). The density of Marenzelleria spp. in the sediment is water depth dependent with maximum densities at water depths 25 between 30 and 50 meters 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 4 ) sampling and the third core for sediment dating using 210 Pb. In August, two cores were collected per site for porewater analyses and for CH 4 . At site NB8, a third and a fourth core were collected in August for 30 mineral analyses and additional 210 Pb measurements.
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 analyser (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 4 and HNO 3 (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 5 was subsequently dissolved in 25 ml of 4.5% HNO 3 at 90 • C overnight. Total elemental concentrations of Al, Fe, Mn, P and S in the HNO 3 solution 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 4-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 hours with 1 10 M Na-acetate (pH brought to 4.5 with acetic acid); (ii) easily reducible Fe oxides (including ferrihydrite and lepidocrocite), extracted for 24 hours with 1 M hydroxylamine-HCl in 25% v/v acetic acid; (iii) reducible (crystalline) oxides (including goethite, hematite and akagenéite), extracted for two hours with Na-dithionite (pH 4.8); (iv) recalcitrant Fe oxides (mainly magnetite), extracted for two hours 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 15 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 3-step S speciation procedure of (Burton et al., 2006(Burton et al., , 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 2 S into a tube filled with 7 ml Zn-acetate; (ii) elemental sulfur (S 0 ) by extracting overnight with 25 ml methanol; (iii) chromium reducible sulfur (CRS, i.e. mostly FeS 2 ) by addition of acidic 20 chromium chloride solution and trapping of the released H 2 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 25 for 30 minutes with 1 M MgCl 2 (pH 8); (ii) Fe-bound P fraction (including Fe oxide bound P and vivianite; Nembrini et al. (1983); Dijkstra et al. (2014)) extracted for 8 hours with citrate-dithionite-bicarbonate (CDB) (pH 7.5) followed by extraction for 30 minutes with 1 M MgCl 2 (pH 8); (iii) authigenic Ca-P (including carbonate fluorapatite, biogenic hydroxyapatite and carbonate-bound P), extracted for 6 hours with Na-acetate buffer (pH 4) followed by extraction for 30 minutes with 1 M MgCl 2 (pH 8); (iv) detrital Ca-P, extracted for 24 hours with 1 M HCl; (v) organic P, extracted for 24 hours with 1 M HCl after ashing 30 at 550 • C for 2 hours. Steps (i)-(iv) were performed under an argon atmosphere to avoid oxidation artefacts (Kraal et al., 2009;Kraal and Slomp, 2014). 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.

Sedimentation rate and P burial
Sediment accumulation rates at all 5 sites were determined from depth profiles of 210 Pb for April (NB1,N6,N10 and N7) or August (NB8). 210 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 210 Pb was calculated as 5 the difference between the measured total 210 Pb and the estimate of the supported 210 Pb activity as given by 214 Pb ( 210 Pb exc = 210 Pb total − 214 Pb). Sediment accumulation rates at each site were estimated by fitting a reactive transport model (Soetaert and Herman, 2008) to the 210 Pb depth profiles assuming depth dependent changes in porosity and bioturbation ( Fig. S.2).
Total P burial (mol m −2 yr −1 ) for all sites was calculated as follows: (1) 10 where P total is the averaged concentration of total P (mol g −1 ) over the deepest 10 cm of the sediment sampled, φ is the porosity in the same interval (vol vol −1 ), sed. rate is the sedimentation rate (cm yr −1 ) and ρ is the density of dry sediment, 2.65 g cm −3 (Burdige, 2006).

Scanning electron microscopy (SEM)
Sediment samples from site NB8 were analyzed by SEM to determine whether large vivianite crystals (≥38 µm size) were Al, silica (Si), P, calcium (Ca), Mn, Fe). Samples were detected and photographed in secondary electron imaging (SEI) mode. 25 Measurements were performed with a 1 µm beam in backscatter mode imaging (BEI).
3 Reactive transport modeling

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. The model describes the mass balance of 11 dissolved and 18 particulate species (  (Soetaert et al., 1996;Wang and Van Cappellen, 1996;Boudreau, 1997). Bioirrigation is modeled as a non-local exchange process (Boudreau, 1984;Emerson et al., 1984).

Model Parameterization
The model was parameterized using data for the field site, information from the literature and by fitting modeled porewater and 15 solid phase depth profiles to the measured data (model constrained;  Fig. S.3A). The density of dry sediment and the C:N ratio for organic matter were taken from literature (Redfield, 1958;Burdige, 2006). Bioturbation and bioirrigation coefficients were used as fitting parameters ( (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 supplements.

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 2 S (Table S.1). In both cases, a highly reactive fraction (α) and a less reactive fraction 25 (β) 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 OM. 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 30 (Jørgensen, 1978;Westrich and Berner, 1984;Middelburg, 1989). 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 (Stockenberg and Johnstone, 1997;Algesten et al., 2006;Leipe et al., 2011). The three fractions of Fe oxides are assumed to have different P:Fe ratio's, 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 S.3; Gunnars et al. (2002)). In the bioirrigated zone, porewater HPO 2-4 is assumed to bind to the Fe oxide β fraction present ( Fig. S.3C&D) to a maximum P:Fe ratio of 0.28 5 (mol mol −1 ; Table S.3). This allows for enhanced formation of Fe oxide bound P due to bioirrigation (Norkko et al., 2012).
Non-reactive P is deposited at the sediment-water interface as detrital P, authigenic Ca-P, P bound to non-reactive Fe oxides (γ fraction) and P bound to non-reactive organic matter (γ fraction).
In the sediment, the rate of vivianite formation (R viv ) is modeled by means of Michaelis-Menten kinetics for dissolved Fe 2+ and HPO 2-4 , which implies that the rate of vivianite formation depends on both porewater species and that there is a maximum 10 rate of formation (Reed 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 HPO 2-4 is not high enough to maintain a high 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 2 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.  Table S.2), were only implemented in the last 12 years of the run, from 2003 onwards.

10
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 first 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 −1 . The input of organic matter was varied between -70% and +25%; (3) input of Fe oxides was varied between -25% and +25% from the baseline scenario. Subsequently, a run was performed in 15 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 (Björkvald et al., 2008;Algesten et al., 2006).

Porewater profiles 20
At all sites, depth profiles of porewater constituents showed relatively little difference between April and August 2015 (Fig. 4).
Porewater SO 2-4 decreased with depth, with a distinct SMTZ only being present in the sampled depth interval at sites N7 and NB8 (Fig. 4)

Solid phase profiles and sedimentation rate
Organic carbon contents in the surface sediment ranged between 2-4 wt% ( Table 1)  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.
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 5 (Fe oxides) or varying (Fe-bound P) concentrations at greater depth (Figs. 6 and S.5). Maxima in AVS coincided with maxima in total S, total Fe and Fe/Al (Figs. 5&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 −1 ( Table 1). Rates of P burial ranged from 0.026 to 10 0.145 mol m −2 yr −1 ( Table 1). The rates of sedimentation and P burial were highest at site NB8. is the release from Fe oxides, with H 2 Sdriven reductive dissolution accounting for 61% of the HPO 2-4 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. S.6; R14&R15).

SEM-EDS
Approximately 30% of the HPO 2-4 removal from the porewater takes place through release to the overlying water by bioir-30 rigation and diffusion, with the former being the most important process. Binding of HPO 2-4 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 Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-327 Manuscript under review for journal Biogeosciences Discussion started: 6 August 2018 c Author(s) 2018. CC BY 4.0 License.
for HPO 2-4 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).
In the model, total P burial rates increased from 0.12 to 0.34 mol m −2 yr −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 5 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 mol m −2 yr −1 when assuming sedimentation rates between 0.25 and 2 cm yr −1 (Fig. 9). At a sedimentation rate of 0.5 cm yr −1 , 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 10 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; Fig. S.7A&B). The response to an increase in the rate of sedimentation is non-linear, with the greatest increase in P burial between 0.5 and 0.85 cm yr −1 . For organic matter, the model runs reveal an optimum input rate beyond which P and vivianite burial rates decrease ( Fig. 10C; Fig. S.7C).

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 mol m −2 yr −1 (Asmala et al., 2017). We find substantially higher rates of P burial in the Öre estuary ranging from 0.026 to 0.145 mol m −2 yr −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 20 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 25 retention. As discussed in detail by Egger et al. (2015a), vivianite formation is generally most pronounced in sediments below the SMTZ, where both Fe 2+ and HPO 2-4 accumulate in the porewater. 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 30 temperature. This implies that any vivianite formed below the SMTZ likely will 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 toward an authigenic origin. The average Fe:P ratio of the crystals of 3.3 mol mol −1 (Table S.5) 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 5 are known to be included in the structure of the mineral (Dijkstra et al., 2016(Dijkstra et al., , 2018Egger et al., 2015a). The presence of Al and Si (Table S.5) likely reflects the presence of clay particles that were not removed prior to the SEM-EDS analysis (Egger et al., 2015a;Dijkstra et al., 2018).
The Fe-P below the SMTZ (Fig. 6) likely is 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). The results of the reactive transport model suggest that 10 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. S.5). The remainder of the P is buried as non-reactive P.

Vivianite formation in coastal sediments in the northern Baltic Sea
Vivianite formation in sediments strongly depends on the balance between the formation of H 2 S and the input of Fe oxides 15 (Ruttenberg, 2003). When there is an excess of Fe oxides over H 2 S, not all Fe oxides will be converted to Fe sulfides, and more Fe 2+ will be available to precipitate with HPO 2-4 as vivianite (Rozan et al., 2002;Gächter and Müller, 2003). 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 20 Fe oxides below the SMTZ. In our model scenario, sedimentation rates greater than 0.5 cm yr −1 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 2+ to the porewater due to Fe oxide reduction decreases to such an extent that formation of vivianite is limited by the availability of Fe 2+ . Consequently, P burial decreases.
An increased input of organic matter relative to the baseline scenario enhances the formation of H 2 S and leads to a decline 25 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. 30

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 (Rember and Trefry, Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-327 Manuscript under review for journal Biogeosciences Discussion started: 6 August 2018 c Author(s) 2018. CC BY 4.0 License. 2004;Hölemann et al., 2005;Sarkkola et al., 2013). The variations in sediment Fe, S and Fe/Al with depth, with distinct maxima at depths of 21, 42 and 60 cm at site NB8, are in accordance with such a variation in input. 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). We suggest that these high inputs of Fe and organic matter are directly linked to variations in river discharge. In 1996, a dry period affected the entire Bothnian Sea region (Fig. 5 S.8; Marttila et al. (2016)). In that year, the Öre River discharge was unusually low (Fig. 2). In the following wet 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 10 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 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 15 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).

Implications
In most marine systems, the balance between H 2 S formation and Fe oxide input is such that most Fe oxides are converted to 20 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 (Bouwman et al., 2013;Asmala et al., 2017).
Our study also highlights the role of bottom water salinity in vivianite formation. Because many studies focus on the rela- 25 tively high salinity parts of estuaries (e.g. as compiled for the Baltic Sea in Asmala et al. (2017)), the role of vivianite as a sink for P 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 2+ input from rivers in Europe and North America increased over 30 the past few decades, possibly due to changes in land use and expanded forestry (Kritzberg and Ekström, 2012;Björnerås et al., 2017). 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 13 Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-327 Manuscript under review for journal Biogeosciences Discussion started: 6 August 2018 c Author(s) 2018. CC BY 4.0 License. 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.

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 mol P m˘2 yr˘1). These sediments act as a highly efficient sink for P from rivers and the open sea. We demonstrate 5 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 2 S versus Fe oxides. We show that vivianite formation is promoted when bottom water salinity is low and 10 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 wet period directly following a dry period on land in 1976 and 1996, respectively. The enhanced input of Fe and organic 15 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 2+ 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.

20
Data availability. All data will be made available in the database Pangea upon acceptance of the manuscript.
Code and data availability. The model code is available from WL and CS upon request.
Author contributions. CS, WL, NH and DC designed the research. WL, ME, NH and DC performed the analyses. WL performed the model simulations. WL, ME, NH, EK and CS interpreted the data. WL and CS wrote the paper with comments provided by ME, NH, EK and DC.
Competing interests. The authors declare that they have no conflict of interest.       9. Sensitivity of the rates 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.  Fig. 10. Sensitivity analysis of P burial as a function of A: changing sedimentation rate between 0.25 and 2 cm yr −1 ; B: organic matter input ranging between -70% and +25% C: changing input of Fe oxides ranging between -25% and 25%. The corresponding burial forms of P are given in Figure S.7.