Dynamics of oxygen sources and sinks in the Baltic Sea under different nutrient inputs

The Baltic Sea is one of the marine systems suffering from pronounced man-made hypoxia due to the elevated nutrient loads from land. To mitigate hypoxia expansion and to return the Baltic Sea to a good environmental state, the Baltic Sea Action Plan (BSAP), regulating the waterborne and airborne nutrient input, was adopted by all states surrounding the Baltic Sea. However, at the moment, no signi ﬁ cant shrinking of the hypoxic area is observed. In this study, two scenario simulations of the future state of the deep parts of the central Baltic Sea (deeper than 70 meters) were carried out, utilizing a 3-dimensional numerical model. Climate change effects on meteorology, hydrology, and oceanic state were not included. We focused on O 2 and H 2 S sources and sinks under different nutrient input scenarios. We found that under the BSAP scenario, all subbasins in the central Baltic Sea, especially the northern and western Gotland Basin, show signi ﬁ cant improvement, namely, oxygenation and oxidation of the deposited reduced material, ceasing its advection to the upper layers and neighboring basins. We found that the nutrient loads are responsible for more than 60% and 80% of the O 2 and H 2 S sources and sinks variability, respectively, at the interannual time scale. We showed that the Baltic Sea could return to the initial state in 1948, but under the more rigorous 0.5 BSAP scenario (nutrient input is halved compared to the BSAP). However, since we observed no hysteresis effect, the system would probably reach the initial state but over a timeframe longer than the 71-year future simulation period.


Introduction
Hypoxia, or dissolved oxygen concentrations in the water column below a certain threshold (usually 2 ml O 2 /l) (Conley et al., 2009;Savchuk, 2018;Stoicescu et al., 2019), is a pronounced problem in numerous marine systems worldwide, which is currently deteriorating (Diaz and Rosenberg, 2008;Breitburg et al., 2018).The main factor favoring the low oxygen concentrations is elevated anthropogenic nutrient loads from land promoting eutrophication (Heathwaite et al., 1996;Knuuttila et al., 2017;Ator et al., 2020).However, climate change also affects global deoxygenation (Whitney, 2022).One example of a marine system suffering from elevated hypoxia is the Baltic Sea -a semienclosed sea in Northern Europe.Due to some natural properties, such as limited water exchange with the North Sea, which lead to a long residence time (Leppäranta and Myrberg, 2009) and a pronounced permanent halocline, located around 60 meters depth (Väli et al., 2013;Uurasjärvi et al., 2021), which limits the exchange between the upper and lower layers, the Baltic Sea is naturally prone to hypoxia (Gustafsson et al., 2012).Lenz et al. (2015) studied the sediment cores from the Baltic Sea and concluded that hypoxic conditions occurred in the Baltic Sea after the establishment of the halocline.However, hypoxic area in the Baltic Sea has been dramatically increasing since the 1960s (Almroth-Rosell et al., 2021;Kõuts et al., 2021;Krapf et al., 2022).It was attributed to the elevated nutrient (P and N) loads from land (Larsson et al., 1985).To mitigate the ongoing eutrophication, the Baltic Sea Action Plan (BSAP) was developed by the Helsinki Commission (HELCOM) in 2007 (HELCOM, 2007;Backer et al., 2010).Later, the BSAP was updated a few times, with the last version applied in 2021 (HELCOM, 2021).One function of the BSAP is to provide information about the maximum allowable nutrient input (MAI) to the Baltic Sea.Adherence to the MAI should guarantee hypoxic area reduction and transition to a better state of the Baltic Sea.However, despite the nutrient loads reduction policy, no significant improvement has been observed yet (Hansson and Viktorsson, 2020).Vahtera et al. (2007) coined the term "vicious circle" for the Baltic Sea, which describes the possible damping mechanism for oxygenation.A model study (Neumann et al., 2002) supported the "vicious circle" mechanism, namely, the stable cyanobacteria biomass supported by the release of the sedimentary phosphorus under anoxic conditions.The "vicious circle" mechanism was further discussed and developed in several studies, for example (Rydin et al., 2017;Meier et al., 2018;Savchuk, 2018).Despite the overall awareness of the oxygen dynamics in the Baltic Sea, only a few studies disentangled the oxygen sources and sinks in the Baltic Sea (Gustafsson and Stigebrandt, 2007;Schneider et al., 2010;Naumov et al., 2023).Sources and sinks of oxygen and hydrogen sulfide in the central Baltic Sea demonstrated substantial changes during the last 70 years.Meier et al. (2018) highlighted the switch between oxygen consumption in sediments to more water column consumption.Naumov et al. (2023) conducted a trend analysis of oxygen and hydrogen sulfide sources and sinks and came to similar conclusions for oxygen.As for hydrogen sulfide, they found an increasing spread between its production in the sediments and consumption in the water column leading to elevated advection to the upper layers.There are also studies projecting oxygen concentrations in the future (e.g., Meier et al., 2011;Friedland et al., 2012;Neumann et al., 2012;Meier et al., 2022).Those projections include both climate and nutrient forcings.Our aim is to study only the nutrient forcing, focusing on the oxygen and hydrogen sulfide sources and sinks under changing nutrient input.It will help to disentangle the two forcing factors and to understand the possible future state of the Baltic Sea and might be useful in future adjustments of the BSAP.

