Fate of legacy ammonium in the coastal Baltic Sea

Eutrophication has enriched coastal sediments globally with organic matter (OM), fuelling internal nutrient loading once benthic hypoxia occurs. Internal ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) loading may be particularly prominent in the coastal Baltic Sea due to the existence of a substantial ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool in its sediments, accumulated via OM deposition and mineralization over long-term eutrophication. However, despite its potential to exacerbate eutrophication, internal ${\mathrm{NH}}_{\mathrm{4}}^{+}$ loading has so far received little attention in the coastal Baltic Sea and it remains poorly understood for how long the sediment ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool could sustain internal ${\mathrm{NH}}_{\mathrm{4}}^{+}$ loading via sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes, especially under varying oxygen availabilities. To reconstruct past and predict future ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes in response to different OM loading and bottom water oxygen conditions, we developed a reactive transport model for muddy, OM-rich sediments of the coastal Baltic Sea. Our model results suggest a strong legacy effect in the sediment ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool, driven by mineralization of OM that accumulated during past and ongoing eutrophication, as expected. The influence of the sediment ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool on the sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes in our model is, however, smaller than the direct influence of the external OM loading. The sediment ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool keeps the sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes 5 %–10 % elevated compared to pre-industrial baseline values for at least 140 years after OM loading has reached constant baseline levels, albeit at a small efflux magnitude. While this magnitude can increase under increasing OM loading and oxygen deficiency, the overall legacy effect on the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes from coastal sediments is mild and the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes react well to external OM load reductions. If strict OM load reductions following the direction of the Baltic Sea Action Plan are implemented, sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes can naturally return to a level similar to the 1950s, the start of widespread eutrophication, within 110 years, given that sediments stay normoxic and bioturbated during the reduction period. Without external OM load reductions, sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes stay indefinitely elevated and fuel continuous internal ${\mathrm{NH}}_{\mathrm{4}}^{+}$ loading to the coastal system. Our findings highlight the complexity of eutrophication legacy effects in regard to the sediment pool and sediment-water efflux of ${\mathrm{NH}}_{\mathrm{4}}^{+}$, and their varying impact on the coastal Baltic Sea ecosystem. The knowledge obtained is beneficial also for other anthropogenically impacted coastal seas with similar geomorphology as the Baltic Sea.

1 Introduction

Excessive land-derived nutrient loading has led to eutrophication of the coastal oceans globally, thereby strongly enriching these systems with organic matter (OM). The combination of high production rate and fast sedimentation of OM in coastal systems can quickly deplete oxygen (O₂) in the benthic habitat, leading not only to ecosystem degradation (Breitburg et al., 2018) but also to the rapid accumulation of OM in sediments, creating a eutrophication-induced OM legacy pool (Gustafsson et al., 2012). Over time, aerobic and anaerobic mineralization of this OM generates nutrient pools in the sediment, including ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) and dissolved inorganic phosphorous (DIP; Aller, 2014). The retention (DIP) and turnover (${\mathrm{NH}}_{\mathrm{4}}^{+}$) of these nutrients within the sediment is usually high under oxic bottom water conditions; yet, when benthic hypoxia (O₂ concentration <63 µM or 2 mg L⁻¹, Middelburg and Levin, 2009) develops, both nutrients are preferentially released into the water column (Koop et al., 1990; Cowan and Boynton, 1996; Norkko et al., 2015), leading to internal nutrient loading with high potential to reinforce eutrophication (Pitkänen et al., 2001; Conley et al., 2007). Research on internal nutrient loading has often focused on DIP due to its crucial role in the repeated occurrences of massive cyanobacteria blooms and close coupling to metal cycling (Conley et al., 2002; Vahtera et al., 2007; Viktorsson et al., 2013; Matisoff et al., 2016), whereas internal ${\mathrm{NH}}_{\mathrm{4}}^{+}$ loading has so far received little attention. This omission is a critical oversight, considering that excess bioavailable nitrogen (N) is a key driver of coastal eutrophication (Ryther and Dunstan, 1971).

Figure 1 (a) Sediment-water ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) efflux at different bottom water oxygen (BW O₂) concentrations from coastal sites around the world, and (b) statistical significance (p≤0.05) of each site's relationship between ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux and BW O₂ tested with linear regression analysis (n=number of measurements, $*=\text{statistical}$ significance; 0 µM O₂ was changed to 0.1 µM O₂ to conform with the analysis). Sites include the Oregon Shelf (Fuchsman et al., 2015), northern Gulf of Mexico (McCarthy et al., 2008; Nunnally et al., 2014; Roberts and Doty, 2015), Chesapeake Bay (Cowan and Boynton, 1996; Testa and Kemp, 2012), Waquoit Bay (Foster and Fulweiler, 2019), Orust Fjord (Norkko et al., 2019), Baltic Sea (Conley et al., 2007; Bonaglia et al., 2014; Ekeroth et al., 2016; Gammal et al., 2017), Mandovi Estuary (Pratihary et al., 2009), Bohai Bay (Mu et al., 2017), coastal East China Sea (Song et al., 2021), Jinhae Bay (Huang and An, 2022), and Derwent Estuary (Banks et al., 2012). All data originate from core- or chamber-incubations and were included in this overview if they supplied ${\mathrm{NH}}_{\mathrm{4}}^{+}$ flux measurements under varying O₂ concentrations; we do not claim an exhaustive literature review, such as can be found by Boynton et al. (2017). ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes are given as average of replicate measurements per investigation without standard-deviation or -error; ${\mathrm{NH}}_{\mathrm{4}}^{+}$ uptake into the sediment is shown as zero efflux. The red dashed line gives the limit of hypoxia (O₂ concentration<63 µM, $\text{resp.}<\mathrm{2}\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}$, Middelburg and Levin, 2009). Map is drawn with Ocean Data View (Schlitzer, 2023).

