Deep-water inflow event increases sedimentary phosphorus release on a multi-year scale

Phosphorus fertilisation (eutrophication) is expanding oxygen depletion in coastal systems worldwide. Under low-oxygen bottom water conditions, phosphorus release from the sediment is elevated, which further stimulates primary production. It is commonly assumed that re-oxygenation could break this “vicious cycle” by increasing the sedimentary phosphorus retention. Recently, a deep-water inflow into the Baltic Sea created a natural in situ experiment that allowed us to investigate if temporary re-oxygenation stimulates sedimentary retention of dissolved inorganic phosphorus (DIP). Surprisingly, during this 3-year study, we observed a transient but considerable increase, rather than a decrease, in the sediment efflux of DIP and other dissolved biogenic compounds. This suggested that the oxygenated inflow elevated the organic matter degradation in the sediment, likely due to an increase in organic matter supply to the deeper basins, potentially combined with a transient stimulation of the mineralisation efficiency. As a result, the net sedimentary DIP release per m² was 56 %–112 % higher over the years following the re-oxygenation than before. In contrast to previous assumptions, our results show that inflows of oxygenated water to anoxic bottom waters can increase the sedimentary phosphorus release.

1 Introduction

Eutrophication is one of the main causes of oxygen depletion in coastal systems worldwide (Breitburg et al., 2018; Diaz and Rosenberg, 2008). Excess input of nutrients from land stimulates primary production, resulting in a higher delivery of organic matter to deeper coastal waters and the seafloor (Breitburg et al., 2018; Middelburg and Levin, 2009). Oxygen is consumed as this organic matter is degraded, and in severe cases anoxia can develop. Depletion of oxygen can lead to decreasing biodiversity and thereby a loss of important ecosystem functions and services (Breitburg et al., 2018; Middelburg and Levin, 2009). One important example is food supply, since many coastal areas are key sites for fishing and aquaculture (Diaz and Rosenberg, 2008).

The Baltic Sea is strongly affected by eutrophication and oxygen depletion. The area has been argued to be a “time machine” for future coastal areas worldwide, due to its long history of anthropogenic disturbances and impacts of climate change (Reusch et al., 2018). As such, understanding the feedback mechanisms between eutrophication and oxygen depletion in the Baltic Sea can help us understand the future of our coastal zones. The central parts of the Baltic Sea have been naturally oxygen-depleted for thousands of years due to a strong water column stratification, a long water residence time and limited deep-water renewal (Zillén et al., 2008). However, high inputs of nutrients from land have severely worsened the situation and the areal extent of benthic oxygen depletion has increased 10 times in size over the last 115 years (Carstensen et al., 2014). Despite improvements in wastewater treatment and decreased nutrient input from agriculture, the oxygen-depleted area in the Baltic Sea has not decreased in size (Reusch et al., 2018).

The low-oxygen conditions in the Baltic Sea are largely sustained by enhanced recycling of dissolved inorganic phosphorus (DIP, mainly phosphate) from sediments with oxygen-depleted overlying water (Vahtera et al., 2007). This “vicious cycle” has been suggested to drive enhanced recycling of DIP. High sedimentary release of DIP fuels extensive growth of cyanobacteria and other phytoplankton, which fuel microbial respiration and oxygen consumption as they are degraded, resulting in increased oxygen depletion and further release of DIP (Vahtera et al., 2007). The elevated release of DIP from oxygen-depleted sediments is mainly a consequence of two factors. Firstly, iron (Fe) and (to some extent) manganese (Mn) oxides, which adsorb and retain DIP in the sediment under oxic conditions, are reduced and solubilised in anoxic environments (Jilbert and Slomp, 2013; Yao and Millero, 1996). Secondly, P is preferentially regenerated from organic matter relative to nitrogen and carbon, which causes an elevated sedimentary DIP release in environments lacking Fe and Mn oxides (Jilbert et al., 2011; Steenbergh et al., 2011, 2013). In anoxic environments, DIP is thus released back to the water column instead of being retained in the sediment where it eventually could be buried as authigenic P minerals (Slomp et al., 1996). Due to high sedimentary DIP concentrations there is also some formation of authigenic P phases in the anoxic areas of the Baltic Sea (Jilbert and Slomp, 2013), but the overall burial efficiency of P in this environment is low (0.2 %–12 %, compared to >14 % in oxic parts of the Baltic Sea; Viktorsson et al., 2013).

It is commonly assumed that re-oxygenation of anoxic sediments would reduce the release of DIP (Ruttenberg, 2003; Stigebrandt, 2018), which could break the vicious cycle. Intrusions of large volumes of salty, oxygenated water from the North Sea (major Baltic inflows; MBIs), can temporarily bring new oxygen to the anoxic areas of the Baltic Sea (Liblik et al., 2018). In the winter of 2014–2015, the largest MBI since 1951 occurred and previously long-term anoxic areas were oxygenated (Liblik et al., 2018). Studies conducted shortly after the inflow did indeed show a lowered DIP release from the sediment, but the total DIP release only decreased by 5 %–23 % on a basin-wide scale (Hall et al., 2017; Sommer et al., 2017), and sequestration of P in the surface sediment was low (Hermans et al., 2019b). It has so far been unknown how the sedimentary P dynamics evolved after anoxia was re-established.

Here we present a time series study from the Eastern Gotland Basin (EGB, Fig. 1a) in the Baltic Sea, a long-term anoxic basin that was temporarily oxygenated by the 2014–2015 MBI. To investigate the effect of the oxygenation event on the sedimentary P cycle, we conducted in situ measurements of sediment–water fluxes and collected pore water and solid-phase samples during 3 years (2016–2018). By combining our results with data from 2008 (Nilsson et al., 2019; Viktorsson et al., 2013), 2010 (Nilsson et al., 2019; Viktorsson et al., 2013), and 2015 (Hall et al., 2017), we construct a unique time series from prior to until after the inflow and show that, contrary to previous assumptions, the inflow of oxygenated water to this anoxic area increased the sedimentary P release on a multi-year scale.

