Understanding Environmental Changes in Temperate Coastal Seas: Linking Models of Benthic Fauna to Carbon and Nutrient Fluxes

Coastal seas are highly productive systems, providing an array of ecosystem services to humankind, such as processing of nutrient effluents from land and climate regulation. However, coastal ecosystems are threatened by human-induced pressures such as climate change and eutrophication. In the coastal zone, the fluxes and transformations of nutrients and carbon sustaining coastal ecosystem functions and services are strongly regulated by benthic biological and chemical processes. Thus, to understand and quantify how coastal ecosystems respond to environmental change, mechanistic modeling of benthic biogeochemical processes is required. Here, we discuss the present model capabilities to quantitatively describe how benthic fauna drives nutrient and carbon processing in the coastal zone. There are a multitude of modeling approaches of different complexity, but a thorough mechanistic description of benthic-pelagic processes is still hampered by a fundamental lack of scientific understanding of the diverse interactions between the physical, chemical and biological processes that drive biogeochemical fluxes in the coastal zone. Especially shallow systems with long water residence times are sensitive to the activities of benthic organisms. Hence, including and improving the description of benthic biomass and metabolism in sediment diagenetic as well as ecosystem models for such systems is essential to increase our understanding of their response to environmental changes and the role of coastal sediments in nutrient and carbon cycling. Major challenges and research priorities are (1) to couple the dynamics of zoobenthic biomass and metabolism to sediment reactive-transport in models, (2) to test and validate model formulations against real-world data to better incorporate the context-dependency of processes in heterogeneous coastal areas in models and (3) to capture the role of stochastic events.