Internal ${\mathrm{NH}}_{\mathrm{4}}^{+}$ loading may be particularly prominent in the long-term eutrophied coastal Baltic Sea (defined here as the coastal areas of the Baltic Sea extending one nautical mile seaward from territorial water baselines, following the Water Framework Directive; European Communities, 2003). Compared to coastal regions worldwide, the coastal Baltic Sea displays the highest sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes of all studied sites and a highly pronounced relationship between ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux and bottom water O₂ condition (Fig. 1A). ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes increase significantly with decreasing O₂ concentrations across all investigated sites in the coastal Baltic Sea, despite their dispersed geographic locations (Denmark, Sweden, and Finland) and only moderate number of replicates in the cited studies (n=6–14; Fig. 1B).

This clear pattern might be explained by Baltic Sea-specific environmental factors that reduce the influence of natural variability in measurements, such as low species diversity, negligible tidal influence, and a long water residence time, as well as by the presence of a significant benthic legacy pool of ${\mathrm{NH}}_{\mathrm{4}}^{+}$. Since the acceleration of Baltic Sea eutrophication in the 1950s (Zillen and Conley, 2010; Gustafsson et al., 2012), OM and eventually ${\mathrm{NH}}_{\mathrm{4}}^{+}$ have accumulated in the porewaters of muddy accumulation bottoms, visible today as characteristic ${\mathrm{NH}}_{\mathrm{4}}^{+}$ profiles of strongly increasing concentrations with depth that are consistent across seasons (Fig. 2). In a first-order estimate, using coastal boundaries as defined above (European Communities, 2003), coastal sediment classifications from Al-Hamdani and Reker (2007), and example porewater ${\mathrm{NH}}_{\mathrm{4}}^{+}$ profiles from representative sandy (site VE15, southern Bay of Gdansk, Thoms et al., 2018) and muddy (site Storfjärden, north-western Gulf of Finland, D. Hellemann Fig. 2) sediments, we suggest that >100 kt ${\mathrm{NH}}_{\mathrm{4}}^{+}$ could be stored in the top 20 cm of Baltic Sea coastal sediments. This benthic pool of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ is not only larger than the pelagic pool of dissolved inorganic nitrogen (DIN) in the coastal Baltic Sea, estimated at ∼61 kt N using a 21 % share of DIN in total nitrogen (TN; C. Lønborg, pers. comm.) and a coastal pelagic TN pool estimate from Lønborg and Markager (2021), but, furthermore, it is only a minimum estimate of the entire ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool in coastal Baltic Sea sediments considering that concentrations continue to increase with depth below 20 cm.

Figure 2 Example porewater profiles of ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) in muddy sediments of the coastal Baltic Sea, taken from (a) site Storfjärden (33 m deep, April–October 2016) and (b) site Långholmsbranten (21 m deep, December 2022–November 2023), both north-western Gulf of Finland; data: D. Hellemann. The pattern of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations increasing with sediment depth and not showing significant variations between seasons is in line with profiles from other muddy sites of the coastal Baltic Sea, such as the Himmerfjärden estuary, western Baltic Proper (Bonaglia et al., 2014), and the coastal Luebeck Bight and Tromper Wiek, south-western Arkona Basin (Gogina et al., 2018).

The impact of this sediment ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool on pelagic processes is largely regulated by O₂ availability, although factors such as OM availability and benthic infauna also affect the magnitude of sediment-water fluxes (Boynton et al., 2017). Commonly, under normoxic conditions, a substantial portion of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ that diffuses upwards from porewater is oxidized to nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) via nitrification in the oxic surface sediment and subsequently is reduced to bio-unavailable di-nitrogen gas (N₂) via denitrification in the underlying anoxic sediment (Fig. 3A). The same also happens in the numerous 3-dimensional oxic-anoxic interfaces created by benthic infauna (Andersen and Kristensen, 1991). Further, microphytobenthos, if present on the sediment surface, can intercept and assimilate nutrients effluxing from the underlying sediments (Sundbäck and Miles, 2000). Together, these factors commonly buffer against large ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes from the sediment, leaving only residual fluxes originating from, for example, faunal activity (Andersen and Kristensen, 1991; Karlson et al., 2007). When benthic hypoxia or anoxia develops, the oxic sediment surface is lost, benthic infauna and microphytobenthos are extirpated, and the benthic N turnover processes change from N removal to N recycling (Kemp et al., 1990; Jensen et al., 1990; Childs et al., 2002; Karlson et al., 2007; Jäntti and Hietanen, 2012). The extent of these effects depends on frequency and duration of O₂ deficiency (Conley et al., 2007; Jäntti and Hietanen, 2012; Villnäs et al., 2012), but they all fuel ${\mathrm{NH}}_{\mathrm{4}}^{+}$ release into the water column compared to normoxic conditions, enriching the pelagic pool of bioavailable N (Fig. 3B). While exceptions to this pattern exist (Foster and Fulweiler, 2019; Testa et al., 2025), the available field data from the coastal Baltic Sea show a significant trend following this paradigm, resulting in internal loading of bioavailable N under hypoxia (Fig. 1). In the open Baltic Sea, such internal N load would be efficiently removed via denitrification at the redoxcline of the permanently hypoxic deep waters (Hietanen et al., 2012; Dalsgaard et al., 2013). However, in the coastal Baltic Sea, pelagic denitrification is not a common feature as density- and wind-driven water mixing usually prevent large-scale permanent hypoxia (Leppäranta and Myrberg, 2009; Conley et al., 2011) and thus the establishment of a pelagic denitrifier community (Falk et al., 2007; Kim et al., 2011). Consequently, in the coastal Baltic Sea, the internal N load from sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes remains in the pelagic system and is thus available for coastal primary production, potentially enhancing coastal eutrophication.

