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₂) and hydrogen sulfide (H₂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.

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₂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₂ or NO 3 − . For the description of all processes contributing to the O₂ and H₂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 70-meter 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 .

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.4 Kton/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.

Following the approach of Naumov et al. (2023), we aggregated all O₂ and H₂S budget terms into three categories based on their origin (physical or biological) and domain (water column or sediments): physical processes encompassing all oxygen fluxes with physical origin, mainly lateral and vertical advection (phy), oxygen fluxes with biological origin situated in the water column, e.g., remineralization of OM and zooplankton respiration (bio), and oxygen fluxes located in the sediments, e.g., remineralization of the sedimentary detritus (sed). Results are shown in Figure 2 (oxygen) and Supplementary Figure 3 (hydrogen sulfide). Figure 2 shows significant changes in the oxygen consumption pattern under both BSAP and 0.5 BSAP scenarios. Especially significant changes happened in the nGB, where oxygen consumption shifted back to the sediments at the end of the study period in both scenarios (with 0.5 BSAP amplifying the reported changes). The wGB demonstrates a similar pattern. However, it achieved sediment dominance in the consumption only under the 0.5 BSAP scenario. Under the BSAP scenario, a significant negative trend was only observed in sedimentary consumption. Elevated oxygen consumption in the sediments under BSAP and 0.5 BSAP scenarios, most visible in the nGB, can be interpreted as a positive sign indicating reoxygenation of the sediments. As was concluded by Naumov et al. (2023), less oxygen consumption in the sediments indicates the absence of oxygen, which means no electron acceptor available for the reactions. As sediments contact with oxygen again, the reduced material is getting oxidized promoting oxygen consumption in the sediments. This mechanism is more pronounced under the more rigorous 0.5 BSAP scenario. Noticeably, a significant positive trend in physical fluxes was observed under the 0.5 BSAP scenario, which is related to better ventilation due to the improved oxygen conditions in the neighboring nGB. Unlike nGB and wGB, BB and eGB are closer to the Baltic Sea’s entrance and therefore receive more oxygen via advection (Liblik et al., 2018; Naumov et al., 2023), do not demonstrate striking trends in consumption terms. The only significant trend in the BB is the positive trend in sedimentary consumption under the 0.5 BSAP scenario, indicating the best oxygen conditions among the studied subbasins. In the eGB, significant positive trends in the water column oxygen consumption were observed under both BSAP and 0.5 BSAP scenarios. It also points out the improvement because less consumption in the water column indicates less reduced material stored there. H₂S sources and sinks ( Supplementary Figure 3 ) also demonstrated substantial changes. By the end of the study period, the eGB came to the dynamic balance between H₂S sources and sinks, indicating that H₂S is not deposited. The same changes by the end of the period are observed even in the nGB under the 0.5 BSAP scenario. Under the BSAP scenario, a significant reduction in H₂S production and consumption is observed in the nGB. However, it stabilizes in the 2070s. Despite that, advection to the upper layers in the nGB stops by the end of the study period in both scenarios. The wGB demonstrated the most noticeable difference between the BSAP and 0.5 BSAP scenarios. Under the BSAP scenario, wGB is still exporting H₂S to the upper layer by the end of the study period, although to a much lower extent than in the beginning. Sedimentary production exhibits a significant negative trend but stabilizes in the 2070s, and water column consumption shows no trend under the BSAP scenario. Under the 0.5 BSAP scenario, both consumption and production of H₂S converge to zero in the wGB, removing any export to the upper layers via advection.

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₂S budget is generally more complex (both regionally and in different scenarios). It points to the strong connection between H₂S dynamics and nutrient forcing (see Supplementary Figure 4 for additional information). Overall, it can be stated that H₂S dynamics is more affected by the nutrient loads reduction than oxygen dynamics, with some regional differences existing in both cases.

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 (1948-2018). 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₂, the first EOF explains around 60% of the total variability, the second EOF about 20%, and the third about 10%. For H₂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₂ and H₂S data can be attributed to the same process based on its loadings (Panel E in both figures). Both for O₂ and H₂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₂S production due to the 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₂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.

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 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₂ and H₂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).