INTRODUCTION
Coastal marine ecosystems have high ecological and societal values. They provide a range of services to society, such as fish and shellfish as livelihoods and recreational opportunities. Coastal systems also contribute to the regulation of climate and nutrient cycles, by efficiently processing anthropogenic emissions from land before they reach the ocean (Tappin, 2002;Laruelle et al., 2009;Regnier et al., 2013a,b;Ramesh et al., 2015). The high value of these ecosystem services is obvious considering that a large proportion of the world population lives close to the coast (Costanza et al., 1997(Costanza et al., , 2014. Currently, coastal seas around the world are undergoing major ecological changes driven by human-induced pressures, such as climate change, anthropogenic nutrient inputs, overfishing and the spread of invasive species (Halpern et al., 2008;Cloern et al., 2016). In many cases, the changes alter underlying ecological functions to such an extent that new states are achieved and that baselines are shifted (Lotze et al., 2006;Duarte et al., 2009;Kemp et al., 2009). Globally, eutrophication is one of the major environmental problems in coastal ecosystems. Over the last century riverine inputs of nitrogen (N) and phosphorus (P) to the oceans have increased from 19 to 37 Tg N yr −1 and from 2 to 4 Tg P yr −1 (Beusen et al., 2016). Regionally, these increases were even more substantial as observed in the United States, Europe and China; in the Baltic Sea N and P loads increased by roughly a factor of three and six, respectively . The riverine N flux has increased by an order of magnitude to coastal waters of China within thirty years, while P export has tripled between 1970 and 2000 (Qu and Kroeze, 2012;Cui et al., 2013).
Efforts to mitigate eutrophication through nutrient load reductions are hampered by the effects of climate change (Cloern et al., 2016). Changes in precipitation increase the runoff of N, P and carbon (C) from land, which together with warming and increased CO 2 dissolution alter the coupled marine nutrient and carbon cycles (Falkowski et al., 2000;Gruber and Galloway, 2008).
In contrast to the open ocean where biogeochemical cycling is largely dominated by pelagic processes driven primarily by ocean circulation, in the coastal zone, pelagic and benthic processes interact strongly and are driven by a complex and dynamic physical environment (Statham, 2012). Eutrophication in coastal areas leads to shifts toward rapidly growing opportunistic algae, and generally to a decline in benthic macrovegetation because of decreased light penetration, substrate change and more reducing sediments (Sand-Jensen and Borum, 1991;Cloern, 2001;Grall and Chauvaud, 2002). Increased production and warming waters have caused expanding hypoxia at the seafloor with a consequent loss of benthic fauna (Diaz and Rosenberg, 2008;Breitburg et al., 2018). Hypoxic systems tend to lose many long-lived higher organisms and biogeochemical cycles typically become dominated by benthic bacterial processes and rapid pelagic turnover (Griffiths et al., 2017). However, if hypoxia does not occur, benthic fauna tends to increase in biomass with eutrophication (Cederwall and Elmgren, 1980;Pearson and Rosenberg, 1987;Josefson and Rasmussen, 2000).
Changes in benthic biota have far-reaching impacts on biogeochemical cycles in the coastal zone and beyond. In the illuminated zone, benthic micro-and macrophytes mediate biogeochemical fluxes through primary production, nutrient storage and sediment stabilization and act as a habitat and food source for a variety of animals (Figure 1). Benthic animals contribute to biogeochemical transformations and fluxes between water and sediments both directly through their metabolism and indirectly by physically reworking the sediments and their porewaters and stimulating bacterial processes (Figure 1). Grazing on pelagic organic matter and biodeposition of feces and pseudofeces by suspension-feeding fauna increases organic matter sedimentation rates (Kautsky and Evans, 1987;Newell, 2004). In addition, nutrients and carbon are retained in biomass and transformed from organic to inorganic forms through metabolic processes (Herman et al., 1999;Josefson and Rasmussen, 2000;Ehrnsten et al., 2019b). Bioturbation, including sediment reworking and burrow ventilation activities (bioirrigation), redistributes particles and solutes within the sediment and enhances sediment-water fluxes of solutes (Kristensen et al., 2012;Villnäs et al., 2013). Bioturbation can also enhance resuspension of particles, a phenomenon termed 'bioresuspension' (Graf and Rosenberg, 1997;Le Hir et al., 2007). Together, all these processes affect physical and chemical conditions at the sediment-water interface (Rhoads, 1974;Aller, 1982;Lohrer et al., 2004;Stief, 2013), and strongly influence organic matter degradation (Aller, 1982;Furukawa, 2001;Canuel and Hardison, 2016). When up-scaled to the ecosystem level, such modified conditions can significantly alter the functioning of coastal ecosystems and ultimately, the role of the coastal zone in filtering and transforming nutrients and carbon.
Accordingly, quantitative insight in the presence and activities of benthic biota and the tight and complex mechanisms controlling benthic ecosystem functioning at various spatial and temporal scales is of primary importance for unraveling nutrient and carbon cycles in coastal areas. Globally, numerous studies have measured bioturbation processes (Wheatcroft et al., 1990;Maire et al., 2010) and assessed the relation between benthic fauna and biogeochemical cycling through manipulative experiments (Michaud et al., 2005;Gilbertson et al., 2012;Forster, 2013, 2014;Braeckman et al., 2014). However, the results of such observations and experiments are generally not fully up-scaled to field situations. Especially in dynamic and heterogeneous coastal seas, lab-controlled experiments only consider limited parts of natural systems (Snelgrove et al., 2014). One way to circumvent lack of system knowledge is to use empirical models (i.e., models that are built solely on empirical observations without explicit descriptions of the underlying mechanisms) to analyze the response of a system to environmental changes. These may provide realistic system-level integrated information based on large-scale in situ observations of e.g., biogeochemical fluxes. However, as empirical models do not describe the underlying mechanisms, their predictive power beyond the observations is uncertain. Therefore, mechanistic models, i.e., models that mathematically describe the individual processes within the modeled system, FIGURE 1 | Major processes related to vegetation and fauna controlling benthic biogeochemical fluxes. White arrows: solute fluxes, black arrows: particulate fluxes. Primary production: nutrient and CO 2 uptake and oxygen release (1), enhanced sedimentation and sediment stabilization by benthic primary producers (2), food uptake (3), egestion/biodeposition of feces (4), nutrient excretion and respiration (5), and bioturbation, including bioirrigation (6) and mixing of sediments (7). are needed to elucidate the effects of environmental change on ecosystem functioning.
However, in existing mechanistic models of marine ecosystems, the sediment is often simplified to a reactive boundary layer (Soetaert et al., 2000;Savchuk, 2002;Lessin et al., 2018), while models focusing on sediment biogeochemistry (early diagenesis) have traditionally considered the role of fauna in mixing and transport processes (Aller, 1980;Boudreau, 1984), but not its metabolism and biomass dynamics (Middelburg, 2018;Snelgrove et al., 2018). While these approaches may be adequate in deep basins where biogeochemistry is dominated by pelagic processes and sediment biological dynamics are relatively stable, these assumptions do not hold in shallow coastal ecosystems where the coupling between benthic and pelagic systems is strong and highly dynamic over time. For example, studies from a long-term monitoring station in the North Sea suggest that the majority of sinking organic matter is processed by benthic macrofauna before it reaches the microbial community (Tait et al., 2015;Zhang et al., 2015;Lessin et al., 2019). Moreover, macrofaunal respiration has been estimated to account for up to 45-70% of benthic carbon mineralization in coastal areas of the Baltic Sea and Black Sea, with seasonal variations spanning one order of magnitude (Wenzhofer et al., 2002;Rodil et al., 2019).
Evidently, in order to represent both the current state of coastal ecosystems and, especially possible future developments, there is a clear need to include benthic biological processes in mechanistic models. Several authors have reviewed the effects of benthic primary producers (Nielsen et al., 2004;Duarte et al., 2005;McGlathery et al., 2007;Krause-Jensen and Duarte, 2016) and bioturbating fauna (Meysman et al., 2006b;Mermillod-Blondin, 2011;Norkko and Shumway, 2011;Kristensen et al., 2012;Stief, 2013) on nutrient and carbon processing, but the effects of biomass production and metabolism of benthic fauna, and especially mechanistic models thereof, have received relatively little attention.
In this review, we discuss current approaches that are used to model benthic fauna and their role as drivers of biogeochemical processes in the coastal zone, with a focus on fluxes of carbon (C), nitrogen (N) and phosphorus (P). First, we provide an overview of existing mechanistic models of biomass and metabolism that can be used to simulate the quantity of benthic fauna and its direct effects on biogeochemical fluxes in the coastal zone. We then discuss how these models can be used to quantify physical effects on the sediment, including bioturbation and bioresuspension, which require a coupling between benthic biology and physical transport models. In the third section, we provide examples of ecosystem models that allow these processes to be combined and up-scaled. We also examine the capabilities and limitations of current models for simulating changes in the benthos and associated nutrient fluxes in coastal areas under environmental change, and give suggestions for future research directions. The opportunities and challenges for modeling nutrient and carbon cycling in estuaries and marine sediments in general have recently been reviewed (Arndt et al., 2013;Regnier et al., 2013a;Ganju et al., 2016;Lessin et al., 2018). To avoid repetition, we only focus on the opportunities and challenges when integrating models of benthic fauna and the biogeochemistry of coastal ecosystems. We do not attempt an extensive review of all the relevant literature, but rather discuss current capabilities and future challenges through representative examples. Most of the modeling studies to date focus on temperate systems, and this region is also the focus of the current review.

MODELING THE EFFECTS OF BENTHIC FAUNA ON BIOGEOCHEMICAL FLUXES
The impact of benthic fauna on nutrient and carbon cycling depends on e.g., density, biomass and the functional traits of individual species in the benthic community Norkko et al., 2013;Gammal et al., 2019), as well as on sediment structure (Braeckman et al., 2014;Bernard et al., 2019). Benthic animals process organic matter by consuming particulate organic matter (POM) and transforming a fraction of the ingested food into inorganic components that are respired and excreted (Figure 2). Organic matter is also incorporated in the biomass, transferred in the food web through predation and reproduction, and released back to the inanimate POM pool by egestion of feces and mortality. Detailed knowledge on the biomass and activity of benthic fauna is thus needed in order to quantify their effects on benthic nutrient fluxes. Below we go through different modeling approaches that couple biomass and metabolism of benthic fauna to biogeochemical processes in marine ecosystems.

Biomass Model Types and State Variables
The main sources and sinks of faunal biomass can be seen in Figure 2. The level of detail or aggregation of these source and sink processes varies among models, as seen in the examples in Table 1. The simplest models describe biomass of benthic fauna as a function of growth rate, mortality rate and environmental constraints, similar to conventional zooplankton formulations FIGURE 2 | Biomass source and sink processes of a typical benthic animal with links to C, N, P and O cycles. POM, particulate organic matter; DIN, DIP, dissolved inorganic N and P, respectively. in pelagic biogeochemical models (Murray and Parslow, 1999;Fulton et al., 2004a;Kim and Montagna, 2009;Montagna and Li, 2010). For example, Parslow (1997, 1999) model the change in biomass of benthic suspension-feeders (B, in mg N m −2 ) over time (t) as: where growth is a function of food uptake and a growth efficiency constant, and mortality is a function of biomass. While these kinds of models may include carbon and nutrient fluxes out of the organisms (cf. Figure 2) as derived variables, a more detailed and mechanistic approach is to model the metabolic fluxes explicitly. For example, in the benthic biological submodule of the European Regional Seas Ecosystem Model (ERSEM) (Baretta et al., 1995;Ebenhöh et al., 1995), the dynamics of a standard benthic organism representing a functional group (in units of C, N and P concentration) is expressed as a function of food uptake (U), egestion of feces (F), respiration of carbon (R) or excretion of nutrients (E), predation (P) and other mortality (M): This formulation is closely related to the ecological bioenergetics concepts developed in the International Biological Programme (Holme and McIntyre, 1971;Grodzinski et al., 1975), which also produced practical advice for the measurement of the component processes in benthic invertebrates (Holme and McIntyre, 1971;van der Meer et al., 2013). Similar formulations can also be found in Soetaert and Herman (1995), Cerco and Noel (2007), Hansen (2011), Timmermann et al. (2012) and Ehrnsten et al. (2019b).
In large-scale models like ERSEM, the macrofaunal biomass is often described by functional groups, representing feeding mode (e.g., suspension-feeders, deposit-feeders, predators), size (e.g., meiofauna, macrofauna) and/or position in relation to the sediment (e.g., infauna, epibenthos, hyperbenthos) (Soetaert and Herman, 1995;Murray and Parslow, 1999;HydroQual, 2000;Sohma et al., 2000Sohma et al., , 2018Fulton et al., 2004a,b;Maar and Hansen, 2011;Heath, 2012). Even though the dominant species belonging to a functional group may change over time and space, the rate of functional processes are generally more stable, making the functional group approach attractive when modeling systems on large scales (Vichi, 2002). However, the parameterization of such a model is not straightforward, as process rates are generally species and life-stage specific Norkko et al., 2013).
More detailed models follow the biomass or energy dynamics of different life-stages of a population. For example, the Scope for Growth-model (Winberg, 1960) describes physiological processes of a standard organism similar to Eq. 2, but applies an allometric scaling function to the processes for different size classes. Spillman et al. (2008) use this approach to model the biomass  Parslow, 1997, 1999;Fulton et al., 2004a Functional group C, N, P components Food uptake Respiration, excretion, egestion, predation, mortality Ebenhöh et al., 1995;Butenschön et al., 2016 Clam Tapes philippinarum size class C, N, P components Food uptake, seeding, recruitment Respiration, excretion, egestion, mortality, harvesting, recruitment Spillman et al., 2008 Blue mussel Mytilus edulis individual reserve energy, structural volume, reproductive tissue Food uptake, assimilation Respiration, excretion, egestion, mortality, spawning Maar et al., 2009 dynamics of three size classes j of the infaunal clam Tapes philippinarum as: where H is harvesting, Sd is seeding and Tr is recruitment between the size classes. Such structured population models avoid the error introduced by aggregating species and life-stages, and often include further processes, such as reproduction, which is usually omitted in functional group-based models. Structured population models are often used to study growth of farmed bivalves, but can also include feedbacks to biogeochemical cycles through ingestion of phytoplankton and detritus, excretion of dissolved nutrients and carbon, biodeposition of feces, and CO 2 production (Cranford et al., 2007;Grant et al., 2008;Spillman et al., 2008). A simpler model type based on the universal allometric relationship between organism size and metabolic rate (Gillooly et al., 2001;Brown et al., 2004b) are size-based models, where diverse communities can be aggregated into size classes. This type of model has been used to simulate present and future C biomass of benthic meio-and macrofauna on a global scale (Yool et al., 2017).
Probably the most detailed and mechanistic approach to simulate faunal biomass is that used in Dynamic Energy Budget (DEB) models (Kooijman, 2000;van der Meer, 2006;Augustine and Kooijman, 2019). DEB models describe the processes that distribute energy through individual organisms from the uptake and assimilation of food to the utilization for maintenance, growth, development and reproduction. The models include energy reserves and structural tissue (e.g., shells) as separate state variables, allowing for a mechanistic description of energy distribution and conservation and varying C:N:P ratios within the individual. There are a few examples of DEB models integrated into in biogeochemical-hydrodynamic models used to study the effects of bivalves on nutrient cycles in coastal areas (Maar et al., 2009;Grangeré et al., 2010;Ren et al., 2010;Saraiva et al., 2017), but their more extensive use in this context is probably hampered by their complexity (Brown et al., 2004a;Filgueira et al., 2011), and because they require parameters that are difficult to derive from commonly measured rates (van der Meer, 2006).
A part of biomass dynamics rarely resolved is the migration and recruitment of pelagic life-stages. This requires a coupling between benthic biology and pelagic hydrodynamics, as well resolving different life-stages, including continued post-larval dispersal that may dominate dispersal and hence the maintenance of several benthic species (Pilditch et al., 2015). Supply and successful establishment of pelagic larvae may be an important factor limiting species occurrence and biomass especially in enclosed coastal areas (Barnes, 1994;Palmer et al., 1996). Saraiva et al. (2017) showed the importance of larval recruitment and survival for blue mussels Mytilus edulis in the Wadden Sea using a DEB model coupled to a 3D hydrodynamic-biogeochemical model.

