Iron-Phosphorus Feedbacks Drive Multidecadal Oscillations in Baltic Sea Hypoxia.

Hypoxia has occurred intermittently in the Baltic Sea since the establishment of brackish-water conditions at ∼8,000 years B.P., principally as recurrent hypoxic events during the Holocene Thermal Maximum (HTM) and the Medieval Climate Anomaly (MCA). Sedimentary phosphorus release has been implicated as a key driver of these events, but previous paleoenvironmental reconstructions have lacked the sampling resolution to investigate feedbacks in past iron-phosphorus cycling on short timescales. Here we employ Laser Ablation (LA)-ICP-MS scanning of sediment cores to generate ultra-high resolution geochemical records of past hypoxic events. We show that in-phase multidecadal oscillations in hypoxia intensity and iron-phosphorus cycling occurred throughout these events. Using a box model, we demonstrate that such oscillations were likely driven by instabilities in the dynamics of iron-phosphorus cycling under preindustrial phosphorus loads, and modulated by external climate forcing. Oscillatory behavior could complicate the recovery from hypoxia during future trajectories of external loading reductions.

Introduction

Marginal marine basins such as the Baltic Sea, Black Sea and many estuaries are naturally susceptible to hypoxia due to stratification and isolation of deep water masses (Conley et al., 2009; Eckert et al., 2013; Lefort et al., 2012). Hypoxia has been further exacerbated in modern coastal systems by anthropogenic nutrient inputs, leading to an expansion of marine ‘dead zones’ worldwide (Rabalais et al., 2010). Positive feedbacks in the phosphorus (P) cycle play a key role in intensifying hypoxic conditions in aquatic systems. Phosphorus release from sediments accelerates under low‐oxygen conditions, due to the dissolution of iron oxide‐bound P (Fe‐P; Mortimer, 1941) and the preferential regeneration of P from organic matter (Ingall et al., 1993). Enhanced sedimentary P release, in turn, sustains high productivity and oxygen demand (Vahtera et al., 2007; Van Cappellen & Ingall, 1994).

An important caveat to the positive feedback between P cycling and hypoxia is that in spatially finite systems, the release of Fe‐P from sediments cannot proceed indefinitely. Low‐oxygen conditions at any given location may rapidly exhaust the local sedimentary inventory of Fe‐P (Jilbert et al., 2011; Reed et al., 2011), while the maximum area of potential hypoxia may be limited by thermohaline stratification (Carstensen et al., 2014). Theoretically, if external P loads are insufficient to maintain high productivity after the exhaustion of the Fe‐P pool, a reverse of the positive feedback may be induced, leading to sequestration of P and abatement of hypoxia. This dual directionality of the Fe‐P feedback has been proposed to explain the rapid transitions at the onset and termination of centennial‐timescale hypoxic events in the Baltic Sea (Jilbert & Slomp, 2013).

The hypoxic events defined in Jilbert and Slomp (2013) are subdivisions of longer hypoxic intervals (HI) in the sedimentary record of the Baltic Sea (Zillén et al., 2008), which occurred during the Holocene Thermal Maximum (HTMHI, ∼8,000–4,000 cal. yr. B.P.), the Medieval Climate Anomaly (MCAHI, 1,200–750 cal. yr. B.P.) and the modern period (Mod.HI., ∼1980–present). Hypoxic conditions during these intervals are thought to have been induced by a combination of physical oceanographic forcing and external nutrient inputs, and sustained by P cycle feedbacks. Due to constraints on the sampling resolution of sediment cores, investigation of multiannual to multidecadal variability in past hypoxia during these intervals is challenging. However, such variability is essential to study because of the relevance for managing the recovery from modern hypoxia, one major goal of the Baltic Sea Action Plan (HELCOM, 2007, 2021). Here, we present new, ultra‐high resolution records of past Baltic Sea hypoxia and Fe‐P dynamics using LA‐ICP‐MS line scanning of resin‐embedded sediments. These records cover the entire HTMHI and MCAHI intervals and thus significantly extend the initial LA‐ICP‐MS data presented in Jilbert and Slomp (2013). Using time‐dependent box‐model, we show that multidecadal oscillations observed in the sediment proxies could plausibly have been sustained by instabilities in the feedbacks between Fe‐P dynamics and hypoxia on these timescales.