Figure 3 Conceptualized ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) turnover in coastal sediments under (a) normoxic and (b) hypoxic conditions, based on Kemp et al. (1990), Jensen et al. (1990), Conley et al. (2007), Jäntti and Hietanen (2012) and others. Bioavailable nitrogen (N)-compounds shown in green, non-bioavailable N-compounds shown in violet. Arrow thickness symbolizes relative flux magnitude; font size symbolizes relative compound abundance. Anaerobic ammonium oxidation is not shown as it is often negligible in organic-rich coastal sediments (Dalsgaard et al., 2005). (a) Normoxic conditions: ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux across the sediment-water interface is small, as most of the upwards diffusing ${\mathrm{NH}}_{\mathrm{4}}^{+}$ is oxidized to nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$; nitrification) when passing the oxic sediment layer. Due to typically strong coupling of nitrification and denitrification in coastal sediments, most of the resulting ${\mathrm{NO}}_{\mathrm{3}}^{-}$ is subsequently reduced to non-bioavailable di-nitrogen gas (N₂; denitrification) in the anoxic sediment layer. Hence, the addition of bioavailable N compounds to the water column is relatively small. (b) Hypoxic conditions: ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux across the sediment-water interface is large, as the oxic sediment layer is lost, disabling nitrification in the sediment and thus the transformation of upwards diffusing ${\mathrm{NH}}_{\mathrm{4}}^{+}$ to ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and its subsequent removal in denitrification. Once in the water column, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ accumulates under hypoxic conditions or is oxidized to ${\mathrm{NO}}_{\mathrm{3}}^{-}$ under oxic conditions, either way increasing the pelagic pool of bioavailable N. Some of this water column ${\mathrm{NO}}_{\mathrm{3}}^{-}$ can be expected to diffuse back into the sediment where it can be removed in denitrification or transformed to ${\mathrm{NH}}_{\mathrm{4}}^{+}$ in dissimilatory nitrate reduction to ammonium (DNRA); however, being subject to a dynamic coastal system, the additional N might also be quickly exported to other water layers and, for instance, taken up in primary production. If hydrogen sulphide (H₂S) builds up in the surface sediment, DNRA commonly dominates over denitrification (Jäntti et al., 2021), further enriching the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool.

The timeline of internal ${\mathrm{NH}}_{\mathrm{4}}^{+}$ loading affecting the Baltic Sea coastal ecosystem is uncertain, as the long-term fate and persistence of the sediment ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool is largely unknown. For the TN pool in Baltic Sea sediments, Lønborg and Markager (2021) suggest a legacy effect of several centuries. This gap in knowledge limits the capacity of eutrophication mitigation measures, such as the Baltic Sea Action Plan (BSAP), an agreement between the Baltic states on specific actions to bring the Baltic Sea back to a good environmental status (HELCOM, 2007, 2021). Coastal hypoxia is likely to increase in the future due to climate change driven water warming (Caballero-Alfonso et al., 2015; Breitburg et al., 2018), increasing the potential for internal ${\mathrm{NH}}_{\mathrm{4}}^{+}$ loading. Furthermore, present-day OM loading to coastal Baltic Sea sediments remains at $\sim \mathrm{10}\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (estimated from the physical-biogeochemical BALTSEM model, Gustafsson et al., 2012), which is lower than the peak loading of $\sim \mathrm{14}\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in the 1990s, but still several-fold higher than before the onset of eutrophication (Fig. 4a). As a result, new OM is continuously added to sediments still containing the OM legacy from past eutrophication, which is likely to extend the lifetime of the existing sediment ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and its potential to impact the coastal ecosystem.

Interestingly, while data from the Chesapeake Bay and the Derwent Estuary also show a significant relationship between sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux and O₂ deficiency (Fig. 1B), the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ porewater profiles of the Chesapeake Bay do not indicate the existence of a ${\mathrm{NH}}_{\mathrm{4}}^{+}$ legacy pool (Kemp et al., 1990; Boynton et al., 2017; no porewater data available for the Derwent Estuary). Despite a eutrophication history similar to the Baltic Sea, the OM deposition environment of the Chesapeake Bay does not seem to support the build-up of benthic legacy pools, owing to higher energy of physical processes, including tidal flushing, and shorter water residence time due to a more open geomorphology. Hence, eutrophication-induced legacy effects similar to the ones from the coastal Baltic Sea might be expected more commonly in eutrophic systems of a semi- to fully enclosed geomorphology and long water residence time, such as the Black Sea (BSC, 2019), Sea of Marmara (Demirel et al., 2023; Akçay et al., 2025) or even inland waters such as Lake Erie (Watson et al., 2016). Owing to its unique data-richness, the coastal Baltic Sea can serve as a key example for these systems and their potential future evolution (Reusch et al., 2018).

2 Reactive transport modelling to investigate benthic ${\mathrm{NH}}_{\mathrm{4}}^{+}$ turnover and eutrophication legacy effects

To shed light on the long-term fate and persistence of the legacy ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool in sediments of the coastal Baltic Sea, we investigate future scenarios of sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes in response to OM loading and bottom water O₂ availability, using a 1-dimensional reactive transport model (RTM).

2.1 Model set-up

The 1-dimensional RTM is developed using the R package “ReacTran” (Soetaert and Meysman, 2012) and comprises a multi-component diagenetic setting to examine sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux in response to two different OM loading scenarios and three different bottom water O₂ conditions (Table 1). The RTM assumes a 1 m long sediment column, as the top meter of sediment is usually highly reactive and determines the mass exchange between water and sediment. Into this sediment column, the redox ladder of marine OM degradation is fully integrated, so that OM is sequentially oxidized with O₂, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, iron and manganese oxides, and sulphate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$), followed by methanogenesis when porewater ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations are low. These six primary reactions of OM oxidation are complemented by eight secondary reactions (Appendix A, Table A1), following Reed et al. (2011) with some modifications. All reaction and transport parameters and reaction rates are given in Appendix A (Tables A2 and A3). Nitrification and denitrification, the two main processes controlling benthic N cycling, are also included, whereas anaerobic ammonium oxidation is not considered due to its commonly negligible occurrence in organic-rich coastal sediments (Dalsgaard et al., 2005).

