Modelling Marine Sediment Biogeochemistry: Current Knowledge Gaps, Challenges, and Some Methodological Advice for Advancement

The benthic environment is a crucial component of marine systems in the provision of ecosystem services, sustaining biodiversity and in climate regulation, and therefore important to human society. With the contemporary increase in computational power, model resolution and technological improvements in quality and quantity of benthic data, it is necessary to ensure that benthic systems are appropriately represented in coupled benthic-pelagic biogeochemical and ecological modelling studies. In this paper we focus on five topical challenges related to various aspects of modelling benthic environments: organic matter reactivity, dynamics of benthic-pelagic boundary layer, microphytobenthos, biological transport and small-scale heterogeneity, and impacts of episodic events. We discuss current gaps in their understanding and indicate plausible ways ahead. Further, we propose a three-pronged approach for the advancement of benthic and benthic-pelagic modelling, essential for improved understanding, management and prediction of the marine environment. This includes: (A) development of a traceable and hierarchical framework for benthic-pelagic models, which will facilitate integration among models, reduce risk of bias, and clarify model limitations; (B) extended cross-disciplinary approach to promote effective collaboration between modelling and empirical scientists of various backgrounds and better involvement of stakeholders and end-users; (C) a common vocabulary for terminology used in benthic modelling, to promote model development and integration, and also to enhance mutual understanding.


INTRODUCTION
The benthic environment is a key component of marine systems: it represents the transfer zone between the biosphere and geosphere; it provides a habitat for organisms from all domains of life; and it modulates biogeochemical cycles of carbon, macroand micronutrients and trace elements. The benthos can be a source of material for the water column (e.g., releasing nutrients from mineralisation of organic matter) as well as a temporary or permanent sink (e.g., C-sequestration). Furthermore, from an ecological point of view, benthic systems play a major role in sustaining biodiversity both directly, with large number of species/orders/phyla, and indirectly, by providing a large number of diverse habitats via ecosystem engineering (Meysman et al., 2006).
Despite this, benthic systems are generally under-represented in marine ecosystem models, in both regional and global applications, and are mostly present as a simple closure term for mass conservation (Soetaert et al., 2000;Hülse et al., 2017). Historically, the parallel development of marine hydrodynamic models and pelagic biogeochemical or plankton models, starting with NPZ formulations (Steele and Henderson, 1981;Fasham et al., 1990), coalesced because their spatial and temporal scales were compatible. Benthic modelling, however, often started from a different conceptual basis such as biodiversity (Hughes, 1984) or diagenesis (Berner, 1980;Boudreau, 1997;Paraska et al., 2014) considering different spatial and temporal scales, sometimes different currencies, and different mechanisms of particle and solute transport, which in sediments is enhanced by biological activity. Consequently, this divergence in development has led to marked differences in the characteristics of the two modelling approaches. Benthic models tend to be applied over smaller spatial scales (e.g., mm to m), typically in 1D (vertical), and solved under steady-state assumptions, while pelagic models usually operate over larger (m to km) scales in 2D or 3D, and are generally dynamic. Additionally, computational constraints have historically limited the ability to dynamically solve advective and diffusive transport in pore waters at large spatial scales; benthic data are relatively sparse and communication between the benthic and the pelagic scientific community has often been hindered by a lack of common vocabulary and objectives (Queirós et al., 2015). Nevertheless, fully coupled regional 3D pelagic-benthic models have been successfully developed and applied (Baretta et al., 1995;Wakelin et al., 2012) illustrating that there is no fundamental barrier between the two domains.
With the contemporary increase in computational power, model resolution and technology-led improvements in benthic data availability and synthesis, many of these constraints are now surmountable. Improved understanding of system complexity and feedbacks are providing renewed impetus for a holistic approach. For instance, understanding long-term evolution of the climate system requires a consideration of feedbacks from the sediments (Soetaert et al., 2000;Hülse et al., 2017); similarly, assessments of the impact of policies on marine environmental status can be worthless without accounting for the sediments that act as the memory of the system (Artioli et al., 2008;Soetaert and Middelburg, 2009). Conversely, investigating the evolution of benthic systems without considering pelagic dynamics may lead to imperfect answers.
Here, we present five priorities for further developing benthic process models that, due to their relevance and intrinsic interdisciplinary value, would ensure that the benthic system is given an appropriate consideration in coupled pelagic-benthic biogeochemical and ecological modelling studies (Figure 1).