Materials and Methods

Sediment Sampling

A gravity core (0–406 cm) was collected from site F80 (58°00.00°N, 19°53.81°E, water depth 191 m) in the Fårö Deep of the Baltic Sea during the HYPER/COMBINE cruise of R/V Aranda (May/June 2009), and sectioned in a nitrogen‐filled glovebox at 0.5–2.0 cm intervals. Undisturbed U‐channel subsamples were taken from laminated intervals and embedded with Spurr's epoxy resin by the fluid displacive method under nitrogen (Jilbert et al., 2008).

Sample Preparation and Geochemical Analysis

Discrete sediment samples were freeze‐dried and ground under nitrogen. Aliquots of 0.125 g dry sediment were digested in 2.5 ml HF (40%) and 2.5 ml of HClO4 (70%)/HNO3 (65%) mixture at 90°C overnight, before evaporation at 160°C and redissolution in 1M HNO3. Solutions were analyzed by ICP‐OES for total Al, Fe and Mo and P (precision and accuracy <10%).

Resin‐embedded sediment trays were cut to expose the internal surface of the sediment, and sectioned into shorter blocks to fit inside the Laser Ablation (LA) ICP‐MS sampling chamber. LA‐ICP‐MS line scanning was performed at 0.04 mm s−1 using an Excimer laser (193 nm, spot size 120 μm, repetition rate 10 Hz, fluence 8 J cm−2) coupled to an Element 2 ICP‐MS (Hennekam et al., 2015). Count intensities of the following isotopes were measured for use in this study: 27Al, 31P, 57Fe, 79Br, and 98Mo. Calibration of the LA‐ICP‐MS data was performed according to a modified version of the two‐step procedure described in the supplement of Jilbert and Slomp (2013). Full details of the calibration procedure are given in the Supporting Information S1.

Rationale of the Geochemical Proxies

Sedimentary Mo/Al and Fe/Al are established proxies for hypoxia intensity and Fe shuttling, respectively (Jilbert & Slomp, 2013; Lenz et al., 2015). The ratio of bromine to phosphorus (Br/P) is here used as a surrogate for C/P, an indicator of the intensity of preferential P regeneration (Jilbert & Slomp, 2013). The rationale for using Br/P is that C cannot be measured in the epoxy‐embedded samples by LA‐ICP‐MS, but Br and C are strongly correlated in sediments due to the bromination of organic matter by seawater bromide during degradation (Ziegler et al., 2008), hence LA‐ICP‐MS‐derived Br provides a surrogate for C. Absolute Br contents were also determined on selected samples by repeated digestion in H2O2 followed by ICP‐MS analysis (Ziegler et al., 2008). These values were used to generate discrete sample Br/Al ratios for the final step of the calibration of the LA‐ICP‐MS data.

Sediment Dating and Time Series Analysis

The age‐depth model used for the sediments in this study is described in the supplement of Jilbert and Slomp (2013). Briefly, the discrete‐sample C org profile was tuned to the Loss on Ignition (LOI) profile of reference core 372740‐3 from the Gotland Deep (Lougheed et al., 2012). 372740‐3 is dated using known paleomagnetic secular variation (PSV) features and lead (Pb) isochrons (12 features total in the 406 cm core covering ∼7,000 years, giving an average spacing of ∼600 years). Linear sedimentation rate is assumed between dating points. In the early part of the record, this age‐depth model differs by up to 500 years from a recently published model for Baltic Sea sediments based on benthic foraminiferal 14C dating (Warden et al., 2017). However, this offset is not expected to impact on the conclusions of the present study, which focuses on frequencies of variability rather than absolute ages of events.