Boundary conditions at the top and bottom of the sediment model domain are defined for both solid and liquid phases. For solid phases, an OM deposition flux ($\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) derived from the simulation of the BALTSEM model (Gustafsson et al., 2012) is prescribed at the sediment-water interface and as the model runs, the OM accumulates over time in the sediment column or is removed from the domain bottom by deep burial. To describe the deposited OM, the model uses three OM pools: a reactive pool with a molar carbon/nitrogen ($\mathrm{C}/\mathrm{N}$) ratio of 6.6 representing marine-derived fresh OM (62 % of deposited OM), a less reactive pool with a molar $\mathrm{C}/\mathrm{N}$ ratio of 13 representing terrestrial-derived OM (15 % of deposited OM), and a non-reactive pool that gets buried directly without degradation (23 % of deposited OM), according to model calibration data from Jilbert et al. (2021; see also Appendix A). For liquid phases, the upper boundary condition is set to be a constant bottom water concentration, while the lower boundary is open so that the RTM determines how much ions could be “leaked out” of the model domain via diffusion and deep burial. The ion concentrations of the bottom water are assumed to be constant and zero, except for ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (5 mmol L⁻¹), ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (0.007 mmol L⁻¹) and O₂ (0, 0.063, and 0.2 mmol L⁻¹, equivalent to anoxia, hypoxia, and normoxia; Table 1). Together, the set-up allows investigation of how the sediment system responds to changes in OM loading, which in turn regulate the sediment ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and the sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux.

The RTM runs on a 1 cm depth resolution and uses sediment porosities representative of muddy, OM-rich accumulation bottoms, which, with ∼41 %, is the predominant bottom type of the coastal Baltic Sea (calculated using Al-Hamdani and Reker, 2007 and coastal boundaries according to European Communities, 2003). The RTM output is given in 1 year time step. To calibrate and verify the model, we used the geochemical porewater profiles and environmental conditions, such as sedimentation rate, nutrient concentrations, and salinity, of the coastal site Storfjärden in the north-western Gulf of Finland (Jilbert et al., 2021, Appendix A, Fig. A1–A3, Table 1; https://doi.org/10.5281/zenodo.22120110).

2.2 Scenario modelling

The OM loading used in the RTM is generated from the BALTSEM model (Gustafsson et al., 2012). Initially, the RTM is run into a steady state using pre-industrial levels of OM flux to the sediment surface as a baseline prior to eutrophication (period 1800–1850). Then, hindcast OM sedimentation rates are used as a forcing from 1850–2020, followed by two future scenarios from 2020–2400: In the “Present scenario”, OM loading is held constant at the 2020 level, and in the “Pre-industrial scenario”, OM loading is reduced in accordance with the BSAP target until 2100, and subsequently continues to decline until pre-industrial levels are reached again around 2160 (Table 1). By adopting this approach, we are able to compare the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and fluxes in 2100, when BSAP targets of “good environmental status” are considered to be achieved (Murray et al., 2019), with the later situation when true pre-industrial conditions are re-established. Each scenario is run for three bottom water O₂ availabilities (normoxia, hypoxia, and anoxia; Table 1), representing best and worst case O₂ conditions. To account for the effects of bioturbating infauna on OM burial and O₂ transport during normoxia, we use the empirical parameter D_b (Boudreau, 1997; Appendix A, Table A2) and furthermore force O₂ to penetrate 2 cm deep into the sediment and subsequently to decrease exponentially to 0 µM at 5 cm depth. Initial tests showed no significant difference in ${\mathrm{NH}}_{\mathrm{4}}^{+}$ fluxes between a forced O₂ penetration depth of 1 or 2 cm. The model code for all simulations can be found at https://doi.org/10.5281/zenodo.22120110.

Table 1 Key input variables and scenario assumptions used in the RTM. Salinity, sulphate, and sedimentation rate are adopted from the coastal site Storfjärden, north-western Gulf of Finland. OM = organic matter, C = carbon, O₂ = oxygen, BSAP = Baltic Sea Action Plan.

Due to ongoing climate change, bottom water temperatures will likely rise over the modelled future period and affect reaction rates and resulting ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes (Boynton et al., 2017), as heterotrophic metabolic processes are expected to increase at twice the rate of autotrophic production (Harris et al., 2006). To provide a first-order assessment of temperature effects on ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes, we carried out a sensitivity test by increasing all used rate constants in the model by 20 % (Gudasz et al., 2015). The resulting effect on the modelled ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes was negligible, hence, we do not further discuss the effects of climate change on our simulations. However, since our model is a generalization that does not account for variable temperature effects, for example along water depth gradients or seasons, our modelled ${\mathrm{NH}}_{\mathrm{4}}^{+}$effluxes may be considered to represent a lower limit of potential values if bottom waters continue to warm in the future.

3 Coastal ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes under future scenarios of OM loading and O₂ availability

The BALTSEM-derived OM loading in both model scenarios can be divided into four phases (Fig. 4). Phase I gives the pre-eutrophication OM level (including pre-industrial baseline levels) representing natural loading. Phase II is characterized by a sharp increase in OM loading from the onset of eutrophication in the 1950s to the early 2000s, during which the loading reaches its peak of $\sim \mathrm{14}\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, tripling the OM loading of phase I. Following this peak, phase III covers the decrease in OM loading to its present-day value of $\sim \mathrm{10}\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Present scenario), or to the projected pre-industrial baseline levels (Pre-industrial scenario). Phase IV represents the far-future return to pre-industrial OM loading and thus can only be seen in the Pre-industrial scenario.

Figure 4 Sediment organic matter (OM) loading ($\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), sediment ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) pool (mmol m⁻³ bulk sediments), and sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) in the northern coastal Baltic Sea under anoxic (red), hypoxic (green), and normoxic (blue) bottom water oxygen (O₂) conditions, modelled for the OM loading scenarios (a) “Present scenario” (constant OM loading) and (b) “Pre-industrial scenario” (reduced OM loading down to pre-industrial levels). Normoxic O₂ conditions include bioturbation, which is absent under O₂-deficient conditions as the benthic infauna dies. Modelling period 1800–2400. The ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool is the amount of porewater ${\mathrm{NH}}_{\mathrm{4}}^{+}$ integrated over a 1 m sediment model domain. Grey shadings indicate the different phases of OM loading with (I) pre-eutrophication, (II) strongly increasing, (III) after-peak, and (VI) return to pre-industrial. Horizontal dashed lines give the pre-industrial baselines of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and efflux at normoxic and O₂-deficient conditions. The vertical dashed line in (b) gives the end point of OM load reductions if following the target of the Baltic Sea Action Plan (BSAP, see Table 1). The green arrow indicates the onset of widespread eutrophication around the 1950s. Note that the hypoxic ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool differentiates very little from the anoxic ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and therefore is only marginally visible.