Figure 1 Description of sampling sites and oxygen conditions. (a) Map of the sampling sites in the Eastern Gotland Basin, Baltic Sea. Bathymetric data from the Baltic Sea Bathymetry Database (v. 0.9.3), the map was created using the software ArcMap™ (v. 10.6) by Esri. (b) Water column profiles of oxygen (O₂). Data from SMHI (SMHI, 2021; 2014, station BY15), Hall et al. (2017; 2015, station F) and van de Velde et al. (2020b; 2016–2018, station F). (c) Bottom water oxygen (O₂) concentrations at stations D–F measured by optodes on the benthic chamber lander (average ± standard deviation) ∼ 0.2 m above the sediment. The detection limit was 0.5 µM. Data from Hall et al. (2017) and van de Velde et al. (2020b). ^* No measurement.

2 Materials and methods

2.1 Site characteristics

The EGB is the largest sub-basin of the Baltic Sea, situated between Latvia and the Swedish island of Gotland (Fig. 1a). The water column is normally hypoxic below the halocline at 60–80 m depth and anoxic at the seafloor (maximum depth 249 m). The MBI that occurred in the winter of 2014–2015 reached the EGB in March 2015, resulting in a more oxygenated but still hypoxic (40–60 µM) water layer below 130–140 m depth (Fig. 1b, c; Hall et al., 2017; Liblik et al., 2018). The oxygen had been consumed by the end of 2015, but a new inflow oxygenated the area again between January and May 2016 (Liblik et al., 2018). In 2017, several smaller inflows also brought oxygen to the 80–170 m depth layer (Liblik et al., 2018).

During yearly expeditions to the EGB in April 2016–2018, we visited three stations along a depth transect (Fig. 2a, Table A1 in Appendix A). The stations are situated below the halocline and were oxygenated by inflows at different points during the sampling period. The sediment is rich in organic matter (organic carbon 10 %–15 % dry weight) and has a high surface porosity (0.95–0.98; Nilsson et al., 2019; van de Velde et al., 2020b). Station D (130 m) is situated at the top boundary of the oxygenated layer. It was anoxic during the sampling in 2016 (Fig. 2c), but monitoring data from the Swedish Meteorological and Hydrological Institute (SMHI, Fig. A1 in Appendix A; SMHI, 2021) and previous studies (Dellwig et al., 2018; Hermans et al., 2019b; Sommer et al., 2017) from the eastern part of the EGB suggest that this station might have been temporarily oxygenated by the MBI prior to our visit. Station D received some oxygen with the smaller inflow in 2017 but was anoxic again in 2018. Stations E (170 m) and F (210 m) were oxygenated by the MBI (Fig. 2c). The oxygen at station F had already been consumed in 2017, while it decreased stepwise at station E. Hydrogen sulfide was only detected in the bottom water at station D in 2016 (Broman et al., 2020; Marzocchi et al., 2018).

Figure 2 Sediment–water fluxes measured with the benthic chamber lander. Flux data are from before the inflow (2008, 2010; Hall et al., 2017; Nilsson et al., 2019; Viktorsson et al., 2013), right after the inflow (2015; Hall et al., 2017), and post-inflow (2016–2018; this study). (a) Changes in fluxes with time. Small circles mark individual chamber measurements, crosses mark station averages, and lines represent station standard deviations. The black line in the DIP graph for station F in 2016–2017 marks the expected flux for those years, based on the decrease in solid phase inventory. N.D. stands for no data. (b) Fluxes of DIC vs. DIP and DIN. The grey line represents the Redfield ratios of 106 : 1 for DIC : DIP and 106 : 16 for DIC : DIN, the shaded area marks the DIC : DIP ratios between 20 : 1 and 70 : 1. Colour-coding is the same as in (a).

For clarity, the remainder of this paper will focus on the results from station D and F. These stations represent the least- and most-affected sites, respectively. The results from station E, which are consistent with the results from station F, are presented in Figs. A2–A4.

2.2 Water column and benthic in situ flux measurements

Upon arrival at the stations, a conductivity–temperature–depth instrument (CTD; SBE 911, Sea-Bird Scientific) with a high-accuracy oxygen sensor (SBE 43, Sea-Bird Scientific) was deployed to record water column profiles of salinity, temperature, and oxygen (van de Velde et al., 2020b). Sediment–water fluxes of nutrients and dissolved inorganic carbon (DIC) were subsequently measured using the Gothenburg benthic chamber lander (Kononets et al., 2021). Deployment of the lander during these samplings has been described in detail elsewhere (Hall et al., 2017; van de Velde et al., 2020b). Briefly, the lander carried four incubation chambers incubating 400 cm² of sediment with overlying water each. At preset times, nine syringes per chamber sampled the overlying water. Sensors monitored conductivity, temperature, and oxygen inside and outside the chambers. Chambers that initially contained oxygen did not become anoxic over the course of the incubations. Chamber volumes were calculated from the decrease in salinity following the injection of a small, known volume of Milli-Q water. Sediment and water were incubated for 14 h. Immediately after lander recovery, the syringes were emptied and samples were filtered through pre-rinsed 0.45 µm cellulose acetate filters (Sartorius). The lander was deployed twice per station except for two occasions (F 2017, E 2018) when time constraints only allowed one deployment.

Calculation of benthic nutrient and DIC fluxes are described in van de Velde et al. (2020b). The fluxes from 2008, 2010, and 2015 (Hall et al., 2017; Nilsson et al., 2019; Viktorsson et al., 2013) were measured using the same lander systems as in this study and have been re-evaluated to assure that the same evaluation method is used and thus the results are comparable. In short, a linear regression model, either simple or quadratic, was fitted to data of the concentration change over time in the chamber water. Models were selected based on their Akaike information criterion score and corrected for small sample sizes. After verifying model and data quality, the flux was calculated by multiplying the chamber water height with the initial slope of the regression line. About 5 % of the regression models gave slopes for which p>0.05. After visual inspection of the data, all fluxes were kept in the data set in order to not overestimate the average fluxes by excluding low flux values.

The effects of year and station (fixed factors) on sediment–water fluxes were tested with a type III ANOVA, with lander deployment as a nested factor within year and station. Before analysis, flux data were inspected visually for homogeneity of variance and were log₁₀ transformed. Significant effects of year and station were investigated using Student–Newman–Keuls (SNK) tests. All tests were two-sided with a significance level of 0.05. Analyses were conducted in the statistical software R (R Core Team, 2020), using the functions “lm” and “anova” from the standard CRAN:stats package and “SNK.test” from the CRAN:agrocolae package (de Mendibur, 2020).