Model description
To investigate oxygen and hydrogen sulfide dynamics in the central Baltic Sea, we used the 3-dimensional coupled regional MOM-ERGOM model.Modular Ocean Model (MOM) version 5 (Griffies, 2012) served as a hydrodynamical model reproducing ocean motion and dynamics of the two principal tracers (temperature and salinity) by solving the set of primitive equations.K-profile parametrization (KPP, Large et al., 1994) was used as the turbulence closure scheme.Our MOM setup utilizes regular orthogonal Arakawa B grid with z* coordinate vertical scheme and a predictor-corrector scheme (Griffies, 2012).Ecological Regional Ocean Model (ERGOM) (Radtke et al., 2019;Neumann et al., 2021;Neumann et al., 2022) was used as a biogeochemical model.ERGOM reproduces cycles of the main nutrients (C, N, P) as well as oxygen (O 2 ) and hydrogen sulfide (H 2 S), which is represented by the separate active tracer in the model.The complete model description can be found in (Neumann et al., 2022).An open boundary was placed in Skagerrak permitting exchange with the North Sea.Our model setup has three nautical miles horizontal resolution, while the vertical resolution varies from 0.5 to 2 meters.This setup has been used multiple times to model the Baltic Sea dynamics.It includes a recent study by Naumov et al. (2023) and a few others (e.g., Neumann, 2010;Voss et al., 2011;Kuznetsov and Neumann, 2013).Recently the model was thoroughly validated (see Naumov et al., 2023).As a short validation summary, the model reasonably reproduced the central Baltic Sea dynamics.However, possibly due to the overestimated exchange with the North Sea and, therefore, too strong halocline, hypoxic and anoxic areas were larger in the model compared to the reanalysis data.

Formulation of O 2 and H 2 S budgets
A budget concept implies the balance between sources and sinks of a given substance in a certain box.If the sources exceed the sinks, the total amount of an arbitrary substance within the box increases, and vice versa (Yurkovskis et al., 1993;Fennel and Testa, 2019).Mathematically it could be summarized by the following equation: Where Tr stands for any tracer.The first term in the equation represents the total change of the amount of tracer in time, which, by definition, equals the internal change of the tracer within the box, for instance, due to the biochemical processes (the third term) plus the total flux across all boundaries of the box (the second term).For oxygen, the budget consists of advective and diffusive supply across the box's boundaries, as well as the supply due to the photosynthesis, which is situated in the photic level during the vegetation period, and biochemical consumption in the water column and sediments expressed as higher trophic level organisms' respiration, mineralization of organic matter, nitrification, and oxidation of H 2 S. For hydrogen sulfide, the budget is formulated as advection and diffusion fluxes across the boundaries and mineralization of organic matter (sulfide reduction) in the sediments and water column.Hydrogen sulfide is consumed via the oxidation by O 2 or NO − 3 .For the description of all processes contributing to the O 2 and H 2 S budgets in ERGOM and their aggregations, see Supplementary Table 1.Four boxes situated in the central Baltic Sea, each representing a specific subbasin, were investigated: the Bornholm Basin (BB), the eastern Gotland Basin (eGB), the northern Gotland Basin (nGB), and the western Gotland Basin (wGB).The upper boundary for each box was set to a 70meter depth, excluding the upper oxygen-reach layer.The boxes spanned the water column down to the bottom.Their locations can be found in Figure 1.