REACTIVITY OF ORGANIC MATERIAL
One of the most important yet enigmatic aspects of benthic biogeochemistry is the reactivity of the particulate organic material (POM) degraded within the sediment (Berner, 1980). Constraining reactivity rates is crucial as it determines the sediment burial and recycling efficiencies (Martin and Bender, 1988;Soetaert et al., 1996), which allow the quantification of benthic-pelagic solute fluxes and, ultimately, feedbacks on ocean chemistry. Organic matter reactivity controls the nature and magnitude of the electron acceptor sinks from the water column: close to the sediment surface underlying oxic waters POM will be mainly respired aerobically and through denitrification, whereas deeper within the sediment it will mainly be mineralised through anaerobic pathways such as sulphate-reduction and methanogenesis (Regnier et al., 2011).
Parameterisation of POM reactivity in models requires information on external (environmental) and internal (intrinsic to the POM composition) regulating factors. The role each of these factors play in POM mineralisation is still very poorly understood (Burdige, 2007;Arndt et al., 2013). Attempts have been made to mathematically describe the vertical distribution of degradation by compartmentalising POM into discrete reactive fractions. Yet, in areas of limited data availability different empirical relationships are used to constrain first-order rate constants as a function of a master variable such as the POM rain rate or the sedimentation rate (Toth and Lerman, 1977;Middelburg et al., 1997;Martin and Sayles, 2006;Thullner et al., 2009). However, their large-scale spatial patterns remain poorly known, and POM reactivity constants do not strongly correlate with these global master variables (Arndt et al., 2013). Nonetheless, POM rain rate may be used for predicting its reactivity when export production is reasonably well described (Dunne et al., 2007).
Less common approaches, in particular the power (Middelburg, 1989) and reactive continuum (Boudreau and Ruddick, 1991) models, assume a continuous distribution of reactive types that capture the full range of degradation time scales. They require less parameterisation than discrete models, but are difficult to apply to bioturbated sediments where age and reactivity of organic matter depend on the burial velocity and bioturbation rate (Meile and Van Cappellen, 2005). Parameterisation of the continuum model to bioturbated sediments better captures the degradation of labile material within the top few millimetres of sediment, and more accurately predicts the drawdown of oxygen and nitrate from the water column (Stolpovsky et al., 2015). Continuum models can also be used to define the reactivity of discrete POM fractions, which can Frontiers in Marine Science | www.frontiersin.org 2 February 2018 | Volume 5 | Article 19 FIGURE 1 | Traditional view of benthic-pelagic coupling and approach to its representation in models (Left), the five cross-disciplinary research priorities that challenge this traditional approach (Middle), and lead to the new paradigm (Right), which is crucial for the improvement of biogeochemical and ecological modelling studies of the marine environment.
then be coupled to the total rain of POM to the seafloor . However, such approaches are generally not well suited to coastal areas and seas where seasonality is especially important. A better understanding of the spatiotemporal variability of lifetime and composition of POM is crucial for an improved representation of the benthic-pelagic coupling at regional to global scales.