Deployment of benthic chamber landers is a commonly used technique to measure in situ fluxes and avoid potential manipulation artefacts during ex situ incubations. However, some caution should be taken when interpreting the absolute values of the measured fluxes. For example, during incubation of sediment from sites with stagnant water, erosion of the diffusive boundary layer by stirring of bottom water can increase the sedimentary oxygen uptake by up to 30 % (Glud et al., 2007; Hall et al., 1989). It is likely that changes to the thickness of the benthic boundary layer also affect benthic releases of DIC and nutrients, but this is yet to be confirmed. Nevertheless, the fluxes in this study may be slightly overestimated. However, during deployment the stirring mechanism was made to resemble the natural bottom currents as much as possible, and we expect the impact on our measured fluxes to be smaller than 30 %. Furthermore, such an impact on the fluxes would only affect absolute values, not patterns between years.

2.3 Sediment sampling

Pore water was extracted from sediment cores (9.9 cm inner diameter) collected with a multiple corer. At each station, one core was sectioned in a portable glove bag (N₂ atmosphere; Aldrich AtmosBag, Sigma-Aldrich). The cores were sliced at 0.5 cm resolution from 0 to 2 cm depth, at 1 cm resolution between 2 and 6 cm depth, and in 2 cm slices from 6 to 20 cm depth. The sediment was transferred into centrifuge tubes into which Rhizon samplers (Rhizosphere research products) were inserted for collection of pore water.

Samples for solid-phase speciation of P were collected at station D and F in 2016 and 2017. Due to logistical constraints, samples were not collected on the other occasions. Sediment was collected using a modified box corer (Blomqvist et al., 2015) from which 6 cores were taken using transparent PVC core liners (6 cm inner diameter; 30 cm long). Two cores per station were sectioned in a portable glove bag (N₂ atmosphere; Captair Field pyramid, Erlab) under constant oxygen monitoring (Portable oxygen analyser 3110 with trace oxygen sensor, Teledyne). The cores were sliced at 0.5 cm resolution from 0 to 3 cm depth, at 1 cm resolution between 3 and 10 cm depth, and in 2 cm slices from 10 to 14 cm depth. The sediment was transferred into polypropylene centrifuge tubes (TPP), centrifuged at 2500 g for 10 min (Sigma 3-16L, Sigma) and was returned to the anoxic glove bag. After removal of the overlying pore water, the solid phase was freeze-dried and stored in an aluminium bag under N₂ atmosphere for later speciation of P. In 2016, the pore water was filtered (0.42 µm cellulose acetate filters, Chromafil Xtra) and acidified with 100 µL mL⁻¹ H₂SO₄ (1 M) for later analysis as a control of our pore water sampling procedure.

2.4 Bottom water and pore water analysis

Nutrient samples from the lander and multiple corer pore water samples were determined by segmented flow colorimetric analysis (Alpkem Flow Solution IV, OI Analytical; Koreleff, 1983). The analytical precisions were 3 % for DIP, 3 % for dissolved inorganic nitrogen (DIN) and 2.5 % for dissolved silica (DSi). In 2016, pore water DIP samples from the box corer were determined by segmented flow colorimetric analysis (QuAAtro, SEAL Analytical; Aminot et al., 2009). The precision was <4 %.

DIC concentrations in the lander samples were determined on-board the ship by non-dispersive infrared spectrometry (LI-COR 6262) after acidification with phosphoric acid (Nilsson et al., 2019). Repeated measurements of certified reference material (CRM from Dickson Laboratory, Scripps Inst. of Oceanography) were conducted to obtain a two-point calibration and correction for potential drift in the system. The analytical precision was 0.2 % (relative standard deviation).

2.5 Sequential phosphorus extractions

A subsample of ∼ 100 mg sediment from each depth was used for sequential P extraction (SEDEX; Ruttenberg, 1992; Slomp et al., 1996). The SEDEX procedure separates total sedimentary P into five operationally defined fractions: exchangeable or loosely sorbed P (P_ex), easily reduced or reactive Fe bound P (P_Fe), authigenic P (P_auth), detrital P (P_det), and P associated with organic matter (P_org). All extractions were performed at room temperature and under constantly agitated conditions; extraction solutions and times are shown in Table 1. After each extraction, the sample was centrifuged (2500 g for 10 min) and the supernatant was filtered (0.45 µm cellulose acetate) and subsequently analysed spectrophotometrically. The first two extraction steps were performed under oxygen-free conditions. The P_Fe fraction was analysed by inductively coupled plasma mass spectrometry (ICP-MS) since the reagent interferes with the spectrophotometric colour reagent. Performance of the method was monitored by analysing one sample in triplicate during each batch of extractions (to assess the repeatability by calculating the relative standard deviation) and by analysing subsamples of one sample with every batch of extractions (to assess reproducibility by calculating the relative standard deviation). For each extraction step, repeatability was around 5 % and reproducibility was around 10 %. Calculations of sedimentary inventories of different P fractions are described in Appendix B.

Table 1 Specifics for the sequential extractions of phosphorus. A brief summary of the reagents, methods, and times used for SEDEX. Each step was followed by a wash with 1 M MgCl₂ (pH 8) for 30 min to avoid re-adsorption of DIP.

3 Results and discussion

3.1 Enhanced sedimentary P release after the inflow

Sediment–water fluxes of DIP measured in situ showed significant changes with time following the MBI (Fig. 2a, Tables A2, A3). At station F, the DIP efflux decreased markedly in 2015 (by 81 %) compared to the 2008 and 2010 (pre-inflow) values and even turned into an uptake in about half of the chambers (Hall et al., 2017). In 2016, however, the DIP flux was about 4.5 times higher than before the inflow, after which it decreased stepwise in 2017 and 2018, returning to pre-inflow values. Although the DIP fluxes were highest at station F, similar patterns were observed at all stations (Fig. A2, Tables A2, A3). Intriguingly, the DIN fluxes followed the same pattern as the DIP fluxes, and similar trends were seen for DIC and DSi.