Nutrient input scenarios
This study features three different nutrient input scenarios.The reference scenario is based on the simulation analyzed in Naumov et al. (2023).It employs HELCOM's actual nutrient loads from 1948 to 2018 (Svendsen and Gustafsson, 2020).Data gaps were filled with linear interpolation.The BSAP scenario imposes constant nutrient loads based on the level of the Maximum Allowable Input (MAI) (HELCOM, 2021).The total input of P and N into the Baltic Sea at levels no more than the BSAP MAI should guarantee the Baltic Sea's transition into a "good environmental status" (Borja et al., 2015).Since BSAP specifies only total loads without separating airborne and waterborne loads, the climatology seasonal cycle of the atmospheric nutrient input was applied.The halved BSAP MAI scenario (0.5 BSAP) assumes that the loads would be constant on the level of the half from the BSAP MAI (both waterborne and airborne input).The two reduction scenarios span the timeframe from 2019 to 2089 (71 years) and start from initial conditions taken from the last year of the reference scenario.Total nitrogen loads to the Baltic Sea under the BSAP scenario are set to 792.2 Kton/a (97% of the actual nitrogen loads averaged for the last ten years -815.4Kton/a).Under the 0.5 BSAP scenario, the total nitrogen loads equal 396.1 Kton/a (48% of the average actual loads for the last ten years).Total phosphorus loads under the BSAP and 0.5 BSAP scenarios equal 21.72 and 10.86 Kton/a, respectively.This constitutes 82 and 41% of the current total phosphorus loads -26.33 Kton/a, respectively.So the BSAP goals have not been fully achieved yet, but the actual loads are very close to them.For more information, see Supplementary Figures 1, 2. In reference and nutrient reduction scenarios, we applied identical atmospheric forcing, namely CoastDat2 atmospheric fields from 1948 to 2018 (Geyer, 2014).Therefore, we neglect the impact of future climate change and natural variability, focusing only on the nutrient loads effect on the Baltic Sea's eutrophication.

Budgets' composition
To evaluate different sources and sinks' contribution to the oxygen and hydrogen sulfide variability, we employed a linear regression framework proposed by Naumov et al. (2023).This linear model allows us to estimate the explained interannual variability by each group of processes.We used the same groups (phy, bio, and sed) as in the previous section.The results are shown in Figure 3.It can be seen that the general pattern is similar compared to the one by Naumov et al. (2023), which is advection dominance for oxygen in both considered scenarios.Another pattern for oxygen is less explained variability by the water column processes (especially in the BB and eGB) moving from the reference scenario to the 0.5 BSAP.This pattern complements the conclusions from the previous section, namely, the elevated variability of the sedimentary processes towards the end of the simulation.The composition of the processes contributing to the H 2 S budget is generally more complex (both regionally and in different scenarios).It points to the strong connection between H 2 S dynamics and nutrient forcing (see Supplementary Figure 4 for additional information).Overall, it can be stated that H 2 S dynamics is more affected by the nutrient loads reduction than oxygen dynamics, with some regional differences existing in both cases.