Time series analysis of the LA‐ICP‐MS Mo/Al, Fe/Al, and Br/P profiles was performed in Analyseries 1.1.1 (Paillard et al., 1996). Raw elemental ratios were detrended and normalized to unit variance, resampled to 1 year resolution and processed by Fast Fourier Transform (FFT, Blackman‐Tukey method with Bartlett window) to yield a power spectrum of cyclic components in each profile. Cross‐spectral analysis (same settings as above) was also applied to investigate coherence and phase relations between the elemental ratio data. A 20–100‐year bandpass filter was subsequently applied to the detrended, normalized data for visual comparison of the time series.

Box Modeling

A box model was created to represent the coupled cycling of P and Fe in the Baltic Sea. The basic formulation is the same for the HTMHI, MCAHI, and Mod.HI (Figure 1a), with six state variables related by a system of equations. A full model description, including calibration and validation against present day data, is provided in the Supporting Information S1. Production of organic matter in surface waters is limited by the phosphorus concentration P. Organic matter sinks vertically and accumulates in the sediments, where remineralization occurs. Phosphorus may be regenerated from sediments due to remineralization or release from sedimentary Fe‐P, or buried. Oxygen is constant in the surface layer but allowed to vary in the deep‐water box, simulating the observed oxygen depletion in the sub‐halocline water column of the Baltic Sea (Carstensen et al., 2014). A key feature of the model is a sigmoid function describing the relationship between the sedimentary Fe‐P inventory and the oxygen supply‐demand ratio, the dimensionless ratio between oxygen flux to sediments and oxygen demand by organic matter mineralization (Figure 1b). The oxygen flux to sediments is directly proportional to oxygen concentration. Advective exchange of dissolved constituents is permitted, both vertically between the two layers and laterally between these layers and a hypothetical adjacent sea (simulating the real exchange between the Baltic and the North Sea). External loading of P is prescribed for each time interval based on literature estimates. The sensitivity and stability of steady state solutions to changes in external loads and vertical exchange were assessed. Simulations with variable frequency and amplitude of vertical mixing were conducted to demonstrate the potential of variable climate forcing to trigger and control oscillatory behavior. Spectral analysis was performed on 10,000‐year long annual oxygen time‐series using FFT with Hamming window of the autocovariance limited to 512 time steps.

Figure 1

Key features of the box model used in this study. (a): Schematic of the model design. Surface water and deep water boxes (blue) are each underlain by a reactive sediment layer (gray). Six time‐dependent state variables are included: P = water column phosphorus concentration, O2 = water column oxygen concentration (deep water only), C org = sediment organic carbon inventory, P org = sediment organic phosphorus inventory, Fe‐P = sediment iron‐bound phosphorus inventory. Equations governing exchange fluxes and their relationship to state variables are given in the Supporting Information (Text S3, Figure S4). (b): Principle of the sigmoid function used in the model. The attainable Fe‐P inventory of the sediments is controlled by the dimensionless ratio between oxygen flux to sediments and sedimentation rate of organic matter (O).

Results and Discussion

Evidence for Oscillations in Past Hypoxia and Phosphorus Regeneration

During HTMHI and MCAHI, multidecadal oscillations are observed in the LA‐ICP‐MS line scan profiles of Mo/Al and Br/P (Figure 2, left). Spectral analysis of the full data series shows multiple peaks across the decadal to multidecadal bands. A clear peak is observed in the 60–100‐year band in both parameters (Figure 2, center). This band is also characterized by high coherence and close‐to‐zero phase relation between Mo/Al and Br/P (Figure 2, right), indicating in‐phase oscillatory signals. Especially in the intervals of lower overall Mo/Al (e.g., 1.6–1.2 ky BP; 5.4–5.0 ky BP), the oscillation is clearly visible in the raw and bandpass‐filtered data (Figure 2, left and Figures S2 and S3). This result implies that a close coupling existed between hypoxia and phosphorus regeneration on multidecadal timescales during past hypoxic events, and that both were highly variable on these timescales.

Figure 2