The decrease in DIP flux at station F in 2015 temporarily raised the DIC:DIP flux ratio above the Redfield ratio (Redfield, 1958) of 106 : 1 (Fig. 2b; Hall et al., 2017). As the DIC : DIN ratio did not change, the lowered DIP flux was most likely due to trapping of P by Fe and Mn oxides, formed at the sediment surface as a result of the oxygenation. When anoxia returned in 2016, P_Fe that had been temporarily trapped during oxic conditions was expected to be released and increase the sedimentary DIP efflux, which would lower the DIC : DIP flux ratio compared to pre-MBI conditions. However, the DIC : DIP ratio only decreased to 30–40 in 2016–2018, which is within the range of what is normally observed in the anoxic parts of the EGB (∼ 20–70; Sommer et al., 2017; Viktorsson et al., 2013; Fig. 2b). These low DIC : DIP flux ratios are the result of preferential regeneration of P from organic matter, causing an elevated sedimentary DIP release in the absence of Fe oxides (Jilbert et al., 2011; Steenbergh et al., 2011, 2013), and release of DIP adsorbed to Fe oxides that are shuttled from shallower areas (Dellwig et al., 2010). The consistent DIC : DIP ratio in all years apart from 2015 (right after the MBI) implies that factors other than release of temporarily trapped P_Fe contributed to the increased benthic DIP efflux at station F in 2016–2017. An intriguing clue to what could have caused the elevated DIP fluxes was given by the fluxes of other biogenic solutes, which followed the same pattern of temporarily increased fluxes (Fig. 2a, Tables A2, A3). The co-variation suggests that a common factor affected all solute fluxes, hinting towards an increase in the organic matter mineralisation rate and/or organic matter supply. The pattern in fluxes over time observed at station F was also significant at station D, but the fluxes were lower and the pattern was less pronounced. A sedimentary uptake of DIP was measured in two out of seven chambers in 2017, suggesting that the oxygenation of the sediment was weak and heterogeneous. Significant differences in fluxes between lander deployments indeed suggest that the sediment environment in the entire area was heterogeneous (Table A2). Station D is situated at the edge of the area that was affected by the MBI and the oxygenation there was less strong (Fig. 1c). The conditions at station D can thereby be assumed to represent a lower boundary for effects from the inflow. At both stations D and F, the DIC : DIN flux ratios were around 10–12 (Fig. 2b) rather than the Redfield ratio of 6.6 (106 : 16). The flux ratios in this study agree well with carbon : nitrogen ratios at the sediment surface and in sinking particles in the EGB (Cisternas-Novoa et al., 2019; Nilsson et al., 2021).

3.2 Sedimentary phosphorus dynamics after the inflow

Sediment profiles from all stations showed clear signs of oxygenation after inflows (Figs. 3a, 4a). At station F, the pore water DIP displayed distinct subsurface peaks in 2016 (Fig. 3a). This type of profile indicates production of DIP in the upper sediment layers (Ruttenberg, 2003), probably caused by a sudden influx of P in the form of reactive organic matter or P adsorbed to Fe and Mn oxides. In contrast, pore water DIP gradually accumulated with depth in 2017 and 2018. Both the P_ex and P_Fe concentrations were higher in 2016 than in 2017 (Fig. 3a; other P fractions in Fig. A4). This suggests that there was a transient enrichment of these fractions in 2016 that had disappeared in 2017, which would explain the DIP peak.

Figure 3 Pore water and sediment solid-phase data from station F. (a) Profiles of pore water DIP, solid-phase exchangeable or loosely-sorbed P (P_ex), easily reduced or reactive iron-bound P (P_Fe), and organic P (P_org). Profiles with light blue markers are from big sediment cores (inner diameter 9.9 cm), and profiles with dark blue and orange markers are from replicate small sediment cores (inner diameter 6 cm). Markers with black edges are missing data that have been replaced with average concentrations from the slices directly above and below (see Appendix B). (b) Excess inventories (min–max range) of P_ex and P_Fe and inventories of P_org in the top 5 cm.

Figure 4 Pore water and sediment solid-phase data from station D. (a) Profiles of pore water DIP, solid-phase exchangeable or loosely sorbed P (P_ex), easily reduced or reactive iron-bound P (P_Fe), and organic P (P_org). Profiles with light blue markers are from big sediment cores (inner diameter 9.9 cm), and profiles with dark blue and orange markers are from replicate small sediment cores (inner diameter 6 cm). Markers with black edges are missing data that have been replaced with average concentrations from the slices directly above and below (see Appendix B). (b) Excess inventories (min–max range) of P_ex and P_Fe and inventories of P_org in the top 5 cm.

The sequestration of P after oxygenation was calculated as the excessive inventories (Fig. B1 in Appendix B) of P_ex and P_Fe, i.e. the integrated enrichment of these fractions compared to their background concentrations. The excess inventory of P_ex and P_Fe suggest that the net P retention at station F by the sampling occasion in 2016 was 4.6–6.4 mmol P m⁻² (Fig. 3b), consistent with a study from the eastern part of the EGB where 12 mmol P m⁻² was estimated to have been retained due to the MBI (Hermans et al., 2019b). By the sampling occasion in 2017, the net excess inventory had decreased to 3.1–3.5 mmol P m⁻². This means that 1.1–3.3 mmol P m⁻² was released from temporary storage in the sediment between the 2016 and 2017 samplings. As no increase was seen in the P_auth fraction (Fig. A4), it appears that most of this P was lost to the overlying water. A release of 3.3 mmol P m⁻² in 1 year would lead to an average sedimentary efflux of 0.01 mmol m⁻² d⁻¹. However, the benthic DIP flux measured in 2016–2017 was 1.2–2 mmol m⁻² d⁻¹ higher than before the inflow (Fig. 2a). As suggested by the DIC : DIP ratio, another source of P must have contributed to the DIP flux at station F in 2016–2018.

Signs of oxygenation were also seen at station D. In 2016 and 2018, the pore water DIP accumulated with sediment depth. However, the profile was very different in 2017 and displayed a peak at 4 cm depth, suggesting an increased release of reactive P (Ruttenberg, 2003). Subsurface peaks of P_ex and P_Fe at station D were visible in one of the cores from 2016. The variation between replicate cores in 2016 raises the question whether changes in pore water DIP, P_ex, and P_Fe at station D were caused by spatial variability or if there was a change over time in response to oxygenation. However, the changes in sedimentary P profiles and fractionation follow what would be expected if station D had been partially affected by the MBI; a likely scenario based on oxygen monitoring data (Fig. A1; SMHI, 2021). The 2017 inflow had a smaller effect on the sediment geochemistry at station D, since no peaks were present in the cores collected that year. At station D, the P_ex + P_Fe excess inventory change between the 2016 and 2017 samplings ranged from +1.3 to −21.5 mmol m⁻² (Fig. 4b). As there was no increase in the authigenic P fraction (Fig. A4), it is likely that this P was released to the water column as DIP. An excess P release at station D of 21.5 mmol m⁻² would only have increased the DIP flux by 0.06 mmol m⁻² d⁻¹ during 2016–2017. The DIP flux at station D did indeed not change substantially over the sampling (Fig. 2).