Particulate Fluxes: Food Uptake, Egestion and Mortality
Benthic fauna redistributes and transforms POM in the sediment through its feeding and metabolic processes. Additionally, suspension-feeders increase pelagic-benthic fluxes of POM by grazing on plankton and detritus in the water column and biodepositing unassimilated particles on the sediment surface (Graf and Rosenberg, 1997). Biodeposition can have substantial effects on benthic nutrient cycles, especially in bivalve farming areas (Chapelle et al., 2000;Cranford et al., 2007;Mesnage et al., 2007). The high organic carbon and nitrogen content in biodeposits stimulate deposit-feeding and bacterial remineralization in the sediment, leading to accelerated removal of nitrogen through denitrification and recycling due to dissimilatory nitrate reduction to ammonium (DNRA) (Henriksen et al., 1983;Kautsky and Evans, 1987;Norkko et al., 2001).
When linking models of faunal metabolism to organic matter fluxes, it is important not only to quantify the amount of organic matter that is ingested, retained in biomass and released back through egestion and mortality, but also to determine how the composition of this organic matter changes along the way. Organic matter degradability (lability/bio-availability) is one of the most important parameters when quantifying nutrient and carbon cycling, but constraining this parameter in sediment models is challenging Holstein and Wirtz, 2009;Lessin et al., 2018).
In models, food uptake is generally related to food concentration, and can be limited by various environmental variables such as temperature, salinity (Cerco and Noel, 2007;Spillman et al., 2008) and bottom water oxygen concentration (Blackford, 1997;Sohma et al., 2001;Maar et al., 2010;Timmermann et al., 2012) and, especially in the case of suspension-feeders, by crowding (Blackford, 1997;Fulton et al., 2004a). The Michaelis-Menten or type II functional response of asymptotic saturation (Holling, 1966) is a common expression for food uptake (or growth rate), but for suspension-feeders a linear type I response is often used (Dowd, 2005; see Jeschke et al., 2004 for a review). There are also various other functions in the literature, including type III sigmoid expressions (Fulton et al., 2003b), that are less frequently used.
The choice of response function is especially important to consider when models are used to extrapolate beyond present states. What may look like a linear or exponentially increasing response in the present ecosystem state (e.g., food uptake in relation to food availability or respiration in relation to temperature), may in fact reach saturation or even a decline at higher levels (e.g., as a consequence of eutrophication or climate change). Thus, choosing the appropriate parameterization requires knowledge of the shape of response curves for the whole range of simulated environmental conditions (Sharpe and DeMichele, 1977;van der Veer et al., 2006).
Only part of ingested food is assimilated, while the rest is egested as feces, adding to the pool of detrital POM. The production of feces (F) is often inversely related to an assimilation factor (AF), that reflects the assumed nutritional quality of the food (Soetaert et al., 1992;Ebenhöh et al., 1995;Timmermann et al., 2012): Although crucial, the assimilability of dead organic matter (detrital POM) is poorly restrained, with values for AF spanning from 0.15 (Blackford, 1997) to 0.8 (Heymans et al., 2016) used in different models. In a model of a deep-sea benthic food-web, van Oevelen et al. (2011) divided detrital organic carbon (POC) into labile, semi-labile and refractory fractions, with AFs depending both on the lability class of detritus and the group of fauna, based on an extensive literature review. They also assumed that egested POC 'moves down' one category, e.g., unassimilated labile detritus becomes semi-labile. A similar approach with two sediment food banks of different quality was used by Ehrnsten et al. (2019b), while Butenschön et al. (2016) used a dilution coefficient to account for the reduced nutrient content of feces compared to ingested material. In some DEB models, the assimilated food is assumed to have the same chemical composition as the reserve tissues in the animal, and thus the amount and composition of feces is dynamically modeled as the difference between ingested food composition and reserve composition (Saraiva et al., 2011).
Some suspension-feeding bivalves, primarily epifaunal mussels and oysters (Hawkins et al., 1998;Gallardi, 2014), have the capability to sort out nutrient-poor particles before ingestion and deposit them as mucus-bound pseudofeces. Pseudofeces production can be formulated as a simple threshold function: if food uptake exceeds a threshold value, the excess is released as pseudofeces (Klepper, 1989;Soetaert et al., 1992;Soetaert and Herman, 1995). More mechanistic approaches relate pseudofeces production to the amount of inorganic matter (Thomas and Bacher, 2018) or size and chemical composition of filtered particles (Saraiva et al., 2011).
In models where food assimilation is not explicitly modeled, the production of POM waste, including feces, pseudeofeces and loss through 'messy feeding' can be derived from uptake (U) and growth efficiency (GE) assuming a constant fraction (W POM ). Thus, for example, Murray and Parslow (1997) model the production of POM waste as: The formulation of biomass removal terms depend on the aims of the study, and can include predation from benthos or a coupled pelagic food-web Fulton et al., 2004b), harvest by humans (Fulton et al., 2004a;Spillman et al., 2007;Bossier et al., 2018) and 'natural mortality.' Additionally, cannibalism is not uncommon among the benthos, and especially in models where several species are aggregated into a functional group, a cannibalism term may be appropriate to represent predation within the group. This term also stabilizes the biomass dynamics in simple systems (Blackford, 1997), as does wide-spread omnivory. While predation and harvest generally remove biomass from a system, natural mortality, like feces, adds to the pool of POM that is readily available for bacterial decomposition. The different forces causing natural mortality are seldom modeled explicitly, but generally represented as a fixed fraction (M nat ) of biomass (B). Thus, the flux from mortality to POM (POM mort ) can be expressed as: A fraction of the dead tissue can be assumed to be directed to a pool of dissolved organic matter (Butenschön et al., 2016). In seasonally hypoxic environments, hypoxia-induced mortality can be an important structuring factor of benthic populations and included as a separate mortality term (Blackford, 1997;Bossier et al., 2018;HydroQual, 2000;Maar et al., 2010;Soetaert and Herman, 1995;Sohma et al., 2001Sohma et al., , 2018Timmermann et al., 2012). In estuarine areas, mortality due to salinity variations can also be important (Spillman et al., 2008). DEB models can also include starvation mortality when catabolic needs exceed reserve density (Maar et al., 2009;Saraiva et al., 2017).
Some predators or other mortality inducing factors may act outside the model boundaries. In such cases, model closure by mortality is commonly achieved through either a linear or quadratic loss term (Fulton et al., 2003b). In cases where different kinds of biomass removal are aggregated into one mortality term (M unresolved ), it is not straightforward to link this to production of detrital POM, as part of the flow may be directed to sinks outside the model system. In such cases, an additional term is needed, describing the fraction of mortality directed to detrital POM (M POM ):

