Large-scale effects of benthic fauna on carbon, nitrogen, and phosphorus dynamics in the Baltic Sea

. Even though the effects of benthic fauna on aquatic biogeochemistry have been long recognized, few studies have addressed the combined effects of animal bioturbation and metabolism on ecosystem–level carbon and nutrient dynamics. Here we merge a model of benthic fauna (BMM) into a physical–biogeochemical ecosystem model (BALTSEM) to study the long–term and large–scale effects of benthic fauna on nutrient and carbon cycling in the Baltic Sea. We include both the direct 10 effects of faunal growth and metabolism and the indirect effects of its bioturbating activities on biogeochemical fluxes of and transformations between organic and inorganic forms of carbon (C), nitrogen (N), phosphorus (P) and oxygen (O). Analyses of simulation results from the Baltic Proper and Gulf of Riga indicate that benthic fauna makes up a small portion of seafloor organic stocks, but contributes considerably to benthic–pelagic fluxes of inorganic C, N and P through its metabolism. Results also suggest that the relative contribution of fauna to mineralisation of sediment organic matter increases with increasing 15 nutrient loads. Further, bioturbation decreases benthic denitrification and increases P retention in sediments, the latter having far–reaching consequences throughout the ecosystem. Reduced benthic–pelagic P fluxes lead to a reduction of N fixation and primary production, lower organic matter sedimentation fluxes and thereby generally lower benthic stocks and fluxes of C, N and P. This chain of indirect effects overrides the direct effects of faunal respiration, excretion and bioturbation. Due to large uncertainties related to parameterization of benthic processes, we consider this modelling study a first step towards 20 disentangling the complex large–scale effects of benthic fauna on biogeochemical represent hydrogen sulphide), and a vNOXY is a fitting constant. the N is oxidized. of the whole benthic fauna community over the whole Baltic Sea. Further, these studies only included the effects of bioturbation, but excluded metabolic fluxes. Here, we estimated that even though the excretion of DIP by benthic fauna constituted a significant proportion of benthic–pelagic fluxes, about ¼ to ⅓, it did not reverse the effect of bioturbation on P retention. These results are supported by Berezina et al. (2019), who found a positive correlation between macrofaunal biomass and P excretion, but a negative correlation between biomass and total sediment to water DIP flux in sediment cores with natural macrofauna communities in the eastern Gulf of Finland. redoxclines denitrified red oxcline we case