The in situ measured DIP fluxes were considerably higher than diffusive fluxes calculated from pore water profiles (Table C1 in Appendix C). Moreover, the P_org inventories in the top 5 cm of the sediment (Figs. 3b, 4b) could only support the DIP fluxes measured in 2016–2017 for 59–159 d at station D and 15–24 d at station F. This discrepancy between measured fluxes and calculated diffusive fluxes and sedimentary P_org inventories may partly be explained by potential overestimation of measured fluxes (see Sect. 2.2). However, even in the unlikely case that the measured fluxes overestimated the in situ fluxes by 50 %, the P_org inventories would only support the DIP fluxes for 118–318 d at station D and 30–48 d at station F. Alternatively, it is well known that chemical gradients at the sediment–water interface may not be accurately captured by the limited vertical resolution of pore water profiles, resulting in underestimates when calculating diffusive fluxes (Nilsson et al., 2019; Noffke et al., 2012; Sundby et al., 1986; van de Velde et al., 2020b). Since the sediment at stations D and F had such high porosities (∼ 0.98; van de Velde et al., 2020b; Fig. A5), it is also possible that the topmost fluffy and organic-rich layer was not captured during sampling and processing of sediment cores. Solute fluxes caused by processes in this sediment layer would have been detected in the benthic chamber lander incubations. The high in situ measured fluxes relative to calculated diffusive fluxes and sediment P_org inventories thus suggest that processes occurring at the sediment–water interface were the main driver behind the benthic fluxes, whereas sedimentary inventories played a minor role. Previous studies have indeed suggested that organic matter mineralisation is rapid and efficient in the deep basins of the EGB (Nilsson et al., 2019, 2021), indicating that it mostly takes place at the sediment surface. In addition to vertical deposition of phytoplankton material, there is an extensive horizontal shuttling of material from shallow to deeper parts of the EGB (Christiansen et al., 2002; Jonsson et al., 1990; Leipe et al., 2011), which is thought to support about half of the benthic organic matter mineralisation in the deep basin on an annual basis (Nilsson et al., 2021).

If the topmost sediment layer was not accurately captured during sampling, the calculated P_ex and P_Fe inventories are likely also too low. This does not necessarily mean that the change in inventories between years is underestimated, if equally large fractions of P_ex and P_Fe were missed both years. Much of the change in P_ex and P_Fe between 2016 and 2017 is expected to have taken place at the sediment–water interface, however, suggesting that the inventory change was in fact larger than estimated here. Nevertheless, the constant DIC : DIP flux ratio in all years but 2015, together with similar flux patterns for multiple biogenic solutes, suggests that the contribution from P_ex and P_Fe to the total DIP flux was in fact small.

3.3 Increased mineralisation after the inflow

Our data suggest that the elevated fluxes of DIP and other dissolved biogenic elements were caused by increased rates of sedimentary organic matter mineralisation. Changes in sedimentary mineralisation rates can be caused by an enhanced input of reactive organic matter to the sediment (for example after a spring bloom or by enhanced organic matter export out of the water column) or by an increase of the intrinsic mineralisation rate of the organic matter already present in the sediment. We will discuss these different possibilities below.

3.3.1 Seasonality, primary production, and benthic fauna

Data for this study were collected during different seasons (2008, 2010: August/September; 2015: early July; 2016–2018: April), and thus some variation in fluxes due to seasonality is expected. However, even though DIC and nutrient effluxes should have been highest in 2008, 2010, and 2015 due to sedimentation of fresh phytoplankton material (Berelson et al., 2003), similar flux values were measured in 2018. Instead, the highest fluxes were recorded during the spring samplings (2016–2018, Fig. 2a, Table A3). The differences in timing of the sampling relative to the spring bloom could still be a potential cause. However, both the lack of clear surficial peaks in the chlorophyll a sediment profiles (Fig. D1 in Appendix D) and monthly water column monitoring data of phytoplankton and nutrient concentrations from monitoring data (Fig. A5; SMHI, 2021) suggest that all samplings in 2016–2018 were conducted before the spring bloom. In fact, the year when the sampling was most likely to coincide with the spring bloom according to the SMHI data was in 2018, which was the year with the lowest fluxes of all spring samplings. It thus seems highly unlikely that seasonality was driving the observed flux pattern.

The MBI could also have pushed up deep water to the surface, thereby providing nutrients to phytoplankton blooms. No increase in primary production was observed, however, during either the spring bloom or summer bloom in the years after the MBI (Johansen and Skjevik, 2016–2017). Instead, the variation of fluxes closely follows the pattern of oxygenation, with higher fluxes just after the oxygenation and a gradual decrease to pre-MBI values once all oxygen had been consumed in the bottom waters (Figs. 1c, 2).

Alternatively, organic matter mineralisation rates are known to increase with temperature (Arnosti et al., 1998). The water temperature below the halocline in the EGB does not follow seasonal patterns (Fig. A6; SMHI, 2021); different sampling months (April vs. August/September) should thus not have affected the benthic fluxes. The MBI increased the bottom water temperature by ∼ 1 ^∘C, however (Fig. A6; Liblik et al., 2018; SMHI, 2021). Despite lower temperatures in 2008 and 2010 compared to 2018, the benthic fluxes were similar during these years (Fig. 2a). Temperature differences are thus unlikely to alone have driven the benthic flux pattern.

Oxygenation can allow re-colonisation of previously anoxic sediments by benthic animals. Their presence would result in sediment reworking (bio-mixing) and bio-irrigation, which are believed to stimulate organic matter degradation (Ekeroth et al., 2016; van de Velde et al., 2020a). Two individuals of Bylgides sarsi were found near our sampling sites in 2015 (Hall et al., 2017; Stigebrandt et al., 2018). However, all sediment retrieved from the newly oxygenated area during our sampling was laminated and there were no signs of colonisation by animals, which is consistent with other studies conducted in the EGB in 2016 (Hermans et al., 2019a). It is thus likely that the two specimens of Bylgides sarsi had arrived with the inflowing water (Stigebrandt et al., 2018). The anoxic water layer at 100–140 m depth would further have created a barrier that prevented re-colonization of the newly oxygenated sediments. We therefore exclude animal activity as an explanation for the enhanced sediment–water fluxes after the MBI.