THE CONTINUUM OF THE BENTHIC-PELAGIC BOUNDARY LAYER
The benthic-pelagic boundary is often simplistically regarded as a discontinuity, the interface between a completely fluid, dynamic water column and a rigid, porous benthic system. The reality is a spatially and temporally varying dynamic bottom boundary layer (BBL), characterised by strong physico-chemical gradients. While the importance of this layer on the mediation and transformation of biogeochemical fluxes, the transfer of organic material and sustaining heterogeneity of benthic habitats has been widely recognised, many uncertainties remain (Boudreau and Jorgensen, 2001). The benthic-pelagic transition is often blurred by the presence of easily resuspendable unconsolidated "fluff " and fluid muds, which can form layers extending from a few centimetres to several meters into the pelagic. The flocculent "fluff " usually consists of newly deposited POM, and can sustain rich microbial communities, also being a valuable food resource for benthic fauna (Laima et al., 2002). Fluid muds containing high concentrations of fine particular material are distinctive for the major delta-forming rivers, but also common in other coastal systems (McAnally et al., 2007). They may be subject to repetitive redox successions and mixing of refractory riverine material with more labile estuarine organic matter (Aller, 1998;Aller and Blair, 2006). To better account for the fate of organic material within the BBL, better knowledge of its composition, degradation and consumption as well as transport dynamics, is necessary. Supporting this development requires data collected with equipment that preserves the sediment-water interface structure.
At the opposite side of benthic-pelagic continuum, properties of pore water in non-cohesive sandy sediments are similar to the water column due to enhanced exchange of solutes, dominated by advective rather than diffusive transport driven by physical forcing, and enhanced by burrowing organisms (Huettel and Webster, 2001;Volkenborn et al., 2012). This affects redox structure, biological activity and ultimately the potential for carbon sequestration (Reckhardt et al., 2015). Better characterisation of sediment types and their regional distribution is necessary to more accurately quantify material fluxes and transformations in coupled models (Almroth Rosell, 2011).
Representation of a continuous benthic-pelagic interface in models requires adequate resolution: layers of few millimetres in thickness can be important for the bulk of biogeochemical transformations. A compromise usually lies between implementing high vertical resolution within the BBL and consolidated sediments (Yakushev et al., 2017) or choosing implicit parameterisations of vertical distributions (Ruardij and Van Raaphorst, 1995). To account for smooth benthic-pelagic transition in models, decreasing turbulent diffusivity with depth can be prescribed in the BBL, or potentially parameterised based on gradients of organic matter concentration and properties.

MICROPHYTOBENTHIC PRIMARY PRODUCTION
A large area of the coastal sediment receives enough light to support photosynthesis and sustain benthic primary production (Gattuso et al., 2006), whether seagrasses, macroalgae or microphytobenthos (MPB), which are generally not yet considered in holistic model systems. For example the MPB-typically dominated by pennate diatoms-can be biogeochemically significant, particularly in oligotrophic settings, supported by nutrients regenerated during benthic mineralisation (Glud et al., 2009b). Their diel rhythm of respiration and photosynthesis has implications for benthic redox conditions and affects the availability of oxygen and labile organic material in the upper centimetres of the sediment (Cook et al., 2007). This dynamic directly affects key processes in the nitrogen cycle (nitrification, denitrification, DNRA), the redox cycling of metal-oxides and metal-sulphides and thereby the mobility of metals, trace-metals, phosphorus, sulphur, and pollutants (Risgaard-Petersen et al., 1993;Fenchel and Glud, 2000;Dalsgaard, 2003). Recent investigations suggest that MPB metabolism also may directly influence nitrogen turnover through DNRA and intracellular storage of nitrate (Kamp et al., 2011). Furthermore, extensive production of exopolymeric substances by MPB can affect the permeability and erosion thresholds of sediments (Hanlon et al., 2006;Pierre et al., 2014). However, only a few theoretical models on element cycling, food web structure or biogeochemistry in the coastal zone illustrate the potential of MPB (Blackford, 2002;Baird et al., 2016). Inclusion of MPB in biogeochemical and trophic models requires high quality data on distribution and processes mediated by MPB. New field approaches for mapping biomass (Glud et al., 2002;Kazemipour et al., 2012), trophic coupling (Evrard et al., 2008;Oakes et al., 2012) and performance of MPB (Berg et al., 2013;Attard et al., 2014) may better facilitate model parameterisation, although basic understanding of in situ growth and grazing on MPB communities is still rudimentary. However, to fully appreciate the importance of MPB, field efforts have to be closely linked to investigations on physiology, metabolism and behaviour of MPB-dominated communities by, for instance, taking advantage of new imaging approaches (Grunwald and Kühl, 2004;Ralph et al., 2005;Hancke et al., 2014).