During phase I, the model reaches a steady state of a constant ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool in the sediment and low ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux to the water column. In phase II, the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and efflux respond to the fast increase in OM loading. In phase III of the Present scenario, the OM loading becomes constant in 2020, while the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool continues to increase steadily for over two centuries before approximating a slight plateau around 2250; in contrast, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes reach an actual plateau around 2100 and thereafter stay high indefinitely under all O₂ conditions (Figs. 4a, 5a and c). In phase III of the Pre-industrial scenario, the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool under normoxic conditions decreases in accordance with the drastic decrease in OM loading, albeit at a comparatively much slower rate: when OM loading reaches pre-industrial baseline values around 2160, the normoxic ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool is still as large as in the 1970s. In contrast, the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool under O₂-deficient conditions continues to increase for about 100 years despite the reductions in OM loading, before eventually starting to slowly decrease (Figs. 4b and 5b). Under neither O₂ condition does the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool of the Pre-industrial scenario reach pre-industrial baseline levels by year 2400, the end of our modelling period (Fig. 4b phase IV). Despite the difference in pool development, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes in phase III of the Pre-industrial scenario follow the reduction of OM loading closely under all O₂ conditions (Figs. 4b and 5d). Nevertheless, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes do not reach pre-industrial baseline levels by the end of phase III, but show remaining efflux-“tails” as can be seen when zooming into the model-output (Appendix B, Fig. B1). True pre-industrial baseline levels of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux are reached in phase IV first after 2300 (normoxic conditions) or even around 2400 (anoxic conditions; Appendix B, Fig. B1); the flux magnitude of these “tails” under all O₂ conditions is however small and comparable to the pre-eutrophication period (Appendix B, Fig. B1). In both scenarios, the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool is smaller and increases at a lower rate under O₂-deficient than under normoxic conditions, whereas the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes are larger under O₂-deficient than under normoxic conditions (Figs. 4 and 5). The latter one is in line with field data from the coastal Baltic Sea (Fig. 1) and the common understanding of more efficient internal ${\mathrm{NH}}_{\mathrm{4}}^{+}$ turnover under normoxic conditions (Fig. 3).

Figure 5 Sediment ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) pool (a, b) and efflux (c, d) plotted against organic matter (OM) loading under two conditions of bottom water oxygen (O₂) and its consequential absence/presence of benthic fauna (triangles=0 µM O₂; circles=200 µM O₂ and bioturbation), shown for the “Present scenario” (a, c) and the “Pre-industrial scenario” (b, d). Arrows and color-coding indicate the temporal evolution, with dark blue=period 1800–1900, light blue=period 1900–2000, orange=period 2000–2100, and red=period 2100–2400.

Depending on the O₂ condition, the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes of the Pre-industrial scenario remain elevated compared to pre-industrial baseline levels for 280–380 years after the start of drastic OM load reductions (Appendix B, Fig. B1). However, about 110 years following the reductions, i.e. in 2130, normoxic ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes have shrunk to a level comparable to the 1950s, the starting point of widespread eutrophication; under O₂-deficient conditions, this level will be reached about 140 years following the reductions, i.e. in 2160 (Fig. 4b). This modelled timeframe of decline in ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux is 30–60 years longer than the estimate of when BSAP-driven good environmental status is assumed to be reached (year 2100, Fig. 4; Murray et al., 2019), which could suggest that it might take longer than previously thought for the BSAP to achieve its environmental targets.

In both scenarios, the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes seem to respond more rapidly to the external OM loading than to the internal ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool (Figs. 4 and 5). This result is interesting, as the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool shows a strong legacy effect and thus creates high internal ${\mathrm{NH}}_{\mathrm{4}}^{+}$ availability, even long after strong external OM loading has ceased, such as in the Pre-industrial scenario (Fig. 5b). This legacy effect is caused by ${\mathrm{NH}}_{\mathrm{4}}^{+}$ being continuously regenerated from the OM that accumulated in the sediment over subsequent years of strong eutrophication (Fig. 4 phase II), and results in (i) the increase of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool in the Present scenario, despite constant OM loading (Figs. 4a phase III and 5a), as well as (ii) the slow decline (normoxic condition) and initial increase (O₂-deficient conditions) of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool in the Pre-industrial scenario, despite rapid OM load reductions (Figs. 4b phase III–IV and 5b). The different evolution of the normoxic and O₂-deficient ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pools can be attributed to the more efficient mineralization and internal ${\mathrm{NH}}_{\mathrm{4}}^{+}$ turnover under normoxic than O₂-deficient conditions (Fig. 3). The influence of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool on the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes eventually becomes visible in the Pre-industrial scenario, where the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes remain 5 %–10 % above pre-industrial baseline values after OM loading has returned to constant baseline levels (Appendix B, Fig. B1, Fig. 4b), creating ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux-“tails”. While these tails are small in magnitude, they last for 140–240 years, depending on the O₂ condition (Appendix B, Fig. B1), clearly indicating the influence of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and its strong legacy effect. In consequence, the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes are subject to a eutrophication legacy effect, which is highly persistent, but only small in effect size. Hence, the majority of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes remain controlled by the external OM loading, which suggests that ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes and their subsequent impact on the coastal environment can be mitigated in a straightforward manner by reducing the external OM loading, such as in our Pre-industrial scenario.

