C-GEM (v 1.0): a new, cost-efficient biogeochemical model for estuaries and its application to a funnel-shaped system

Reactive transport models (RTMs) are powerful tools for disentangling the complex process interplay that drives estuarine biogeochemical dynamics, for assessing the quantitative role of estuaries in global biogeochemical cycles and for predicting their response to anthropogenic disturbances (land-use change, climate change and water management). Nevertheless, the application of RTMs for a regional or global estimation of estuarine biogeochemical transformations and fluxes is generally compromised by their high computational and data demands. Here, we describe C-GEM (Carbon-Generic Estuary Model), a new one-dimensional, computationally efficient RTM that reduces data requirements by using a generic, theoretical framework based on the direct relationship between estuarine geometry and hydrodynamics. Despite its efficiency, it provides an accurate description of estuarine hydrodynamics, salt transport and biogeochemistry on the appropriate spatio–temporal scales. We provide a detailed description of the model, as well as a protocol for its set-up. The new model is then applied to the funnelshaped Scheldt estuary (BE /NL), one of the best-surveyed estuarine systems in the world. Its performance is evaluated through comprehensive model–data and model–model comparisons. Model results show that C-GEM captures the dominant features of the biogeochemical cycling in the Scheldt estuary. Longitudinal steady-state profiles of oxygen, ammonium, nitrate and silica are generally in good agreement with measured data. In addition, simulated, system-wide integrated reaction rates of the main pelagic biogeochemical processes are comparable with those obtained using a highresolved, two-dimensional RTM. A comparison of fully transient simulations results with those of a two-dimensional model shows that the estuarine net ecosystem metabolism (NEM) only differs by about 10 %, while system-wide estimates of individual biogeochemical processes never diverge by more than 40 %. A sensitivity analysis is carried out to assess the sensitivity of biogeochemical processes to uncertainties in parameter values. Results reveal that the geometric parameters LC (estuarine convergence length) and H (water depth), as well as the rate constant of organic matter degradation (kox) exert an important influence on the biogeochemical functioning of the estuary. The sensitivity results also show that, currently, the most important hurdle towards regionalor global-scale applications arises from the lack of an objective framework for sediment and biogeochemical process parameterization. They, therefore, emphasize the need for a global compilation of biogeochemical parameter values that can help identify common trends and possible relationships between parameters and controlling factors, such as climate, catchment characteristics and anthropic pressure.