BIOLOGICAL TRANSPORT AND SMALL-SCALE HETEROGENEITY
Benthic systems are characterised by spatial heterogeneity at all scales, making sediments a complex mosaic of redox conditions and habitat niches, often not sufficiently represented by average rates or attributes (Glud et al., 2009a). Without understanding this variability, benthic systems may exhibit apparently paradoxical behaviour. For instance, nitrification in deep nominally anoxic sediment layers (Satoh et al., 2007), occurring as a consequence of burrowing fauna pumping oxygenated waters into the deeper sediment, locally modifying redox conditions (Wenzhöfer and Glud, 2004;Volkenborn et al., 2010). Biogeochemical models mostly consider bio-irrigation (biological enhancement of solute transfer) as a factor enhancing the diffusion coefficient (Blackford, 1997;Reed et al., 2011), resulting in a thicker oxic layer near the sediment surface. However, the transport of oxygen (and other solutes) will likely not occur homogenously, rather being concentrated along the network of burrows (Glud et al., 2016), which become an extension of the oxic-anoxic interface. It has been shown that burrowing activity of benthic macrofauna can lead to as much as 400% increase in denitrification rates (Gilbert et al., 1998;Webb and Eyre, 2004), mostly due to the high rates of nitrification occurring within the burrows (Howe et al., 2004). Burrows thus offer an ideal environment for diverse microbial communities to work very closely and efficiently: indeed, genetic analysis has shown that in-burrow bacterial communities are more similar to the surface sediment community than the ambient sediment community at similar depth (Satoh et al., 2007;Laverock et al., 2010). Conversely, anaerobic microniches may be formed in surface aerobic sediments when the consumption rate due to the mineralisation of POM by bacteria is higher than the oxygen diffusion into the particle. Under these conditions, anaerobic processes such as sulphate reduction and denitrification (Jørgensen, 1977;Lehto et al., 2014) may occur at the sediment surface.
Recently, sulphide oxidation has been observed even within anoxic sediments, at the expense of the free oxygen in the distant oxic layer, thanks to the presence of filamentous cable bacteria that physically connect the two layers and act as micro-cables transferring electrons (Nielsen et al., 2010;Pfeffer et al., 2012). In a seasonally hypoxic basin in the North Sea, cable bacteria are very abundant and active during winter when they build a large buffer of iron (hydro-)oxides before the onset of summer hypoxia that is able to prevent the formation of euxinic toxic waters by chemically oxidising the sulphides formed during summer (Seitaj et al., 2015). Sharp discontinuities in the redox structure of the benthic environment are rarely represented in models despite their potential significance. Explicit representation of these processes would require a significant increase in the vertical resolution of benthic models, which may be incompatible with operational modelling. In order to properly upscale these processes, detailed models implemented at a small scale, or a sub-grid scale approach could be used to derive optimal parameterisations.