While our modelled ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes in both scenarios show distinctly different responses between normoxic and O₂-deficient conditions, confirming the common assumption of O₂ as a key control of sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes (Fig. 3), the magnitude of our modelled O₂-deficient fluxes is below the range of field-observations from the coastal Baltic Sea (Fig. 1). This underestimation of O₂-deficient ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes can be explained by the resolution of our RTM, which runs on an annual basis and thus does not capture seasonal fluctuations and anoxic extreme events that can drastically increase the efflux. Such events are generally challenging to incorporate in a model, as they originate from individual combinations of OM loading and loading history, anoxia duration, temperature, and more. Without the addition of extreme values to the RTM's annual average, the modelled O₂-deficient ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes are lower than they can periodically be in nature, and our model consequently provides an estimate of O₂-deficient NH₄ effluxes at the lower range of their annual value. This underrepresentation of seasonal extreme values in our model affects only the O₂-deficient, but not the normoxic fluxes, which fit well to the range observed in the coastal Baltic Sea (Fig. 1).

4 Sensitivity of benthic ${\mathrm{NH}}_{\mathrm{4}}^{+}$ turnover to changed boundary conditions

Our model results represent one site in the northern coastal Baltic Sea, which is considered representative for a muddy, OM-rich accumulation bottom in low-saline, eutrophic conditions. To test the sensitivity of our modelled ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and efflux to different environmental conditions, we applied a range of OM loadings, salinities (representative of different ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations in bottom water), sedimentation rates, bioturbation coefficients, $\mathrm{C}/\mathrm{N}$ ratios of OM, and different shares of reactive material in OM to our RTM, both under normoxic (Fig. 6) and anoxic (Appendix B, Fig. B2) bottom water conditions. As a baseline for comparison, we used the conditions of our model setting under the Pre-industrial OM loading scenario, i.e. a salinity of 6 with a ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$concentration of 5 mmol L⁻¹, a sedimentation rate of 6 mm yr⁻¹, strong bioturbation under normoxic conditions, and a $\mathrm{C}/\mathrm{N}$ ratio of 6.6 in the reactive fraction and 13 in the less reactive fraction. The OM loadings were generated by the BALTSEM simulations.

Figure 6 Effect of changed organic matter (OM) loading (I), salinity (II), sedimentation rate (III), bioturbation coefficient (IV), carbon to nitrogen ($\mathrm{C}/\mathrm{N}$) ratio in OM (V), and different shares of reactive material in OM (VI) compared to baseline conditions representative of the northern coastal Baltic Sea (Pre-industrial OM loading scenario, salinity of 6, sedimentation rate of 6 mm yr⁻¹, strong bioturbation, and OM with a $\mathrm{C}/\mathrm{N}$ ratio of 6.6 (reactive fraction) and 13 (less reactive fraction); black line) on (a) sediment ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) pool and efflux, and (b) its percent change relative to baseline condition. Modelled for normoxic bottom water conditions (200 µM O₂). Dashed lines give the pre-industrial baselines of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and efflux. Note the different y-axis scales between some parameters.

Under normoxic conditions and compared to the baseline, an increase in OM loading would strongly increase the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and efflux and both their legacy effects (Fig. 6I). While the increase in legacy effect in the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool is very visible, its increase in the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux is more subtle in larger magnitudes of the efflux-“tails” (Fig. 6a I), proportionally to the level of OM loading: the highest OM loading (5-times the BSAP OM loading) shows the strongest effect on the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes (Fig. 6b I). An increase in salinity does not change the size of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and only marginally increases the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux (Fig. 6a II), with the largest change observed when salinity changes from 6.6 to 10 (Fig. 6b II). This result is caused by the comparatively low availability of OM that limits the reduction of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$; hence, higher salinity can only increase the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and efflux, if also the OM availability is higher, such as can be expected under the Present scenario. A change in sedimentation rate affects the size of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool but not the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux (Fig. 6 III). The presence or absence of bioturbation has a strong effect on both ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and efflux (Fig. 6 IV). In the presence of bioturbation, the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool increases, likely because bioturbation promotes the degradation of OM (Andersen and Kristensen, 1991) and thus the release of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ to the sediment porewater. However, at the same time, the ${\mathrm{NH}}_{\mathrm{4}}^{+}$efflux decreases, as bioturbation generally increases the availability of O₂ in surface sediments, thereby facilitating sediment internal ${\mathrm{NH}}_{\mathrm{4}}^{+}$ turnover (Fig. 3). If bioturbation is absent, the opposite happens. This result highlights the importance of bioturbating infauna for sediment O₂ and N-turnover. An increase in $\mathrm{C}/\mathrm{N}$ ratio towards less bioavailable and more refractory OM will decrease both the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and efflux (Fig. 6V), as the OM becomes less rich in N and thus less ${\mathrm{NH}}_{\mathrm{4}}^{+}$ is produced in mineralization. The $\mathrm{C}/\mathrm{N}$ ratio is also affected by different shares of reactive OM, thereby affecting the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and efflux (Fig. 6VI). In this test, the share of refractory OM is kept at 23 %, and the share of reactive OM is 47 % and 10 %. Lower shares of reactive OM lead to a larger ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and a lower ${\mathrm{NH}}_{\mathrm{4}}^{+}$efflux (Fig. 6VI), which is caused by in turn higher shares of less reactive OM that degrade slowly but last longer in the sediments. Under anoxic bottom water conditions, similar patterns as the ones described can be observed but with slightly stronger effects on the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and efflux (Appendix B, Fig. B2).

While these results point to a general direction for ${\mathrm{NH}}_{\mathrm{4}}^{+}$ cycling in eutrophied coastal systems of different environmental conditions, they need to be interpreted while considering the geomorphological context for each system and the understanding that parameters occur in combination not in isolation. For instance, high salinity might often occur together with a low sedimentation rate, given that many coastal systems of high salinity are connected to the open ocean and thus experience high water energy from waves and tidal force, which enhances dispersion of OM and ultimately affects the sedimentation rate. In our sensitivity analysis, we do not consider such parameter combinations. Hence, understanding a habitat's geomorphological setting and environmental conditions are highly important when aiming to realistically model ${\mathrm{NH}}_{\mathrm{4}}^{+}$ flux predications of coastal sediments.