3.3.2 Enhanced organic matter input

The oxygenation of the bottom waters following a MBI could potentially lead to enhanced reactive particulate organic matter input to the sediment. In oxygen minimum zones, particle fluxes are normally higher than in surrounding oxic waters, possibly because of lower particle mineralisation rates as a result of less energetic electron acceptors or loss of zooplankton that consume particles (Le Moigne et al., 2017). However, a study conducted in the EGB in the summer of 2015 suggested that the particle sedimentation rate and freshness of the organic material reaching the sediment increased due to the MBI (Cisternas-Novoa et al., 2019). In the newly oxygenated EGB, the flux of particulate organic carbon (POC) increased by 18 % between 40 and 180 m depth, while it decreased by 28 % at a permanently anoxic site nearby (Cisternas-Novoa et al., 2019). The cause of the increased POC transport in the EGB was thought to be the formation of Mn oxides in the oxygenated water column, which scavenged organic matter while sinking (Cisternas-Novoa et al., 2019). Several other studies have confirmed that the MBI led to a considerable formation and sedimentation of Mn oxides in the water column of the EGB (Dellwig et al., 2018; Hermans et al., 2019b; van de Velde et al., 2020b). Enhanced transport of organic matter to the sediment via Mn oxides thus have led to a higher deposition and increased benthic organic matter degradation rates. According to the numbers presented by Cisternas-Novoa et al. (2019), however, the increased vertical POC flux after oxygenation would still have remained relatively low.

In contrast, lateral transport could have been a substantial source of POC and associated elements after the inflow. Measurements in the south-eastern part of the EGB in March 2015 showed increased turbidity close to the seafloor below 140 m depth, as a gravity current generated by the inflow moved along the slope (Schmale et al., 2016). As the inflow moved north, it could thus have enhanced the constantly ongoing lateral transport of relatively fresh organic matter from shallower to deeper areas (Nilsson et al., 2021). Furthermore, the speed of the gravity current generated by the inflow was higher than current speeds generally observed in the area during stagnation periods (Hagen and Feistel, 2004; Schmale et al., 2016). The higher current speed could have resulted in a more intense resuspension and of larger particles (Danielsson et al., 2007), including material that would otherwise have been buried. This mixing of old and new organic matter could additionally have induced a “priming” effect, where older material is broken down more efficiently in the presence of the fresher material (Bianchi, 2011; van Nugteren et al., 2009).

3.3.3 The effect of oxygen

Oxygen can directly contribute to an enhancement of the organic matter mineralisation rate. Increased mineralisation rates under oxic conditions compared to anoxic conditions have been suggested for refractory (or aged) and some types of fresh organic matter (Arndt et al., 2013; Bianchi et al., 2018; Hulthe et al., 1998). The reason for this is debated, but factors such as production of different enzymes and higher energy yield during degradation of organic matter under oxic conditions have been proposed (LaRowe et al., 2020). Although oxygen can affect the mineralisation of organic matter, the connection between oxygen and increased DSi fluxes is less clear. It is possible that oxygen promoted the degradation of the organic coating covering diatom frustules, which in turn could enhance the dissolution of biogenic silica (Bidle and Azam, 1999). The elevated DSi fluxes observed here could also indicate that an increased input of organic matter (with which biogenic silica is associated, such as diatom material), rather than the presence of oxygen, was causing the change in fluxes.

Interestingly, a study from the EGB conducted shortly after the MBI found that while the intrusion of oxygen in itself did not enhance the heterotrophic bacterial production in the water column, the removal of hydrogen sulfide did (Piontek et al., 2019). The same mechanism could be relevant for the surface layer of the sediment, since the MBI also resulted in the removal of hydrogen sulfide right below the sediment–water interface (Marzocchi et al., 2018).

3.4 Environmental implications

Following the MBI, there was a substantial elevation of the sedimentary DIP release in the EGB despite oxygenation of the sediment. Although Fe and Mn oxides were present at the sediment surface (Figs. 3a, 4a; Dellwig et al., 2018; Hermans et al., 2019b), they did not prevent DIP from escaping to the water column. The oxygen penetration was restricted to the top 2 mm of the sediment, and fluxes of hydrogen sulfide and other reduced compounds from deeper sediment layers appear to have prevented extensive oxidation of the sediment (Hermans et al., 2019b; Marzocchi et al., 2018). In the eastern part of the EGB, the ratio between Fe oxides and P_ex + P_Fe and was ∼ 1 (Hermans et al., 2019b), suggesting that the adsorption sites on the Fe oxides were saturated since the lowest ratio Fe oxide : P ratio that has been observed in marine environments is ∼ 2 (Gunnars et al., 2002). Other minerals, such as Mn oxides and carbonates, could also act as sorption sites for phosphate (Kraal et al., 2017; Yao and Millero, 1996). After the MBI, large quantities of Mn oxides were deposited on the sediment (Dellwig et al., 2018; Hermans et al., 2019b; van de Velde et al., 2020b). Hermans et al. (2019b) found clear associations between Mn and P enrichments in the sediment and stipulated that Mn minerals had an important function in the retention of phosphate, which could explain the relatively large amount of P_ex + P_Fe relative to Fe oxides. Nevertheless, the fact that DIP escaped from the sediment to the water suggests either that all adsorption sites were saturated or that the adsorption process was not efficient (e.g. due to kinetic limitations) relative to the diffusion of DIP toward the water column. Overall, it appears that the sediment of our study site had limited capacity to trap DIP during this short oxygenation event.

Figure 5 Upscaling of DIP release. (a) DIP release at station F according to three different scenarios: background, MBI min, and MBI max. The lines show the DIP flux according to the different scenarios; the shaded area represents the total DIP release during the period. (b) Sediment areas (western flank of the EGB) used for upscaling calculations. The DIP flux per square metre in the light grey area (150–190 m) is assumed to be the same as at station E, while the dark grey area (> 190 m) is represented by station F. Bathymetric data from the Baltic Sea Bathymetry Database (version 0.9.3); the map was created using the software ArcMap™ (v. 10.6) by Esri.

To estimate the potential impact on the water column P inventory in the EGB, we propose three scenarios for the total, time-integrated DIP release from the sediment between 1 March 2015 and 16 April 2018 (3.13 years, Fig. 5a). The background scenario describes a situation without the MBI, assuming that the flux equals the average of the fluxes in 2008, 2010, and 2018. The MBI min scenario describes a situation where the flux decreased to the value measured in 2015 as soon as oxygen reached the EGB, stayed at that level until oxygen disappeared and then changed linearly between the measured values. The MBI max scenario describes a situation where the flux begins at the background value and then changes linearly between the measured values. At station F, the background scenario gives a total DIP release of 0.91 mol m⁻², while the MBI min scenario results in a total release of 1.42 mol DIP m⁻² (56 % increase) and the MBI max scenario gives a total release of 1.80 mol DIP m⁻² (98 % increase). At station E, the total release of DIP is increased by 72 % and 112 % (Table E1 in Appendix E).

Upscaling of these results (Appendix E) further shows the importance of the increased DIP release. The calculation is limited to the western flank of the EGB, since previous studies have shown that the sediment–water fluxes are generally lower on the eastern side (Hall et al., 2017; Nilsson et al., 2019; Sommer et al., 2017; Viktorsson et al., 2013). We assume that the sedimentary DIP release per square metre at 150–190 m depth equalled that at station E, while the release at sediments deeper than 190 m were represented by station F (Fig. 5b). Since measurements are lacking from station D in 2015, only stations E and F are included in the upscaling. Since MBI had a relatively small effect at the fluxes at station D, however, these calculations should still give an acceptable estimate of the total DIP release. The MBI scenarios predict an extra release of 9.57–15.39 kton yr⁻¹. These numbers can be put in relation to the total release of DIP from oxygen-depleted sediments in the central Baltic Sea (Baltic Proper), which has been estimated to 60–150 kton yr⁻¹ (Almroth-Rosell et al., 2015; Hall et al., 2017; Sommer et al., 2017; Viktorsson et al., 2013). Even our conservative estimates suggest a substantial increase in sedimentary DIP release in the area following the MBI.

The sedimentary release of DIP led to a steady increase of the deep water DIP concentration in the years following the MBI (Fig. A8; SMHI, 2021). Yet the effect of the excess DIP release on an ecosystem scale depends on the cause of the elevated benthic fluxes. If the MBI enhanced the degradation of the organic matter already present in the sediment, the excess DIP release would increase the water column DIP pool. As discussed in Sect. 3.2, however, the sedimentary P_org pools were not large enough to support the in situ measured fluxes. Enhanced degradation of organic matter due to oxygenation is thus unlikely to have been the main cause of the elevated benthic fluxes. The data instead suggested that the primary driver behind the benthic fluxes was input of organic matter. If the MBI enhanced the organic matter input by accelerating vertical and horizontal shuttling, material would have been rearranged in the basin. The excess DIP release would then not have increased the DIP pool in the water column, since the same material would have been degraded elsewhere in the basin had the MBI not occurred. However, as discussed above, the events following the MBI (transient oxygenation, resuspension caused by the gravity current) could also stimulate the degradation of materials that was already present in the sediment. Hence, even if most of the excess DIP release was caused by accelerated shuttling, the MBI may have led to an increase in the total DIP release (and the water column DIP inventory) on a basin-wide scale over a 3-year period.

4 Conclusions

In this paper, we show that an oxygenated inflow to long-term anoxic sediments temporarily elevated organic matter degradation, thereby increasing the sedimentary release of DIP. Oxygen is generally believed to be able to break the positive feedback between anoxia and eutrophication (Ruttenberg, 2003; Stigebrandt, 2018), but that was not the case with the short and moderate oxygenation in this study. Instead, our results suggests that the inflow accelerated vertical and horizontal shuttling of particulate matter which, combined with potential stimulation of organic matter breakdown in the sediment, increased the sedimentary nutrient and DIC effluxes. The enhanced mineralisation of organic matter after oxygenated inflows observed in this study is a previously unrecognised feedback mechanism in eutrophic coastal systems (Breitburg et al., 2018; Middelburg and Levin, 2009).

In estuaries and coastal seas around the world, climate change and land use are predicted to both strengthen and expand oxygen-depleted zones (Altieri and Gedan, 2015). While oxygenated inflows have been thought to mitigate eutrophication feedbacks by enhancing P burial, our results show that they also can lead to increased mineralisation of organic matter. In long-term anoxic areas, irregular inflows might not bring enough oxygen to oxidise reduced chemical species and promote removal of nutrients. Instead, inflows may cause an internal redistribution of organic matter, leading to a temporal stimulation of organic matter degradation in the sediment. This stimulation of benthic nutrient recycling can lead to increased availability of P and other nutrients, thereby strengthening eutrophication.

Appendix A: Supplementary figures and tables

Figure A1 Oxygen monitoring data. Water column oxygen concentrations at depths corresponding to stations D–F at the SMHI (Swedish Meteorological and Hydrological Institute) monitoring station BY15 (SMHI, 2021). The grey lines indicate when samplings for this study were conducted; the dotted line shows the limit for hypoxia (63 µM). Colours on the lines show the bottom water concentration measured by the benthic chamber lander at each station. Following van de Velde et al. (2020b).

Figure A2 Sediment–water fluxes measured with the benthic chamber lander. Flux data are from before the inflow (2008, 2010; Hall et al., 2017; Nilsson et al., 2019; Viktorsson et al., 2013), right after the inflow (2015; Hall et al., 2017), and post-inflow (2016–2018; this study). There are no data from station D in 2015.

Figure A3 Pore water profiles of dissolved inorganic phosphorus (DIP). Profiles from 2008 are from Viktorsson et al. (2013). Different colours show replicate cores (same as Table C1): light blue markers show core 1, orange markers show core 2, dark blue markers show core 3, and red markers show core 4.

Figure A4 Solid-phase phosphorus speciation. P_ex stands for exchangeable or loosely sorbed P, P_Fe stands for easily reduced or reactive iron-bound P, P_aut stands for authigenic P, P_det stands for detrital P, and P_org stands for P associated with organic matter. The blue and orange markers represent replicate small sediment cores (inner diameter 6 cm).

Figure A5 Porosity profiles from stations D and F in 2016 and 2017.