LA‐ICP‐MS line scan data of Baltic Sea sediments from (a) HIHTM and (b) HIMCA. (left) High‐resolution geochemical profiles of sediments from site F80. Horizontal gray bars indicate the subdivision of HTMHI and MCAHI into numbered hypoxic events, as given in Jilbert and Slomp (2013). Upper panels = calibrated LA‐ICP‐MS profiles of Mo/Al. Lower panels = 20–100 years bandpass‐filtered profiles of Mo/Al, Br/P and Fe/Al (detrended and normalized to unit variance prior to filtering). (center) Blackman‐Tukey spectral analysis of Mo/Al, Br/P, and Fe/Al data for the entire time intervals shown on the left. Gray field indicates the period 60–100 years. (right) Coherence and phase analysis of Mo/Al and Br/P (solid lines) and Mo/Al and Fe/Al (dashed lines) for the entire time intervals shown on the left. Phase in radians (0 = in phase, П = antiphase).

Evidence for Oscillations in Past Fe Shuttling

In the intervals where the oscillations in Mo/Al and Br/P are most pronounced (e.g., 1.6–1.2 ky BP; 5.3–4.9 ky BP), the LA‐ICP‐MS data also show similarly paced variability in Fe/Al (Figure 2, left and Figures S2 and S3). Accordingly, a similar peak in the 60–100‐year band is observed in the power spectrum of the Fe/Al data (Figure 2, center), as well as high‐coherence and close‐to‐zero phase relation between Mo/Al and Fe/Al (Figure 2, right). This indicates that shelf‐to‐basin shuttling of Fe in the Baltic Sea was also sensitive to multidecadal variability in hypoxia. Namely, during periods of low deep‐water oxygen, more Fe was transported laterally downslope into the deep basins via cycles of dissolution and reprecipitation (Lyons & Severmann, 2006).

Oscillations in Box Model Simulations

Model simulations suggest that multidecadal oscillations in hypoxia and phosphorus regeneration may have been an intrinsic feature of biogeochemical cycles in the Baltic Sea under the forcing conditions of the HTMHI and MCAHI (Figure 3). Steady‐state solutions of the model are unstable when parameters including basin geometry and external loading of P are set to constant realistic values for these intervals, indicating the presence of unforced oscillations. The periodicity of these unforced oscillations in the model is typically 130–170 years, slightly longer than observed in the sediment records (Figure 3). Within each oscillation, during the period of low deep water oxygen, the sediment Fe‐P inventory is at a minimum, whereas deep water phosphorus (P), and sediment organic carbon (C org) and phosphorus (P org), show maximum values (Figure 3). Conversely, when deep water oxygen is high, the opposite trends are observed.

Figure 3

State variables in box model simulations of biogeochemical cycles during past hypoxic intervals in the Baltic Sea. (a and b): sedimentary P inventories, (c and d): oxygen concentration in deep water box, (e): sensitivity analysis showing mean values, at steady state, of state variables in response to changes in external P loading. (f): sensitivity analysis showing mean values, at steady state, of state variables in response to changes in vertical exchange between shallow and deep‐water boxes. Blue shaded areas in (e‒f) indicate the range in P input and vertical exchange, where the steady state solutions are unstable, according to linear stability analysis, causing oscillations in time‐dependent solutions.

Mechanism and Frequency of the Oscillations

The presence of relatively large amounts of Fe‐P in Baltic Sea sediments under low P loading conditions is a prerequisite for the observed instability. The quantitative representation of the sigmoid function used in the model (Figure S5) shows that the sensitivity of the Fe‐P inventory to the oxygen supply‐demand ratio is high at intermediate values (Fe‐P = 50–150 mmol m−2). During the recharge phase (Figure 1b), storage of Fe‐P in sediments provides a reinforcing feedback toward a higher supply‐demand ratio. Conversely, the discharge phase is characterized by release of Fe‐P and a reinforcing feedback toward lower supply‐demand ratio. However when oxygen is plentiful and productivity low (high values of the supply‐demand ratio), the sediment Fe‐P pool becomes increasingly saturated (close to 200 mmol m−2 Fe‐P), leading to a leveling‐off in the Fe‐P inventory. Similarly, when oxygen is scarce and productivity high, the Fe‐P concentration levels off due to the approaching exhaustion of the sedimentary Fe‐P pool. This insensitivity at high and low values makes the system vulnerable to a switch in directionality.