O 2 and H 2 S dynamics induced by the nutrient forcing
Nutrient loads dynamics determine a significant fraction of oxygen variability in the Baltic Sea (Meier et al., 2019).Here, we quantify its contribution to the oxygen and hydrogen sulfide sources and sinks variability under the considered scenarios.We employed Empirical Orthogonal Functions (EOFs), based on eigenvectors and eigenvalues (Hannachi et al., 2007).It allows a decomposition of the complex multidimensional data into a set of orthogonal functions with a reduced number of dimensions by calculating eigenvalues and eigenvectors of the data's correlation/ covariation matrix and multiplying the latter with the original matrix.We applied such an EOF analysis to the matrixes of oxygen and hydrogen sulfide budget terms' annual consumption/ production anomalies during the reference period .Then we projected the resulting EOFs into future projections.Resulting EOFs are guaranteed to represent the same process during the reference period and future scenarios.The results are presented in Figure 4 (for oxygen) and Supplementary Figure 4 (for hydrogen sulfide).Based on the eigenvalues, the first three EOFs were considered significant for both elements (Panel A in Figure 4; Supplementary Figure 4).For O 2 , the first EOF explains around 60% of the total variability, the second EOF about 20%, and the third about 10%.For H 2 S, the eigenvalues converge quicker, with the first EOF already explaining 85% of the total variability and the last two EOFs around 5% each.The leading EOF in both O 2 and H 2 S data can be attributed to the same process based on its loadings (Panel E in both figures).Both for O 2 and H 2 S data, the leading EOF demonstrates a strong positive connection with anomalies of sedimentary detritus' oxidation (either by oxygen or sulfate reduction).This means that when the score of the first EOF is positive (from the 1970s to the 2030s, according to Panel B in both figures), there is an intensified H 2 S production due to the Location of the four studied sub-basins (boxes).Note that the upper boundary is located at 70 meters depth, so the depth counting starts at 70 meters.
mineralization of sedimentary detritus in the nGB.The same is true (but to a lesser extent) for the other regions of the Gotland Basin.At the same time, detritus is less mineralized by oxygen in the nGB and wGB but more in the eGB, indicating reduced material deposition in the nGB and wGB.The negative score of the first EOF, which is observed at the beginning of the reference period and after the 2030s, indicates the overall improvement of the oxygen conditions in the central Baltic Sea (no long-term deposition of the reduced material and oxidation of existing H 2 S, especially in the remote basins).This allowed us to attribute the leading EOF to the nutrient forcing.The last two EOFs were attributed to the inflow activity and natural variability.They do not exhibit any significant difference between the reference period and the future projections (for both BSAP and 0.5 BSAP scenarios).The first EOFs for BSAP and 0.5 BSAP are highly correlated (>0.9), which suggests that there are no significant differences in how the nutrient load reduction affects the marine system for the BSAP and more rigorous 0.5 BSAP scenarios; however, the first EOF for 0.5 BSAP scenario shows a more negative score compared to the first EOF of the BSAP scenario, which means the more resilient state and more improvement.

Discussion
Our results suggest that oxygen conditions in the central Baltic Sea will substantially improve both under BSAP and 0.5 BSAP scenarios, especially in the remote nGB and wGB, which, however, did not significantly diminish the hypoxic area, only the anoxic area significantly shrank at the end of the study period (see Supplementary Figure 5), which is related to the significant oxygen debt in the deep central Baltic Sea (Rolff et al., 2022).However, our study uses atmospheric forcing from 1948-2018, which completely ignores possible changes, such as changes in stratification and internal nutrient cycle (Meier et al., 2011;Hordoir and Meier, 2012), that could worsen the results to some extent, but, based on the results by Saraiva et al. (2019); Meier et al. (2021) and Bartosova et al. (2019), the effect of nutrient load reduction policy should dominate climate change impacts, at least in the near future.Elevated halocline strength observed in the model (Naumov et al., 2023) also worsens oxygen conditions in the deep sea, which could mimic the negative effects in the real system related to climate change.Still, the dynamics of cyanobacteria blooms might be underestimated in the model since Temporal dynamics of the total annual oxygen fluxes by the three categories (physical fluxes -Phy (blue line), consumption in the water column -Bio (red line), and consumption in the sediments -Sed (dark brown line)).The categories are explained in more detail in the Section 3.1.For a list of processes comprising each category see Supplementary Table 1.Zero is marked by the horizontal black translucent line.Positive fluxes mean oxygen supply (only fluxes by the physical processes can be positive), and negativeoxygen consumption (by the water column and sedimentary processes).For negative oxygen fluxes (oxygen consumption), a negative linear trend means amplified oxygen consumption, and a positivereduced consumption.The reversed logic is applicable to the positive fluxes (oxygen supply), where a positive trend means increasing supply and a negativedecreasing.Grey lines represent the reference scenario.Only significant linear trends (p< 0.05) are shown in the figure.Mt/a stands for 10 9 kg per year.
their blooms may get amplified by climate change.Since significant positive trends in nitrates were observed everywhere across the Gotland Basin (see Supplementary Figures 6-9), it could potentially lead to a serious deterioration in case of a sudden short-time anoxia and release of the sedimentary phosphorus (Mort et al., 2010).Another possible shortcoming of the study is related to the difficulties in validating the sources and sinks of O 2 and H 2 S due to the lack of proper observational data (long-term monitoring at the specific station in the deep Baltic Sea).Only studies by Schneider et al. (2010) and Gustafsson and Stigebrandt (2007) reconstructed the oxygen sources and sinks based on observational data.Unfortunately, their results cannot be qualitatively compared to our model findings due to the different spatiotemporal scales and different formulations of processes in the model.Still, they exhibit similar patterns (elevated oxygen consumption by nitrification after more oxygen enters the system, for example).