Solute Fluxes: Respiration and Nutrient Excretion
Respiration by benthic fauna consumes oxygen and thereby affects the redox conditions of sediments and bottom water. This, in turn, can alter the pathways of microbial degradation and recycling of carbon and nutrients. The mineralization of organic carbon through food uptake and subsequent respiration can also have major effects on the carbon cycle. Macrofaunal respiration has been estimated to account for ca 10-60% of annual benthic carbon mineralization in the Baltic Sea area (Elmgren, 1984;Hansen and Bendtsen, 2013;Ehrnsten et al., 2019b;Rodil et al., 2019). In a mussel-dominated area on the northwestern shelf of the Black Sea, up to 70% of carbon mineralization was attributed to macrofaunal respiration (Wenzhofer et al., 2002). Respiration is commonly modeled as release of dissolved inorganic carbon, and oxygen consumption may then be calculated using a fixed stochiometric ratio Fulton et al., 2004b;Spillman et al., 2008;Butenschön et al., 2016). However, oxygen consumption can also be modeled explicitly (Ichikawa et al., 2007;Maar et al., 2009). Respiration is commonly divided into a basal rate (R b ) dependent on biomass and active respiration as a fraction of food uptake (R a ), as feeding can be expected to correlate to some degree with activity. For example, Maar and Hansen (2011) model respiration (R) as a function of functional group biomass (B) and assimilated uptake (U x AF): In structured population models, respiration is usually allometrically scaled to individual biomass (Grant et al., 2008;Spillman et al., 2008), while in DEB models respiration is related to structural volume and reserve density (van der Meer, 2006). DEB models also specifically resolve the respiration related to the breakdown and assimilation of food. In all types of models, (basal) respiration is generally dependent on temperature, and might also be related to other environmental constraints such as oxygen concentration (Blackford, 1997) and salinity (Spillman et al., 2008).
The metabolic processes of animals lead to waste products that are released through excretion. Excretion of dissolved nutrients by dense bivalve populations can increase pelagic concentrations of nitrogen and promote phytoplankton production (Asmus and Asmus, 1991;Dame and Libes, 1993;Prins and Smaal, 1994), but the balance between effects of grazing, denitrification and DNRA in biodeposits and nitrogen excretion on the total N budget is highly dependent on the local hydrological and biological setting (Smyth et al., 2018). A depleting effect on phytoplankton can generally be expected when the clearance time is approximately equal to phytoplankton turnover time and shorter than water residence time (Dame and Prins, 1998;Grall and Chauvaud, 2002).
The effect of infaunal excretion on benthic nutrient fluxes has been suggested to be negligible compared to microbial processes (Aller, 1978;Tuominen et al., 1999) and is often omitted in sediment diagenetic models. It is difficult to accurately measure infaunal excretion in situ, but several studies estimate considerable contributions to sediment-water fluxes.
For example, Kristensen (1984) estimated that a third of sediment-water NH 4 flux originated from excretion by the polychaete worm Nereis virens in a fjord system in the Baltic Sea, while Magni et al. (2000) estimated excretion rates of NH 4 and PO 4 by infaunal bivalves one order of magnitude higher than fluxes from defaunated sediments in a sandy tidal flat environment in Japan. They also estimated 10-to 100-fold seasonal variations in excretion fluxes, highlighting the importance of temporal variations in biologically mediated fluxes in coastal areas. In another study from the Baltic Sea, Norkko et al. (2013) showed that the outflux of NH 4 and PO 4 from shallow sandy sediments was directly proportional to the biomass of infaunal bivalves.
The chemical composition of an animal is generally fixed within certain limits. In models including both carbon and nutrient components in the biomass, a fixed carbon to nutrient ratio is often assumed and maintained through excretion of excess nutrients as dissolved ammonium NH 4 Maar and Hansen, 2011) or, more rarely, phosphate PO 4 (Vichi, 2002;Spillman et al., 2008). When food is poor in nutrients, excess carbon can instead be respired . Some models also assume that a proportion of ingested food is excreted as dissolved organic carbon (Spillman et al., 2008). In simpler model formulations where stoichiometry is not resolved, nutrient release cannot be dynamically modeled in this way. However, nutrient release can be formulated as a constant proportion of uptake or growth rate similar to Eq. 5 (Murray and Parslow, 1997;Fulton et al., 2004a).