Abstract. Reactive transport models (RTMs) are powerful tools for disentangling the complex process interplay that drives estuarine biogeochemical dynamics, for assessing the quantitative role of estuaries in global biogeochemical cycles and for predicting their response to anthropogenic disturbances (land-use change, climate change and water management). Nevertheless, the application of RTMs for a regional or global estimation of estuarine biogeochemical transformations and fluxes is generally compromised by their high computational and data demands. Here, we describe C-GEM (Carbon-Generic Estuary Model), a new one-dimensional, computationally efficient RTM that reduces data requirements by using a generic, theoretical framework based on the direct relationship between estuarine geometry and hydrodynamics. Despite its efficiency, it provides an accurate description of estuarine hydrodynamics, salt transport and biogeochemistry on the appropriate spatio-temporal scales. We provide a detailed description of the model, as well as a protocol for its set-up. The new model is then applied to the funnelshaped Scheldt estuary (BE/NL), one of the best-surveyed estuarine systems in the world. Its performance is evaluated through comprehensive model-data and model-model comparisons. Model results show that C-GEM captures the dominant features of the biogeochemical cycling in the Scheldt estuary. Longitudinal steady-state profiles of oxygen, ammonium, nitrate and silica are generally in good agreement with measured data. In addition, simulated, system-wide integrated reaction rates of the main pelagic biogeochemical processes are comparable with those obtained using a highresolved, two-dimensional RTM. A comparison of fully transient simulations results with those of a two-dimensional model shows that the estuarine net ecosystem metabolism (NEM) only differs by about 10 %, while system-wide estimates of individual biogeochemical processes never diverge by more than 40 %. A sensitivity analysis is carried out to assess the sensitivity of biogeochemical processes to uncertainties in parameter values. Results reveal that the geometric parameters LC (estuarine convergence length) and H (water depth), as well as the rate constant of organic matter degradation (k ox ) exert an important influence on the biogeochemical functioning of the estuary. The sensitivity results also show that, currently, the most important hurdle towards regionalor global-scale applications arises from the lack of an objective framework for sediment and biogeochemical process parameterization. They, therefore, emphasize the need for a global compilation of biogeochemical parameter values that can help identify common trends and possible relationships between parameters and controlling factors, such as climate, catchment characteristics and anthropic pressure. and biologically modified along the estuarine gradient, with likely consequences for the coastal biogeochemical dynamics and, ultimately, for global biogeochemical cycles (e.g., Jahnke, 1996;Gattuso et al., 1998;Rabouille et al., 2001;Laruelle et al., 2009;Liu et al., 2010;Arndt et al., 2011;Jiao et al., 2011;Regnier et al., 2013a;Bauer et al., 2013).
The limited number of comparative studies covering a large range of estuarine systems hampers the identification of global patterns and precludes a robust assessment of the quantitative role of estuaries in global element cycles (Borges and Abril, 2011). In addition, individual estuarine systems reveal tremendous internal spatial and temporal heterogeneity, making it difficult to quantify the net carbon balance for a single estuary and even more for a set of representative systems upon which regional and global estimates could rely (Bauer et al., 2013). In this context, the long tradition of research in estuarine physics provides a suitable framework for addressing the large-scale estuarine biogeochemical dynamics. Dominant features of the estuarine transport can be constrained from hydrodynamic parameters (e.g., Stommel and Farmer, 1952;Hansen and Rattray, 1966;Prandle, 1985;Jay et al., 2000) or geometrical parameters (e.g., Pritchard, 1955;Davies, 1964;Dyer, 1973;Pethick, 1984;Dalrymple et al., 1992;Dürr et al., 2011), two seemingly distinct approaches that can be related to one another through the interdependence between estuarine geometry and hydrodynamics (Savenije, 1992). Hence, important transport and mixing properties can be directly deduced from readily available geometric data (Savenije, 2005(Savenije, , 2012. Taking into account that the hydrodynamics also exerts a first-order control on the estuarine biogeochemistry (e.g., Alpine and Cloern, 1992;Friedrichs and Hofmann, 2001;, a logical step is to use these interdependencies to predict the biogeochemical dynamics from the main geometrical features of estuaries. The tight hydrodynamic-biogeochemical coupling has already been partly recognized in the past, for instance by correlating the biogeochemical behavior of an estuary with given hydrodynamic characteristics such as residence time or tidal forcing (Monbet, 1992;Nixon et al., 1996;Laruelle, 2009), yet these correlations are based on a limited number of data sets (< 40) that do not cover the diversity of estuarine systems and do not resolve their seasonal and inter-annual variability (e.g., Brion et al., 2008;Arndt et al., 2009). Such a correlative approach also does not provide fundamental insights into the complex interplay of multiple reaction and transport processes in estuarine systems (Nielsen et al., 1995;Geyer et al., 2000;Arndt et al., 2009). The aim is thus to extend the approach and to develop generalized methods for up-scaling that resolve the strong spatio-temporal variability of the estuarine environment and explicitly account for the process interplay that controls the biogeochemical cycling of carbon and nutrients along the estuarine gradient.
Over the last three decades, increasingly complex processbased models have been applied to unravel the organic and inorganic carbon and nutrient cycles on the scale of individual estuaries (e.g., O'Kane, 1980;Soetaert and Herman, 1995;Lin et al., 2007;Arndt et al., 2009;Cerco et al., 2010;Baklouti et al., 2011), yet none of these models are currently suitable for regional or global applications (Bauer et al., 2013). In particular, model applications remain limited by data requirements, calibration and validation procedures as well as by the high computational demand required to address important physical, biogeochemical and geological processes on the relevant temporal and spatial scales (Regnier et al., 2013b). Therefore, applications on scales larger than individual, well-constrained systems require simplifications to afford the treatment of a large number of estuaries, including those for which morphological, hydrodynamic and biogeochemical data are incomplete or absent. A generalization of simulation results from a representative set of systems covering contrasting climate, hydromorphology and catchment properties will ultimately provide better estimates of the quantitative contribution of estuaries to global biogeochemical cycles.
Here, we propose the Carbon-Generic Estuary Model (C-GEM), a new, one-dimensional, generic reactive transport model (RTM) for the biogeochemical dynamics of carbon and associated bio-elements (N, P, Si) in estuaries. RTMs are well-established quantitative tools for disentangling the complex biogeochemical dynamics of estuaries (Thouvenin et al., 1994;Regnier et al., 1997Regnier et al., , 2003Vanderborght et al., , 2007Arndt et al., 2009), including their response to anthropogenic perturbations (Paerl et al., 2006;Thieu et al., 2010) and the complex process interplay that underlies system-wide key biogeochemical indicators, such as net ecosystem metabolism (NEM), an integrative measure of the whole system biogeochemical dynamics defined as the difference between net primary production (NPP), aerobic degradation and denitrification on a system scale (Odum, 1956;Andersson and Mackenzie, 2004). C-GEM is not only computationally efficient, but also reduces data requirements by using an idealized representation of the estuarine geometry to support hydrodynamic calculations and, subsequently, transport and biogeochemical reaction processes. The C-GEM modeling platform is thus compatible with hundreds to thousands of stationary or fully transient simulations (including daily to seasonal fluctuations) on a time span of years to decades, using geometric information readily available through maps or remote sensing images. Moreover, unlike simpler box model approaches, which are still widely used to assess global estuarine dynamics (e.g., Andersson et al., 2005;Slomp and Van Cappellen, 2007;Laruelle, 2009;Mackenzie et al., 2012), C-GEM resolves the most important temporal and spatial scales and provides an accurate description of the estuarine hydrodynamics and transport. It may thus represent a promising avenue towards the development of a generalized method for exploring and quantifying biogeochemical transformations and fluxes in alluvial estuaries on the regional and/or global scale.
In the first part of this paper, the general structure of C-GEM is described. This includes detailed descriptions of the model support, of the fundamental equations for the hydrodynamics and transport and their parameterization and of the biogeochemical reaction network. In addition, a generic protocol for the set-up of C-GEM for an estuarine system is illustrated and different strategies will be proposed depending on the availability of data to constrain model parameters. The second part of this paper presents, as a proof of concept, the application of C-GEM to the funnel-shaped Scheldt estuary (Belgium-Netherlands). The macro-tidal Scheldt estuary is among the best-surveyed estuarine systems worldwide and has been the subject of intense modeling efforts (e.g., Wollast and Peters, 1978;Soetaert and Herman, 1995;Regnier et al., 1997;Vanderborght et al., , 2007Billen et al., 2005;Desmit et al., 2005;Hofmann et al., 2008;Arndt et al., 2009Arndt et al., , 2011Gypens et al., 2013). In order to test the performance of C-GEM in predicting the estuarine hydrodynamics and biogeochemical dynamics, both steady-state simulations for average summer conditions as well as transient simulations for an entire year (2003) are carried out. Steady-state simulations are compared with a comprehensive set of field observations, while mass budget results, as well as NEM, derived from the transient simulation, are compared with results from a highly resolved 2D-RTM for the same period (Arndt et al., 2009). This model-data, model-model comparison allows one to assess the model's performance on different temporal and spatial scales. In addition, a sensitivity analysis is carried out to identify model parameters that exert the most important control on biogeochemical processes and to assess the sensitivity of estimated process rates to uncertainties in these parameter values. Finally, current model limitations with respect to local, regional and, ultimately, global-scale applications are critically analyzed.

Model support
Alluvial estuaries are commonly defined as systems that are characterized by a movable bed, consisting of sediments of both marine and terrestrial origin, and a measurable influence of freshwater discharge (Savenije, 2005(Savenije, , 2012. In such estuaries, the amount of water flow entering or leaving the estuarine channel is entirely controlled by the shape of the estuary (Pethick, 1984). In turn, the water movement, driven by tides and freshwater discharge, leads to a redistribution of the unconsolidated sediments and determines the shape of the estuary. Alluvial estuaries display a wide variety of shapes ranging from funnel-shaped estuaries with a dominant tidal influence to prismatic estuaries with a large fluvial influence. Nevertheless, they bear common geometric characteristics that are compatible with an idealized representation of an estuary (Savenije, 1992(Savenije, , 2005(Savenije, , 2012. For tidally averaged conditions, their cross-sectional area A or width B can be described by decreasing exponential functions with distance, x, from the mouth (Savenije, 1986(Savenije, , 2005(Savenije, , 2012: where A 0 and B 0 are the cross-sectional area and the width at the estuarine mouth (x = 0), respectively, a is the crosssectional convergence length and b is the width convergence length. Combining Eqs. (1) and (2) leads to an expression for the mean longitudinal variation in estuarine depth, h (Savenije, 2005): (3) Savenije (1992) showed that alluvial estuaries can be classified according to the Canter-Cremers number, N, and the estuarine shape number, S. The dimensionless hydrodynamic Canter-Cremers number for flood discharge is defined as the ratio between the volume of the river discharge and the volume of saline water flowing into the estuary during a tidal period (Savenije, 2012): where Q b is the bankfull discharge, defined as the momentary maximum flow, which has an average recurrence interval of 1.5 years, associated with a state of maximum velocity in the channel and, therefore, with the maximum ability to govern the shape and the size of the channel. T is the tidal period, which corresponds to the interval between successive high (or low) tides, and P is the tidal prism that represents the amount of water that flows in and out an estuary between high and low tide. The dimensionless estuarine shape number is a geometric parameter defined as the ratio between the convergence length a and the tidally averaged depth at the estuarine mouth (h 0 ): These two numbers provide a theoretical framework to analyze the tight link between the geometry and the hydrodynamics of estuaries ( Fig. 1). We can see that estuaries with a large Q b are more riverine and have a long convergence length. On the other hand, estuaries with a large tidal prism are generally deep and have a short convergence length. Based on Fig Savenije, 1992). The Scheldt estuary, where C-GEM has been tested, is highlighted in red.
N (> 15) and S (> 15000) and mixed-type estuaries fall in between these two end-member cases. For instance, estuaries such as the Limpopo estuary ( Fig. 2a) have a long convergence length and a dominant fluvial influence and show a longitudinal salt intrusion distribution that exponentially declines towards the land. At the opposite end of the shape spectrum, the Scheldt estuary has a short convergence length and a marine character, with a dome-shaped salt intrusion curve (Fig. 2c). The Incomati estuary is a good representation of the mixed category, showing a half-Gaussian shaped salt intrusion curve (Fig. 2b). The recognition of this tight link between estuarine geometry, hydrodynamics and transport (Fig. 2) and the identification of three main estuarine types ( Fig. 1) becomes important when thinking about estuarine biogeochemical dynamics and its significance for global biogeochemical cycles. Because estuarine hydrodynamics exert a first-order control on transport and biogeochemical processes (Fig. 3), estuarine biogeochemical characteristics, such as NEM, carbon and nutrient filtering capacities or CO 2 exchange fluxes can potentially be directly linked to hydrodynamic and thus geometrical characteristics. Such direct relationships between biogeochemical and readily available geometric characteristics would not only serve as a promising basis for a biogeochemical classification scheme, but would also significantly facilitate a quantitative assessment of the role of estuaries in global biogeochemical cycles and their response to anthropogenic perturbations including land-use and climate change (Regnier et al., 2013b).

Hydrodynamics
Estuaries are subject to tidal forcing and freshwater inflow. At the estuarine mouth, tidal variations in water level induce a tidal wave. This wave travels upstream and is progressively distorted due to the combined influence of the estuarine geometry and river discharge. The tidal range is, to a first order, determined by the balance between energy gain through channel convergence and energy loss through friction on the estuarine bed. As a result, fundamental hydrodynamic characteristics, such as tidal range, tidal excursion and the phase lag of the tidal wave vary along the estuarine gradient and can be related to key geometric characteristics, such as convergence lengths or depth. For weakly stratified or wellmixed estuaries whose depth is much smaller than width, the hydrodynamics can be described by the one-dimensional barotropic, cross-sectionally integrated mass and momentum conservation equations for a channel with arbitrary geometry (Nihoul and Ronday, 1976;Regnier et al., 1998;Regnier and Steefel, 1999): where t = time (in s), x = distance along the longitudinal axis (in m), A = cross-sectional area (A = H · B) (in m 2 ), Q = cross-sectional discharge (Q = A · U ) (in m 3 s −1 ), U = flow velocity (in m 2 s −1 ), r s = storage ratio (r s = B s /B) (-), B s = storage width (in m), C = Chézy coefficient (in m 1/2 m −1 ), and H = water depth (H = h + ξ (x, t)) (in m). The coupled partial differential equations (Eqs. 6 and 7) are solved by specifying the elevation ξ 0 at the estuarine mouth and the river discharge Q r (t) at the upstream limit of the model domain. Bed friction exerted on the moving water is described by means of a roughness formulation following Manning-Strickler (Savenije, 2012): where C is the Chézy coefficient, n is the channel roughness coefficient or the dimensionless Manning number and H is the water depth. The bed roughness, which depends on the bottom material and on the depth of the flow, is a notoriously difficult parameter to measure and is generally constrained via model calibration by fitting simulated water elevations, tidal wave propagation and current velocities to observations. In the absence of data, realistic C values range between 40 and 60 m 1/2 s −1 (Savenije, 2001(Savenije, , 2012. Lower values can typically be applied in the shallow tidal river where bottom friction is significant, while higher values can be applied in the saline estuary.

Mass conservation for solutes
The one-dimensional, tidally resolved, advection-dispersion equation for a solute C(x, t) in an estuary can be written as (e.g., Pritchard, 1958) In Eq. (9), Q and A are provided by the hydrodynamic model and P is the sum of all production and consumption process rates for the solute C. The effective dispersion coefficient D (m 2 s −1 ) implicitly accounts for dispersion mechanisms associated with sub-grid scale processes (Fischer, 1976;Regnier et al., 1998). In general, D is maximal near the sea, decreases upstream and becomes virtually zero near the tail of the salt intrusion curve (Preddy, 1954;Kent, 1958;Ippen and Harleman, 1961;Stigter and Siemons, 1967). The effective dispersion at the estuarine mouth can be quantified by the following relation ( Van der Burgh, 1972): where h 0 (m) is the tidally averaged depth at the estuarine mouth, N is the dimensionless Canter Cremers estuary number defined as the ratio of the freshwater entering the estuary during a tidal cycle to the volume of salt water entering the estuary over a tidal cycle (Eq. 4) and g (m s −2 ) is the gravitational acceleration. The variation in D along the estuarine gradient can be described by Van der Burgh's equation (Savenije, 1986): where K is the dimensionless Van der Burgh coefficient and the minus sign indicates that D increases in the downstream direction (Savenije, 2012). The Van der Burgh coefficient is a shape factor that can be shown to have values between 0 and 1 (Savenije, 2012), which depends on geometry and tidally average conditions. Therefore, each estuarine system has its own characteristic K value, which correlates with geometric and hydraulic scales (Savenije, 2005). It has thus been proposed, based on a regression analysis covering a set of 15 estuaries, that K can be constrained from the estuarine geometry (Savenije, 1992):

Biogeochemical reactions
The reaction network for the water column estuarine biogeochemistry includes total (particulate and dissolved) organic carbon (TOC), oxygen (O 2 ), ammonium (NH 4 ), nitrate (NO 3 ), phosphate (PO 4 ), dissolved silica (dSi) and phytoplankton biomass (PHY) as state variables. The reaction network considers the essential biogeochemical processes that affect carbon and associated bio-elements: primary production, phytoplankton mortality, aerobic degradation, denitrification, nitrification and O 2 exchange across the air-water interface. Variables and process rates included in C-GEM are schematized in Fig. 4 and their formulations and stoichiometric equations are summarized in Table 1. Despite its limited set of reaction processes, the simplicity of the biogeochemical network warrants application in data-poor systems. The gross primary production rate, GPP, is controlled by the underwater light regime that explicitly accounts for the effect of the suspended particulate matter (see below) and neglects phytoplankton self-shadowing, an effect that is generally weak in turbid estuarine systems (Desmit et al., 2005). In addition, macronutrient concentrations (dSi, DIN = NO 3 + NH 4 and PO 4 ) limit phytoplankton growth through a succession of Michaelis-Menten terms, each with their corresponding half-saturation constant, K MM . Net primary production, NPP, is calculated as the difference between GPP and autotrophic phytoplankton respiration, which accounts for biosynthesis, maintenance and excretion. Biosynthesis and excretion terms are assumed to be linearly proportional to GPP (Weger et al., 1989;Langdon, 1993;Lancelot et al., 2000), while the maintenance term is a direct function of the total phytoplankton concentration . The gradual switch between ammonium and nitrate utilization pathways for NPP is controlled by the availability of ammonium. Phytoplankton mortality is linearly proportional to the phytoplankton concentration through a mortality rate constant, k mort , which integrates the combined effects of cell lysis and grazing by higher trophic levels. Upon death, phytoplankton contributes to the total organic matter pool. The latter is represented as a single pool including only the fraction of the organic carbon, which actively contributes to the short-term supply of inorganic nutrients (Regnier and Steefel, 1999). Thus, the model does not account for burial of (refractory) particulate organic carbon in estuarine sediments (Abril et al., 2002;. Organic matter is degraded by aerobic degradation, aer_deg, and denitrification, denit. If oxygen concentrations are sufficient, aer_deg is the most energetically favorable pathway, and thus dominates the other metabolic processes (e.g., Stumm and Morgan, 1996). denit becomes important in polluted estuaries where oxygen levels drop to limiting concentrations. The heterotrophic degradation processes are described by Michaelis-Menten terms for both organic carbon and electron acceptor concentration . By oxidizing NH 4 to NO 3 , nitrification (nit) consumes large amounts of O 2 in polluted estuaries (Soetaert and Herman, 1995;Regnier and Steefel, 1999;Andersson et al., 2006;Hofmann et al., 2008). It is formulated as a one-step process including two Michaelis-Menten terms with respect to O 2 and NH 4 . The temperature dependence of maximum degradation rates, k ox and k denit , and maximum nitrification rate, k nit , is expressed via a function with a Q 10 value. Oxygen transfer through the air-water interface, O 2,ex , exerts an important influence on the oxygen concentration in the water column. The exchange rate is expressed by the product of the piston velocity (vp) and the difference between oxygen concentration and oxygen saturation. The latter is calculated as a function of temperature and salinity (Benson and Krause, 1984), while the piston velocity is calculated as the sum of two terms attributed to the current velocity and the wind speed at 10 m above the air-water interface (Regnier et Arndt et al. (2009), c Garnier et al. (1995). * If PHY = DIA, nlim needs to account for the silica limitation in the phytoplankton growth.

660
Switch between NH 4 and NO 3 utilization a f NH 4 = NH 4 , 2002). At this stage, the benthic-pelagic exchange is not included in the model, although cost-efficient numerical approaches are available for carbon and nutrients (e.g., Jahnke et al., 1982;Ruardij and van Raaphorst, 1995;Soetaert et al., 1996;Gypens et al., 2008). Hence, the application of C-GEM to shallow, pristine estuarine systems subject to intense element recycling within the sediments is not recommended at this stage.

Suspended particulate matter
The simulation of the suspended particulate matter (SPM) dynamics is required for the prediction of the light availability within the water column that exerts an important control on primary production in turbid estuaries, mainly. The onedimensional, tidally resolved, advection-dispersion equation for suspended particulate matter (SPM) dynamics follows an equation similar to that of solutes (Eq. 9) with the addition of two extra terms describing the mass exchange with the material surfaces of the estuarine bed: where R ero and R dep denote the erosion and deposition rates, respectively. In the theory of cohesive sediment transport, they are often considered to be mutually exclusive (Sanford and Halka, 1993) and expressed according to the wellestablished formulation of Partheniades (1962) and Einstein and Krone (1962): where H denotes the water depth and p ero and p dep (-) are the probabilities for erosion and deposition, respectively. E (mg m −2 s −1 ) is the erosion coefficient, while w s (m s −1 ) is the settling velocity of particles. p ero and p dep are given by (Einstein and Krone, 1962;Dyer, 1986;Mehta et al., 1989) where τ cr (N m −2 ) is the critical shear stress for erosion and deposition. The bottom shear stress, τ b (N m −2 ), is calculated dynamically using the quadratic friction law where ρ w (kg m −3 ) is the pure water density. All SPM parameters (τ cr , τ b , E, w s ) implicitly account for geomorphological and biological processes, such as sediment composition or biological stabilization mechanisms that are not explicitly resolved (e.g., Wolanski et al., 1992;Cancino and Neves, 1999;van Ledden et al., 2004). SPM parameter values are generally derived by model calibration against locally observed SPM data and their transferability to other estuarine systems may thus be limited.

Numerical solution
The non-linear partial differential equations are solved by a finite difference scheme on a regular grid, with a grid size x = 2000 m and using a time step t =150 s. If required, both spatial and temporal resolution can easily be modified. Transport and reaction terms are solved in sequence within a single time step using an operator-splitting approach . The advective term in the transport equation is integrated using a third-order accurate total variation diminishing algorithm with flux limiters, ensuring monotonicity (Leonard, 1984), while a semi-implicit Crank-Nicholson algorithm is used for the dispersive term (Press et al., 1992).
The schemes have been extensively tested using the CON-TRASTE estuarine model (e.g., Regnier et al., 1998;Regnier and Steefel, 1999; and guarantee mass conservation to within < 1 %. The erosion-deposition terms, as well as the reaction network, are numerically integrated using the Euler method (Press et al., 1992). The primary production dynamics, which requires vertical resolution of the photic depth, is calculated according to the method described in Vanderborght et al. (2007).

Protocol for the set-up of C-GEM
The following section is a step-by-step protocol describing how to set up C-GEM and specifying data requirements at each step. Each step of the set-up is described using direct references to the corresponding source code file of C-GEM provided as supplementary material (refer to the end of the manuscript for more details).

Step 1: construction of the idealized geometry
The idealized estuarine geometry is defined by the estuarine length (EL) and the depth (DEPTH), as well as the width (B). The depth and the width are specified in define.h for both upper (B_ub and DEPTH_ub) and lower (B_lb and DEPTH_lb) boundaries. In general, and especially for navigable channels, estuarine bathymetric data are available or can be derived from navigation charts. If no data are available, the depth can be approximated using remote sensing data (Gao, 2009) or assumed to be about 7 m for alluvial estuaries (e.g., Savenije, 1992). The estuarine width at both boundaries of the model domain can be easily derived from local maps. The width convergence length, LC, is then calculated in init.c using Eq. (2). The cross-sectional area is then calculated at every grid point by the product of water depth and estuarine width (see Eq. 6).

3.2
Step 2: set-up of the hydrodynamic module

Step 2.1: parameters
The Chézy coefficient (C) is the only control parameter in the equation of motion. Its value is defined at the two boundaries of the model domain (define.h) and its variation in space is specified in init.c. The Chézy coefficient is rarely measured and, thus, generally calibrated (Savenije, 1992). If observations for model calibration are missing, typical values reported in the literature for alluvial estuaries are 60 m 1/2 s −1 in the saline zone and 40 m 1/2 s −1 in the freshwater reaches (Savenije, 1992(Savenije, , 2001.

Step 2.2: boundary conditions
The boundary conditions for the hydrodynamic module are specified in define.h and consist of the freshwater discharge (Q r ) at the upstream boundary and the tidal elevation at the estuarine mouth, which requires specification of the amplitude (AMPL) and the frequency (pfun). Tidal elevation can be deduced from water level data obtained from gauging stations or estimated theoretically using an astronomical model (e.g., Regnier et al., 1998). The freshwater discharge is often monitored in rivers, but when missing, it can be derived from local or global watershed model outputs Fekete et al., 2002).

Step 2.3: validation
Hydrodynamics can be validated by comparing simulated and observed tidal amplitude profiles. If water level time series are not available, remote sensing data, such as laser altimetry, can be used to validate tidal wave amplitude and propagation (Cazenave and Savenije, 2008). Although promising, this method remains currently limited to a few locations (e.g., Syed et al., 2008).

3.3
Step 3: set-up of the salt transport module

Step 3.1: parameters
The dispersion coefficient at the estuarine mouth, D 0 , and its longitudinal variation are the only controlling parameters of the transport module. They are calculated in init.c. according to Eqs. (10), (11) and (12).

Step 3.2: boundary conditions for salinity
Boundary conditions for salinity are specified in init.c. In general, the upper boundary condition is set to 0, while the lower boundary condition can be extracted from local measurements or regional or global databases such as the World Ocean Atlas (http://www.nodc.noaa.gov/OC5/indprod.html).

Step 3.3: validation
The validation of the transport module is typically performed by comparing simulated longitudinal salinity profiles with observed data collected along the estuarine gradient or by comparing simulated and measured time series at a given location (e.g., Regnier et al., 1998). Note that the transport module is based on a predictive model, which only requires geometrical information. Hence, it can also be applied in estuaries for which salinity data are not available.

3.4
Step 4: set-up of the SPM module

Step 4.1: parameters
The sediment settling velocity, w s , the critical shear stress for erosion and deposition, τ ero and τ dep , and the erosion coefficient, Mero, are specified in define.h. τ ero , τ dep and Mero need to be defined at both the upper and lower boundaries. If longitudinal variations in sediment parameters need to be implemented, their formulations are defined in sed.c. These parameters generally require calibration. However, since the bottom material of the wider part of alluvial estuaries consists of mud or fine sediments (Savenije, 1986), w s rarely exceeds 1 mm s −1 (Winterwerp, 2002). Other parameters such as τ ero , τ dep and Mero are calibrated on the basis of observed SPM profiles. The latter is an important step where observations still remain essential.

Step 4.2: boundary conditions
Boundary conditions for SPM are specified in init.c. SPM concentrations are usually available for navigable channels, in particular those where dredging works are carried out. In the case of data-poor systems, the upper boundary condition can be derived from global statistical models, such as Glob-alNEWS2 . When no observations or models are available to constrain lower boundary conditions, SPM values can be deduced from remote sensing data (e.g., Bowers et al., 1998;Fettweis and Nechad, 2011).

Step 4.3: validation
SPM dynamics may be validated by comparing simulated longitudinal profiles along the estuarine axis and/or time series modeled at a given location with observed sediment concentrations. Otherwise, simulated concentrations can be validated using remote sensing and satellite data (e.g., Stumpf, 1988;Moore et al., 1999;Robinson et al., 1999;Doxaran et al., 2002Doxaran et al., , 2009van der Wal et al., 2010). The C-GEM biogeochemical module is implemented in biogeo.c by defining all biogeochemical reaction equations and by implementing all stoichiometric coefficients for each variable of the model. This structure allows for a flexible implementation and a rapid extension of the network by, for instance, different phytoplankton groups or additional transformation processes, such as adsorption-desorption or benthicpelagic exchange processes.

Step 5.2: parameters
All parameter values for the biogeochemistry are specified in define.h. In most estuaries, system-specific values for all required parameters are not available, but a literature survey can provide reasonable ranges within which a calibration can be performed (e.g., Cerco and Cole, 1994;Garnier et al., 1995;Le Pape et al., 1999;Desmit et al., 2005 for the phytoplankton mortality rate constant or Regnier et al., 1997Regnier et al., , 1999Park et al., 2005; for the nitrification rate constant). Unfortunately, estuarine parameter values for the biogeochemistry remain to be assembled in a global database (Regnier et al., 2013b).

Step 5.3: boundary conditions
The boundary conditions required for the biogeochemical module are assigned a numerical value in init.c. If direct observations are not available, boundary conditions for the riverine inputs of organic carbon and nutrients can be extracted from the GlobalNEWS2 global watershed statistical model , while boundary conditions at the downstream limit can be obtained from the World Ocean Atlas (http://www.nodc.noaa.gov/OC5/indprod.html).

Step 5.4: external forcings
The biogeochemical module requires specification of a number of external forcings depending on the formulation used to describe biogeochemical processes. For instance, in this study, phytoplankton growth depends on irradiance, photoperiod and temperature. The latter also influences other biogeochemical transformations, such as heterotrophic degradation and nitrification, while wind speed is required to constrain the exchange rate at the air-water interface. In C-GEM, photoperiod, temperature and wind speed are specified in define.h, while irradiance is calculated in fun.c. All external forcings should preferably be derived from observations, but, if direct observations are not available, irradiance and photoperiod can be constrained using radiation models (e.g., van der Goot, 1997) or may be extrapolated as a function of time, year and latitude using the astronomical equation of Brock (1981). Other external forcings can be obtained from global databases, such as the World Ocean Atlas (http://www.nodc.noaa.gov/OC5/indprod.html) for the water temperature and the CCMP data set (Atlas et al., 2011) for the wind velocity.

Step 6: sensitivity analysis
A sensitivity analysis is a crucial part of the iterative revision process of the model set-up. Depending on the results of each model validation and sensitivity analysis, the user may be required to repeat a step or even return to a previous step. The sensitivity analysis also provides useful information regarding the uncertainty in model predictions.
4 Application to the funnel-shaped Scheldt estuary: a test case

The Scheldt estuary
The Scheldt river and its tributaries drain an area of 21 580 km 2 in northern France, western Belgium and southwestern Netherlands before discharging into the southern North Sea (Fig. 5a). Its hydrographical basin includes one of the most populated regions of Europe, heavily affected by human activities (e.g., Wollast and Peters, 1978;Billen et al., 1985;Soetaert et al., 2006). The part of the river that is influenced by the tide is referred to as the Scheldt estuary extending 160 km from the estuarine mouth at Vlissingen (Netherlands) to Gent (Belgium), where a sluice blocks the tidal wave. The tide is semi-diurnal with an amplitude of about 4 m (Regnier et al., 1998). Salt intrudes as far as 100 km from the estuarine mouth. Upstream of 100 km, the estuary is characterized by a complex network of six tributaries (Dender, Durme, Grote Nete, Kleine Nete, Zenne and Dijle). The latter four form the Rupel, a single stream, which rejoins the main channel of the Scheldt at the salt intrusion limit.

Geometry
The Scheldt estuary is characterized by a large tidal range inducing a short convergence length (Table 2) and can thus be classified as a funnel-shaped system ( Fig. 1) (Savenije, 2005). Figure 5 compares the geometry of the Scheldt estuary (Fig. 5a) to its idealized geometry ( Fig. 5b and c) derived from the width convergence length, water depth and tidal amplitude. A variable depth (h) is applied here to account for a small, constant bottom slope over the total estuarine length. This idealized geometry (Fig. 5b and c) forms the support for C-GEM and illustrates the typical features of a funnel-shaped estuary: wide and deep at the mouth with a short convergence length, which induces a rapid upstream decrease in width.

Boundary conditions
Both steady-state and transient model simulations are conducted to test the performance of C-GEM. For both cases, a spin-up period of two months is imposed. In addition, a constant tidal amplitude is applied at the estuarine mouth. The tidal amplitude only accounts for the dominant semi-diurnal component M2, characterized by a period of 12.42 h and a frequency of 0.080 cycles h −1 (Regnier et al., 1998). For the steady-state simulations, a constant river discharge is specified at the inland limit of the Scheldt and its tributaries. In addition, constant biogeochemical boundary conditions and physical forcings (e.g., temperature and light intensity), representative of the summer conditions during the 1990s (Table 3; for further details see Vanderborght et al.,  Arndt et al., 2009 for details) are performed to test the performance of C-GEM in quantifying integrative, system-scale biogeochemical indicators, such as NEM. These integrative indicators cannot easily be quantified on the basis of observations alone and its quantitative assessment thus requires the application of model approaches (e.g., Arndt et al., 2009Arndt et al., , 2011Regnier et al., 2013b). Here, C-GEM results are compared to the outputs from a carefully calibrated and validated, highly resolved horizontal 2-D reactive transport model (Arndt et al., 2009). The latter uses a total of 56 000 computational points and provides a very detailed representation of the estuarine morphology. Both models are forced with identical boundary conditions and physical forcings (see Arndt et al., 2009 for a detailed description).

Suspended particulate matter and biogeochemistry
For the sake of comparison, all biogeochemical parameters and the biogeochemical reaction network, described in Sect. 2.4, are identical to those used in Arndt et al. (2009), with the exception of the Michaelis-Menten constant for phosphate (K PO 4 ), a variable not included in Arndt et al. (2009), and the maximum specific photosynthetic rate (P B max ), which is constant in the stationary simulation and varies with temperature in the transient simulation (see Table 1). A complete list of biogeochemical parameters is presented in Table 4. In the Scheldt estuary, diatoms are the dominant phytoplankton species (e.g., Mulyaert and Sabbe, 1999). Hence, GGP is assumed to be carried out by diatoms only (PHY = DIA). Because of the large anthropogenic influence on the Scheldt estuary, which favors net heterotrophy, nitrogen and phosphorous levels are typically well above  Vanderborght et al., 2007) and silica can be assumed to be the only limiting nutrient for diatom growth . Sediment parameters are calibrated on the basis of SPM observations and by comparing the simulated annual evolution of NPP and sediment concentration with results obtained from the 2-D model. SPM parameter values are provided in Table 5.

Lateral loads and the Rupel network
Lateral inputs from domestic, industrial and agricultural activities are accounted for in the model and are applied in all runs as constant point sources of organic matter, ammonium and nitrate distributed along the estuarine gradient Arndt et al., 2009). Their values and their input locations are given in Table 6. Differences between lateral loads use for stationary and transient simulations mainly reflect the improvement in wastewater treatment in the Scheldt catchment at the end of the 20th century ( Vanderborght et al., 2007).
In addition, C-GEM also accounts for the river network of the Rupel, the most important tributary of the Scheldt (Hellings and Dehairs, 2001) in the form of a simple box model with a volume of about 1.5 × 10 7 m 3 that discharges unilaterally into the main channel at 102 km ( Fig. 5b and c). This approach allows for a better comparison between simulation results and field data. Rupel boundary conditions are listed in Table 3.

Sensitivity study
A sensitivity analysis, using a one factor at a time (OFAT) method, was conducted to assess the influence of model parameter variations on net primary production (NPP), aerobic degradation (R), denitrification (D), nitrification (N), O 2 exchange across the air-water interface (O 2 ex) and net ecosystem metabolism (NEM). The original parameter set adopted by the 2-D model (Arndt et al., 2009) serves as a reference case for the sensitivity study. The sensitivity of spatially and temporally integrated rates to parameter changes is investigated. Table 7 provides an overview of the model parameters, their baseline values, as well as the tested parameter range. Note that the Chézy coefficient is considered as a sediment parameter despite its dual role in hydrodynamics and sediment erosion/deposition dynamics (see Eqs. 7 and 18). Although sediment and biogeochemical parameters, such as for instance the rate constant of organic matter degradation (e.g., Arndt et al., 2013), can vary over orders of magnitude, here they are varied arbitrarily over a range of ±50 % of their baseline value because our aim is to test the relative sensitivity of the model response and establish priorities for future research rather than to assess the variability arising from different ranges in parameter values reported in the literature. On the other hand, geometric parameters (convergence length and depth) are varied over a smaller range (±10 % and ±20 %, respectively) since they can be constrained on the basis of observations.  Vanderborght et al. (2007). b from Billen and Garnier (1997). All other values are from Arndt et al. (2009 Table 5. Calibrated sediment parameters used in C-GEM for stationary and transient simulations. Note that a linear variation is applied to the Chézy coefficient (C) and the critical shear stress for erosion and deposition (τ cr ) between 100 km and 158 km is applied. Numerical values assigned to C 158 km and τ cr,158 km correspond to their value imposed at the estuarine upper boundary.

Hydrodynamics and transport
The simulated longitudinal profile of the tidal amplitude (Fig. 6) reveals the characteristic features of a funnel-shaped, macro-tidal estuary (Savenije and Veling, 2005;Nguyen, 2008). In the lower, tidally dominated part of the estuary, channel convergence results in the amplification of the tidal wave. However, the influence of fluvial energy progressively increases as the tidal wave moves upstream. It acts primarily through bottom friction and induces a dampening of the tidal amplitude (Fig. 6). High water levels are less influenced by friction than low water levels and thus contribute less to the decrease in tidal range. Figure 6 shows that the model slightly underestimates the tidal amplitude in the saline estuary (km < 100), while it overestimates the tidal amplitude in the tidal river. In particular, mean relative Table 6. Lateral loads (mmol s −1 ). For more information, refer to Vanderborght et al. (2007) and Arndt et al. (2009 Sediment parameters E = erosion coefficient (mg m 2 s −1 ) Variable ±50 τ cr = critical shear stress for erosion and deposition (N m −2 ) Variable ±50 C = Chézy coefficient (m 1/2 s −1 ) Variable ±50 W s = settling velocity (m s −1 ) 1 × 10 −3 ±50 Primary production parameters α = photosynthesis efficiency (m 2 s s −1 µE −1 ) 5.8 × 10 −7 ±50 k excr = excretion constant (-) 0.03 ±50 k growth = growth constant (-) 0.3 ±50 k maint = maintenance rate constant (s −1 ) 9.26 × 10 −7 ±50 k mort = mortality rate constant (s −1 ) 7.1 × 10 −7 ±50 Biogeochemical reaction rates k nit = nitrification rate constant (µM N s −1 ) 1.5 × 10 −4 ±50 k ox = aerobic degradation rate constant (µM C s −1 ) 2.0 × 10 −4 ±50 k denit = denitrification rate constant (µM C s −1 ) 1.0 × 10 −4 ±50 O 2 air exchange parameter k flow = current component for piston velocity (m s −1 ) Variable ±50 differences between observed and simulated tidal amplitudes are smaller than 5 % and 22 % in the saline estuary and in the tidal river, respectively. Discrepancies between model results and observations are mainly related to the seasonal and interannual variability in freshwater discharge, which cannot be captured by the steady-state simulation. Part of the deviation may also arise from the use of an idealized geometry, which does not resolve the complex bathymetry of the Scheldt estuary that is characterized by deep tidal channels and shallow tidal flats.
The dispersion coefficient D is quantified according to Eq. (11) using the idealized geometry of the Scheldt estuary (shown in Fig. 5b and c and summarized in Table 2) and assuming a constant freshwater discharge of 39 m 3 s −1 corresponding to the mean value for which observations were available. These assumptions yield a Van der Burgh coefficient K of 0.39. Figure 7 illustrates the evolution of the dispersion coefficient D along the estuarine gradient and reveals a dome-shaped profile with a maximum value of about 124 m 2 s −1 near the estuarine mouth that reduces to 0 in the tidal river.  The longitudinal distribution of salinity is controlled by the balance between upstream dispersion and downstream advection (Savenije, 2005(Savenije, , 2012. The steady-state salinity profile (Fig. 8) also follows a dome-shaped distribution characterized by a small salinity gradient at the estuarine mouth. This shape is typical of funnel-shaped estuaries (e.g., Savenije, 2005). Simulation results (Fig. 8) agree well with salinity distributions observed under similar hydrodynamic conditions (Regnier et al., 1998).

SPM and biogeochemistry
The estuarine SPM distribution is mainly controlled by the total dissipation of tidal and fluvial energies (Chen, 2003;Chen et al., 2005;. Although SPM concentrations in the Scheldt estuary show a very patchy pattern in time and space due to their high sensitivity to changes in physical forcing conditions , a typical trend, which relates to three well-defined energy regimes along the longitudinal axis of the estuary, can be identified (e.g., Jay et al., 1990;Dalrymple et al., 1992;Arndt Figure 8. Comparison between salinity measurements (Regnier et al., 1998) and simulated longitudinal distribution of the tidally averaged salinity for a mean tidal amplitude of 3.7 m, modeled using a constant freshwater discharge Q = 39 m 3 s −1 . et al., 2007). In the lower estuary, where mechanical energy is almost exclusively provided by the tide, observed SPM concentrations are generally low and range between 0 and 150 mg L −1 . Moving upstream, channel convergence induces an upstream increase in energy dissipation and the associated intensification in tidal amplitude (e.g., Fig. 6) triggers an increase in SPM concentrations from the mouth to the turbidity maximum zone (TMZ), where maximum values of up to 600 mg L −1 can be observed . The exact location of the TMZ shifts in response to the tidal excursion and the river discharge and is generally found between 60 km and 100 km (e.g., Wollast and Marijns, 1981;Chen et al., 2005). Beyond the TMZ, friction progressively reduces the tidal influence (Horrevoets et al., 2004) and energy dissipation becomes progressively controlled by the seaward flux of fluvial energy. At the so-called balance point, where both contributions are of similar but low magnitude, low SPM concentrations are observed (0-250 mg L −1 , Van . Upstream of the balance point, close to the estuarine upper limit, the magnitude of the riverine input flux controls the SPM concentration (Chen et al., 2005). The simulated steady-state longitudinal SPM profile (Fig. 9) is in agreement with this general pattern. Direct comparison with an observed SPM profile is however not possible because the simulated steadystate conditions do not reproduce a situation observed in the field. SPM concentrations are strongly controlled by local exchange processes with the estuarine bed. Hence, already small changes in the physical forcing, as well as their history, exert a large impact on local SPM concentrations and result in large local fluctuation, rendering a direct comparison of simulation results and the range of observed SPM values not very informative.
Longitudinal steady-state profiles of oxygen, ammonium, nitrate and silica generally show a good agreement with measured data (Fig. 10). These profiles are discussed in detail in  Vanderborght et al. (2007) and some key features are briefly summarized here. In the tidal river, high riverine loads of carbon and reduced nitrogen drive intense heterotrophic processes rates and, thus, trigger low oxygen concentrations (Fig. 10a). Further downstream, the decrease in consumption rates and the increase in air-water exchange fluxes result in a progressive increase in O 2 levels. In contrast, nutrient concentrations are generally high in the upper tidal reaches, but decrease along the estuarine gradient due to the progressive dilution and the decrease in autotrophic process rates (Fig. 10b-d). A short increase in NH 4 (Fig. 10b) and a concomitant decrease in O 2 and NO 3 concentrations (Fig. 10a and c) around 100 km reflect an increase in heterotrophic process rates that is mainly driven by the influence of the Rupel tributary.
Despite the overall agreement between model results and observations, Fig. 10 also reveals some discrepancies. For instance, the simulated O 2 , NH 4 and dSi gradients are steeper than in the observed profiles and simulated concentration minima are located further downstream. Part of this discrepancy can be explained by the highly dynamic nature of the estuarine environment and the strong inter-annual variability (e.g., . Steady-state simulations forced with average summer conditions do not resolve such complex dynamics (e.g., Regnier et al., 1997;Arndt et al., 2009). Nevertheless, steady-state simulations results show that, despite numerous simplifying assumptions during model set-up, C-GEM is able to capture the general features of the biogeochemical dynamics in the Scheldt estuary.

Biogeochemical functioning
Long-term seasonal to decadal biogeochemical dynamics or system-wide biogeochemical indicators, such as the NEM, are difficult to assess through observations only. Their quantification requires the application of fully transient RTMs to complement field measurements (Regnier et al., 2013b). The quantification of such system-wide biogeochemical Figure 10. Comparison between longitudinal distributions of field data averaged over the period May-September for the years 1990-1995 (dots; vertical bars correspond to the standard deviation) and steady-state maximum and minimum O 2 , NH 4 , NO 3 and dSi concentrations over a tidal cycle (solid line). Physical conditions are summarized in Table 2, boundary conditions and external forcings are summarized in Table 3 and parameters are listed in Tables 4 and  5. indicators provides an important integrative measure for the overall performance of C-GEM.
Therefore, the simulated annual evolution of spatially integrated NPP, aerobic degradation, denitrification, total heterotrophic degradation (denitrification and aerobic degradation), nitrification rates and NEM are compared to those obtained with the highly resolved 2D-RTM by Arndt et al. (2009). The integration is performed over the entire estuarine domain. Figure 11 shows that C-GEM captures the main seasonal evolution of biogeochemical process rates. Autotrophic process rates are low during winter and autumn, but increase to a maximum in early summer (Fig. 11a), when Figure 11. Comparison between annual evolution of biogeochemical rates modeled by C-GEM (solid line) and the 2D-RTM (dashed line) by Arndt et al. (2009). favorable temperature and light conditions, large nutrient inventories and low turbidities drive high in situ NPP rates. Heterotrophic process rates and nitrification are high during both winter and summer months (Fig. 11b-e). These high rates are sustained by high riverine inputs in winter and elevated ambient temperatures in summer (Fig. 11b-e).
In addition, Fig. 11b-e show that nitrification, denitrification and aerobic degradation are tightly coupled. For instance, high nitrification rates (Fig. 11e) are supported by the ammonium supplied by high aerobic degradation rates (Fig. 11b). Moreover, during summer, high nitrification and aerobic degradation rates result in a depletion of oxygen and thus contribute to the increase in denitrification rates (Fig. 11c). Furthermore, heterotrophic degradation processes are enhanced by the supply of organic matter derived from dead phytoplankton in the aftermath of the summer algae bloom (Fig. 11d). Model results indicate that the heterotrophic degradation in the Scheldt is largely dominated by the aerobic organic matter degradation. The simulated NEM profile (Fig. 11f) closely follows the total heterotrophic degradation rate profile (Fig. 11d). During summer, the influence of heterotrophic processes on NEM is partly compensated for by primary production rates (Fig. 11a), but the simulated NEM remains negative throughout the year, reflecting the heterotrophic nature of the estuary.
Although the idealized simulation performed with C-GEM captures the general seasonal pattern of system-wide process rates, Fig. 11 also reveals discrepancies between C-GEM and 2-D simulation results. Whole-estuarine aerobic degradation rates are lower than those obtained with the 2-D model during the first period of the year (day < 60), while differences in NPP rates are more pronounced during the summer months. Moreover, C-GEM simulates lower nitrification and denitrification rates. These discrepancies can be traced back to differences in simulated water depth, estuarine circulation, residence times and/or turbidity. The idealized geometry provides a highly simplified representation of the complex estuarine bathymetry with deep tidal channels and extensive intertidal mud flats. As a consequence, C-GEM ignores the cross-sectional variability in water depth, circulation and, thus, residence times. For instance, C-GEM underestimates residence times in the upper reaches and, therefore, simulates lower biogeochemical rates. These cross-sectional variabilities in residence time, turbidity and residual circulation also exert an important influence on summer NPP rates. Twodimensional simulation results highlight the pronounced differences between NPP rates in tidal channels and intertidal flats (e.g., , a feature that cannot be resolved by the idealized bathymetry of C-GEM. The simplification of the estuarine bathymetry may thus also explain the observed differences in simulated NPP rates. In addition, C-GEM simulates lower nitrification rates but slightly higher aerobic degradation rates during the summer months. These discrepancies probably arise from different estimates of the transient overlap in TOC and O 2 for aerobic degradation and in NH 4 and O 2 for nitrification, which induce different values of the Michaelis-Menten terms involved in these two processes. Despite these discrepancies, integrated biogeochemical reaction rates estimated with C-GEM concur well with the 2-D results. Annually integrated biogeochemical process rates are compared in Fig. 12. C-GEM slightly underestimates nitrification, denitrification and aerobic degradation rates, as well as the oxygen exchange with atmosphere with a relative error of 36 %, 24 %, 4 % and 17 %, respectively. Simulated NPP rates are slightly higher with a relative error of 23 %, while the simulated NEM value is slightly lower by about 10 %. Thus, all integrated measures fall within the same order of magnitude. Figure 13 illustrates the sensitivity of biogeochemical process rates to parameter variations (Table 7). Geometrical parameters generally exert an important influence on all integrated process rates (Fig. 13a). For instance, a 10 % variation in convergence length (LC) triggers large changes (> 15 %) in NPP, aerobic degradation and nitrification rates and also exerts a somewhat smaller influence (∼ 10 %) on denitrification and air/water exchange rates. This difference is system specific and can be explained by the effect of convergence length on estuarine volume and residence time (Eqs. 1 and 2). Fixing the estuarine width, B, at the inland limit, as done during this sensitivity test and following Eq. (2), a shorter convergence length increases the volume and the residence time in the estuarine system, a central parameter that in turn promotes all processes and increases their biogeochemical rates (Fig. 13a). A larger convergence length has the opposite effect on the rates. Denitrification is the most sensitive process to variations in water depth, H (Fig. 13a). The volumetric reduction of the estuary induced by a shallower water depth translates into a decrease in aerobic degradation, denitrification and nitrification rates. The large reduction in denitrification may also be related to the positive effect of shallow water depth on oxygen exchange rate, which, inducing an increase in O 2 levels in the water column, strongly inhibits denitrification. The increase in NPP rates to both positive and negative relative variations in water depth highlights the strong dependence of this process on the underwater light field. Shallow waters increase the photic depth to water depth ratio, while deep waters decrease light attenuation through a dilution effect on suspended sediment concentrations (results not shown; Chen et al., 2005;Desmit et al., 2005). Despite their strong influence on biogeochemical processes, estuarine geometric features do not limit the application of C-GEM to data-poor estuarine systems, since they can be readily extracted from nautical charts or maps.

Sensitivity analysis
Integrated NPP rates are also highly sensitive to variations in primary production and SPM parameters ( Fig. 13 and e), while they are not affected by variations in gas exchange parameters and biogeochemical rate constants ( Fig. 13c and d). This reflects the fact that underwater light field rather than nutrient availability controls NPP. As a consequence, NPP is also sensitive (> 66 %) to changes in the Chézy coefficient, C, which affects SPM dynamics and thus the light availability, and in phytoplankton parameters (Fig. 13e). Variations in the maintenance rate constant exert the largest influence on system-wide integrated NPP (> 77 %) because the maintenance term is directly proportional to the total phytoplankton concentration (see Table 1). Although both growth and excretion are linearly proportional to gross primary production, the integrated NPP only responds to variations in the growth constant because its value is one order of magnitude larger than that of the excretion constant. Photosynthesis efficiency also has a significant effect on NPP variations as shown in Fig. 13b and integrated rates vary by as much as 53 %. Overall, simulation results indicate that NPP rates are most sensitive to uncertainties in the Chézy coefficient and the rate constant for maintenance. These parameter values are difficult to determine and are generally obtained from model calibration. In particular, the Chézy coefficient is never measured directly, while the maintenance term generally varies across different phytoplankton groups. Heterotrophic and oxygen exchange rates are most sensitive to variations in biogeochemical reaction rate constants (Fig. 13d) and to a lesser degree variations in the current component for the piston velocity (Fig. 13c) and in the Chézy coefficient (Fig. 13e). On the other hand, NPP parameters exert virtually no effect (Fig. 13b), emphasizing the strongly heterotrophic character of the estuarine system (Figs. 11 and 12). While aerobic degradation and nitrification show only small variations (< 10 %) associated with changes in the current contribution to the piston velocity, Fig. 13c confirms the sensitivity of denitrification to the O 2 exchange process at the air-water interface and to the O 2 level in water. To a variation in the gas exchange rate corresponds an opposite variation in denitrification. Hence, estimates of these two processes require a good resolution of the flow velocity field and the water depth in order to constrain well the flow component for the piston velocity.  Simulation results emphasize the fact that a robust quantitative estimation of the estuarine biogeochemical functioning calls for well-constrained biogeochemical rate constants. However, these constants are difficult to constrain as they implicitly account for factors that are not resolved in C-GEM, such as the structure and the abundance of the microbial community or a complete description of the environmental conditions within the estuarine systems. The lack of an objective framework for model parameterization and the limited transferability of system-specific parameter values potentially may limit the generic approach of C-GEM. Hence, a sensitivity study should be an integral part of the model application and can help to estimate uncertainties in predicted rates.
Despite the relatively large variations applied in the sensitive runs, the estuary never becomes net autotrophic and NEM always remains negative within the range −6235 and −10461 kmol C d −1 . Figure 14 identifies the parameters that lead to a NEM variation larger than 5 %. Since the NEM is always negative, a positive relative variation in its value implies a more heterotrophic status of the system. These results again highlight the fact that an increase in volume and, thus, in residence time (induced by a decrease in LC and by an increase in depth; see above) and in the aerobic degradation rate constant induce a more negative NEM, while an increase in LC and a decrease in depth and aerobic degradation constant rate have the inverse effect. A comparison of Fig. 14 with Fig. 13a and d shows that variations in NEM closely follow the variations in aerobic degradation, induced by these three parameters (LC, H , k ox ), reflecting the overall dominance of this process in the NEM estimates.
Note that, while the general pattern emerging from this sensitivity study is valid across systems, the quantitative influence of parameter variations is highly system dependent. For instance, prismatic systems with a longer convergence length and, thus, a stronger fluvial influence are characterized by much shorter residence times. Therefore, integrated biogeochemical reaction rates in prismatic systems will, likely, reveal a much weaker response to variations in biogeochemical parameters than in funnel-shaped systems.

Scope of applicability and model limitations
Site-specific, multi-dimensional models generally perform satisfactorily at reproducing the biogeochemical dynamics of estuarine systems, but are highly demanding in terms of data and numerical requirements. At the other end of the model spectrum, box models are very efficient, but generally fail to resolve the spatial and temporal variability of estuarine systems and are not well suited for model-data comparison. However, our ability to assess the role of the estuarine environment for global biogeochemical cycles and greenhouse gas budgets, as well as their response to ongoing global change requires tools that are computationally efficient and can extrapolate knowledge from well-studied to data-poor systems, while at the same time resolving the most important hydrodynamic and biogeochemical processes and scales. The new C-GEM one-dimensional model proposed here is such a computational tool. It represents a valid compromise between performance and computational efficiency and reduces data requirements by using an idealized representation of the estuarine geometry. Its scope of applicability covers the entire range of alluvial estuaries, from tidally dominated systems with a large tidal range and low river discharge to fluvial-dominated systems characterized by significant freshwater input (Regnier et al., 2013b). It can be used to resolve the complex process interplay that drives the estuarine biogeochemical dynamic and to quantify estuarine carbon and nutrient budgets. In addition, the computational efficiency of C-GEM offers the possibility of simulating simultaneously the biogeochemical dynamics of a large number of estuaries and the contiguous coastal ocean. Although not considered so far, C-GEM could theoretically be applied to the tidally influenced, inland sections of very large river systems (e.g., the Amazon). The value of such an application is however questionable because very large rivers contribute disproportionally to the overall land to ocean carbon fluxes and might thus deserve a dedicated model. In addition, their tight estuarine-continental shelf coupling and importance, as well as the complex multi-dimensional dynamics of their coastal plumes, require a multi-dimensional model representation. Numerous models have already been developed for these systems (e.g., Gallo and Vinzon, 2005;Denamiel et al., 2013), and in the future they could be explicitly represented in high-resolution Earth system models (Bauer et al., 2013). In contrast, for the smaller alluvial estuarine systems, mechanistically rooted upscaling strategies need to be designed to constrain their roles in the global carbon cycle better (Bauer et al., 2013), and C-GEM is a tool of choice in this context. However, C-GEM is associated with a certain degree of simplification and, therefore, is characterized by some limitations. Currently, the model does not include a benthicpelagic exchange module. Hence, its application is not recommended for estuaries that are subject to an intense benthic-pelagic coupling. The resulting lack of a representation of particulate organic carbon burial might result in an overestimation of estuarine organic carbon export fluxes to the coastal ocean. The most important hurdle towards generalization arises from the lack of an objective, global framework for SPM and biogeochemical process parameterization. These parameters implicitly account for a large number of controlling factors that are usually not explicitly resolved in estuarine models. They are typically derived by model calibration on the basis of observations, and their transferability to other systems is thus limited. Comprehensive sets of model parameters are now available for some estuaries of the world, such as those in Europe, North America and Australia, but are essentially missing in the tropical and polar regions (Regnier et al., 2013b). The limited transferability of model parameters and the lack of observational data call for the creation of a global data set of estuarine sediment and biogeochemical parameters on which a statistical analysis is strongly desirable in order to identify common trends and possible relationships between parameters and control factors, such as latitude, catchment characteristics and anthropic pressure.

Conclusion and perspectives
The model developed in this study represents a first attempt to quantify the biogeochemical dynamics in estuaries using a one-dimensional reactive transport model that relies on idealized geometries to support the estuarine hydrodynamics and transport. Despite its highly simplified geometric support, C-GEM captures the dominant features of the biogeochemical behavior along a complex system as the Scheldt estuary (BE/NL) and the system-wide integrated reaction rates for the main biogeochemical pelagic processes are comparable with those obtained using a high-resolved site-specific 2D-RTM. A sensitivity analysis, based on the OFAT method, has been performed in order to assess the importance of the internal parameters on the estuarine biogeochemistry. It reveals that geometry and hydrodynamics exert a strong first-order control on the biogeochemical functioning and therefore supports our hypothesis that the estuarine response is a systemspecific attribute that cannot be reduced to a simple and direct signal response, such as the nutrient filtering capacity and the residence time relationship proposed, for instance, by Nixon et al. (1996). Results also provide a rational support to identify the model parameters that are the most sensitive with respect to integrative measures, such as the NEM, and emphasize the need for a global compilation of estuarine sediment and biogeochemical parameters. In addition, such a compilation could help identify trends between parameter values and control factors, such as climate, catchment properties and anthropic pressure, and compensate for the current lack of an objective, global framework for parameterization in data-poor areas.
The structure of C-GEM, which optimizes the ratio between the number of parameters and the availability of data, provides an easy and cost-efficient tool that can be used to quantify the biogeochemical dynamics of estuaries and to forecast their response to combined climate and environmental changes over the coming century. In the future, C-GEM could be applied in combination with, for example, Global-NEWS2 models  for riverine inputs, to a wide range of estuarine systems characterized by different climatic regimes, geometries and chemical loadings. This, together with the compilation of a global data set for sediment and biogeochemical parameters, could help in the quantification of estuarine biogeochemical cycles on regional and global scales.

Code availability
The C-GEM source code related to this article is provided as a supplementary package together with a Read Me file, where hardware and software requirements, source code files and model output file management are fully described.
The Supplement related to this article is available online at doi:10.5194/gmd-7-1271-2014-supplement.