5 Implications for restoring eutrophied coastal ecosystems

Our model results highlight the presence of legacy effects from past eutrophication in present-day and future ${\mathrm{NH}}_{\mathrm{4}}^{+}$ cycling in coastal sediments of the Baltic Sea. While the legacy effects are highly persistent in both the sediment ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and the sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes, the legacy's effect size is much smaller in the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes than in the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool. This means, that the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool affects the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes for a long time, but potentially plays a minor role in its ecological impact. The impact of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux legacy on the coastal ecosystem is ultimately determined by the external OM loading, which not only strongly drives the effluxes (Fig. 4b) but also determines the magnitude and thus ecological significance of the legacy effluxes (Fig. 6I). The strong connection between external OM loading and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes gives a clear direction for the restoration of eutrophied coastal ecosystems. Our model results show that if strict OM load reductions are implemented, sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes will decrease promptly and nearly proportionally to the load reductions. About 110 years after the start of OM load reductions, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes could reach the same level as at the start of eutrophication in the 1950s. This is a relatively short time-frame for natural ecosystem recovery of a heavily eutrophied system, such as the coastal Baltic Sea. The time-frame can be met when bottom waters are normoxic and sediments are strongly bioturbated, meaning that the benthic system stays as healthy as possible during the duration of OM load reductions. If bottom waters become O₂-deficient and the bioturbating infauna dies, the duration to reach 1950s levels of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes will be extended by approximately 30 years, during which also the magnitude of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes will be higher than under normoxic conditions (Fig. 4b), thereby sustaining the strain on the ecosystem. Hence, the key factors affecting internal N loading via sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes in the coastal Baltic Sea are the external OM loading and the bottom water O₂ conditions. Consequently, internal N loading and its impact on the coastal environment cannot be mitigated if OM loading remains on its present-day level (Fig. 4a).

The limited direct impact of the sediment ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool on the sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes in our model is unexpected, but important knowledge for ecosystem management, as it suggests that ecosystem recovery in regard to N overloading is achievable within realistic time-frames. A stronger influence of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool on the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes than what we see in our modelling results would increase both the magnitude of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes and their resistance to external management measures, potentially extending their impact on the coastal ecosystem by many years. Further, it is worthwhile to keep in mind that, despite its limited direct impact, our modelled sediment ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool nevertheless maintains the long-lasting legacy effect of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes, whose magnitude and thus ecological impact can increase under O₂-deficient (Appendix B, Fig. B2) and OM enriched (Fig. 6a I) conditions. The estimated duration of legacy effects in the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes is similar to Lønborg and Markager (2021), who suggested a TN legacy in Baltic Sea sediments of several centuries.

Considering that our model is 1-dimensional and not directly coupled to the water column, we cannot estimate any direct feedback between sediment and water column. In a natural system, a large part of the N released to the water column via sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux would be taken up in primary production and eventually settle back to the sediment to become new OM, which adds complexity to the internal nutrient loops in both modelling scenarios and is not included in the current model. It is likely that the coastal current of the Baltic Sea would keep OM and nutrients for a certain time in the coastal zone before both end up in the open Baltic Sea (Radtke et al., 2012), which would increase the chance that the coastal N load is ultimately removed from the aquatic system via denitrification at the permanent redoxcline. How long it would take to remove the internal N loading from the coastal zone and thus from the potential to favour coastal primary production, remains, however, unclear. Furthermore, the impact of internal N loading on the coastal ecosystem may be localized, as despite a predicted rise in hypoxia it is highly unlikely that the entire coastal zone of the Baltic Sea will eventually become O₂-deficient. Localized hypoxic or anoxic events can lead to increases in N availability and subsequent primary production, yet, such increases are likely masked by the overall coastal eutrophication. Nevertheless, for management efforts to mitigate coastal eutrophication and restore a healthy ecosystem, every addition of bioavailable nutrients that is not accounted for in mass balance budgets should be noteworthy. In the big picture, internal loading adds nutrients to a system and even if these nutrients are transient, they can cause ecological harm and impede eutrophication mitigation measures.

Following OM load reductions, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes in our model decrease slightly faster and with a smaller flux magnitude under normoxic than under O₂-deficient conditions (Fig. 4b), which might raise the question whether artificial aeration methods could be a justifiable option to aid the ecosystem recovery. However, the effects of artificial aeration on benthic N-cycling in the coastal Baltic Sea are not straight forward and sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes both did (site Lännerstasundet, Lehtoranta et al., 2022) and did not (site By Fjord, De Brabandere et al., 2015) significantly decrease under aeration, or even increased (site Sandöfjärden, Lehtoranta et al., 2022) due to aeration induced sediment warming and a return of anoxia despite aeration. Similar conflicting results were found also outside the coastal Baltic Sea (e.g. Harris et al., 2015, Testa et al., 2025), including increased air-water methane exchange under conditions of aeration (Lapham et al., 2022), which suggests that engineered aeration approaches often have unexpected and unintentional impacts on both nutrient cycling and greenhouse gas emissions. Given these variable responses, we do not consider artificial aeration to be a clear-cut solution to shorten the duration of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes in the coastal Baltic Sea.

6 Conclusions