Even in our relatively simple model formulation, the frequency of the intrinsic oscillations is a complex function of the prescribed parameters. Frequency is influenced partly by the geometry of the Baltic Sea basins, which is a major factor controlling the residence time of P in both the shallow and deep layers in the model. Indeed, the periodicity of simulated oscillations during the HTMHI, when the subhalocline hypoxic area of the Baltic was expanded due to glacio‐isostatic effects (Jilbert et al., 2015) is slightly longer than that during the MCAHI (Figures 3a‒3d). However, the frequency of real‐world oscillations may also be externally driven. Indeed, the inherent instability of the model can be triggered by variability in external driving forces as discussed below.

Influence of Climatic Drivers on Oscillations

The observed oscillations in the sediment records indicate dominant periods of 60–100 years, which is consistent with periods of the Atlantic Multidecadal Oscillation (AMO), or Atlantic Multidecadal Variability (AMV), that may have been present throughout the Holocene (Knudsen et al., 2011). The AMO in turn modulates the influence of the North Atlantic Oscillation (NAO) on the conditions in the Baltic Sea region (Börgel et al., 2020). Thus, it is reasonable to assume that the observed oscillations in Baltic Sea sediment records may to some extent have been externally driven by climatic variability as well as the internal instability of biogeochemical cycles.

To explore the impact of climatic variability, we imposed a sinusoidal perturbation to the vertical exchange relative to the parameter settings for the model simulations of the MCAHI. Vertical exchange is chosen because it has a direct and influential impact on oxygen conditions and subsequent effects on the Fe‐P inventory. In Figure 4, oxygen power spectra are shown for a number of simulations with different amplitudes and periods of the vertical exchange perturbation. It is evident that for certain periods of climatic oscillations, the inherent model instability is triggered and the oxygen oscillations are strongly amplified compared to the unforced case (e.g., for amplitude = 3 m y−1 and T = 80 and 100 years). In such cases, the period of oscillations is largely determined by the climate. This external influence could partly explain the generally shorter periods observed in the sediment records (60–100 years, similar to AMO/AMV forcing) compared to the inherent oscillations of the model (130–170 years). In terms of actual historical variability in AMO/AMV forcing, we note that Wang et al. (2017) showed the magnitude of AMV in the 60–100‐year band to be higher in the centuries following the MCAHI than during it. However, since the Baltic Sea was in a stable oxic state after 0.8 ky BP (Jilbert & Slomp, 2013) there was no amplification of the climatic signal to be observed in the sediment records.

Figure 4

Power spectra of deep‐water oxygen concentrations for different cases of imposed cyclic variation of the vertical exchange rate (amplitude in m y−1 as deviation from 11 m y−1, the MCAHI default value), representing variability of climatic forcing. Each row represents a given period (T) and the columns different amplitude of the oscillation. The dashed line graphs show the power spectra of the inherent oscillation with constant vertical exchange rate (cf. Figure 3c).

Sensitivity to External Phosphorus Loading and Vertical Exchange

The presence of multidecadal oscillations in the model simulations is highly sensitive to external loading of P. Using a parameterization for the modern Baltic Sea, external loading of <8,000 tonnes P/year yields solutions in which no oscillations are observed because the system remains in a quasipermanent oxic state (Figure 3e). As external loading increases, oscillations are observed within a range of ∼8,000–12,000 tonnes P/year. However at high external loading (>12,000 tonnes P/year), oscillations disappear again because the system remains in a quasipermanent hypoxic state.