EPISODIC EVENTS AND RECOVERY
Seasonal variability is an important driver of benthic-pelagic coupling in the marine environment (Martens, 1976;Aller, 1994). While it is relatively well described in current models, the role of episodic events has been long neglected because of the lack of observations and their unpredictability (Durrieu de Madron et al., 2011). Along the continental slope, earthquakes may generate landslides that affect deposition and, consequently, the benthic ecosystem (Oguri et al., 2013). Other events such as storms and river floods have been shown to dominate particle transport in key areas of the coastal zone which implies a large effect on community structuration and biogeochemical function.
In the North-Western Mediterranean, for example, the Rhône River shows very large floods that represent 70-80% of the annual particulate delivery (Antonelli et al., 2008) including POM. These very large deposition rates generate disturbances of the benthic communities and their functions, and recovery is often terminated by another event (Cathalot et al., 2010). At the same time, river floods deliver large amounts of nutrients that induce phytoplankton blooms which feed the pelagic and benthic biota (Auger et al., 2011). Waves and storm-induced currents can resuspend fine particulates from centimetres in the sediment (Toussaint et al., 2014;Bourrin et al., 2015) and transport them away from the nearshore coastal zone, thus feeding deeper benthic communities spread on the shelf and deep slope (Jahnke et al., 1990). Export from the shelf to the deep-sea is also dominated by intense events (Sanchez-Vidal et al., 2008) through cascading or downwelling in canyons. In the North-Western Mediterranean, cascading arises from a peak of cold winds which generate cold and dense waters on the shelf that are transported to the shelf break and sink, whereas downwelling events are linked to storms and cyclonic circulation. Both these phenomena contribute to major export of particles (Ulses et al., 2008) which impoverish the shelf and have profound impacts on the deep-sea ecosystem (Sanchez-Vidal et al., 2012) including beneficial feeding of the shrimp populations (Canals et al., 2009). Post-disturbance recovery dynamics of benthic systems will differ depending on nature of disturbance: hypoxia (Rosenberg et al., 2002), benthic trawling (Collie et al., 2000), dredge dumping (Bolam et al., 2006) or deep sea mining (Jones et al., 2017) will have different impacts and dynamics of recovery. One of the modelling challenges is the representation of lateral recolonization in benthic communities (Dittmann et al., 1999;Fowler, 2002), since microbial functions such as denitrification only return fully after recovery of deep-digging fauna (Van Colen et al., 2012). Modelling the impacts of intense events on benthic systems requires consideration of both the physical reorganisation of the sediments and the associated response of flora and fauna. Synthesis of current understanding of effects of these events on organisms' physiology and ecology will facilitate better parameterisation of their impacts.

CONCLUSIONS
Benthic environments are a crucial component in the functioning of marine systems and their services, particularly under the increasing pressure from anthropogenic activities and climate change. Therefore, development and application of numerical models that integrate our understanding of benthic systems is essential. However, benthic models often require complex, specific structures and scales, very different from faster-developing pelagic counterparts. In this paper, we have discussed some topical challenges related to different aspects of benthic environments, gaps in their understanding, and research priorities that will lead to tackling these gaps using a modelling approach (Figure 1). We believe that these challenges can be addressed properly with a different and broader approach than has been used in the past, one that maximises the utility of benthic modelling within a whole systems approach. In particular: • Building on (Soetaert et al., 2000) approach, a framework of traceable and hierarchical complexity for benthic-pelagic models should be developed. Benthic systems exert different effects and functions depending on the spatial and temporal scale being considered. Depending on the scientific question, some aspects can be neglected, others can be simplified. The function of the framework would be to assure that these choices are fully transparent both to scientists and end-users. This approach will facilitate integration among models, reduce risk of bias, and clarify any model limitations.
• A cross-disciplinary approach is needed: the complexity of benthic systems requires an integrated approach toward research, which necessitates effective communication and collaboration between modelling and empirical scientists of various backgrounds, involving stakeholders and end-users (Queirós et al., 2015). This approach has to be fully integrated rather than sequential or in parallel, including the design of the conceptual framework of the modelling and experimental approach.
• A common vocabulary for terminology used in benthic modelling community needs to be developed. With many strains and flavours of benthic models developed to date, an effective interdisciplinary dialogue requires a common language not only to promote model development and integration, but also to enhance mutual understanding and avoid ambiguity.
These three pillars form the foundation for a robust advancement in multi-disciplinary collaborations in benthic and benthicpelagic research necessary to better understand, manage and predict the marine environment.

AUTHOR CONTRIBUTIONS
All authors contributed to EuroMarine workshop "A multidisciplinary consortium approach for advancement of conceptual frameworks, modelling and experimental approaches to benthicpelagic coupling" held at Plymouth Marine Laboratory on May 18-19, 2015, and to consequent discussions leading to this contribution.