Muddy sediments in the coastal Baltic Sea commonly contain large ${\mathrm{NH}}_{\mathrm{4}}^{+}$ accumulations from past and ongoing eutrophication, which are favourably released to the water column under O₂-deficient conditions, resulting in internal N loading with the potential to facilitate coastal primary production. To better understand the fate and persistency of this ${\mathrm{NH}}_{\mathrm{4}}^{+}$ legacy and its impact on the coastal ecosystem, we modelled the past to future development of the sediment ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes in the coastal Baltic Sea under two contrasting OM loading scenarios and three O₂ conditions. Our model results clearly show a strong eutrophication legacy effect in the development of the sediment ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool, which affects the sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes for a long time, but with only a mild effect size. The key driver of the sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes is the external OM loading, which continuously adds to the OM- and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ accumulations in coastal surface sediments. This result strongly indicates, that present-day elevated sediment-water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes can be effectively reduced if strict OM load reductions are realized, such as proposed by the BSAP. Oxic bottom water and strong bioturbation aid the reduction of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes, as both aspects increase the sediment internal turnover of ${\mathrm{NH}}_{\mathrm{4}}^{+}$. However, these factors are secondary behind the reduction in OM loading: following our model results, internal N loading to the coastal Baltic ecosystem cannot be mitigated if OM loading continues as at present, which has also been suggested in other studies (Pitkänen et al., 2001, Saraiva et al., 2019). Our results are beneficial also for other coastal areas of similar geomorphology as the coastal Baltic Sea that face similar challenges in nutrient reduction and ecological management.

Appendix A:  

Table A1 Chemical reactions in the RTM; organic matter (OM) is of the form $\left({\mathrm{CH}}_{\mathrm{2}}\mathrm{O}{)}_a\right({\mathrm{NH}}_{\mathrm{4}}^{+}{)}_b({\mathrm{H}}_{\mathrm{3}}{\mathrm{PO}}_{\mathrm{4}}{)}_c$ with a=106, b=16, and c=1.

Partial differential equations

We established a 1-dimensional RTM to simulate the variations of OM and porewater ions given in Table A1. The RTM was achieved by the ReacTran package in R language (Soetaert and Meysman, 2012) via the following partial differential equations:

$$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }\mathit{\varphi }\cdot {\displaystyle \frac{\partial C}{\partial t}}={\displaystyle \frac{\partial \left(\mathit{\varphi }\cdot (D_s+D_b)\cdot \frac{\partial C}{\partial x}\right)}{\partial x}}-{\displaystyle \frac{\partial (\mathit{\varphi }\cdot v\cdot C)}{\partial x}}\\ {\displaystyle }& {\displaystyle }+\mathit{\varphi }\cdot \sum R_{\mathrm{1}}\end{array}\text{(A2)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }(\mathrm{1}-\mathit{\varphi })\cdot {\displaystyle \frac{\partial S}{\partial t}}={\displaystyle \frac{\partial \left((\mathrm{1}-\mathit{\varphi })\cdot D_b\cdot \frac{\partial S}{\partial x}\right)}{\partial x}}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\partial \left(\right(\mathrm{1}-\mathit{\varphi })\cdot w\cdot S)}{\partial x}}+(\mathrm{1}-\mathit{\varphi })\cdot \sum R_{\mathrm{2}}\end{array}\end{array}$$

where C denotes the dissolved ion concentration in porewater (mmol L⁻¹), S represents OM contents, x is depth (cm below sea floor), and ∑R₁ and ∑R₂ define the net reaction rates in porewater and dry sediments, respectively. More parameters and definitions are given in Table A2.

Table A2 Reaction and transport parameters used in the RTM, following (a) Boudreau, 1997, (b) Jilbert et al., 2021, and (c) Reed et al., 2011. Limiting concentrations double as inhibition coefficients in the Monod scheme (Boudreau, 1996).

Table A3 Rates of reactions given in Table A1; i stands for the reactive or less reactive OM pools.

RTM calibration and validation

The RTM is calibrated using the sediment property and geochemical porewater data of the coastal site Storfjärden, north-western Gulf of Finland, Baltic Sea (Jilbert et al., 2021). The included reactions and parameters are given in Tables A1–A3. The following Figs. A1–A3 show the model outputs fitted on the measured sediment profiles.

Figure A1 The fit of downcore porosity (blue circles = measured, red line = fitted), the calculation of burial rates of organic matter (OM) and porewater, and the effective ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) diffusion coefficient corrected by tortuosity.

Figure A2 The prescribed organic matter (OM) loading onto surface sediments (3–4 times BALTSEM output) in order to fit the measured OM profile of the coastal site Storfjärden by Jilbert et al. (2021), as shown in Fig. A3.

Figure A3 Sediment and porewater profiles of major model parameters, measured at the coastal site Storfjärden, north-western Gulf of Finland (blue circles; Jilbert et al., 2021), and modelled results of those parameters (red line); with OC = organic carbon content, TOC$/$TON = total organic carbon to nitrogen ratio, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ = ammonium, Fe²⁺ = iron(II), and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ = sulphate.

Appendix B:  

Figure B1 Zoom into the normoxic (blue; 200 µM O₂), hypoxic (green; 63 µM O₂), and anoxic (red; 0 µM O₂) ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes along time of the Pre-industrial scenario, showing efflux-“tails” starting from circa year 2160 when organic matter (OM) loading reaches pre-industrial baseline values. These ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux-“tails” are indicative of a eutrophication legacy effect. Due to the effect, true baseline values of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux are first reached after 2300, depending on the O₂ conditions. Grey shadings indicate the different phases of OM loading with (I) pre-eutrophication, (II) strongly increasing, (III) after-peak, and (VI) return to pre-industrial, in line with Fig. 4. The vertical dashed line gives the end point of OM load reductions if following the target of the Baltic Sea Action Plan (BSAP, Table 1).

Figure B2 Effect of changed organic matter (OM) loading (I), salinity (II), sedimentation rate (III), carbon to nitrogen ($\mathrm{C}/\mathrm{N}$) ratio in OM (IV), and different shares of reactive material in OM (V) compared to baseline conditions representative of the northern coastal Baltic Sea (salinity of 6, Pre-industrial OM loading scenario, sedimentation rate of 6 mm yr⁻¹, and OM with a $\mathrm{C}/\mathrm{N}$ ratio of 6.6 (reactive fraction) and 13 (less reactive fraction); black line) on (a) sediment ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) pool and efflux, and (b) its percent change relative to baseline condition. Modelled for anoxic bottom water conditions (0 µM O₂). Dashed lines give the pre-industrial baselines of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool and efflux. Note the different y-axis scales between some of the parameters.