The modern P loading to the Baltic Proper is ∼25,000 tonnes P/year, and loading during the late 20th century was significantly in excess of this value (Gustafsson et al., 2012). Such high loads are sufficient to suppress oscillatory behavior in the model, by sustaining high productivity and holding the system in a quasipermanent hypoxic state. A simulation of the modern development of hypoxia since 1900 supports this observation (Figures S6 and S7). Therefore, under constant future loading scenarios, we would not expect to observe multidecadal oscillations in Baltic Sea hypoxia during the coming centuries.

Although enhanced nutrient loading to the Baltic during the MCA has been suggested (Zillén & Conley, 2010), little evidence for eutrophication has been found in coastal paleoenvironmental data from this period (Jokinen et al., 2018; Ning et al., 2016; Norbäck Ivarsson et al., 2019). Therefore, we imposed an external loading of 8,000 tonnes P/year for both the MCAHI and HTMHI simulations, as derived from estimates of the preindustrial loading to the Baltic Sea (Gustafsson et al., 2012). Because this value is at the lower limit of the range yielding oscillations in most parameterizations (Figures 3c and 3d), we adjusted the vertical exchange parameter D in MCAHI and HTMHI simulations to ∼20% lower than the modern value (Table S1). Lower values of D favor stagnation and the development of deep water hypoxia, and tip the system into an oscillatory state at external loading of <8,000 tonnes P/year (Figures S8 and S9). Our simulations therefore suggest that climatic factors controlling vertical mixing may be capable of triggering hypoxia under preanthropogenic P loading, despite a recent study indicating otherwise (Meier et al., 2019).

Implications for Future Management of the Baltic Sea

Under the BSAP (HELCOM, 2007, 2021), riparian countries agreed to limit external P loading to the Baltic in order to reduce the occurrence of hypoxia and cyanobacterial blooms in the future. As part of the plan, maximum allowable inputs of P to the whole Baltic Sea were set at 21,716 tonnes P/year (about 18,500 tonnes P/year in the region covered by the model). Our findings suggest that such targets, if met, will improve water quality but are unlikely to force the system to a fully oxic state (Figure 3e). Rather, deep water oxygen concentrations are likely to remain close to the hypoxia threshold of 63 μmol/L. This finding is consistent with existing simulations from more complex coupled physical‐biogeochemical models for the Baltic Sea (e.g., Meier et al., 2018). To achieve a transition beyond hypoxia to quasipermanent oxic conditions would require a further load reduction and our model shows that this will likely proceed via an interval during which the system is vulnerable to multidecadal oscillations in hypoxia, creating a challenge for assessments of the success of management efforts. These results highlight the need to include more details of Fe‐P interactions in models used for management decisions in the Baltic Sea (e.g., BALTSEM; Savchuk et al., 2012). We suggest that the extent to which oscillations may impact upon the actual recovery trajectory in the coming centuries will depend on the speed with which loading reduction targets are achieved. More rapid reductions are more likely to starve the system of P and therefore to promote a more linear recovery trajectory.

Conclusions

Our study shows that feedbacks in the coupled cycling of Fe and P drove multidecadal oscillations in the intensity of oxygen depletion during past hypoxic intervals in the 8,000 years history of the Baltic Sea. In‐phase oscillatory profiles of Mo/Al, Br/P and Fe/Al in sediments suggest large‐scale synchronicity in deep water hypoxia, P regeneration and Fe shuttling, respectively. A box model of coupled Fe and P cycling suggests that internal multidecadal instability in hypoxia intensity is controlled by nonlinear recharging and discharging of the sedimentary Fe‐P reservoir in deeper areas. Crucially, our model simulations show that oscillatory behavior is sensitive to external P loading. While low external loads favor stable oxic conditions, and high loads sustain stable hypoxia, intermediate loads lead to instability causing oscillatory behavior on multidecadal timescales. This observation has implications for predicting the trajectory of recovery from modern hypoxia in the Baltic Sea, since external P loading is expected to decline in accordance with the BSAP. Variable external climate forcing may influence the frequency and amplitude of oscillations through its impact on vertical water mass exchange.

Supporting Information S1

Click here for additional data file.