Conclusions
1.The central Baltic Sea under the MAIs of the BSAP and the halved BSAP showed an improvement and transition to a more oxic state after 2018 within the simulated 71 years.However, the hypoxic area did not reduce dramatically in both scenarios, pointing out the high oxygen debt in the deep central Baltic Sea. 2. Positive changes were observed in all regions of the Gotland Basin, especially in the remote northern and western Gotland basins.In the nGB, a shift from consumption in the water column to consumption in the sediments was observed in both BSAP and 0.5 BSAP.In the wGB, the same changes were observed under 0.5 BSAP.Positive trends in the water column oxygen consumption are mostly explained by the reduced oxidation of hydrogen sulfide and nitrification towards the end of the simulation (the year 2089).In the sediments, the increased consumption is mostly attributed to the elevated oxidation of organic matter towards the end of the simulation.Faster improvement in the remote basins was attributed to the positive feedback related to the reduction of H 2 S concentration in eGB and, therefore, its less advection to the nGB and wGB.Elevated oxygen concentrations in the eGB could also lead to more effective ventilation of the remote basins by inflowing oxygen due to reduced consumption on the way to the remote basins.3. The general trend to the less explained variance by water column oxygen consumption and more by the sedimentary oxygen consumption was found across all studied subbasins.Physical fluxes (mainly advection) explain most oxygen variability in all scenarios.Hydrogen sulfide dynamics was found to be more region dependent and influenced by the different nutrient input scenarios.4. Three EOFs were identified as patterns governing the dynamics of oxygen and hydrogen sulfide in the reference scenario.For oxygen budget terms, the first EOF explained approx.60% of the variability, the secondapprox.20%, and the thirdapprox.10%.For the hydrogen sulfide dynamics, the first EOF explained more than 80% of the variability, and all three EOFs together explained more than 95%.The last two EOFs were attributed to the inflows' activity and the natural variability, and the first EOFto the change in the nutrient forcing with positive (deposition of reduced material and organic matter, deoxygenation) and negative (oxidation of the reduced material, oxygenation) phases.The first EOF Maps showing different processes' contributions to the O 2 (left) and H 2 S (right) budgets in the water column of the four studying sub-basins.
Processes are aggregated into the same three groups (phy, bio, and sed) as in Figure 2.More information about processes in each group can be found in Supplementary Tables 1 (oxygen) and 2 (hydrogen sulfide).In bar plots, the y-axis represents a fraction of the total variability explained by a certain group of processes, and the x-axis represents scenarios.Note that oxygen charts have two y-axes with different scales.The left and right axes represent the fraction of variability explained by the group phy and the groups bio as well as sed, respectively.
went into the constant negative phase (no deposition of reduced material and oxidation of already existing one) in the 2050s (O 2 , BSAP), the 2030s (O 2 , 0.5 BSAP), and the 2030s (H 2 S, both BSAP and 0.5 BSAP), indicating resilient transformation to the oxic regime.5.According to our model study, it is possible for the Baltic Sea to return to its initial state (the year 1948) within 71 years under the 0.5 BSAP scenario (see Figures 2, 4).Under the BSAP scenario, the Baltic Sea state comes close to the initial state but did not reach it within the simulation time.
We thank the Open Access Fund of Leibniz Association for funding the publication of this article.EOF decomposition of the spatial-temporal matrix of oxygen consumption terms aggregated into specific groups.Group names [x-axis labels in E)] should be read as follows:<domain>_<name of processes>.There are two domains: biowater column, and sedsediments; and six processes: nitr/nitdenitnitrification, photphotosynthesis, resprespiration of biota, min_sulfnineralization of sulfur, min_ommineralization of organic matter, detmineralization of detritus only.To get more information about individual processes in the group, see Supplementary Table

FIGURE 2
FIGURE 2 FIGURE 4 1 . (A) Shows the fraction of variability explained by all calculated EOFs.Only the first three EOFs were considered significant and used later in the analysis [colored red in (A)].(B-D) Depict the temporal variability of the first three EOFs correspondingly.Here, the black vertical line demarcates the reference scenario (grey curve) and BSAP, 0.5 BSAP scenarios (vivid and translucent red lines, correspondingly).(E) Demonstrates the loadings of the first three EOFs.Loadings vary from minus one to one and show the direction and magnitude of the connection between an EOF and a variable.Black lines highlight the spatial structure of the matrix by separating the sub-basins.