Figure A6 Primary production and nutrient monitoring data. Data from SMHI station BY15 at top 10 m of the water column (SMHI, 2021). The graphs show (from a to d) phytoplankton biovolume, chlorophyll a concentrations, DIN (nitrate + nitrite + ammonium) concentrations, and dissolved inorganic phosphorus (DIP) concentrations. The grey lines indicate when samplings for this study were conducted.

Figure A7 Temperature in the surface layer (a) and deep layer (b) at the SMHI station BY15 (SMHI, 2021). The grey lines indicate when samplings for this study were conducted (2008, 2010: August/September; 2015: early July; 2016–2018: April).

Figure A8 Dissolved inorganic phosphorus concentrations in the deep layer of the SMHI station BY15 (SMHI, 2021). The grey lines indicate when samplings for this study were conducted (2008, 2010: August/September; 2015: early July; 2016–2018: April).

Table A1 Description of stations. Bottom water salinity, temperature, and oxygen concentration as measured by the benthic chamber lander or CTD. Limit of detection for oxygen was 0.5 µM.

Table A2 Results from the type III ANOVA of sediment–water solute fluxes. Degrees of freedom (df) values for DIN within are given brackets. p values < 0.05 marked in bold and italic font.

Table A3 Results from Student–Newman–Keuls (SNK) tests. For years or stations with different group letters, the average fluxes differed significantly. Fluxes were highest in groups named a, with decreasing fluxes in groups following alphabetical order.

Appendix B: Calculation of solid phosphorus inventories

The solid-phase inventory between the sediment surface and sediment depth z was calculated by multiplying the solid-phase density (σ_sed) with the depth-integrated porosity (ϕ) and concentration of the P fraction (C_P) of interest:

$$\begin{array}{}\text{(B1)}& \mathrm{Inventory}={\mathit{\sigma }}_{\mathrm{sed}}\underset{\mathrm{0}}{\overset{z_{\mathrm{end}}}{\int }}\left(\mathrm{1}-\mathit{\phi }\left(z\right)\right)C_{\mathrm{P}}\left(z\right)\mathrm{d}z.\end{array}$$

The value for σ_sed (1.9 g cm⁻³) was taken from van de Velde et al. (2020b). In a few cases (see Figs. 4a and 5a), P concentrations from certain sediment depths were missing. To not underestimate the sedimentary inventories, the missing values were replaced by average concentrations calculated from the sediment depths just above and below.

We were interested in the “excess” inventories of P_ex and P_Fe (Fig. B1), i.e. is the transiently increased P fraction. For P_ex, we consider the excess P inventory to be the inventory as calculated by Eq. (1), as the concentration of P_ex decreased to 0 below 5 cm in both 2016 and 2017. The P_Fe profile had a background concentration of ∼ 3 µmol g⁻¹ below 5 cm, which was consistent between sampling years. The concentration at the water interface in 2017 (i.e. the concentration in the settling particles that had not been exposed to the oxygen-rich water of the MBI) was also close to 3 µmol g⁻¹. Therefore, we determined the excess P_Fe concentration as the inventory as calculated by Eq. (1), subtracted by a background inventory calculated by replacing C_P in Eq. (1) with the average concentration below 5 cm.

Figure B1 Definition of “excess” phosphorus (P). Schematic figure showing the definition of excess phosphorus (P), which is the integrated enrichment of a P fraction compared to its background concentration.

Appendix C: Calculation of diffusive fluxes

Diffusive fluxes (J) were calculated from the pore water profiles of DIP (Fig. A3) according to Eq. (C1):

$$\begin{array}{}\text{(C1)}& J=-\mathit{\phi }{D}_{\mathrm{S}}{\displaystyle \frac{\partial C}{\partial z}},\end{array}$$

where D_S is the sediment diffusion coefficient and C is the concentration (Boudreau, 1997). D_S was calculated according to Eq. (C2):

$$\begin{array}{}\text{(C2)}& D_{\mathrm{S}}={\displaystyle \frac{D}{{\mathit{\theta }}^{\mathrm{2}}}}={\displaystyle \frac{D}{\mathrm{1}-\ln \left({\mathit{\phi }}^{\mathrm{2}}\right)}},\end{array}$$

where D is the diffusion coefficient in free solution and θ is the tortuosity (Boudreau, 1997). The R package “marelac” (Soetaert et al., 2020) was used to calculate D.

Table C1 Diffusive fluxes of dissolved inorganic phosphorus. The fluxes are calculated from pore water profiles, values are given in mmol m⁻² d⁻¹. The 2008 values are calculated from profiles from Viktorsson et al. (2013). When there are results from multiple cores, the core number is given in brackets and italic font. For profiles that displayed production peaks, both the upward flux (J_up) and the downward flux (J_down) from that peak are presented.

Appendix D: Pigment extractions

Pigment samples were collected in 2008, 2017, and 2018. Sediment cores were collected with a multiple corer and were sliced at 0.5 cm resolution from 0 to 2 cm depth. From each slice, 1 mL of sediment was transferred to 15 mL centrifuge tubes and was subsequently frozen in liquid nitrogen. The samples were extracted in 10 mL of an acetone : methanol (80 : 20) solution at −20 ^∘C for 24 h before they were sonicated for 60 s and filtered through 0.45 µm nylon filters. Extracts were separated and identified by high-performance liquid chromatography (HPCL, Shimadzu LCsolutions System) with photodiode array detection (SPD-M20A). The relative concentration of sample chlorophyll a was determined from calibration curves acquired from analysis of an external standard (DHI LAB).

Figure D1 Profiles of chlorophyll a (Chl a) in the top 2 cm of the sediment.

Appendix E: Upscaling calculations

Bathymetry data from the Baltic Sea Bathymetry Database (version 0.9.3) were used for the upscaling calculations. Sediment surface areas were extracted using the tool “Surface Volume” in the software ArcMap™ (v. 10.6) from Esri. The area at 150–190 m depth was represented by station E, while sediments deeper than 190 m were represented by station F. The extracted areas were multiplied with the total DIP release calculated for the background, MBI min, and MBI max scenarios and were then divided by the time used in the total release scenarios (3.13 years) to obtain an areal release per year.

Table E1 Total release of dissolved inorganic phosphorus (DIP). The total DIP release represents the entire time period from the arrival of oxygen until the end of the study (3.13 years). The total DIP release was not calculated for station D since no flux measurements were conducted at the site in 2015.