Physical Effects on the Sediment
Bioturbation Benthic animals physically disturb the sediments in different ways (Figure 1). Their activities can alter the distribution of chemical substances and microorganisms and change redox conditions, and thus play a dominant role in biogeochemical fluxes across the sediment-water interface. The physical mixing of sediments by burrowing and by ingestion and defecation is generally referred to 'biodiffusion' , a random, diffusive-like process (Goldberg and Koide, 1962;Aller, 1982;Boudreau, 2000;Meysman et al., 2003). Mixing between non-adjacent parts of the sediments is called 'non-local mixing, ' but in the case of a diverse benthic community, this is commonly overshadowed by biodiffusion. Another activity is water exchange between the sediment and the overlying water through burrows or tubes built by benthic macrofauna, which is termed 'bioirrigation.' Bioirrigation can alter benthic fluxes by expanding the area of the sediment-water interface (Rhoads, 1974;Aller, 1988Aller, , 2001Kristensen, 1988Kristensen, , 2001Robbins, 1986;Boudreau, 1998). Bioturbation is defined as an umbrella term covering the physical effects and transport processes carried out by faunal activities, including biodiffusion and bioirrigation (Meysman et al., 2006b;Kristensen et al., 2012). However, it should be noted that the definitions of these terms vary between studies and fields.
Biodiffusion is widely modeled by assuming a single mixed surface layer of sediment with a biodiffusion coefficient (D b ) describing the intensity of biological mixing, and that below this layer, bioturbation is absent (D b = 0) (Sanford, 2008;Rooze et al., 2016). The mixed layer depth varies from site to site in marine sediments and the average depth is around 5 -10 cm, while the shallowest mixed layers in oxic coastal systems are found in the Baltic Sea, about 0.2 -2 cm (Boudreau, 1998;Teal et al., 2008). Also D b varies across marine sediments, with much higher values in coastal sediments (often on the order of 10-100 cm 2 yr −1 ), while in deep-sea sediments D b is as low as 0.01-1 cm 2 yr −1 (Boudreau, 1994;Tromp et al., 1995;Middelburg et al., 1997). Furthermore, Berg et al. (2001) concluded that the effects of bioturbation on solutes are up to 20 times stronger than that on solids in coastal sediments, suggesting that the D b for solutes and solids should be different. In general, biodiffusion is defined in reactive-transport models (RTMs) for simulating distributions of solutes and solids over time and depth in sediments as follows, Eq. 9 for solids and Eq. 10 for solutes: where C solid and C solute stand for the concentration of a solid and a solute, respectively, t is time, z is depth, ϕ is porosity, ω is sedimentation rate, D b is the biodiffusion coefficient, D s is the molecular diffusion coefficient and R solid and R solute are the reactions for a solid and a solute, respectively. The biodiffusion coefficient (D b ) is the predominant parameter in RTMs to represent biological activities in marine sediments. In most cases, the magnitude of D b is difficult to constrain, as the abundance and activity of benthic fauna are variable in time and space, particularly in coastal areas. To account for this variation without explicitly including fauna in the model, D b has been scaled with functions of various natural conditions as proxies for faunal activities. One common way is to scale D b with oxygen concentration in the sediment (Wang et al., 2003;Reed et al., 2011) or bottom water (Fossing et al., 2004), as benthic animals are aerobic. D b has also been scaled with temperature and with labile organic matter concentration as an indicator of food availability (Boudreau, 1998;Wang et al., 2003;Reed et al., 2011).
There are some examples of an explicit coupling of models of benthic fauna and RTMs. In the ERSEM model, the benthic system is modeled in three submodules dealing with (1) benthic biology, (2) bioturbation and (3) nutrient profiles in a multi-layer reactive-transport model Ruardij and Van Raaphorst, 1995). The submodules exchange information, for example the biomass of macrofauna affects the distribution of oxygen and detritus in the sediments through bioturbation, which in turn affects the growth of benthic fauna. Each functional group is assigned a coefficient of contribution to bioturbation, and D b is then defined as the weighted sum of the contribution of each functional group according to a Michaelis-Menten function (Supplementary Eq. S1). The same approach has been adopted in the ecosystem models IGBEM, Bay Model 2 and Atlantis (Fulton et al., 2004a,b;Audzijonyte et al., 2017). Several coastal ecosystem models developed for Japanese waters use a similar formulation, but without the weighing coefficients (Sohma et al., 2001(Sohma et al., , 2004(Sohma et al., , 2008(Sohma et al., , 2018. As the bioturbation activity is more strongly related to the activity of fauna than to its biomass, D b has also been related to the net carbon uptake flux in other versions of ERSEM, as feeding can be assumed to correlate with the general activity of the animals (Blackford, 1997;Butenschön et al., 2016). See Supplementary Material for a comparison of different ways to couple fauna to bioturbation in RTMs.
The first model formulation for bioirrigation was provided by Aller's radial-diffusion tube model (Aller, 1980), which approximates the bioirrigated zone in sediments as a series of evenly distributed and closely packed hollow cylinders of equal length, with burrows represented by the hollow portion of the cylinders. For such bioirrigated sediments, Eq. 10 is rewritten as: where r is the radial distance between cylinders. The boundary conditions are time-dependent, as the animals ventilate their burrows periodically, not constantly. This approach has been shown to reproduce porewater profiles in the bioirrigated zone well (Aller, 1982;Boudreau and Marinelli, 1994;Boudreau, 2000;Gilbert et al., 2003). Shull et al. (2009) extended the model to include burrow ventilation by benthic fauna through injection of bottom waters, suggesting that burrow ventilation strongly influences the vertical profiles of porewater and fluxes. In permeable (sandy) sediments where irrigation flows penetrate the surrounding sediments, advection and not diffusion can be the dominant mode of transport. In order to quantify the advective transport induced by bioirrigation, 1D and 2D advection models have been developed and used to study the effects of advective transport on e.g., inert tracers and oxygen (Timmermann et al., 2002(Timmermann et al., , 2006Meysman et al., 2006a). Another way to model bioirrigation was developed by Boudreau (1984) and Emerson et al. (1984). Here, bioirrigation is expressed by a bioirrigation coefficient (α s ), which describes the exchange rate between bottom water and porewater at a certain depth in the irrigated sediments. Eq. 11 can be reformulated as: where C 0 is the concentration at the sediment surface. Although bioirrigation actually occurs by a complex network of sediment burrows and their 3D geometrical sources, Eq. 12 parameterizes the process by assuming exchange between bottom waters and the average porewater at any given sediment depth. Nowadays, this model has become the most common way to express bioirrigation intensity, and thus how to define α s is crucial for modeling bioirrigation properly. Studies have shown that α s depends on the depth distribution of burrows both in laboratory (Furukawa, 2001) and field conditions (Schlüter et al., 2000;Norkko et al., 2012;Dale et al., 2013). Meile et al. (2005) have reported that α s is highly chemical species-dependent, i.e., different solutes should be assigned their own α s rather than sharing the same value. Studies sometimes use 3D models to represent explicitly the burrow networks that are constructed by resident macrofauna populations, but also in these models the bioirrigation intensity is based on α s as a function of depth and the total burrow volume (Koretsky et al., 2002). Isaev et al. (2017) also considered burrows as an expansion of the sediment surface area and made bioirrigation effects dependent on the total surface area and on oxygen concentrations as they decrease over sediment depth. So far, dynamic models of the distribution or activities of bioirrigating animals have not been coupled to these models, and the number and shape of tubes are generally fixed parameters.

Bioresuspension
In addition to effects on processes within the sediment, the activities of benthic fauna can have major impacts on sediment stability and resuspension (Graf and Rosenberg, 1997;Le Hir et al., 2007;Joensuu et al., 2018). Burrows, mounds and tracks produced by benthic fauna affect bottom roughness and erodibility of sediments. The increased sediment resuspension caused directly by faunal activities and indirectly by alteration of sediment properties has been termed 'bioresuspension' (Graf and Rosenberg, 1997). As erosion is fundamentally linked to lateral hydrodynamic processes, it is rarely included in one-dimensional (vertical) sediment reactive-transport models. Additionally, the same animal can have both stabilizing and destabilizing effects on the sediment, making it difficult to include such effects in models (Le Hir et al., 2007;Cozzoli et al., 2019). However, there are several examples of models coupling resuspension of sediment and benthic microalgae to fauna, including effects of infaunal bivalves (Willows et al., 1998;Wood and Widdows, 2002;Orvain et al., 2012;Rakotomalala et al., 2015) and epifaunal snail tracks (Orvain et al., 2003(Orvain et al., , 2012. For example, Wood and Widdows (2002) simulated effects of Limecola balthica on sediment erosion in an intertidal area of the Humber estuary, United Kingdom, and showed that bioresuspension by the clam was of the same order of magnitude as effects of tidal currents. This clearly shows that there is a need to account for biological effects in sediment transport modeling, especially in intertidal areas with dense biological communities.
The studies mentioned above have all used semi-empirical formulations for the relationship between animal density and sediment erosion based on flume studies. van Prooijen et al.
(2011) developed a more process-based model for cohesive (muddy) sediment resuspension by deposit-feeding activities of L. balthica based on the individual rate, area and maximum depth of feeding and density of individuals. A recent study by Cozzoli et al. (2019) proposed another process-based model for cohesive sediment resuspension mediated by a variety of bioturbators and demonstrated metabolic rate as a general indicator of impacts of biological activities on sediment resuspension. They concluded that the overall metabolic rate of the bioturbators was well correlated with changes in hydrodynamic energy gradients of suspended cohesive sediments, and a better predictor than size, density or biomass of the fauna. However, the metabolic rates were not dynamically modeled, but derived from species identity and size of the animals based on an empirical model (Brey, 2001).
The presence of live and dead bivalves increases the structural heterogeneity of the sediment surface (Norkko and Shumway, 2011), affecting the near-bottom flow regime and providing an extended, hard surface area for colonization of nitrifying and denitrifying microbes as well as algae (Commito and Dankers, 2001;Stief, 2013). However, these processes have not been well quantified and to our knowledge there are no models of their effects on biogeochemical fluxes.

Scaling Up to the Ecosystem
As shown by the examples above, there is a large array of models simulating the effects of benthic fauna on biogeochemistry.
The main processes and model formulations discussed are summarized in Table 2. A major challenge is the combination and upscaling of these individual processes to understand system-level effects. Ecosystem models can be used to evaluate the joint dynamics of biological, chemical and physical processes of coastal systems and their response to environmental change. To capture and be able to predict effects of these often non-linear changes (Cloern, 2001;Kemp et al., 2005;Hewitt et al., 2016) is a major challenge in mechanistic modeling. Below we go through illustrative examples of ecosystem models used to simulate interactions between nutrient loading and benthic ecology and geochemistry in coastal areas. We include examples of both natural assemblages of benthic fauna and assemblages influenced by human intervention (e.g., restoration, artificial habitats), but exclude a wealth of models dealing with seafood farming (see section "Biomass Model Types and State Variables" for some examples).
The development of marine ecosystem models stems from Riley's simple nutrient-phytoplankton-zooplankton model (Riley, 1946;Riley et al., 1949). The Ems estuary model was the first to include carbon dynamics of benthic functional groups about four decades later (Baretta and Ruardij, 1988). In the three decades since, ecosystem models have developed to include e.g., nutrient and oxygen dynamics of benthic fauna and vegetation in a range of coastal ecosystems from shelf seas  to shallow lagoons (Plus et al., 2003) and tidal flats (Sohma et al., 2000(Sohma et al., , 2008, and these models vary in complexity from simple box models Biodiffusion (D b ) dynamic function of: • Proxies for faunal activity: oxygen concentration, POC content Wang et al., 2003 • Biomass of fauna Ebenhöh et al., 1995 • Uptake rate of carbon by fauna Blackford, 1997 Bioirrigation, expressed by • A series of evenly distributed and closely packed hollow cylinders of equal length Aller, 1982 • The exchange rate between bottom water and porewater at a certain depth (α s ) Boudreau, 1984 Bioresuspension Related to abundance of fauna Wood and Widdows, 2002 Related to metabolic rate of fauna Cozzoli et al., 2019(Triantafyllou et al., 2000Kim and Montagna, 2009) to 3D-applications with tens of state variables, such as the Chesapeake Bay Environmental Model Package (Cerco et al., 2010; and Atlantis (Fulton et al., 2011;Bossier et al., 2018). A prominent example of a comprehensive model coupling benthic biology to biogeochemistry is ERSEM (Baretta et al., 1995;Butenschön et al., 2016). It explicitly simulates the dynamics of several chemical elements (C, N, P, Si, Fe and O) and living functional groups in the pelagic and benthic subsystems and can be coupled to different hydrodynamic transport models. Recent versions include the options to use either a singleor a multi-layer benthic component. The feeding and feces egestion depth is defined separately for each group of fauna, while faunal mineralization products (CO 2 , NH 4 and PO 4 ) are assumed to be released in the uppermost oxic layer (Ruardij and Van Raaphorst, 1995). The model has been applied to study carbon and nutrient fluxes in coastal areas, including the effects of benthic fauna and bacteria (Petihakis et al., 1999(Petihakis et al., , 2005Triantafyllou et al., 2000;Mussap and Zavatarelli, 2017), benthic microalgae (Blackford, 2002), as well as seagrasses and macroalgae (Aveytua-Alcázar et al., 2008). Triantafyllou et al. (2000) used ERSEM to simulate the effects of increased nutrient input to the shallow, temperate Gialova lagoon in the Mediterranean and showed a significant increase in suspension-feeders leading to increased biodeposition and a switch from nutrient limitation to grazing control for phytoplankton. While the model could capture the main characteristics of the changing system as observed in the field, it failed to match the great variability in measured biomass and nutrient fluxes (Petihakis et al., 1999;Triantafyllou et al., 2000), illustrating the difficulty of modeling shallow, dynamic systems. The authors mention the inclusion of benthic primary producers and effects of hypoxia as major future development needs. Blackford (2002) illustrated the advantages of coupling dynamics of benthic vegetation and fauna for a more realistic representation of shallow productive systems. A simulation of benthic and pelagic processes along a depth gradient in the Adriatic Sea showed significant interactions between benthic primary producers and consumers: the presence of benthic microalgae increased the biomass of bacteria and fauna by increasing sediment oxygenation and providing a food source to grazers, which in turn formed an important limiting factor for microalgal biomass. Together these processes enhanced nitrogen regeneration in the photic zone. The study is also a good example of extensive model validation of several important aspects, including annual means and seasonal cycles of biomass stocks, biogeochemical fluxes as well as measures of net ecosystem metabolism (i.e., the relation between primary production and ecosystem respiration). In general, model outputs were within the range of observed values.
The coupling of physical, chemical and biological processes and detailed descriptions of several benthic functional groups in ERSEM allows for comprehensive testing of the joint effects of benthic biological dynamics on biogeochemical processes. Yet, the model requires estimations of a large number of parameters, which are usually site-specific and laborious to measure. The advantages of model simplification are many including lower computational needs, fewer uncertain parameters and increased interpretability (Fulton et al., 2003a;Fulton, 2010;Ganju et al., 2016). Therefore, many ecosystem models use simpler approaches to describe the biota. For example, Bay Model 2 (Fulton et al., 2004a), a predecessor of Atlantis (Fulton et al., 2011;Bossier et al., 2018), explicitly includes only dynamics of the main limiting nutrient, N, and uses a simple formulation for biomass dynamics of functional groups (see section "Biomass Model Types and State Variables", Table 1). Fulton et al. (2004a) applied Bay Model 2 to Port Phillip Bay, Australia including six functional groups of benthic fauna and three groups of benthic primary producers. Despite the simple description of benthic biology, the model was able to predict emergent complex biological interactions such as competitive exclusion of benthic microalgae by macrophytes and predator-prey interactions, leading to a spatial differentiation of biological communities within the bay. However, the predicted biomasses of benthic microalgae and macrofauna for a range of temperate bays were significantly different from measured values, showing that some constraints were missing or inaccurate in the benthic compartment. The model was also used to study the effects of differing nutrient loads on the ecosystem by applying two nitrogen load scenarios to the same model bay, one corresponding to the eutrophic Chesapeake Bay and the other to the mesotrophic Port Phillip Bay (Fulton et al., 2004a). The model did remarkably well in predicting biological changes due to eutrophication, such as simplification of communities, extinction of seagrasses, a shift from a primary producerbased to a more detritus-based food chain and an increase in small-sized functional groups, corresponding to general theory and field observations from the two bays. Unfortunately, the feedbacks from changing benthic communities on nutrient cycles were not discussed.
Eutrophication and effects of benthic biology on nutrient cycles in Chesapeake Bay were simulated in several studies using the Chesapeake Bay Environmental Model Package, which couples pelagic biogeochemistry and hydrodynamics in 3D with a sediment diagenetic model and several biological components including two functional groups of benthic fauna (HydroQual, 2000;Noel, 2007, 2013;Cerco et al., 2010. Cerco and Noel (2007) explored the potential of oyster restoration as a means to reduce eutrophication effects in the bay. Even though increased suspension-feeder biomass led to increased respiration, DIN excretion and biodeposition, the grazing on plankton simultaneously decreased sedimentation rates, leading to overall lower sediment oxygen consumption. Decreased water turbidity increased biomass of seagrasses and benthic microalgae, further increasing sediment oxygenation. The overall result was decreased sediment-water ammonium fluxes and increased N loss through denitrification. Although this hypothetical simulation cannot be validated against measured data, it shows the potential of complex models to explore largescale, non-linear changes on the coastal ecosystem level.
Also Rasmussen et al. (2009) simulated positive effects of filterfeeding bivalves on restoration of the eutrophied Ghar el Melh lagoon in the southern Mediterranean using a nutrient-planktondetritus model including benthic vegetation and 2D transport of sediments and suspended matter. Simulations with up to 75% reductions in nutrient loads showed a decrease in macroalgae, but only marginal effects on rooted vegetation (Ruppia spp.), which were still limited by shading from resuspended sediments. Introducing filter-feeding by bivalves as a forcing function increased water clarity and allowed for re-establishment of Ruppia spp., but also increased the accumulation of N, P and organic C in the sediments due to biodeposition. The study also nicely illustrates how ecosystem models can be used to build nutrient budgets of coastal ecosystems to elucidate their overall role as a source or sink for nutrient inputs to the open sea.
Similar effects of suspension-feeders were simulated 'the other way around' in the northern Mediterranean Gulf of Trieste using a 1D version of ERSEM (BFM-POM), including five groups of benthic fauna and two groups of bacteria (Mussap and Zavatarelli, 2017). Removal of suspension-feeders increased dissolved P outflux from sediments, fueling increased phytoplankton production and increasing POM content in the water column, while less POM entered the sediments in the absence of suspension-feeders. In another application of ERSEM in the North Sea, Lessin et al. (2016) concluded that increased benthic nutrient outfluxes following a reduction in benthic macrofauna were due to increased activity of bacteria and meiobenthos. Maar et al. (2009) coupled a DEB model of blue mussels to a small-scale 3D hydrodynamic-biogeochemical model to study the effects of mussel beds on a wind farm foundation in a shallow coastal ecosystem in Limfjorden in northern Denmark. Simulations showed that although bivalves increased organic carbon sedimentation rates in the immediate surrounding, the overall effect was reduced sedimentation due to depletion of phyto-and zooplankton by grazing. Another example of a DEB model coupled to a pelagic ecosystem model is given by Saraiva et al. (2017), who simulated the detailed distribution and population dynamics of blue mussels in the shallow Balgzand area of the Wadden Sea. They concluded that the area acts as a sink for phytoplankton due to filtration by the bivalves, and as a source of ammonia, exporting about 40% more than the input flux from land and surrounding waters.
The diverse effects of suspension-feeders on sedimentation fluxes simulated in the different studies may stem from differences in model formulations and assumptions, such as the representation of physical transport in one, two or three dimensions, but it is also highly likely that the results reflect the true complexity and context-dependency of coastal biological and biogeochemical dynamics.
Evidently, mechanistic models are useful tools to quantitatively study the biogeochemistry of coastal ecosystems and to tackle managements questions such as the relationships between nutrient inputs, human interventions (e.g., restoration, constructions, shellfish farming) and the nutrient filtering capacity or ecological state of the ecosystem. Such models also have great potential to investigate the influence of climate change on coastal benthos (Thomas and Bacher, 2018;Ehrnsten et al., 2019a) and the role of coastal ecosystems as sources or sinks for greenhouse gases (Sohma et al., 2018). Isaev et al. (2017) also demonstrate how an ecosystem model can be used to study the combined impact of climate change and invasion of a bioturbating species on coastal nutrient fluxes. In a rapidly changing world, mechanistic models will become increasingly important for management, through their capacity to explore outcomes of different environmental change scenarios. They may also help pin-point variables that can indicate when an ecosystem is nearing a tipping point to an alternate state.

Coupling Benthic Fauna and Benthic Biogeochemistry
Models focusing on benthic biology on the one hand and benthic biogeochemistry on the other, have been developed in parallel in different scientific communities with largely diverging foci and assumptions that may be difficult to reconcile (Meysman et al., 2005;Middelburg, 2018;Snelgrove et al., 2018). However, there is a growing number of model studies coupling benthic fauna and biogeochemistry, although so far primarily focusing on coupling bivalve models to pelagic processes (e.g., nutrient-planktondetritus dynamics) while the coupling of benthic biogeochemistry to biology is less mature.
Many sediment biogeochemical models simplify benthic biology by omitting fauna and vegetation as state variables, and only model their effects implicitly, e.g., as a bioturbation and/or bioirrigation coefficient in reactive-transport models (RTMs). This implicit approach has been successfully applied to systems with relatively simple and stable benthic communities (e.g., Soetaert et al., 1996Soetaert et al., , 2000Paraska et al., 2014), but has been insufficient to explain nutrient fluxes in some shallow, near-shore areas (Murray and Parslow, 1997;Holstein and Wirtz, 2009), suggesting that the temporal variability in fauna and flora and associated processes typical for many shallow coastal areas needs to be accounted for by dynamic state variables. However, linking biodiffusion (D b ) to biomass dynamics of functional groups Sohma et al., 2001Sohma et al., , 2004Sohma et al., , 2008Sohma et al., , 2018Vichi, 2002;Fulton et al., 2004a,b;Audzijonyte et al., 2017) may not be enough to account for the great variability in reworking rates measured between coastal areas (Boudreau, 1994;Tromp et al., 1995;Middelburg et al., 1997), and even between replicates from the same site (Morys et al., 2016).
Ecological studies have developed indices of bioturbation and bioirrigation potential based on abundance and traits such as size, mobility and feeding type of the individual animals observed in natural communities (Solan et al., 2004;Wrede et al., 2018), which correlate to some degree with reworking rates, mixing depth, oxygen penetration depth and ammonium efflux (Braeckman et al., 2014;Gogina et al., 2017). Even though this kind of detail is not feasible in mechanistic modeling, these studies may aid improved model parameterizations. Traitbased descriptions of pelagic groups are emerging in ecosystem models (Follows et al., 2007;Andersen et al., 2015). In general, metabolism seems to be a better predictor of biogenic mixing and resuspension than biomass, advocating the modeling approaches of Blackford (1997) and Cozzoli et al. (2019). As metabolic rate is increasingly recognized as master parameter for the 'pace of life' (Gillooly et al., 2001;Brown et al., 2004b;McClain et al., 2012), this provides a promising avenue for parameterization of the 'pace of bioturbation'. However, there is a clear need to develop and test the generality of metabolism-based model formulations.
Despite decades of bioturbation research, there still seems to be a lack of mechanistic understanding of the complex processes that determine bioturbation rates in natural settings, as well as the diverse effects of bioturbation on physical, chemical and biological processes in the sediment, especially in heterogeneous coastal areas. On top of the effects of individual traits, the environmental surroundings of those individuals also modify their behavior and metabolism, leading to different effects on biogeochemical fluxes, as shown for example in a series of recent studies of a coastal gradient in the Baltic Sea (Joensuu et al., 2018;Bernard et al., 2019;Gammal et al., 2019).
Another long-lasting paradigm in RTMs has been the assumption that animal metabolism is unimportant compared to microbial degradation of organic matter, and can be omitted or implicitly included in a first-order organic matter degradation rate (Soetaert and Middelburg, 2009). As discussed in section "Solute Fluxes: Respiration and Nutrient Excretion", accumulating evidence suggests that faunal metabolisms can contribute significantly to the mineralization of organic carbon (Elmgren, 1984;Wenzhofer et al., 2002;Hansen and Bendtsen, 2013;Ehrnsten et al., 2019b) and nitrogen (Kristensen, 1984;Prins and Smaal, 1994;Magni et al., 2000). How can the metabolism of fauna then be coupled to organic matter degradation in RTMs? In theory, a straightforward way is to include the respiration and excretion of fauna in the sum of reactions affecting solids and solutes in the general diagenetic equations (Eqs. 9-10). In practice, this is not so straightforward. Sediment models generally resolve small-scale 1D processes with a fine depth resolution in mm to cm. Ecological or ecosystem models often focus on large-scale processes on m to km scales, and reduce the sediment depth resolution to a single layer, e.g., because of computational constraints (Soetaert et al., 2000). Additionally, even when sediments are resolved into several layers, the benthic fauna is usually represented as vertically integrated variables (i.e., mass area −1 , e.g., Fulton et al., 2004b;Butenschön et al., 2016), making it challenging to link biological and chemical processes in the vertical dimension. Thus, combining these scales is a major challenge for coupled modeling to understand benthic and pelagic processes, as also pointed out by Lessin et al. (2018).
Another challenge is the representation of organic matter. RTMs usually model the degradation of organic matter as a continuous process over depth, while ecological models generally include sediment organic matter as food for the benthos as a single depth-integrated variable. A good 'compromise' could be to divide the organic matter into a finite number of lability classes with different bio-availability for both faunal assimilation and microbial degradation, as in van Oevelen et al. (2011). Sohma et al. (2004Sohma et al. ( , 2018 took this approach one step further by also dividing the benthic functional groups into lability classes. While this formulation can help keep track of organic matter classes in a coupled RTM, it is hard to validate against observations.

Choosing the 'Right' Processes and Formulations
When modeling benthic nutrient and carbon cycling, there is a need to carefully consider which processes to include and how they are formulated to best answer the specific questions of the study. There are always trade-offs between realism, precision and generality of any given model (Levins, 1966). It is not possible to construct comprehensive guidelines for model combinations and formulations to be used in different circumstances, but here we try to provide some general advice.
The role of benthic fauna in coastal systems could be broadly categorized into three cases: (1) systems dominated by direct effects on pelagic processes, such as the effects of dense bivalve beds on pelagic N cycles and plankton dynamics; (2) systems dominated by direct and/or indirect effects in the sediment, such as effects of metabolism and bioturbation on microbial processes; and (3) systems where benthic fauna has a minor effect on biogeochemical fluxes. Modeling case 1 requires a good representation of the pelagic systems in terms of resolution, hydrodynamic transport and biological dynamics of plankton and suspension-feeders, but may not need a detailed description of sediment reactive-transport. In case 2, there may be a strong need to couple detailed sediment processes with a dynamic description of fauna, but pelagic processes may be simplified, while in case 3, benthic fauna could be omitted from the model.
It has been argued that benthic fauna generally is of minor importance for sediment processes in advection-dominated systems as opposed to diffusion-dominated ones, as hydrological forces override the effects of bioturbation and biodeposition (Mermillod-Blondin, 2011). Also, the effects of suspensionfeeders on pelagic processes (case 1) generally decrease with decreasing water residence time, which is strongly related to hydrodynamic flow speed (Dame and Prins, 1998;Grall and Chauvaud, 2002). Thus, as a general rule, one may expect the importance of benthic fauna to increase with decreasing flow speed and increasing water residence time. However, there are exceptions from this general rule, such as intertidal areas where bioturbation and grazing on microphytobenthos can have large effects on sediment resuspension and transport despite a dominance of advective tidal currents (Wood and Widdows, 2002;Le Hir et al., 2007;Rakotomalala et al., 2015).
In addition to the main types of impact, another important consideration when choosing a modeling approach is the main research question. Understanding responses to major environmental changes that may lead to completely new system states and shifted baselines requires an explicit description of mechanisms that determine the biological and biogeochemical dynamics of the system (e.g., complex mechanistic ecosystem-level models). In contrast, studies focusing on biogeochemical fluxes in systems with little ecological change may work better with empirical formulations for benthic activities that emphasize precision over generality (e.g., as is done in RTMs tuned to a specific site).
A related consideration is in what detail is desired in ecological aspects. Sophisticated individual-level models such as Dynamic Energy Budget models add more realism to the simulation of biomass, reproductive constraints and physiological fluxes compared to simpler functional group approaches, but so far, they have only been applied to single, well-studied species, as they require a large number of species-specific parameter values and several state variables for each individual. If the interest is in the effects of a diverse benthic community on biogeochemical fluxes (case 2), a functional group approach is probably the best option while structured population models and (sub)individual-level models are better suited for systems dominated by single species, e.g., bivalves (case 1). Simplifying, while more physiological detail leads to a more accurate description of biology, there is no clear proof that this automatically improves the description of biogeochemical fluxes at the ecosystem scale. A meta-analysis of pelagic biogeochemical models also concluded that increased complexity did not improve the ability of models to predict spatial and temporal patterns of plankton dynamics (Arhonditsis and Brett, 2004).
Can We Model the Future Biogeochemical Cycles in Coastal Ecosystems?
In principle, mechanistic models can be used to understand the dynamics of a system outside the boundaries of observations as long as the generic mechanisms of the model are validated. However, coastal ecosystems are extremely complex and our understanding of the different mechanisms determining the dynamics of these systems is far from complete. Additionally, accurate spatiotemporal information to force and initialize the models is largely lacking. Therefore, models describing processes in coastal systems all contain a large degree of empiricism. For example, diagenetic models (RTMs) are usually tuned to fit observational datasets of the distribution of sediment components from specific sites (Arndt et al., 2013). Attempts to find globally valid (empirical) relationship for parameterizations have been marginally successful (Boudreau, 1994;Middelburg et al., 1997). Similarly, models of biomass and metabolism of fauna generally contain environment-specific rate constants (such as ingestion half-saturation constants or mortality rates) that are used as tuning parameters. These empirical formulations limit the predictive capacity of models in a changing environment.
However, despite large uncertainties, models are still the most cost-effective tools available to simulate plausible future states of coastal ecosystems and may provide valuable management advice by examining future nutrient and carbon cycling in coastal areas, as long as the uncertainties are properly recognized. The use of model ensembles to quantify uncertainty in model projections has become increasingly popular in climate and marine hydrodynamic-biogeochemical sciences (Murphy et al., 2004;Meier et al., 2012Meier et al., , 2018. When moving up in complexity to biological process modeling, this approach would increase model utility considerably both for scientific and applied purposes.
An important part of model development is the validation against observations. For this purpose, long-term and consistent monitoring of benthic fauna and biogeochemical properties are invaluable. Only when we know that we can simulate past changes with reasonable accuracy, may we endeavor into projections of the future. We also acknowledge that the methods and extent of validation in published studies vary and call for more rigorous validation of benthic model dynamics in future modeling studies.
In addition, to increase the mechanistic understanding of complex coastal dynamics we suggest that model development should concentrate on well-studied coastal systems where process understanding and data availability are good. Increased interdisciplinary collaboration and understanding between e.g., benthic ecologists and biogeochemical modelers is a prerequisite to be able to plan experiments and monitoring that support modeling efforts and vice versa. A prime example is Port Phillip Bay in Australia where large field campaigns supported the development of several coupled hydrodynamic-biogeochemical ecosystem models Parslow, 1997, 1999;Fulton et al., 2004a,b). Another area with large potential for future model development is the Baltic Sea, because of its long research and monitoring tradition focusing on benthic fauna, its responses to environmental change and its role in coupled benthic-pelagic processes (Remane, 1934;Karlson et al., 2007;Griffiths et al., 2017;Ehrnsten et al., 2019b;Rodil et al., 2019). Additionally, the Baltic Sea has been affected by multiple anthropogenic pressures, such as eutrophication, climate change and overfishing, far longer and more severely than most other coastal seas. It also has a long history of management of these pressures, making it a good test case for the study of management challenges in the future coastal ocean under multiple pressures and shifting baselines (Villnäs and Norkko, 2011;Reusch et al., 2018).
To understand the carbon and nutrient cycles and constrain their budgets in coastal systems, it is important to recognize that coastal areas constitute a continuum from land to sea rather than only from the water surface to the bottom. In general, the variability in scales and complexity of coastal areas require models that are carefully built to answer specific questions in the study system. Some ecosystem modeling packages recognize this need for flexibility in their design, enabling the user to combine modules of different complexity, e.g. choosing between a single-or multi-layer benthic component (DHI, 2007;Butenschön et al., 2016).
Another major challenge is to incorporate the role of stochastic events, such as storms, upwelling and heavy precipitation, which may significantly impact biogeochemical cycling and the retention capacity of coastal systems (Canuel and Hardison, 2016;Malta et al., 2017). As some of the examples discussed in this study demonstrate, models are generally better at predicting the average conditions and potential distributions of biota than the large and somewhat stochastic variability often present in real coastal ecosystems (Triantafyllou et al., 2000;Fulton et al., 2004a). As climate change is expected to lead to changing extremes in the future (Easterling et al., 2000;Alexander et al., 2018), being able to incorporate the effects of extreme stochastic events becomes increasingly important when simulating future ecosystem states and biogeochemical fluxes.

CONCLUSIONS
Benthic fauna is an important component of nutrient and carbon cycles in coastal seas. There is a clear need to include and improve the description of benthic fauna in models of coastal nutrient and carbon cycling in order to properly describe the coupled benthic and pelagic processes, especially in the light of environmental change. We have summarized existing model formulations of different complexity to help implement descriptions of relevant zoobenthic processes in a variety of coastal habitats. In general, metabolism appears to be a good descriptor of both direct and indirect effects of benthic fauna on nutrient and carbon fluxes, offering a promising avenue for future developments in model formulations. Although numerical models can capture a large number of processes in a system, a full mechanistic description of the ecosystem dynamics will never be possible; therefore, it is important to clearly delineate the scope of each model with its purpose in terms of scales and complexity.
Major challenges and research priorities identified in this article are (1) to couple the dynamics of zoobenthic biomass and metabolism to sediment reactive-transport in models, (2) to test and validate model formulations against real-world data to better incorporate the context-dependency of processes in heterogeneous coastal areas in models, and (3) to capture the role of stochastic events. To understand coastal carbon and nutrient cycling in the present and future, lateral thinking is needed, recognizing that coastal biogeochemical cycles are heavily influenced by processes acting along the gradient from land to sea, and that the role of physics, chemistry and biology should all be considered.

AUTHOR CONTRIBUTIONS
EE, XS, AN, KT, and BG designed the review. EE, XS, and KT were mainly responsible for collecting the literature. All authors contributed to writing the manuscript, with XS mainly responsible for section "Bioturbation" and EE mainly responsible for all other sections. All listed authors have made significant intellectual contributions to the study.

FUNDING
This study received funding from the EU, the Academy of Finland and Formas through BONUS (Art 185) as part of the BONUS BALTICAPP and BONUS COCOA projects, and from the Swedish Agency for Marine and Water Management through their grant 1.11 -Measure for Marine and Water environment. Further funding was provided by research grants from The Academy of Finland (Project ID 294853 to AN) and the University of Helsinki and Stockholm University strategic fund for collaborative research (the Baltic Bridge initiative; AN, CH, and BG).
Ragnar Elmgren for comments improving the manuscript.