Introduction 25
Coastal ecosystems are highly productive, consist of diverse biological communities and carry out important functions including those supporting a growing world population (Costanza et al., 1997(Costanza et al., , 2014. However, they are facing multiple anthropogenic pressures such as nutrient loading and climate change (Cloern et al., 2016;Halpern et al., 2008). Elucidating the mechanisms of the coupled biogeochemical cycling of carbon (C), nitrogen (N) and phosphorous (P) in these systems is important to understand how they respond to current and future pressures , but also because they contribute to the regulation of 30 global climate and nutrient cycles by processing anthropogenic emissions from land before they reach the ocean (Ramesh et al., 2015;Regnier et al., 2013aRegnier et al., , 2013bSeitzinger, 1988).
In contrast to the deep open ocean, benthic-pelagic coupling plays a large role in biogeochemical cycling in coastal and estuarine ecosystems (Soetaert and Middelburg, 2009). Coastal sediments act as hotspots for organic matter degradation and permanent removal of elements from biological cycling through burial and denitrification (Asmala et al., 2017;Regnier et al., 35 2013a;Seitzinger, 1988). The bioturbating activities of benthic fauna alter the physical and chemical properties of surface sediments, which in turn strongly influence organic matter degradation processes and benthic-pelagic biogeochemical fluxes (Aller, 1982;Rhoads, 1974;Stief, 2013). Additionally, benthic fauna retain carbon and nutrients in its biomass and transform them between organic and inorganic forms through metabolic processes (Ehrnsten et al., 2020b and references therein;Herman et al., 1999;Josefson and Rasmussen, 2000). Together, these direct and indirect effects of benthic fau na have far-reaching 40 consequences for ecosystem functioning in the benthic and pelagic realms (Griffiths et al., 2017;Lohrer et al., 2004).
Ev en though the importance of benthic fauna for sediment biogeochemistry and benthic -pelagic fluxes has long been recognized (Rhoads, 1974), the combined effects of animal bioturbation and metabolism have seldom been studied together (Ehrnsten et al., 2020b;Middelburg, 2018;Snelgrove et al., 2018). Further, empirical studies of faunal effects often focus on temporally and spatially limited parts of the system, omitting important interactions and variability occurring in natural 45 ecosystems (Snelgrove et al., 2014). Here, we extend a physical-biogeochemical model of the Baltic Sea ecosystem (BALTSEM;Gustafsson et al., 2014;Savchuk et al., 2012) with benthic fauna components based on the Benthic Macrofauna Model (BMM; Ehrnsten et al., 2020a). We include both the direct feedbacks from animal growth and metabolism and the indirect effects of their bioturbating activities on biogeochemical cycling to evaluate their relative contributions.
We use the Baltic Sea as a model area for two reasons: (i) the shallow depth (mean depth 57 m) and enclosed geography with 50 a long water residence time (about 33 years) contribute to strong benthic-pelagic coupling (Snoeijs-Leijonmalm et al., 2017;Stigebrandt and Gustafsson, 2003), and (ii) the major features of biogeochemical cycling of C, N and P in the Baltic Sea are well known due to a wealth of oceanographic measurements and studies performed over the past century, making it an ideal system for process-based modelling (Eilola et al., 2011;Gustafsson et al., 2017;Wulff, 2009, 2001). However, the sediment pools and the role of sediment processes in benthic-pelagic exchange are not as well quantified as pelagic pools 55 https://doi.org/10.5194/bg-2022-31 Preprint. Discussion started: 1 February 2022 c Author(s) 2022. CC BY 4.0 License.
The Gulf of Riga is a semi-enclosed coastal bay with mean and maximu m depths of 23 and 51 m, respectively, and a salinity of 4-7 (Snoeijs-Leijonmalm et al., 2017). In contrast to the Baltic Proper, the Gulf of Riga is relatively well mixed and hypoxia only occurs intermittently under the summer thermocline (Kotta et al., 2008), accompanied by increased release of phosphate from the sediments (Eglite et al., 2014). When occurring more often in the recent decade, the intensity and extension of both sporadic hypoxia and phosphate release have somewhat increased (HELCOM, 2018;Stoicescu et al., 2021). Similarly, the 90 summer cyanobacteria blooms sporadically occurring in the Gulf of Riga before the 2010s (Kahru and Elmgren, 2014), regularly and extensively cover the gulf since 2015 (pers. comm. Mati Kahru).

Model description
The biogeochemical cycling in the BALTSEM model was extended to include benthic fauna. BALTSEM simulates physical circulation and biogeochemical transformations of C, N, P, O and Si in the Baltic Sea in response to climatic conditions and 95 nutrient inputs from rivers , point sources and atmospheric deposition. It describes the Baltic Sea as 13 horizontally homogenous boxes with a dynamic depth resolution of generally less than 1 m in the pelagial (Fig. 1). Sediments are represented as terraces at 1 m depth intervals with an area corresponding to the hypsography of each basin. The new benthic components were constructed from the carbon-based Benthic Macrofauna Model (BMM) described in Ehrnsten et al. (2019bEhrnsten et al. ( , 2019aEhrnsten et al. ( , 2020a extended to include nitrogen (N) and phosphorus (P) components. 100 Below we give a short description of the benthic dynamics in the new model version, referred to as BALTSEM-BMM, with a focus on the effects of benthic fauna on biogeochemical processes (Fig. 2). A full mathematical description of the benthic biogeochemical processes of the model is found in Appendix A. For a description of the pelagic biogeochemistry and physics, we refer the reader to Gustafsson et al. ( , 2014 and Savchuk et al. (2012). Additionally, all benthic and pelagic state variables are listed in Table A1. 105

Benthic fauna dynamics
BALTSEM-BMM includes C, N, and P contents in biomass of three functional groups of benthic fauna. The facultative deposit-/suspension-feeding bivalve Limecola balthica is a key species dominating the biomass of benthic communities in large parts of the sea and is therefore represented by its own state variable. The group "deposit-feeders" represents the amphipods Monoporeia affinis and Ponotporeia femorata, the invasive polychaetes Marenzelleria spp. and other macrofaunal 110 species dependent on surface sediment organic matter as their primary food source. The group "predators" represents species feeding on the two former groups, such as the isopod Saduria entomon, the polychaete Bylgides sarsi and the priapulid Halicryptus spinulosus.
The biomass of all groups of fauna are modelled as a dynamic mass balance between fluxes formed by food uptake, assimilation, respiration or excretion, and mortality. The formulations for dynamics of the functional groups and their food 115 banks were kept as in BMM (Ehrnsten et al., 2019a, 2020a as far as possible. The main change is the addition of N and P https://doi.org/10.5194/bg-2022-31 Preprint.  (Sterner and Elser, 2002), therefore the fauna was given a constant C:N:P ratio. The ratio was approximated based on measured ratios for the dominating species in the Baltic Sea (see Appendix A).
Food uptake is modelled as a function of food availability in C units. The uptake of N and P component s of a food source are 120 thereafter calculated proportionally to the C:N:P ratio of the food source. Part of the food is assimilated, with an assimila tion factor depending on the food source. The assimilation factors for C components were applied to N and P as well (Table A3).
The unassimilated part is released as faeces, adding organic matter to the sediment C, N and P pools (Fig. 2).
Respiration and excretion of inorganic C, N and P is divided into three parts: (1) a basal maintenance part related to biomas s; (2) a growth and activity part related to food uptake as a proxy for activity; and (3) excess excretion. As the stoichiometry of 125 assimilated food varies, excretion of excess elements is calculated dynamically to keep the fixed stoichiometry of the benth os.
Formulations are similar to those used for zooplankton in BALTSEM and for benthos in other ecosystem models (Ebenhöh et al., 1995;Spillman et al., 2008). Respiration and excretion fluxes add to the bottom water pools of dissolved inorganic carbon (DIC), nitrogen (NH, representing total ammonia) and phosphorus (PO, representing total phosphate). Respiration also consumes bottom water oxygen with a respiratory quotient of 1 mol O2: mol CO2 (Brey, 2001). 130

Sediment dynamics and bioturbation
As in the standard BALTSEM, sediment bioavailable C, N, P and Si are represented as vertically integrated concentrations in the biogeochemically active surface layer of unspecified thickness. The concentrations are modelled as a dynamic mass balance between fluxes formed by sedimentation, degradation and burial, extended by interactions with benthic fauna. Sediment C, N and P pools are further divided into three banks of different age to resolve the food limitation of benthic fauna (Fig. 2), while 135 benthic Si is represented as a single pool. Oxygen is not a state variable in sediments, but several benthic processes interact with simulated bottom water oxygen.
Bioturbation by benthic fauna, including sediment reworking and burrow ventilation, generally increases the oxygenation of sediments (Michaud et al., 2005;Volkenborn et al., 2012), promoting the binding of phosphate to iron oxides (P sequestration) and stimulating nitrogen oxidation (Ekeroth et al., 2016;Norkko et al., 2012;Renz and Forster, 2014). Similar to Isaev et al. 140 (2017), we use simple formulations for the effects of faunal activities on the oxygen-dependent processes of sediment nitrification/denitrification and P sequestration through a bioturbation enhancement factor Ebio. The formulation for Ebio was taken from Blackford (1997), using the feeding rate of fauna as a proxy of its bioturbation activity: where Emax is the maximu m enhancement, cfBFi is a contribution factor of functional group i. UBFiC is the carbon uptake rate of 145 group i and Kbio is a half saturation constant. As L. balthica is more sedentary than the other groups, a contribution factor of 0.5 was assigned to it and a factor of 1 to the two other groups (Ebenhöh et al., 1995;Gogina et al., 2017 and refences therein).
Within each of the three sediment banks, C, N and P components share the same source and sink processes. Sinking org anic matter is integrated into a bank of fresh organic matter available as food for deposit -feeders. This bank ages into a slightly older bank available as food for L. balthica only. The second bank ages into a third bank considered unavailable as food for 150 benthic fauna, but available for bacterial mineralization.
For each element, mineralization fluxes from the three sediment banks are combined into a total flux (ZSEDXtot, X = C, N, P).
For C, the total sediment mineralization flux directly adds to the pelagic DIC pool. N and P mineralization fluxes are further divided in the same way as in the standard BALTSEM (Gustafsson et al., 2014;Savchuk, 2002) with the addition of bioturbation effects . 155 Depending on oxygen concentrations in the bottom water layer and bioturbation intensity, mineralized sediment N is released to the water column as ammonia (ONH, Eq. 2) or oxidised N (ONO, NO representing NO2 and NO3, Eq. 3) or denitrified to N2 (WDeni, Eq. 4).

=
(2) In anoxic or nearly anoxic conditions, the majority of mineralized N is released to the water column as NH, defined by the proportion νNOXY (Eq. 5), where OXY is bottom water total oxygen concentration (that is allowed to become negative to represent hydrogen sulphide), and avNOXY is a fitting constant. Otherwise the mineralized N is oxidized.
Subsequently, the oxidized portion can be denitrified into N2 or released to the pelagic NO pool. Denitrification is treated as a permanent sink for N. The proportion released as NO (ηN) is positively related to oxygen according to Eq. (7), where qN, aN, bN and cN are fitting constants and Ebio is the bioturbation enhancement factor.
The fraction sequestered is positively related to oxygen and bioturbation and negatively related to salinity according to Eq.
(10), where qP, aP, bP, cP, dP, eP, fP and gP are fitting constants. The fraction has an upper limit of 1 (100% of mineralized P 175 sequestered), but can take on negative values, representing a release of previously sequestered P in severely hypoxic or anoxic conditions. The salinity dependence fSAL in the third term is used as a proxy for the higher availability of the phosphate-binding agents (e.g. iron and humic substances) in the fresher Gulf of Bothnia.

Simulations 180
The model was run over 1970-2020 forced with observed nutrient loads and actual weather conditions as described in Gustafsson et al. ( , 2017 with forcing time-series extended to 2020. Initial conditions in 1970 were based on observations for pelagic variables and hindcast simulations for benthic variables as described in  and Ehrnsten et al. (2020a). As the purpose of this study was to evaluate large-scale dynamics, results were aggregated as means and standard deviations of the last two decades (2000-2020) to capture differences in long-term averages while accounting for interannual 185 variations.
It is difficult to constrain the new parameters as well as to validate the large-scale dynamics against observations from the field or laboratory, which are usually made on much smaller temporal and spatial scales. Instead, we made a sensitivity analysis testing the effects of changing the parameter Emax in the range 0 to 0.6, where 0 represents no bioturbation and 0.6 is the theoretical maximu m value of Ebio, giving 100% P sequestration. Emax = 0.3 was used in the default model run. 190 Additionally, we estimate the contribution of bioturbation to benthic-pelagic nutrient fluxes by calculating the theoretical fluxes without bioturbation enhancement (i.e. with Ebio = 0) for each time-step, while running the default model with bioturbation. In contrast to the sensitivity analysis, this analysis shows the direct effects of bioturbation without accounting for indirect effects through the ecosystem.
Finally, to study the relationship between nutrient loads and the role of benthic fauna in biogeochemical cycling, we ran two 195 future scenarios 2021-2100 with either decreasing loads of N and P according to the Baltic Sea Action Plan (BSAP scenario, total loads to the Baltic Sea 739 kton N year -1 and 21 kton P year -1 ) or increasing loads corresponding to the highest recorded historical loads based on mean N and P loads of 1980-1990(HIGH load scenario, 1235 kton N year -1 and 69 kton P year -1 ).
For comparison, the average nutrient loads in 2000-2020 were 936 kton N year -1 and 36 kton P year -1 . The scenarios were combined with a statistical climate forcing representing no change in climate. Details of the scenarios can be found in Ehrnsten 200 et al. (2020). These results were also aggregated over the last two decades (2080-2100).

Validation
The BALTSEM-BMM behaves very similarly to the standard BALTSEM model, which has been extensively validated (Gustafsson et al., , 2014Savchuk et al., 2012) and shown to perform favourably in relation to similar Baltic Sea models 205 (Eilola et al., 2011;Meier et al., 2018). A comparison of the main pelagic state variables (salinity, temperature and concentrations of oxygen, NH, NO, and PO) to observations over time  and depth shows an overall relative bias of 1.40, while the relative bias of the standard BALTSEM is 1.41. The relative bias index compares model-data difference with variability in the data, giving an estimate of how well the model captures variab ility in nature on seasonal, annual and decadal scales (Savchuk et al., 2012). A detailed description and results of this analysis are found in Appendix B. although it should be noted that the spread of observed data is large. In the high-salinity Kattegat and Danish Straits near the entrance to Baltic Sea, the model is not applicable as the benthic biomass is dominated by groups not included in the present model, such as suspension-feeding bivalves and large echinoderms. Further details on this analysis are presented in Appendix C. 220 The databases used did not contain observations of benthic fauna from the Gulf of Riga, therefore we compared the simulated biomasses to estimates from the literature ( Table 1). The simulated biomass of benthic fauna in the gulf varied substantially over depth and time (29-284 g wwt m -2 ), which is supported by field observations . The simulated mean biomass is within the large range of estimates from the literature. Carman and Cederwall (2001) have estimated the amounts of C, N and P in Baltic Sea sediments based on core samples. It is 225 not straightforward to compare the total amounts to simulations, as the thickness of the simulated sediment is not defined.
However, the estimated C:N:P ratios can be more readily compared. In the Baltic Proper, the molar C:N:P ratio (calculated from their Table 11.4) was estimated to be 116:12:1 in the top centimetre of sediments and 137:14:1 in the top 5 cm, while the simulated ratio was 108:15:1. In the Gulf of Riga, the estimated C:N:P ratio was 73:7:1 in the top 1 cm and 83:8:1 in the top 5 cm, while the simulated ratio was 75:10:1.  -term (2000-2020) average benthic budgets of C, N and P are shown in Fig. 4 for the Baltic Proper and Gulf of Riga.

Long
The results for the Baltic Proper are restricted to the depth interval 0-90 m, as benthic fauna is practically absent in the oxygenpoor waters below this depth (Fig. 5a). Figure 5 also shows the long-term average depth distribution of the bioturbation factor Ebio in the two basins. 235 According to the default simulation, the benthic fauna made up a minor part of the benthic organic C, N and P stocks (1-4 %), but had a proportionally larger share in benthic-pelagic fluxes of DIC (23 % and 31 % in the Baltic Proper and Gulf of Riga, respectively), DIN (43 % and 51 %), and DIP (25 % and 34 %). The budgets also show that input of organic matter to the sediments was higher in the Gulf of Riga compared to the Baltic Proper, resulting in overall higher benthic stocks and benthicpelagic fluxes (Figs. 4 and 5). 240

Bioturbation effects on C, N and P dynamics
When accounting only for the direct effects of bioturbation, it increased NO outflux by 0.41 g N m -2 year -1 (+40 %) and decreased sediment denitrification by the same amount (-14 %) in the Baltic Proper (simulated average of 2000-2020, Fig.   6a). Similarly, P sequestration was increased and PO outflux decreased by 0.09 g P m -2 year -1 (+31 % and -15 %), respectively.
In the Gulf of Riga, both the absolute flux rates and the relative effects of bioturbation on them were larger than in the Baltic 245 Proper (Fig 6b).
When also accounting for the indirect effects of bioturbation in the sensitivity analysis, changes were more complex (Figs. 7-8). For example, comparing the default run (Emax = 0.3) to the run with no bioturbation (Emax = 0) in the Baltic Proper, denitrification was reduced by 0.44 g N m -2 year -1 but NO outflux increased by only 0.38 g N m -2 year -1 (Fig 7c), while P sequestration increased by 0.14 g P m -2 year -1 and PO outflux decreased by 0.03 g P m -2 year -1 (Fig. 7e). 250 In general, increasing bioturbation led to a decrease of most benthic stocks ( Fig. 9) and fluxes (Fig 7). This can be explained by the following chain of effects . Increased P sequestration in the sediments ( Fig. 7e-f) led to less pelagic DIP available for phytoplankton, especially cyanobacteria, growth and thereby lower N fixation and primary production ( Fig. 8) which in turn led to lower organic matter sedimentation rates, lower sediment stocks of organic matter and consequently lower rates of most sediment biogeochemical transformations and fluxes (Fig s. 7-9). Decreasing organic matter s edimentation also led to 255 decreased biomass of benthic fauna (primarily due to a reduction in L. balthica biomass) and excretion of DIN and DIP. An exception is the sediment P stock (Fig 9c), which increased with bioturbation despite decreased sedimentation of organic P ( Fig. 7f). This was due to the increased sequestration adding to the 3 rd sediment P bank.

Nutrient load scenarios
All results below are calculated from means of 2080-2100 for the BSAP and HIGH nutrient load scenarios and 2000-2020 for 260 the default model run in the Baltic Proper (0-90 m depth) and Gulf of Riga. With increasing nutrient loads, primary production and input of particulate organic matter (POM) to the sediments increase d ( Fig. 10a-c), resulting in an increase of most benthic stocks and fluxes. The biomass of benthic fauna responded more strongly to changing nutrient loads than the bioturbation enhancement coefficient linked to the feeding activities of fauna ( Fig. 10d-e).
With changing loads, the relative roles of faunal and microbial processes in the sediment changed (Figure 11). With increasing 265 loads, an increasing proportion of benthic-pelagic fluxes of inorganic nutrients originated from faunal metabolism. Expressed as percent of POM input to the sediments, the respiration and excretion of fauna was 12-13 % in the BSAP scenario and 23-24 % in the HIGH scenario in the Baltic Proper. In the Gulf of Riga, respiration and excretion was 23-24 % of POM input in the BSAP scenario and 35-37 % in the HIGH scenario. Correspondingly, the proportions of POM input released as dissolved inorganic substances resulting from microbial processes in the sediment (DIC, NO+NH and PO outflux in Fig. 11) were lower 270 in the BSAP than in the HIGH scenario.
The relative proportion of POP input sequestered shows a complex pattern in the Gulf of Riga ( Fig. 11f): the proportion of POP input sequestered was lower in the BSAP scenario compared to the default model due to less fauna and thereby less bioturbation. However, the proportion was also lower in the HIGH load scenario, as the increased occurrence of hypoxia ( Fig.   10f) counteracted the effects of increased bioturbation. In the Baltic Proper, relative P sequestration shows a decreasing pattern 275 with increasing loads driven by increasing occurrence of hypo xia (Fig. 11e).
Even though the total amount of POM input to the sediment increased with increasing nutrient loads, it constituted a decreasing proportion of primary production. In the BSAP scenario, almost half of the annual primary production reached the seafloor (48 % and 47 % in the Gulf of Riga and Baltic Proper, respectively) compared to 21 % and 27 % in the HIGH load scenario.
Thus, the proportion of primary production mineralized by fauna varied only slightly with nutrient load scenario because of 280 the opposite responses of sinking organic matter and fauna: 4.7 % (BSAP) to 4.1 % (HIGH) in the Baltic Proper and 10.7 % (HIGH) to 9.6 % (BSAP) in the Gulf of Riga.

Discussion
We have created a new tool to simulate the long-term and large-scale effects of benthic fauna on biogeochemical cycling in a coastal sea by fully merging two existing process-based models. First simulations with the new model indicate that the 285 benthic fauna makes up a small part of benthic organic stocks, but contribute s substantially to organic matter mineralizat io n and benthic-pelagic fluxes of inorganic C, N and P through its metabolism. Further, the stimulation of P binding in sediments by bioturbation significantly reduced N fixation and primary production in the simulations, indicating that benthic fauna can alleviate the 'vicious circle' of eutrophication.

Model performance 290
In general, the BALTSEM-BMM model reproduces the observed Baltic Sea-scale patterns of decreasing biomass of benthic fauna with latitude and depth reasonably well, as also shown for a previous one-way coupled model version (Ehrnsten et al., explain the latter (see Appendix C). Additionally, the model would need an addition of several groups of "megabenthos" (e.g. large echinoderms and suspension-feeding bivalves) to be applicable to the marine areas at the entrance to the Baltic Sea.
The simulated mean biomass of benthic fauna in Gulf of Riga was more than twice as high as the recent estimate by Gogina et al. (2016). Possible reasons for overestimation may be that the model does not take into account the limitations by mobile substrates and low salinity, especially in the southern part of the basin (Carman et al., 1996;Kotta et al., 2008). In this region, 300 a reduction of benthic biomass (e.g. of Monoporeia affinis and Limecola balthica) occurred in the 1990s for reasons not well understood (Kortsch et al., 2021). On the other hand, our biomass estimate is less than half of that by Kotta et al. (2008), assuming reported dry weight is 10 % of wet weight.
The modelled patterns in total benthic biomass are strongly driven by changes in L. balthica biomass in response to changes in food availability, leading to extinction of the group in deep waters and in the oligotrophic Bothnian Bay. These patterns are 305 strongly supported by observations. The BALTSEM model was neither improved nor worsened by the addition of benthic fauna, according to the performance analyses comparing pelagic variables to observations (Appendix B). This shows that increasing model complexity does not necessarily increase accuracy, especially when the functions and/or variables added are not well known (Ehrnsten et al., 2020b ;Levins, 1966). In general, though, the previous assessments of model performance showing that the model is able to reproduce 310 seasonal and long-term variations in biogeochemical variables and performs well in comparison to other Baltic Sea models, remain valid (Eilola et al., 2011;Gustafsson et al., , 2014Meier et al., 2018;Savchuk et al., 2012).
Unfortunately, we cannot properly validate the simulated sediment stocks or fluxes due to a lack of large-scale data and insufficient understanding of the multitude of mechanisms underlying the biogeochemical transformations and fluxes. Thus, the quantitative results should be viewed with caution. The main strength of this study is instead in the "what-if" analysis 315 showing the many interlinkages among C, N, P and O cycles and between benthic and pelagic processes.

Biogeochemical effects of benthic fauna
Similar to estimates made with previous uncoupled versio ns of the model (Ehrnsten et al., 2019a(Ehrnsten et al., , 2020a, the results of this study suggests that respiration by fauna constitutes a significant part of organic matter mineralization in sediments. The fauna mineralized about 8-17 g C m -2 year -1 or 22-31 % of POC input to the sediments in 2000-2020. This agrees well with previous 320 estimates from the Baltic Sea of 22-40 % (Ankar, 1977;Elmgren, 1984;Kuparinen et al., 1984). Similarly, Rodil et al. (2019)  production in shallow estuaries. In these deeper coastal areas, we estimate that the fauna mineralized 3-9 % of annual primary production in the Baltic Proper and 8-15 % in the Gulf of Riga in 1970-2020, with considerable interannual variations. 325 The sensitivity analysis showed a large effect of bioturbation on primary production levels mainly due to increased P retention.
When bioturbation increased P sequestration, this led to a weakening of the 'vicious cycle' in the Baltic Proper where less DIP in the water column led to less N fixation and organic matter production, which in turn led to less organic matter input to sediments, less heterotrophic oxygen consumption, les s hypoxia and thereby further increased P sequestration in the oxygenated sediments. Also in the Gulf of Riga, where hypoxia was rare, the bioturbation -induced reduction of pelagic DIP 330 had large effects on primary production and especially N-fixation. In the two runs with bioturbation there was no or very little pelagic DIP surplus available for the N-fixing phytoplankton group in contrast to the run without bioturbation, where N fixation added on average 0.97 g N m -2 year -1 or ca 17 000 tonnes N year -1 to the basin. It should be noted that these bioturbation effects are only valid where conditions are favourable for P binding to metal oxides, which occur primarily in freshwater and brackish sediments with a high iron and low sulphide content (Van Helmond et al., 2020;Jordan et al., 2008). 335 The significant effect of bioturbation on P retention found here is in line with the results of studies on the effects of the invasive polychaete Marenzelleria spp. in the Gulf of Finland and Stockholm archipelago in the northern Baltic Proper (Isaev et al., 2017;Norkko et al., 2012). However, both of these studies focussed on areas with very high abundances of Marenzelleria spp., and there are some indications that lower abundances may yield an opposite effect, i.e. increase P outflux from the sediment Nyström Sandman et al., 2018). While these studies concentrated on a single taxon in a limited area, we 340 included a dynamic representation of the whole benthic fauna community over the whole Baltic Sea. Further, these studies only included the effects of bioturbation, but excluded metabolic fluxes. Here, we estimated that even though the excretion of DIP by benthic fauna constituted a significant proportion of benthic-pelagic fluxes, about ¼ to ⅓, it did not reverse the effect of bioturbation on P retention. These results are supported by Berezina et al. (2019), who found a positive correlation between macrofaunal biomass and P excretion, but a negative correlation between biomass and total sediment to water DIP flux in 345 sediment cores with natural macrofauna communities in the eastern Gulf of Finla nd.
The effects of bioturbation on sediment N dynamics were less important for eutrophication processes than the effects on P dynamics in this study. We assumed a very simple process formulation, where bioturbation increases oxygen penetration depth in the sediments leading to a larger proportion of organic N mineralization in oxic environments, thus promoting outflux of nitrates over benthic denitrification. In reality, denitrification is a complex process depending on e.g. the 3D-structure of 350 redoxclines in the sediment. If the biogenic structures of tube-dwelling bio-irrigators increase the area of oxic-anoxic interface in the sediment, this can lead to the opposite effect where a larger proportion of nitrate is denitrified at the enlarged red oxclin e (Aller, 1988;Gilbert et al., 2003). However, we believe this to be a special case unlikely to dominate in the Baltic Sea. low irrigation activity (e.g. Arenicola marina), but a decrease in sediments with species common in the Baltic Sea (e.g. 355 Limecola balthica and Mya arenaria).
To better capture alterations in redoxclines, a depth-resolved sediment model with oxygen as a state variable would be needed.
We also recognize that many other possible effects of bioturbation, e.g. on burial (Josefson et al., 2002) and resuspension (Cozzoli et al., 2021) were not included. However, there is always a trade-off between model complexity and generality, with few models to date combining a depth-resolved sediment module together with a full pelagic model (Ehrnsten et al., 2020b ;360 Lessin et al., 2018). One of the main advantages of the BALTSEM model is that its simplicity and fast running time promotes

Benthic-pelagic coupling in a changing environment
The Gulf of Riga had higher simulated benthic stocks and fluxes than the Baltic Proper. This can partly be attributed to the slightly higher primary production (171±19 vs 159±25 g C m -2 year -1 in 2000-2020), but probably more importantly to the shallower mean depth causing a larger proportion of pelagic production to sink to the sediments before it is mineralized in t he water column. In the Gulf of Riga, on average one third (32 %) of primary production reached the bottom as POC during 375 2000-2020, compared to one fifth (21 %) in the Baltic Proper. Besides depth, the amount of organic matter export also depends on e.g. the type of plankton and temperature (Tamelander et al., 2017). Despite large interannual variations, there was a clear decreasing trend in the proportion of primary production exported to the seafloor of about 5 percentage points per decade in the Baltic Proper during 1970-2020 (R 2 = 0.68, F=103, p < 0.0001). This coincides with an increase in water temperature and shift in phytoplankton composition towards an increased 380 proportion of cyanobacteria, seen both in these simulations and in reality (Belkin, 2009;Kahru and Elmgren, 2014). A less clear decreasing trend of 2 percentage points per decade was simulated in the Gulf of Riga (R 2 =0.29, F=20, p < 0.0001).
The proportion of primary production arriving at the seafloor also varied with nutrient loads , constituting almost half of the annual primary production in the reduced load BSAP scenario (47-48 %) compared to 21-27% in the HIGH load scenario.
Simultaneously, the amount of benthic fauna decreased and mineralization of sediment organic matter became more dominated 385 https://doi.org/10.5194/bg-2022-31 Preprint. Discussion started: 1 February 2022 c Author(s) 2022. CC BY 4.0 License. by microbial processes with reduced loads. Thus, we can conclude that both the absolute amount and the relative proportion of POM input mineralized by macrofauna increased with nutrient loads, but in relation to primary production the proportion mineralized by fauna was almost independent of changes in loads because of the opposite responses of sinking organic matter and fauna.

Conclusions and outlook 390
Using a newly developed modelling tool, significant effects of benthic macrofauna on C, N and P cycling were simulated in the semi-enclosed brackish-water Baltic Sea, with impacts on the ecosystem from the extent of hypoxic bottoms to the rates of pelagic nitrogen fixation and primary production. The magnitude of effects was dependent on depth and productivity, as shown by the comparison of two basins and different nutrient load scenarios. Even though these large-scale simulations contain a large degree of uncertainty, they are an important complement to empirical studies, which for practical reasons can only 395 consider temporally and spatially limited parts of the system (Boyd et al., 2018;Snelgrove et al., 2014). Our results suggest that in addition to the much-studied bioturbation, the metabolism of benthic fauna s hould be given more attention in future studies as it may play a significant role in benthic mineralization of organic C, N and P in coastal seas and estuaries.
These simulations confirm the notion that benthic-pelagic coupling is strongest in shallow coas tal areas (Griffiths et al., 2017;Nixon, 1981), but also show that this relationship is modified by multiple physical and biological drivers, which may change 400 over time. Unravelling these interacting drivers and responses on a system scale is important to understand how coastal and global biogeochemical cycles are responding to changes in, e.g. nutrient loads and climate.
Code availability. The model code is available upon request from the corresponding author.

Appendi x A. Mathematical description of benthic model dynamics in BALTSEM-B MM
All state variables are listed in Table A1 and benthic parameters in Tables A2-A4. A graphical overview of benthic state 710 variables and processes is provided in the main manuscript (Fig . 2). Bioavailable surface sediments are represented as terraces at 1 m depth intervals with and area corresponding to the hypsography of each basin. All benthic state variables are calculated in mg m -2 for each terrace.

Sediment dynamics
Sediment concentrations of elements X (X = C, N, P) are divided into three banks that share the same processes (Eq. A1-A3). 715 Deposition of sinking detritus (DDETX) and phytoplankton (DPHYiX, i = phytoplankton group 1, 2, 3) on the sediment is integrated into a bank of fresh organic matter SED1X (Eq. A1). Loss terms of the bank are mineralization (ZSED1X), aging (ASED1X) into the second food bank SED2X and uptake by deposit-feeders (USED1X,BF2). Sources of SED2X include aging from bank 1 (ASED1X) and faeces from deposit-feeders (FBF2X) and predators (FBF3X) (Eq. A2). Loss terms are mineralization (ZSED2X), uptake by Limecola balthica (USED2X,BF1) and aging (ASED2X) into the third sediment bank SED3X, which is considered unavailable as 720 food for the benthic fauna, but available for bacterial mineralization. Mortality (MBFiX) of all benthic groups and faeces of L.
balthica are added to the third bank (Eq. A3). Loss terms are mineralization (ZSED2X) and burial (Bx).
Sediment silica is modelled as a single pool with sinking diatom Si (PHY2Si) as a source and mineralization and burial as sinks (Eq. A4).  For each element X (X = C, N, P), the mineralization fluxes from all three sediment banks are summed into a total flu x (ZSEDXtot), which is further divided in the same way as in the standard BALTSEM with the addition of bioturbation effects . For C, the total sediment mineralization flux goes to the pelagic dissolved inorganic carbon (DIC) pool (ODIC, Eq. A8).
In anoxic or nearly anoxic conditions, the mineralized N is released to the water column as NH, defined by the fraction vNOXY (Eq. A12), otherwise it is oxidized.
Subsequently, the oxidized fraction can be denitrified into N2 (treated as a permanent sink) or released to the pelagic NO pool.
The fraction released as NO (ηN,) is positively related to oxygen concentrations according to Eq. (A13), where qN, aN, bN and 750 cN are fitting constants (Savchuk, 2002) and Ebio is a bioturbation enhancement factor (see Eq. (A37) below).
A fraction ηP of mineralized sediment P is sequestered (KP, Eq. A14), i.e. bound to iron oxides in the sediment, while the rest is released to the pelagic phosphate (PO) pool (OPO, Eq. A15).
The fraction sequestered is positively related to oxygen concentrations (enhanced by bioturbation) and negatively related to salinity according to Eq. (A16), where qP, aP, bP, cP, dP, eP, fP and gP are fitting constants. The fraction has an upper limit of 1 (100% of mineralized P sequestered), but can take on negative values, representing a release of previously sequestered P in severely hypoxic or anoxic conditions. The salinity (SAL) dependence in the third term is used as a proxy for the higher 760 availability of iron in the low-saline Bothnian Bay.

Benthic fauna
The model includes three functional groups of benthic fauna: Limecola balthica (BF1), deposit-feeders (BF2) and predators (BF3) that share the same processes but with group-dependent parameterisations . The change in C biomass of functional group i (i = 1, 2, 3) is the difference between food uptake, faeces production, respiration, mortality and predation (Eq. A18), where 770 UBFiC,BF3 denotes predation on group i (i = 1, 2) by benthic predators.
Since the model does not include recruitment or migration, biomass change is set to 0 when biomass falls below 0.01 mg C m -2 to avoid permanent extinction of any group.
The change in N and P components of a group BFiX (i = 1, 2, 3; X = N, P) share the same processes as above, except that 775 respiratory release of C is replaced by excretion of N or P (Eq. A19). As the fauna has fixed stoichiometry, the dynamics can also be expressed by the change in C biomass and a conversion factor λCX.
Food uptake of element X (X = C, N, P) by group i is the sum of ingestion of food sources j, where IjX,BFi is ingestion rate of food source j by group i (Eq. A20). 780 Predators feed on deposit-feeders and L. balthica, with a strong preference for the former. The ingestion rate of multiple food sources is formulated according to Eq. (A21), where IjC,BF3 is the ingestion rate of food source j by predators, Imax,BF3 is the maximu m specific ingestion rate of predators, prj and prk are preference factors for food sources j and k, Flimj and Flimk and lower feeding limits on j and k and Km,BF3 is a half-saturation constant for predator ingestion rate. 785 Additionally, feeding stops at anoxia for all groups.
The ingestion rate of food source component X (X = N, P) is proportional to the C:X ratio of the food source j (Eq. A23).
Ingested food is divided into assimilated uptake (AUBFiX,, Eq. A24) and feces (FBFiX,, Eq. A25) by an assimilation factor AFj, which depends on the assumed nutritional quality of the food source j. 795 Respiratory release of DIC (RBFiC) is divided into three terms: basal respiration related to biomass and temperature, growth and activity respiration related to food uptake, and possible excess respiration to keep stoichiometry (RexBFiC) according to Eq.

= (A33) 815
Mortality is divided into hypoxia-induced mortality and other mortality. Other mortality rate mother,i is linear for L. balthica (Eq. A34) and quadratic for the other two groups (Eq. A35).
The hypoxia-induced mortality rate mox,BFi is dependent on bottom water oxygen concentration, temperature and the functional 820 group's sensitivity to hypoxia according to Eq. (A36), where m0,BFi is the mortality rate at anoxia and 10°C and Kox,BFi is a hypoxic sensitivity constant . and P sequestration (Eq. A16) through a bioturbation enhancement factor Ebio which mimics increased oxygen penetration into the sediment. Ebio uses the feeding rate of fauna as a proxy of their bioturbation activity (Blackford, 1997) according to Eq.
(A37), where Emax is the maximu m enhancement, cfBFi is a contribution factor of functional group i (i = 1, 2, 3), UBFiC is the carbon uptake rate of group i and Kbio is a half saturation constant. As L. balthica is more sedentary than the other groups, a contribution factor of 0.5 was assigned to it and a factor of 1 to the two other groups (Ebenhöh et al., 1995;Gogina et al., 2017 830 and refences therein).

Appendi x B. Validation of pelagic variables
As a measure of model performance, we calculate the relative bias of simulated and observed long -term monthly means for some of the main pelagic state variables as described in Savchuk et al. (2012), but with data extended to 2015. The index compares model-data difference with variability in the data, giving an estimate of how well the model captures variability in nature on seasonal, annual and decadal scales. 855

Methods
Observations of salinity, temperature, and concentrations of oxygen, ammonium, nitrate, phosphate and silicate were collected from the Baltic Sea Environment Database (BED) and other major data sources around the Baltic Sea such as IOW (German y ), NERI (Denmark), SYKE-FMI (Finland), and SHARK (SMHI, Sweden) databases. The full list of the data contributors can be found at http://nest.su.se/bed/ACKNOWLE.shtml. 860 Basin-wide monthly time-series were prepared from available long-term observations in the following way. All the measurements found in monthly intervals over 1970-2015 for all frequently sampled water layers within every BALTSEM basin, i.e. usually at 5 m intervals for the top 20 m of the water column and 10 -25 m intervals for the deeper parts of basins, were pooled together and averaged. Coastal measurements, defined as being sampled within 12 nautical miles from the shore, were excluded for all basins except the three Danish Straits basins, where the 12 nm coas tal strip covers almost the entire 865 basins. Measurements from several deep and isolated trenches in the northern Baltic Proper were excluded as they often display their own dynamics, asynchronous to that in the domain of the larger basin.
To emphasize both long-term changes and seasonality of variables, time series of a model-data difference of pairwise monthly means were used. Because the seasonal cycle is also reflected in monthly standard deviations, especially in the upper part of the water column these differences were scaled with month-specific standard deviation SDm. SDm was calculated as the standard 870 deviation of data collected in month m during the period 1970 -2015 for each available sampling depth. To remove any remaining outliers, the estimated monthly standard deviations were replaced by a spline smooth fitted by a GAM model. To avoid shifts due to some seasons being over-represented in the field data, in every basin the relative bias RBi at each sampling depth was calculated as an average of the twelve months in the annual cycle: