Linking biogeochemistry to hydro-geometrical variability in tidal estuaries: a generic modeling approach.

. This study applies the Carbon-Generic Estuary Model (C-GEM) modeling platform to simulate the estuarine biogeochemical dynamics – in particular the air– water CO 2 exchange – in three idealized tidal estuaries characterized by increasing riverine inﬂuence, from a so-called “marine estuary” to a “riverine estuary”. An intermediate case called “mixed estuary” is also considered. C-GEM uses a generic biogeochemical reaction network and a unique set of model parameters extracted from a comprehensive literature survey to perform steady-state simulations representing average conditions for temperate estuaries worldwide. Climate and boundary conditions are extracted from published global databases (e.g., World Ocean Atlas, GLORICH) and catchment model outputs (GlobalNEWS2). The whole-system biogeochemical indicators net ecosystem metabolism (NEM), C and N ﬁltering capacities (FC TC and FC TN , respectively) and CO 2 gas exchanges (FCO 2 ) are calculated across the three idealized systems and are related to their main hydrodynamic and transport characteristics. A sensitivity analysis, which propagates the parameter uncertainties, is also carried out, followed by projections of changes in the biogeochemical indicators for the year 2050. Results show that the average C ﬁltering capacities for baseline conditions are 40, 30 and 22 % for the marine, mixed and riverine estuary, respectively, while N ﬁltering capacities, calculated in a similar fashion, range from 22 % for the marine estuary to 18 and 15 % for the mixed and the riverine estuaries. Sensitivity analysis performed by varying the rate constants for aerobic degradation, denitriﬁcation and nitriﬁcation over the range of values reported in the literature signiﬁcantly widens these ranges for both C and N. Simulations for the year 2050 suggest that all estuaries will remain largely heterotrophic, although a slight improvement of the estuarine trophic status is predicted. In addition, our results suggest that, while the riverine and mixed systems will only marginally be affected by an increase in atmospheric p CO 2 , the marine estuary is likely to become a signiﬁcant CO 2 sink in its downstream section. In the decades to come, such a change in behavior might strengthen the overall CO 2 sink of the estuary–coastal ocean continuum.


Introduction
Located at the interface between land and ocean, estuaries are highly dynamic ecosystems, which process variable fractions of land-derived inputs of carbon (C) and nutrients (N, P, Si) through a wide range of chemical and biological processes (Alongi, 1998;Crossland et al., 2005;Bianchi, 2007).Gaseous species, in particular CO 2 , are produced and further exchanged with the atmosphere as a result of these intense biogeochemical dynamics.Recent compilations of observed data reveal that the vast majority of estuarine systems are net CO 2 emitters (Laruelle et al., 2010(Laruelle et al., , 2013;;Borges and Abril, 2011) and are responsible for a global CO 2 outgassing ranging between 0.15 and 0.25 PgC yr −1 (Cai, 2011;Bauer et al., 2013;Regnier et al., 2013a).This flux corresponds to about 20-25 % of the riverine carbon inputs and is of similar magnitude to the global uptake of CO 2 by continental shelves (Laruelle et al., 2014).Estuaries are thus important modulators of the carbon and associated bio-elements fluxes from the land to the open ocean (e.g., Jahnke, 1996;Billen and Garnier, 1997;Gattuso et al., 1998;Lancelot et al., 2005;Mackenzie et al., 2005;Arndt et al., 2009Arndt et al., , 2011a; Laru-C.Volta et al.: Linking biogeochemistry to hydro-geometrical variability in tidal estuaries elle et al., 2009;Borges and Abril, 2011;Cai, 2011;Bauer et al., 2013;Regnier et al., 2013a).However, the extent to which their biogeochemical dynamics and thus their role in the global cycles will change in the future in response to anthropogenically driven changes in land use, climate and atmospheric CO 2 remains poorly known.
Over the past 30 years, highly resolved, process-oriented and often multi-dimensional models have helped disentangle and quantify estuarine biogeochemical dynamics (e.g., Thomann and Fitzpatrick, 1982;Lung and Paerl, 1988;Regnier et al., 1997Regnier et al., , 1999;;Margvelashvili et al., 2003;Baklouti et al., 2011;Cerco, 2000;Arndt et al., 2007Arndt et al., , 2009;;Mateus et al., 2012).Most of these studies have focused on specific estuarine systems, and comparative studies covering the wide range of estuarine systems are limited.Therefore, a quantitative evaluation of the role of estuaries in the global biogeochemical cycles and their potential response to global change are characterized by large uncertainties (Hobbie, 2000;Borges and Abril, 2011;Laruelle et al., 2013;Regnier et al., 2013a).This lack can be partly attributed to the high data requirements for model calibration and validation, as well as the high computational demand, which have prevented their application to regional and/or global scales (Bauer et al., 2013).In addition, the limited availability of comparative studies compromises the identification of common patterns across the wide range of estuarine types (Geyer et al., 2000;Hobbie, 2000;Borges and Abril, 2011;Regnier et al., 2013b).Moreover, the diagnostic modeling of the CO 2 dynamics has so far only been performed for a temperate estuary in Europe (the Scheldt; Vanderborght et al., 2002), although observational data are now available for more than 100 land-ocean transition systems (Chen et al., 2013;Laruelle et al., 2013).To our knowledge, prognostic simulations of greenhouse gas emissions in estuaries are currently lacking.
The objectives of this paper are thus to develop a unified modeling approach to identify similarities and differences in the biogeochemical dynamics across different estuarine systems and to explore quantitative relationships between the estuarine biogeochemical functioning and a set of key hydro-geometrical characteristics.The overarching goal is to enhance our ability to transfer information from wellconstrained estuarine systems to poorly surveyed systems, to improve upscaling strategies and to enable projections.In the first part of this paper, the description of the conceptual framework underlying our generic approach is provided.In particular, we build on the mutual dependency between geometry and hydrodynamics (Savenije, 1992(Savenije, , 2005(Savenije, , 2012) ) and further hypothesize that hydrodynamics exert a strong control on biogeochemistry in alluvial estuaries (e.g., Alpine and Cloern, 1992;Arndt et al., 2007;Volta et al., 2014).Next, three idealized systems, characterized by variable riverine influence and covering the main hydro-geometrical features of tidal alluvial estuaries, are modeled using the recently developed C-GEM modeling platform (Volta et al., 2014).These systems are designed to represent a tidal estuary dominated by marine characteristics, a tidal estuary dominated by its riverine characteristics, and an intermediate case (so-called mixed system).Here, C-GEM uses a generic biogeochemical reaction network and a unique parameter set extracted from a large literature survey of more than 40 local modeling studies.Steady-state simulations representing average conditions for temperate estuaries worldwide are performed for the present decade, as well as for the mid-21st century.The whole-system biogeochemical indicators net ecosystem metabolism (NEM), C and N filtering capacities (FC TC and FC TN , respectively) and CO 2 gas exchanges (FCO 2 ) are calculated across the three idealized systems and are related to their main hydrodynamic and transport characteristics.A sensitivity analysis is also carried out to assess the sensitivity of the estuarine biogeochemical functioning to parameter uncertainties synthesized in the present study.Finally, the main findings are summarized and their significance and limitations are critically analyzed in the context of improving upscaling strategies for regional and global CO 2 emissions of estuaries.

Theoretical support
Most tidal estuaries are alluvial estuaries (Regnier et al., 2013b), which are defined as estuarine systems with movable beds, consisting of material from marine and terrestrial origin, and a measurable freshwater inflow (e.g., Hobbie, 2000;Savenije, 2005Savenije, , 2012)).The global distribution of alluvial estuaries is roughly equivalent to that of the tidal estuaries as defined in the estuarine coastal typology of Dürr et al. (2011).The approach developed here builds on this intercompatibility, already mentioned in Regnier et al. (2013b).In tidal estuaries, two different zones can be identified along their longitudinal gradient: a lower zone, the so-called saline estuary, whose dynamics are essentially controlled by mixing with marine waters, and an upstream zone, referred to as the tidal river, where processes are mainly driven by the freshwater input from rivers (e.g., Jay et al., 1990;Dalrymple et al., 1992;Regnier et al., 2013b).Alluvial estuaries are characterized by a mutual dependency of the estuarine geometry and hydrodynamics.The magnitude of the water flow entering or leaving the estuarine channel is entirely controlled by its shape (Pethick, 1984).In turn, the water movement, mainly driven by tides and freshwater discharge, leads to a redistribution of the unconsolidated sediments that determines the shape of the estuary.This dynamic interplay between hydrodynamics and morphology results in a continuum of estuarine shapes that cover the entire spectrum between two end-member cases: (1) systems with rapidly converging banks towards the land and (2) systems characterized by parallel banks (Savenije, 1992).Although the exact shape of alluvial estuaries can vary, they nonetheless show common geometrical characteristics that are compatible with an idealized representation of the estuarine geometry (Savenije, 2005(Savenije, , 2012)).The tidally averaged estuarine width B (in m) typically shows an exponential decrease in the landward direction (e.g., Pethick, 1992;Lanzoni and Seminara, 1998;Savenije, 1992Savenije, , 2005)), while the tidally averaged estuarine depth H (in m) remains nearly constant along the estuarine gradient (Savenije, 1992(Savenije, , 2005(Savenije, , 2012)): where x (in m) is the distance from the estuarine mouth, B 0 and H 0 denote the estuarine width and depth (in m) at the estuarine mouth (x = 0), respectively, and b is the width convergence length (in m), defined as the distance over which the estuarine width reduces to 37 % (e −1 ) of its value at the mouth.The shape of alluvial estuaries can, thus, be fully defined by the width convergence length, b, and the channel depth, H 0 (Savenije, 2012).The ratio between these two geometrical key parameters is generally defined as the dimensionless estuarine shape number, S (Savenije, 1992): The hydrodynamic characteristics of alluvial estuaries may, in turn, be directly related to their main hydrodynamic forcings, such as tidal influence and freshwater inflow (Wright et al., 1973) by means of the dimensionless hydrodynamic Canter-Cremers estuary number N (Simmons, 1955): where Q b denotes the bankfull discharge or, in other words, the temporary maximum river discharge (in m 3 s −1 ) that is associated with a state of maximum flow velocity and, thus, to the maximum ability to shape the estuary (Savenije, 2012).
T is the tidal period (in s).P is the tidal prism (in m 3 ), which represents the volume of saline water entering the estuary over a tidal period, T , and can be approximated by the product of the cross-sectional area at the estuarine mouth A 0 (in m 2 ) (A 0 = B 0 • H 0 ) and the tidal excursion length E (in m), defined as the maximum distance a water particle travels during a tidal period T .By using a regression analysis based on measurements from 16 alluvial estuaries worldwide, Savenije (1992) showed that the estuarine shape number (S; Eq. 3) is related to the hydrodynamic Canter-Cremers number (N ; Eq. 4) through a power law relationship, resulting in the dimensionless hydro-geometrical relationship: Certain hydro-geometrical characteristics, such as the tidal period T , the tidal excursion E and the estuarine depth at the mouth H 0 can be approximated by characteristic values.
For instance, E is usually close to 10 km for a semi-diurnal (T ≈ 12 h) tidal estuary, while an alluvial estuary flowing in a coastal plain generally reveals a tidally averaged depth (H ) of about 7 m (Savenije, 1992;Volta et al., 2014).Note that, for estuaries in which the freshwater discharge is known, a new formulation proposed by Gisen and Savenije (2015) can be used to estimate the estuarine depth.On the other hand, other characteristics, such as the bankfull discharge Q b , the width convergence length b, and the cross-sectional area at the estuarine mouth A 0 , are system-specific characteristics, depending on the hydrological regime of the estuarine watershed and on the local balance between riverine and marine energies.Therefore, they cannot be easily approximated.
The width convergence length, b, and the bankfull riverine discharge, Q b , can thus be considered key parameters for defining the hydro-geometrical character of alluvial estuaries and for predicting their salt intrusion profiles (Savenije, 2005(Savenije, , 2012)).Important transport and mixing properties can be directly related to these hydro-geometrical characteristics.Hence, Eq. ( 5) not only provides a universal theoretical framework for analyzing the tight link between geometry, hydrodynamics and transport but also offers a theoretical basis for a classification of alluvial estuaries.Savenije (2005Savenije ( , 2012) ) identified two main estuarine types, which differ in terms of geometrical features, hydrodynamics characteristics and salt intrusion patterns: 1. funnel-shaped (or marine-dominated) estuaries that are typically characterized by a short width convergence length, b, and thus rapidly converging banks, a low freshwater discharge, a dome-shaped salinity profile with a small salinity gradient at the estuarine mouth and an intrusion of saltwater far upstream; 2. prismatic (or river-dominated) estuaries that are characterized by a theoretically infinite width convergence length, b, and, thus, a constant channel width, a high river discharge and a steep salt intrusion profile with a strong salinity gradient close to the estuary mouth and a short salt intrusion length.
These estuarine classes represent the extreme ends of the wide range of estuarine hydro-geometrical properties.As a consequence, a series of systems, which show intermediate conditions and fall in between the funnel-shaped and the prismatic end-member cases, can be hypothesized between them (Savenije, 1992).Physical and hydrodynamic characteristics as they relate to the estuary shape are synthesized in nition of the first-order control of hydrodynamics on estuarine biogeochemistry (e.g., Alpine and Cloern, 1992;Nixon et al., 1996;Arndt et al., 2007Arndt et al., , 2009;;Laruelle et al., 2009), on the other hand, allow hypothesizing that each estuarine type might respond in a specific way to the tight coupling between geometry, hydrodynamics, transport and biogeochemistry.Hence, important biogeochemical properties in estuaries might, just as salinity profiles, be predicted on the basis of hydro-geometrical features (Fig. 2; Volta et al., 2014).

Representative estuarine systems
In this study, we explore the link between biogeochemical dynamics and key hydro-geometrical properties in three idealized, tidal alluvial estuaries characterized by variable marine/riverine influence by means of a reactive-transport model.For this purpose, three idealized geometries are defined to be representative of the two extreme classes and the intermediate types as described in Sect.2.1 (marine-and riverine-dominated estuaries and intermediate cases).The width convergence length, b, recognized as a shape and hydrodynamic key parameter (see Sect. 2.1), is used to discriminate between the three estuarine types.First, a reference estuary, characterized by an idealized geometry resembling that tested in Volta et al. (2014), is defined.Then, its width convergence length (b = 30 km) is decreased and increased by 50 % in order to intensify the marine (b = 15 km) and the riverine (b = 45 km) character of the system, respectively.This allows defining two other idealized systems, which can be regarded as representative of the marine and river-dominated estuarine classes and between which the reference estuary can be considered an intermediate case.Henceforth, according to Regnier et al. (2013b) and Volta et al. (2014), the marine-dominated estuary, the reference case and the riverine-dominated estuary will be referred to as the marine, the mixed and the riverine estuary, respectively.The estuarine width at the seaward limit, B 0 , and the estuarine length, EL, of the marine and the riverine estuaries are adjusted so that their total, tidally averaged volume corresponds to that of the mixed estuary (V ≈ 1.5 × 10 9 m 3 ).This allows minimizing the effect of volume variations on the biogeochemical dynamics across the three estuarine types.As a result, the channel width of the marine-dominated estuary reduces, over a distance of 90 km, from 13 830 m at the estuarine mouth to 30 m close to the upper limit, whereas the width of the riverine-dominated estuary decreases, over a distance of 226 km, from 4760 m at the estuary mouth to 30 m at the upper limit.All estuaries are assumed to be coastal plain estuaries.Hence, their tidally averaged water depth is approximated to 7 m (see Sect. 2.1).Furthermore, it is assumed that the three idealized systems are subject to a semidiurnal tidal forcing, thus resulting in an identical tidal excursion length, E, of approximately 10 km (see Sect. 2.1).Based on these geometrical characteristics, Eq. ( 5) can be used to calculate the bankfull discharge, Q b in m 3 s −1 , for each system.New formulations are now available to constrain the bankfull discharge from geometrical parameters (i.e., depth and width; see Gisen and Savenije, 2015) and they could be used in future applications of C-GEM.The geometrical features of the three idealized estuaries are illustrated in Fig. 3 and summarized, together with their hydrodynamic properties, in Table 1.Table 1 also reveals that both the estuarine shape (S, Eq. 3) and the hydrodynamic Canter-Cremers (N, Eq. 4) numbers of the three idealized systems cover the whole range of observed values for temperate tidal estuaries (2500 < S > 6000, N < 0.05; Savenije, 1992).As a consequence, they may be considered as representative of a large range of hydro-geometrical conditions observed in this type of estuaries.Note that higher values of S have been reported for tropical estuaries (Gisen et al., 2015).

Model description
Simulations presented in this study are performed by using the C-GEM modeling platform, fully described and available as supplementary material in Volta et al. (2014).C-GEM dynamically resolves hydrodynamics, transport and pelagic biogeochemistry using the numerical schemes described in Regnier and Steefel (1999) and explicitly accounts for variations induced by tides.Here, we use a version of C-GEM similar to that described in Volta et al. (2014), but its biogeochemical reaction network was extended to include a non-siliceous phytoplanktonic group and an inorganic carbon module.
Table 1.Geometrical and hydrodynamic parameters describing the three idealized estuaries.Following Eq. ( 2), H 0 = H for alluvial estuaries flowing in a coastal plain.V is the total estuarine volume calculated as

Physical model support and hydrodynamic module
The three idealized hydro-geometrical estuarine cases, as described in Sect.2.2, form the physical support and provide the boundary conditions for the hydrodynamic module of C-GEM.The latter resolves the 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 is the time (in s); x is the space (in m); r s is the dimensionless storage water ratio that is typically equal to 1 for idealized representations of estuarine geometries (Davies and Woodroffe, 2010); A is the cross-sectional area (in m 2 ), calculated by the product of the estuarine width B, formulated as in Eq. 1, and the instantaneous water depth H , computed as the sum of the tidally averaged water depth (in m) and the water elevation ζ (x, t) (in m); Q is the cross-sectional discharge (in m 3 s −1 ), calculated by the product of A (in m 2 ) and the flow velocity U (in m s −1 ); g is the gravitational acceleration (in m s −2 ); and C h is the Chézy coefficient (in m 1/2 s −1 ).

Coupling reaction and transport
The coupling of mass transport and chemical reactions is described by using a one-dimensional, tidally resolved advection-dispersion equation for gaseous, solute and solid species in the water column (e.g., Pritchard, 1958): where C i is the concentration of the species i, A is the crosssectional area (in m 2 ) and D is the dispersion coefficient (in m 2 s −1 ), which decreases in the upstream direction according to the Van der Burgh equation (Savenije, 1986) and is dynamically calculated as in Volta et al. (2014).Finally, P i is the sum of volumetric biogeochemical reactions and exchanges through the material surfaces of the estuary (e.g., gas transfer through the air-water interface, erosion and deposition processes), affecting species i.The transport and reaction terms are solved in sequence by applying the operator-splitting approach proposed by Regnier et al. (1998) and a finitedifference scheme on a regular grid ( x = 2000 m) with a time step t = 150 s.In the transport equation, the dispersive and the advective terms are solved by using the semi-implicit Crank-Nicolson algorithm (Press et al., 1992) and the thirdorder Leonard total variation diminishing scheme (Leonard, 1984;Datta-Gupta et al., 1991), respectively.These schemes guarantee mass conservation to within > 1 %.Reaction processes are numerically integrated using the Euler method (Press et al., 1992).A spin-up period of 24 months is imposed.A is calculated as the product of the tidally averaged depth H (H = 7 m along EL) and B and V is the product between A and x.Profiles are obtained by using geometrical parameters reported in Table 1.

Reaction network
The reaction network of C-GEM includes 12 state variables and 10 biogeochemical processes (Table 2).Table 3 summarizes the stoichiometric equation and mathematical formulation of each biogeochemical reaction, while the biogeochemical scheme of the reaction module is shown in Fig. 4. The original version C-GEM 1.0 is described in detail in Volta et al. (2014).Here, it is extended by a second phytoplankton group (nDIA), which represents non-siliceous phytoplanktonic species, such as cyanobacteria, green algae and flagellates, whose growth does not depend on silica concentrations.Furthermore, a new inorganic carbon module, which allows quantifying the estuarine inorganic carbon dynamics, was also implemented.Its mathematical formulation is adapted from Arndt et al. (2011b).pH is considered the master variable for the estimation of dissolution and hydration of CO 2 and is calculated by using the computationally fast approach provided by Follows et al. (2006) by means of an iterative procedure, which accounts for total (TAlk) and carbonate alkalinities.While the model accounts for borate species, contributions from ammonium, fluorine, phosphate, silicate, sulfide and other minor species are neglected because their concentrations are much lower than those of carbonate species (Vanderborght et al., 2002).The apparent equilibrium constants for CO 2 solubility and dissociation of carbonic acid (HCO − 3 ), bicarbonate (CO 2− 3 ), water (H 2 O) and boric acid (B(OH) − 4 ) vary with temperature and salinity according to equations formulated by Cai and Wang (1998) and Dickson (1990).Moreover, new generic temperature-dependent functions for different physiological processes, such as for instance microbial and phytoplankton growth and decay, are implemented in the current version of C-GEM.The implementation of these terms is informed by a large modeling literature survey and further details are provided in the following section (Sect.3).

Sediment parameters
The sediment (SPM) module of C-GEM (Volta et al., 2014) requires specification of six parameters (Table 4).Although assembling a generic data set for sediment parameters is beyond the scope of this study, some generic assumptions are adopted for the Chézy coefficient (C h in m 1/2 s −1 ) and the SPM settling velocity (w s in mm s −1 ).In particular, as suggested by Savenije (2001Savenije ( , 2012)), a C h value of 40 and 60 m 1/2 s −1 is applied in the tidal river and in the saline estuary, respectively, while w s is approximated to 1 mm s −1 (Winterwerp, 2002).On the other hand, sediment parameters, such as the critical shear stress for erosion and deposition (τ cr in N m −2 ) and the erosion coefficient (E ero in kg m −2 s −1 ), are usually calibrated on the basis of local SPM observations and represent system-specific fitting parameters with a limited transferability (Volta et al., 2014).Here, we adapt values reported for an idealized, tidal alluvial estuary by Volta et al. (2014).

Biogeochemical parameters
Biogeochemical parameters and their corresponding numerical values used in the C-GEM simulations are listed in Table 5.In natural waters, the Redfield ratio C : N : Si : P, which is instrumental for estimating carbon and nutrient production/consumption rates (Table 3), can be approximated by the Redfield-Brzezinski ratio 106 : 16 : 15 : 1 (Redfield at al., 1963;Brzezinski, 1985).The background light extinction coefficient (K D1 in m −1 ) and the specific light attenuation of SPM (K D2 in mg −1 m −1 ) are system-specific attributes that are generally derived from local underwater light and SPM observations (Volta et al., 2014) and they are adapted here from values reported for an idealized, tidal alluvial estuary by Volta et al. (2014).All other biogeochemical parameters (n =18, identified in bold in Table 5) are derived from a comprehensive literature review of all published estuarine biogeochemical model applications in temperate regions over the last 30 years.In this study, we define temperate regions as those lying between 30 and 60 • in either hemisphere.As a consequence, the generic biogeochemical parameterization provided and applied in this study should be considered representative of these zones.A more detailed description of the biogeochemical parameter review and analysis is provided in Sect.3.

Climate forcings and boundary conditions
Except for the freshwater discharge, all model simulations are forced with the same set of climate and hydrological forcings, representative of the mean annual conditions in temperate regions during the 2000s (Table 6).Annually averaged values of irradiance and photoperiod are calculated by using the astronomical equation of Brock (1981), while Table 3. Formulations of the biogeochemical and sediment processes with the corresponding stoichiometric equations as implemented in the current C-GEM reaction network.T abs and T denote the absolute and the Celsius temperature, respectively.H is the instantaneous water depth.PHY is the phytoplankton concentration.If PHY = DIA, nlim also accounts for the silica limitation for the phytoplankton growth.Further details can be found in Arndt et al. (2011b) and in Volta et al. (2014).Temperature functions are deduced from a Garcia et al. (2010);  b Cerco (2000), Kim and Cerco (2003), Cerco and Noel (2004); c Chapelle et al. (1994Chapelle et al. ( , 2000)); d calibration; e Solidoro et al. (2005); and f Robson and Hamilton (2004), Zheng et al. (2004).Further details about the temperature functions are given in Sect.3.2.

Gross primary production GPP
Nutrient limitation for phytoplankton growth nlim Switch between NH 4 and NO 3 utilization Light extinction coefficient Wind component for piston velocity

660
Sediment erosion Critical shear stress for erosion and deposition Critical shear stress for erosion and deposition in the saline estuary Critical shear stress for erosion and deposition in the tidal river Erosion coefficient in the tidal river [kg m −2 s −1 ] 6.0 × 10 −8 * 6.0 × 10 −8 * 6.0 × 10 −8 * wind speed and water temperature are extracted from data reported for the coastal temperate zones by the CCMP data set (Atlas et al., 2011) and the World Ocean Atlas global database (http://www.nodc.noaa.gov/OC5/indprod.html),respectively.A constant tidal amplitude (ζ 0 = 3.5 m) is applied at the mouth of all estuaries, whereas a different systemspecific freshwater discharge is imposed at their upper limits (see Sect. 2.2 and Table 1).Model simulations are forced with two different sets of biogeochemical boundary condition (Table 7).The first set, referred to as the baseline set, is representative of present-day conditions, while the second set represents a scenario for the year 2050.The riverine inputs of organic carbon, nutrients and suspended particulate matter of the baseline set are derived from the global statistical model GlobalNEWS2 (Mayorga et al., 2010).Values represent the average calculated over all watersheds in temperate regions that discharge to the sea through a tidal estuary.For these calculations, we use the estuarine coastal typology of Dürr et al. (2011), which identifies four types of estuaries: small deltas, tidal systems, lagoons and fjords.In this typology, tidal systems (type 2) represent a good approximation of the domain of applicability of C-GEM (Regnier et al., 2013b).NO 3 and NH 4 concentrations are derived by applying a NH 4 / DIN ratio equal to 0.2 to DIN concentrations provided by GlobalNEWS2.The latter is calculated as the median of the NH 4 / DIN ratios reported by Meybeck (1982) in more than 40 rivers.Alkalinity is derived from the GLORICH database (Hartmann et al., 2014) and represents the average value for all watersheds in temperate regions.A constant CO 2 concentration typical of temperate rivers worldwide (≈ 90 µM C, from 3-year time series in more than 1000 sampling locations; Lauerwald et al., 2015) is then used to calculate dissolved inorganic carbon (DIC) concentration at the upstream limit.Because of the lack of relevant global database and in order to minimize the influence of boundary conditions on estuarine phytoplankton and oxygen profiles, arbitrary low riverine concentrations are imposed for both phytoplanktonic groups (DIA and nDIA), while the saturation concentration is imposed for O 2 .The downstream boundary is located 50 km away from the estuarine mouth in order to minimize its influence on the estuarine biogeochemical dynamics.Here, salinity, organic carbon, nutrient and oxygen concentrations are extracted from data reported for the temperate regions by the World Ocean Atlas global database (http://www.nodc.noaa.gov/OC5/indprod.html),while phytoplankton concen-  (1963) and Brzezinski (1985).K D1 and K D2 are from Volta et al. (2014).All rates are defined at 20 • C.

Name Description [unit]
Value Phytoplankton maintenance rate constant Phytoplankton mortality rate constant  trations are deduced from SeaWIFS data (http://oceancolor.gsfc.nasa.gov).TAlk and DIC are calculated by assuming a typical seawater pH of 8.2 (Mackenzie et al., 2011) and an average difference between atmospheric and shelf seawater pCO 2 ( pCO 2 ) of 20 µatm for the year 2000 (Cai, 2011) consistent with a CO 2 sink for the coastal ocean under present-day conditions (Bauer et al., 2013;Laruelle et al., 2014).No database is available to constrain average total organic carbon and suspended particulate matter concentrations at the lower boundary and both concentrations are arbitrarily set to 0. In the case of SPM, the implementation in C-GEM of a sediment transport module allows for enough internal production of SPM to generate realistic profiles along the entire estuarine length.In addition, the low concentration reflects the low values typically observed on the shelf compared to those observed in shallow nearshore areas (e.g., Ruddick et al., 2003).On the other hand, although a variable load of organic carbon can be imported from the adjacent coast into the estuary (e.g., Arndt et al., 2011a) and assuming a TOC concentration of 0 at the lower boundary may thus be an approximation, our choice allows focusing on the fate of organic matter brought the estuary from rivers by minimizing the effect of organic carbon produced and imported from the sea into estuaries.The second biogeochemical boundary condition set represents a future scenario  Jr. et al., 2003;Raymond et al., 2008).As a consequence, no future evolution can be confidently attributed to riverine DIC and TAlk at the moment.Because of the increase in atmospheric CO 2 concentrations, ocean DIC concentrations are expected to increase, while ocean alkalinity will stay constant over the next decades due to the buffering capacity of the ocean (e.g., Andersson et al., 2005;Mackenzie et al., 2011;IPCC Report, 2013).However, the exact magnitude of these variations remains unconstrained and the CO 2 uptake rate by the coastal ocean may double or even triple by year 2050 (Andersson and Mackenzie, 2004;IPCC Report, 2013).Here, future marine DIC concentrations are calculated using MatLab csys.m (Zeebe and Wolf-Gladrow, 2001) by assuming that TAlk does not change in the near future and an average atmospheric pCO 2 of 468 µatm, corresponding to the value predicted by the IPCC RCP6 scenario for the year 2050 (IPCC Report, 2013), while a seawater pCO 2 of 408 µatm is imposed in order to test the influence of an hypothetical future 3-fold increase of the atmospheric and marine water pCO 2 (i.e., 60 µatm) with respect to the year 2000.

Biogeochemical indicators
The biogeochemical dynamics of the three representative estuarine systems are investigated and compared by means of four whole-system biogeochemical indicators: the net ecosystem metabolism (NEM), the CO 2 exchange flux with the atmosphere (FCO 2 ) and the total nitrogen and carbon filtering capacities (FC TN and FC TC , respectively).

Net ecosystem metabolism (NEM)
The NEM is defined as the whole-system difference between net primary production (NPP) and heterotrophic degradation (aerobic degradation + denitrification) (Andersson and Mackenzie, 2004).The NEM is, thus, controlled by the input, export, production and decomposition of terrestrial and in situ-produced organic matter (Odum, 1956) and is typically used to assess the trophic status of an estuary (e.g., Caffrey, 2003;Borges and Abril, 2011).Heterotrophic systems are characterized by the dominance of organic matter degradation over production (NEM < 0), which leads to a net regeneration and export of inorganic carbon and nutrients.In contrast, autotrophic estuaries, dominated by photosynthetic production (NEM > 0), are characterized by a net burial and export of organic matter.As a consequence, the estuarine NEM can be used not only in defining the trophic status of system but also as an indicator of carbon and nutrient sources and sinks in an estuary.

CO 2 exchange flux (FCO 2 )
The balance between autotrophic and heterotrophic processes also controls to a large extent the estuarine inorganic carbon dynamics and, thus, the CO 2 exchange across the water-atmosphere interface (FCO 2 ).In general, FCO 2 depends on the overall effect of biogeochemical processes on dissolved inorganic carbon (DIC) and total alkalinity (TAlk).
For instance, the aerobic degradation of organic matter generally releases large amounts of DIC, decreases pH and thus promotes CO 2 outgassing, whereas NPP increases water pH and limits CO 2 evasion.Denitrification, on the other hand, increases both DIC and TAlk and thus exerts a limited influence on the inorganic carbon budget, while nitrification decreases pH and generally sustains CO 2 outgassing.FCO 2 is, thus, an integrative measure of all biogeochemical processes that exert an influence on the carbonate systems in estuaries (Regnier et al., 2013b).

Nitrogen and carbon filtering capacities (FC TN and FC TC )
Total nitrogen and total carbon filtering capacities (FC TN and FC TC , respectively) provide a measure for how much nitrogen or carbon is lost during the estuarine transit.By considering denitrification as the (negative) net process rate affecting estuarine nitrogen, FC TN is calculated as the ratio between denitrification and total riverine nitrogen (TN) flux, defined as the sum of the dissolved inorganic nitrogen and the organic nitrogen bound to living and detrital organic matter (Arndt et al., 2009).Similarly, FC TC is defined as the ratio between carbon loss through the water-atmosphere interface and the total riverine carbon (TC) influx, which accounts for both inorganic and organic carbon (Regnier et al., 2013b).Carbon and nitrogen cycles are typically strongly altered by biogeochemical processes within the estuarine bioreactor and estuarine removal efficiency may thus have relevant implications for the coastal biogeochemistry and for processes producing and releasing greenhouse gases (e.g., N 2 O, CO 2 ) (e.g., Seitzinger, 1988;Soetaert and Hermann, 1995;Voss et al., 2011;Bauer et al., 2013;Regnier et al., 2013a).An analysis of FC TN and FC TC may thus advance the understanding of the role of estuaries in the carbon and nitrogen biogeochemical cycling.

Sensitivity study
A sensitivity analysis is carried out to assess the response of the biogeochemical indicators (FC TN , FC TC , FCO 2 and NEM) in the three idealized estuaries to variations in aerobic degradation, denitrification and nitrification rate constants (k ox , k denit and k nit , respectively).These parameters are selected because of the heterotrophic character of estuaries (Borges and Abril, 2011;Volta et al., 2014) and the significant effect of nitrification rates on water pH and, thus, on DIC speciation, water pCO 2 and FCO 2 (Regnier et al., 2013b).The generic rate constants (Table 5) are used as baseline values and are exponentially increased and decreased over 1 order of magnitude.That range of variation corresponds to the range of k ox values reported in the literature and minimum and maximum values may be regarded as representative of refractory and labile organic carbon loads, respectively (see Sect. 3).Similarly, the outermost k nit and k denit values may be considered as representative of different ammonia-oxidizing and nitrate-reducing microbial communities, respectively.For each estuary, 10 parameter combinations are tested.An additional sensitivity test is also performed in which, k ox and k denit vary over 2 orders of magnitude as in the first test (SA1), but k nit remains constant at its baseline value (Table 5).This second sensitivity analysis (SA2) allows assessing the relative importance of the heterotrophic and nitrification reactions on the estuarine biogeochemical indicators in the three idealized estuaries by comparing results from both series of tests.All the parameter values used in SA1 and SA2 are provided in Table A1.
3 Biogeochemical parameter review and analysis

Literature review
Values for 18 biogeochemical parameters included in the biogeochemical reaction network of the C-GEM modeling platform (Table 3) were compiled by a literature review comprising 49 model applications for tidal estuarine systems in temperate regions (Table 5).The review comprises models of different complexity ranging from 0-to 3-D models, which were developed and applied to investigate different aspect of estuarine biogeochemistry, such as water quality control (e.g., HydroQual Inc., 1987;Lin et al., 2007), bacterial and phytoplankton dynamics (e.g., Robson and Hamilton, 2004;Macedo and Duarte, 2006;Gypens et al., 2013), or estuarine nutrient budgets (e.g., Soetaert and Herman, 1995;Arndt et al., 2009) at different timescales ranging from months to several years.The review covers the following 12 parameters, assumed to be temperature-independent: -Michaelis-Menten half-saturation constants (K DSi , , and K N in µM), which account for the dependency of biogeochemical reaction rates on substrate availability; -the inhibition constant for denitrification (K in,O 2 in µM O 2 ), describing the inhibition of denitrification reaction by oxygen; -the photosynthetic efficiency (α in m 2 s µE −1 s −1 ), which represents the light harvesting efficiency of phytoplankton; -the phytoplankton excretion (k excr ) and growth constants (k growth ), accounting for the fraction (in %) of the C. Volta et al.: Linking biogeochemistry to hydro-geometrical variability in tidal estuaries gross production lost through exudation processes and the maximum growth of phytoplankton, respectively.
In addition, 6 temperature-dependent parameters (as well as their associated temperature functions) are covered: -the maximum specific photosynthetic rate (P B max in s −1 ), which corresponds to the light-saturated, carbon uptake rate by primary production; -the phytoplanktonic maintenance (k maint in s −1 ) and mortality (k mort in s −1 ) rate constants, representing the loss of biomass due to maintenance activity and phytoplankton mortality, respectively; -the rate constants for aerobic degradation (k ox in µM C s −1 ) and denitrification (k denit in µM C s −1 ) describing the reactivity of organic matter towards heterotrophic decay and the nitrification rate constant (k nit in µM N s −1 ), which defines the reactivity of ammonia towards biologically mediated oxidation.
In order to include as many studies as possible, a unit homogenization was performed for the parameters involved in biogeochemical reactions whose mathematical formulations differ from those implemented in C-GEM (Table 3).Firstorder reaction rate constants of aerobic degradation, denitrification and/or nitrification, expressed in s −1 (e.g., Soetaert and Herman, 1995;Robson and Hamilton, 2004;Hofmann et al., 2008), were converted to zero-order constants in µM s −1 by multiplying the first-order reaction constant by the typical watershed concentration of the substrate involved in the specific reaction (i.e., TOC for aerobic degradation and denitrification and NH 4 for nitrification).Both TOC and NH 4 concentrations were extracted from the global statistical model GlobalNEWS2 (Mayorga et al., 2010).Watershed-specific NH 4 concentration are estimated on the basis of reported DIN concentrations by assuming a NH 4 / DIN ratio of 0.2, representing the median of the NH 4 / DIN ratios reported by Meybeck (1982) in more than 40 uncontaminated and polluted rivers worldwide.Similarly, maximum specific photosynthetic rate constants (P B max ), expressed in gC g Chl −1 s −1 (e.g., Cerco, 2000;Kim and Cerco, 2003;Desmit et al., 2005), are converted to first-order carbon-production-based P B max expressed in s −1 , by dividing P B max in gC g Chl −1 s −1 by the associated carbon-to-chlorophyll ratio in gC g Chl −1 .Biogeochemical parameter values reported in the literature are summarized in Fig. 5.The latter shows that, with the exception of P B max , for which very large variability (> 35 orders of magnitude) is found, biogeochemical parameter values typically span over a maximum of 5 orders of magnitude.The Michaelis-Menten constant K O 2 ,ox and the inhibition constant for denitrification (K in,O 2 ), as well as the phytoplankton parameters α, k excr and k growth , vary over 1 order of magnitude.On the other hand, a variability over 2 orders of magnitude is found for the five Michaelis-Menten terms K NH 4 , K NO 3 , K TOC , K O 2 ,nit , K N , as well as for the phytoplankton parameter k maint and the rate constant for aerobic degradation (k ox ).Finally, the Michaelis-Menten terms K DSi , the mortality rate for phytoplankton, k mort , and the nitrification constant rate (k nit ) display a relatively large variability over 3 orders of magnitude, while the Michaelis-Menten parameter K PO 4 and the denitrification constant rate (k denit ) vary over 4 and 5 orders of magnitude, respectively.The large variability range observed for phytoplankton parameters, such as P B max , k maint , k mort and K PO 4 , is likely related to the fact that these parameters typically represent different phytoplanktonic species or groups, with varying traits.However, in the reviewed modeling applications, phytoplanktonic parameter values vary not only from a phytoplankton species to another but also within the same group.For example, different k maint values are reported for the same phytoplankton group (i.e., diatoms) in Soetaert et al. (1994), Garnier et al. (1995) and Kim and Cerco (2003).On the other hand, some modeling applications used the same parameter value for different phytoplanktonic groups (e.g., same mortality rate constant for diatoms, flagellates and Phaeocystis in Blauw et al., 2009).The literature review also reveals that biogeochemical parameter values do not strongly vary from one estuarine system to another, but different values are also used in modeling studies of the same estuary (e.g., four different k ox values used in modeling applications to the Scheldt estuary by Soetaert and Herman, 1995;Regnier et al., 1997;Regnier and Steefel, 1999;Hofmann et al., 2008;Arndt et al., 2009;Volta et al., 2014).
For the temperature-dependent parameters, corresponding temperature functions are also included in the review.The review reveals that the temperature dependence of the aerobic and the denitrification rate constants (k ox and k denit , respectively) is typically expressed as an exponential increase of the rate constants with temperature (e.g., Regnier et al., 1997;Hofmann et al., 2008;Arndt et al., 2009;Volta et al., 2014).On the other hand, the temperature dependence of autotrophic parameters, such as the nitrification rate constant (k nit ), the maximum specific photosynthetic rate (P B max ), and the phytoplankton maintenance and mortality rate constants (k maint and k mort , respectively) can be implemented as exponential functions, where the parameter value increases with temperature (e.g., Peterson and Festa, 1984;Le Pape and Ménesguen, 1997;Guillaud et al., 2000;Laruelle et al., 2009), or as Gaussian functions, where the value increases until an optimum temperature is reached and then progressively decreases (e.g., Garnier et al., 1995;Kim and Cerco, 2003;Gypens et al., 2013;Zheng et al., 2004).
A summary of the reviewed biogeochemical parameters, as well as their values, the location of the respective modeling studies and the corresponding reference, is provided in Tables B1 and B2.For the temperature-dependent parameters (Table B2), temperature functions are also reported.

Parameter analysis
A statistical analysis was performed to determine a generic biogeochemical parameter set representing the mean or the median of the respective published parameter values.Generally, the mean is considered a rigorous estimate of the central tendency of a set of normally distributed numerical scores.However, it is not a well-suited measure for data sets with a skewed distribution since it is largely influenced by outlier values.In this case, a numerical measure able to minimize the outlier influence on the generic trend, such as the median, should be preferred (Mendenhall et al., 2013).
Here, the presence of outliers is used to determine whether a biogeochemical parameter may be approximated with the arithmetical mean or the median of a data set.If no outliers are detected, the literature parameter set can be considered normally distributed and thus the generic biogeochemical parameter can be approximated by the mean.Otherwise, when at least one outlier is identified, the distribution of values in the parameter set is assumed to be skewed and the generic biogeochemical parameter is calculated as the median of the literature values (Fig. 6).Parameter values are classified as outliers if they are larger than q 3 + w • (q 3 − q 1 ) or smaller than q 1 − w • (q 3 − q 1 ), where q 1 and q 3 are the 25th and 75th percentiles, respectively, while w is the maximum whisker length.The latter is set equal to 1.5, which corresponds to approximately 99 % coverage if the data are normally distributed and represents a rational compromise  1 and 4 and baseline boundary conditions reported in Table 7.
between a rigorous parameter value selection, which aims to identify a generic tendency of the literature parameter distributions without being influenced by values too high or too low compared to the rest of the sample, and the conservation of a statistically relevant number of parameter values.
For the sake of generalization, although in some modeling applications phytoplankton parameters were associated with specific groups, no distinctions between different phytoplankton species and/or groups were made during the analysis and any numerical value was assigned the same weight in the parameter estimation.Therefore, the generic set of phytoplankton parameters derived here should represent a generic estuarine phytoplankton group.Moreover, when a modeling study reported different aerobic degradation and denitrification rate constants (k ox and k denit , respectively) for differently reactive organic matter pools (e.g., Soetaert and Herman, 1995;Schroeder, 1997;Hofmann et al., 2008), the value considered in the parameter analysis corresponds to the average value of the two constants reported.Therefore, the generic values for k ox and k denit may be regarded as representative of organic matter decomposition through high-energy-yielding metabolic pathways.
A generic set of temperature functions, associated with temperature-dependent biogeochemical parameters, is derived by analyzing the T functions reported by the reviewed modeling applications.For phytoplankton temperaturedependent parameters, generic temperature functions were preferentially chosen over group-specific temperature formulations.For the phytoplankton maintenance rate constants (k maint ), all temperature functions reported in the literature referred to specific phytoplanktonic groups.In this case, the function associated with the largest number of different phytoplankton groups was considered the most generic and, thus, retained.On the other hand, although all tempera-ture functions reported in the literature for aerobic degradation and denitrification rate constants are exponential, they nevertheless show different trends depending on the initial values and exponential factors applied.In this case, the formulation showing the average trend is selected as generic function (refer to Fig. B1).Finally, for the temperature function associated with the rate constant of nitrification (k nit ), the most frequently used function is chosen (refer to Table B2 for more details).For the generic temperature functions of biogeochemical parameters as implemented in C-GEM, refer to Table 3.

Hydrodynamics and salt transport
In an estuarine system, the interplay between tidal and fluvial influence results in important longitudinal variations in tidal amplitude and salinity (Fig. 7a and b).The comparison of the simulated longitudinal tidal amplitude and the salinity profiles for the three representative estuaries (Table 1) reflects the strong mutual dependency between estuarine geometry and hydrodynamic characteristics.
The distortion of the tidal wave in upstream direction (Fig. 7a) is mainly controlled by the balance between energy gain due to channel convergence and energy loss through friction.In the marine estuary, the strong convergence of the estuarine banks (short convergence length) compensates for the energy loss through friction and the tidal amplitude increases landward to 5.5 m at the upper estuarine boundary.In the mixed estuary, on the other hand, the weak amplification of the tidal amplitude indicates that the energy gain through convergence is almost balanced by the energy lost through friction.Upstream, tidal amplitude increases more rapidly to reach a maximum of about 5 m where convergence and friction effect are both of low magnitude.Beyond this point, close to the inland limit, the relative importance of the friction increases, owing to a relatively higher fluvial energy, and triggers a slight dampening of the tidal amplitude.In the riverine estuary, the weak channel convergence (long convergence length) combined with the dominant fluvial influence leads to a dampening of the tidal amplitude, in particular in the upper reach where frictional energy loss reaches a maximum.In summary, the marine and the mixed estuaries are typically characterized by a net energy gain through convergence and thus by a tidal amplitude amplification along their longitudinal gradients, whereas the large fluvial influence always induces a net reduction of estuarine energy by bottom friction and a consequent tidal damping in the riverine system.
In alluvial estuaries, the salinity intrusion is controlled by the balance between upstream dispersion and downstream advective transport and, thus, by the system's geometric and hydrodynamic properties (Savenije, 2005(Savenije, , 2012)).Simulated salinity profiles for the three idealized systems (Fig. 7b) reflect this dependency.Despite identical lower boundary conditions (S = 34 at 50 km beyond the estuarine mouth; see Sect.2.3.5), the shape of the simulated salinity profile and, in particular, the salinity gradient close to the estuarine mouth, as well as the salt intrusion length, reveal the characteristic differences generally observed across the different estuarine types (Savenije, 1992(Savenije, , 2005(Savenije, , 2012)).These differences can be largely attributed to the relative significance of the tidal versus the fluvial influence in each system.The marine estuary is characterized by a dominant tidal influence and, thus, a small salinity gradient close to the mouth ( S = 7), a concave salinity profile and a long salinity intrusion length of almost the 75 % of the total estuarine length (EL).On the other hand, in the mixed estuary, tidal and fluvial influences are of roughly equal importance and the salinity gradient close to the mouth is thus larger than in the marine estuary ( S = 17).Moreover, the larger fluvial influence, results in a recession shape profile and a shorter salinity intrusion length of 40 % of the EL.Similarly, the riverine estuary displays a recessionshape salinity profile.However, the salinity gradient at the estuarine mouth is much larger ( S = 24) and the salinity intrusion shorter (20 % of the EL) than in the mixed system due to the large riverine discharge and the smaller tidal exchange.

Solid transport
In alluvial estuaries, the main features of the longitudinal distribution of suspended particulate matter (SPM) can be linked to the mechanical energy provided by the tides and the riverine discharge (Jay et al., 1990;Dalrymple et al., 1992).The longitudinal SPM profiles simulated in the marine and the mixed estuary reveal similar SPM trends (Fig. 7c).In both cases, SPM concentrations increase in the lower estuary due to the progressive compression of the incoming flood into a smaller cross-sectional area and to a consequent increase in the flood-tidal current speed and reach locally high values where the total energy (fluvial + tidal) is maximal, while, further upstream, SPM concentrations decrease to a minimum at the so-called balance point, where fluvial and tidal energy contributions are of similar but low magnitude.In the upper reaches, SPM concentrations are largely controlled by the riverine influence.The progressive decrease in fluvial energy from the upper estuarine limit to the sea induces a reduction of erosion and a consequent increase in deposition rates.As a result, decreasing SPM from the upper limit to the energy balance point is simulated.On the other hand, the riverine estuary reveals two turbidity maxima separated by a zone of low SPM concentrations, corresponding to the energetic balance point (Fig. 7c).As for the marine and the mixed cases, the progressive increase in SPM from the downstream limit upwards is essentially related to the progressive compression of the marine inflow into a progressively smaller volume.Nonetheless, in the upper reaches, maximum SPM concen- trations are likely induced by the high river discharges promoting net erosion, as well as by the strong volumetric reduction associated with the large tidal wave dampening (Fig. 7a).
A quantitative comparison of the SPM concentrations simulated in the three representative estuaries (Fig. 7c) reveals that the marine system shows very low SPM levels.The latter are likely related to the strong funnel character of this system that, bearing relatively steeper and larger volumetric variations along the estuarine gradient (Fig. 3), entails a larger dilution effect on SPM concentrations.Overall, modeled SPM distributions agree well with the conceptual, mechanical energy patterns described by Dalrymple et al. (1992), who identified two high-energy zones separated by a low-energy area (the balance point) in tidal estuaries.Moreover, they indicate that the extent of these zones along the three idealized estuaries, as well as their sediment loads, can be strongly affected by different hydro-geometrical characteristics.

Biogeochemical dynamics -baseline simulations
Figure 8 shows the longitudinal distribution of ammonium (NH 4 ), nitrate (NO 3 ), phosphate (PO 4 ), dissolved silica (DSi), oxygen (O 2 ), total organic carbon (TOC), diatoms (DIA), non-diatom phytoplankton (nDIA), dissolved inorganic carbon (DIC), total alkalinity (TAlk), pH and water pCO 2 in the three idealized estuaries simulated with the set of baseline boundary conditions (Table 7).Simulated concentration and biogeochemical rate profiles reveal several characteristic features that are common across the three estuaries.In addition, although a direct quantitative comparisons with observations is difficult, the simulated concentration and rate profiles qualitatively agree with observed chemical distributions from different temperate tidal systems, such as, for instance, the Chesapeake Bay (e.g., Horrigan et al., 1990) and the Delaware (e.g., Sharp et al., 2009), the Scheldt (e.g., Baeyens et al., 1998), the Thames and the Gironde (e.g., Frankignoulle et al., 1998), the Severn (e.g., Jonas and Millward, 2010) and the Tweed (Howland et al., 2000) estuaries.TOC and nutrient (NH 4 , NO 3 , PO 4 and DSi) concentrations show an overall decrease in downstream direction (Fig. 8a) due to dilution and high reaction rates in the upper reaches of the estuarine systems (Fig. 9).However, the estuarine profiles of NH 4 and PO 4 reveal a mid-estuary maximum (Fig. 8a), which results from the slow degradation of TOC in this area (Fig. 8a).Although slightly positive net primary production rates are simulated in the upper estuaries, average annual conditions (e.g., temperature, photoperiod and solar irradiation) prevent the occurrence of phytoplanktonic blooms and phytoplankton concentrations (DIA and nDIA, Fig. 8a) progressively decrease downstream due to dilution and phytoplankton mortality effects.The lower reaches are consistently characterized by high O 2 , DIC and TAlk concentrations and a high pH, which decrease in up-Table 8. Biogeochemical indicators calculated for the three idealized estuaries by using the baseline and the future boundary conditions (columns a and b, respectively).NEM and FCO 2 are in kmol C day −1 , while FC TN and FC TC are expressed as percentage of the total riverine input.FCO 2 is negative when towards the atmosphere.Integrated rates for NEM's constitutive reactions are also provided (R: aerobic degradation in kmol C day −1 ; D: denitrification in kmol C day −1 ; NPP: net primary production in kmol C day −1 ).Relative variations in biogeochemical indicators (and NEM's constitutive reactions) compared to their baseline values, are also reported in column b (values inside brackets).MAR, MIX and RIV correspond to the marine, the mixed and the riverine estuary, respectively.stream direction.Aerobic degradation is by far the dominant pathway of organic carbon degradation in all estuaries (Fig. 9a) and, together with the O 2 transfer across the air-water interface, it is the dominant control on the O 2 longitudinal distribution (Fig. 9b).TAlk profiles are essentially determined by total heterotrophic degradation (aerobic degradation + denitrification) and nitrification (Fig. 9d), while DIC concentrations largely depend on a dynamic balance between production via aerobic degradation and loss through CO 2 outgassing (Fig. 9e).These two processes are also the main drivers of the pH changes along the longitudinal axis of all estuaries (Fig. 9f).In the upper reaches, maximum pCO 2 concentrations (Fig. 8b) correspond to the areas where minimum pH values are simulated.
Although the three idealized estuarine types reveal similar biogeochemical dynamics, they also show distinct features.In particular, the increasing significance of freshwa-  7.
ter discharge from the marine to the riverine estuary results in a stronger residual transport in downstream direction, which systematically pushes the biogeochemical fronts towards the estuarine mouth.Large quantitative differences between the three idealized estuaries are also evident when comparing the values of their integrated biogeochemical indicators (Table 8, column a).Although results show that all estuaries are net heterotrophic (NEM < 0) and act as net sources for atmospheric CO 2 (FCO 2 < 0), the system-scale NEM and FCO 2 increase from the marine (−916 and −2018 kmol C day −1 , respectively) to the mixed (−8161 and −10 940 kmol C day −1 ) and to the riverine system (−21 476 and −25 612 kmol C day −1 ) due to the strong correlation between these indicators and the riverine influx (Regnier et al., 2013b; refer to Tables 1 and 7).Yet, the nitrogen and carbon filtering capacities (FC TN and FC TC , respectively) decrease from the marine to the riverine type (Table 8, column a), reflecting the progressively shorter transit times of the water masses as the discharge increases (e.g., Nixon et al., 1996;Arndt et al., 2009Arndt et al., , 2011a;;Regnier et al., 2013b).The simulated integrative measures of NEM and FCO 2 for the three idealized estuaries, ranging between −2 and −39 mol C m −2 yr −1 and −4 and −47 mol C m −2 yr −1 , respectively, fall within the range of values observed in tidal estuaries in temperate regions (−63 to +25.1 and −76 to +6 mol C m −2 yr −1 for NEM and FCO 2 , respectively; refer to Borges andAbril, 2011, andLaruelle et al., 2013, for a complete list of values).Furthermore, N filtering capacities of the three systems (15-22 %) are comparable to typical values reported in the literature for temperate tidal estuaries, such as the Seine (< 7-40 %; Garnier et al., 2001), the Scheldt (13-78 %; Arndt et al., 2009) and the James (24-40 %; Bukaveckas and Isenberg, 2013) estuary, while their Cremoval efficiencies (22-40 %) are comparable to the range calculated in Regnier et al. (2013b) for three idealized, western European tidal estuaries (25-31 %).O 2 exchange with the atmosphere; M: phytoplankton mortality; FCO 2 : CO 2 exchange with the atmosphere.Negative gas exchange rates indicate gas evasion from water towards the atmosphere, while positive rates represent gas fluxes from atmosphere to water column.Process contributions are calculated as in Regnier et al. (1997Regnier et al. ( , 2013b)).

Sensitivity analysis
Two sensitivity tests (refer to Sect.2.5) are performed to quantify the responses of the biogeochemical indicators (NEM, FCO 2 , FC TN and FC TC ) to different combinations of aerobic degradation, denitrification and nitrification rate constants (k ox , k denit and k nit , respectively) in the three idealized estuaries (Fig. 10).The first sensitivity test (SA1) analyzes the response of the biogeochemical functioning of the three systems to changes in reaction kinetics or, in other words, the characteristic timescales of reactions.In this case, it is  A1.All simulations are forced with the baseline biogeochemical scenario for the year 2000 (Table 7).assumed that the ammonification, induced by the degradation of organic matter, and the nitrification are tightly coupled (e.g., Billen et al., 1985).These results are then compared to those obtained with variable k ox and k denit but constant k nit (SA2) in order to assess the relative importance of heterotrophic degradation and nitrification on the biogeochemical indicators in the three estuarine systems.
In all three idealized estuaries, the NEM becomes obviously more negative as k ox , k denit and k nit increase (Fig. 10a).However, results indicate that this trend is much more pronounced in the riverine than in the marine estuary.Furthermore, despite large variations in parameter values, none of the estuaries becomes autotrophic.Because the strong correlation between the NEM and the system-specific riverine organic carbon influx does not allow for a direct quantitative comparison of the NEM variability across the three systems, the integrated NEM value is normalized with respect  c) denitrification (D, in mmol C m −1 day −1 ) rates in the three idealized estuaries for each of the simulations of the first sensitivity analysis (SA1; see Table A1).S1 to S10 refer to different simulations, each corresponding to a specific parameter combination while BS indicates baseline parameter values (Table 5).Black vertical lines indicate the position of the maximal rates simulated by using the BS biogeochemical parameter sets, while arrows displays the longitudinal shift of the biogeochemical fronts by varying reaction rate constants.
to the total organic carbon input flux of each estuary as a relative measure of the amount of organic carbon that is processed within the estuary (FC org in Fig. 10b).In the marine estuary, FC org is high and weakly sensitive to variations in rate constants.On the other hand, in the mixed and the riverine estuary, FC org is high and fairly constant when parameter values are larger than BS values, but rapidly decreases when reducing k ox , k denit and k nit below their BS values.Above the BS threshold, a rapid degradation of organic matter results in a significant consumption of the organic carbon flux (about 90 % of the riverine flux) irrespective of the estuarine type.These results emphasize the delicate balance between reaction timescales and residence time in determining filtering capacities.In general, increasing the values of the rate constants above BS values shifts the heterotrophic degradation zone further upstream in all estuaries (Fig. 11a), where residence times are, due to smaller volumes, shorter.However, although the shape and the peak value of the total heterotrophic longitudinal profile vary, the integrated degradation rates and thus the FC org are similar across estuaries (Fig. 10a and b).Therefore, results reveal that differences in geometrical characteristics do not exert a strong influence on organic carbon dynamics.On the other hand, a reduction of kinetic rate constants below BS values triggers a downstream shift of the heterotrophic zone (Fig. 11a).Here, larger volumes and thus longer residence times may partly compensate for the low reaction rates.This is particularly true for the marine estuary, where the long residence time strongly compensates for the reduced rate constants and translates in high and fairly constant FC org , while the continuous decrease in FC org values simulated in the mixed and the riverine estuary suggests that the increase in residence time is not sufficient to compensate for low reaction rates.
Figure 10c and d summarizes the FCO 2 and the FC TC simulated in the three idealized estuaries for the set of rate constants SA1.Since the C filtering capacity depends on the estuarine efficiency in scavenging carbon through CO 2 degassing, FC TC 's and FCO 2 's sensitivity results are discussed together.As for the NEM, the strong correlation between FCO 2 and system-specific riverine carbon inputs limits the quantitative comparison of the FCO 2 variability across the three systems.However, a qualitative analysis reveals that CO 2 outgassing (FCO 2 < 0) increases with increasing parameter values in all systems.As a consequence, enhanced C-removal efficiencies, which increase from the riverine to the marine estuary, are also simulated.Closer inspection of the different responses suggests that decreasing k ox , k denit and k nit below BS values only induces small FC TC variations in the marine estuary, whereas a larger sensitivity is simulated in the mixed and the riverine estuaries.Such behavior is very similar to that of the NEM and further confirms that system-specific differences in residence time play an important role in controlling the estuarine biogeochemical functioning when reaction rates are low and biogeochemical reaction zones are pushed further downstream (Fig. 11b).On the other hand, increasing k ox , k denit and k nit above BS values results in a larger sensitivity in the marine estuary, where FC TC can reach values up to 90 %, while the mixed and the riverine systems only reveal a weak sensitivity and fairly constant FC TC values (about 40 %).These differences reflect the direct dependence of FCO 2 on gas exchange and thus the estuarine surface area.As shown in Fig. 11b, increasing rate constants above BS values induces an upstream shift of the zone where maximum CO 2 degassing occurs.Here, smaller surface areas may, however, partly limit the gas exchange.This is particularly true for the mixed and the riverine systems, where long convergence length and thus comparably smaller surface areas limit the potential positive effect of higher reaction rate constants on CO 2 outgassing.In contrast, the large increase in FC TC simulated in the marine estuary by increasing k ox , k denit and k nit above BS values suggests that its strong funnel-shaped character, entailing relatively larger surface, always allows for enhanced C-removal efficiencies even if the FCO 2 front is located in the upper zone.
The FC TN sensitivities in the three idealized estuaries to changes in parameter combinations are presented in Fig. 10e.The close similarity between the results obtained for FC TN and FC TC suggests that N removal is also controlled by residence time when the rate constants decrease below BS values and the denitrification zone is pushed further seawards (Fig. 11c).However, unlike FC TC , the progressive FC TN increase simulated in the three idealized estuaries by increasing the rate constants indicates that denitrification and, thus, FC TN could theoretically be enhanced by a further increase in rate constants in all systems.
The comparison of the results from the two sensitivity tests (SA1 and SA2 in Fig. 10) shows that nitrification rate constants exert small effects on all biogeochemical indicators in the three idealized estuaries.Slightly larger variations are simulated for FCO 2 (and consequently for FC TC ), as well as for FC TN , by varying the k nit value and reflect the biogeochemical coupling between nitrification and inorganic carbon dynamics (higher nitrification promotes CO 2 outgassing) and between nitrification and denitrification (higher nitrification promotes higher denitrification), respectively.On the other hand, although varying k nit may also potentially influence estuarine NEM (and thus FC org ) through the nitrification-denitrification and nitrificationprimary production couplings as implemented in the C-GEM biogeochemical module (refer to Table 3), virtually no effects are simulated, confirming that aerobic degradation is by far the dominant process in controlling the organic carbon dynamics in the three idealized estuaries.

Future scenarios
The potential effects of future environmental changes on the biogeochemical functioning of estuaries are assessed by comparing the biogeochemical indicators (NEM, FCO 2 , FC TN and FC TC ) simulated under baseline conditions (year 2000) with those simulated using a future scenario representing the year 2050 (Table 7).Results and integrated rates for the reactions influencing the NEM are summarized in Table 8.Compared to the biogeochemical boundary conditions of the baseline run, the future scenario using the global orchestration projections (Seitzinger et al., 2010) predicts an increase in PO 4 and dissolved inorganic nitrogen (DIN = NH 4 + NO 3 ) loads from rivers (+57 and +29 %, respectively), as well as a decrease in TOC, DSi and SPM (−6, −1 and −17 %, respectively).Simulation results indicate that, in 2050, NEM will still be dominated by heterotrophic degradation and all estuaries will remain largely net heterotrophic (Table 8), although a slight improvement of the estuarine trophic status (less negative NEM values compared to baseline conditions) is simulated.These results essentially reflect the effect of increasing nutrient and decreasing TOC and SPM concentrations on primary production and heterotrophic degradation reactions.In particular, the larger nutrient and light availability, promoting higher NPP rates, and the smaller organic carbon input from rivers, sustaining lower degradation rates, translate into a less negative NEM in all three idealized systems in the future (Table 8).In quantitative terms, the magnitude of NEM variations in the three estuaries is always smaller than 6 %, because the influence of NPP on the overall carbon balance, as well as the predicted drop in organic carbon inputs from rivers, is marginal.
Table 8 also reveals that the CO 2 outgassing flux from the three idealized estuarine systems will decrease in the future.While such trends may appear to contradict a simulated increase in water pCO 2 (data not shown) and decrease in pH (Fig. 12a), the higher atmospheric pCO 2 expected in 2050 (468 µatm) actually leads to a net decrease in the pCO 2 gradient at the air-water interface and to an acidification of the estuarine water masses.As a consequence, the three estuaries will become less important CO 2 sources for the atmosphere in 2050 because the atmospheric CO 2 increase will largely offset the estuarine CO 2 increase (Fig. 12b, Table 8).Furthermore, Table 8 indicates that the mixed and the riverine estuaries reveal low sensitivities to increased atmospheric CO 2 concentrations (−8 and −4 % in CO 2 evasion for the mixed and the riverine cases, respectively), while the marine estuary shows a larger reduction in CO 2 outgassing (−20 %).Although all estuaries remain net CO 2 sources for atmosphere in 2050 when the exchange rate is integrated over the entire estuarine volume, a reversal trend is simulated in the downstream zone of the marine estuary, which will become a net sink for atmospheric CO 2 (Fig. 12b).On the other hand, the mixed and the riverine systems will remain atmospheric CO 2 sources all along their longitudinal profiles.These re- sults essentially reflect the different influences of the adjacent coastal zone on the estuarine biogeochemistry in the three systems.If simulations are carried out until year 2100, assuming a pCO 2 gradient as high as 234 µatm at the marine boundary and an atmospheric pCO 2 of 660 µatm (Moss et al., 2010), the downstream zone of the marine estuary could become strongly undersaturated with respect to the atmospheric equilibrium and largely compensate for CO 2 outgassing in the upper reaches (Fig. 12b).For instance, volume-integrated CO 2 exchange rate over the entire domain indicates that the marine estuary might become a sink for atmospheric CO 2 in 2100.Furthermore, a sink zone could develop in the mixed estuary in 2100, while the riverine system will remain a net source (Fig. 12b).
By definition, FC TC depends on CO 2 outgassing and on the TC riverine input flux.Therefore, the overall decrease in CO 2 emissions simulated in the future in the three idealized systems translates into consistent reductions in FC TC (Table 8).Finally, simulation results also reveal that the efficiency of the three estuaries in removing nitrogen will decrease in the future (FC TN in Table 8), because the projected increase in riverine nitrogen loads (+29 %) is comparatively larger than the enhanced nitrogen removal through denitrification simulated in each idealized estuary (D in Table 8).The outcomes of our scenarios are in overall agreement with the future increase in carbon and nitrogen export to the coastal ocean predicted at global scale (e.g., Kroeze and Seitzinger, 1998;Bauer et al., 2013).Although these projections can only be regarded as a first-order estimate for future trends because of the great deal of uncertainty in predicting future conditions, our findings highlight that, by the end of the century, changes in carbon and nutrient loads from the catchments may lead to a significantly smaller impact on the overall estuarine C balance than the atmospheric CO 2 increase.Furthermore, the latter may be such that estuaries under strong marine influence could experience a shift from being net CO 2 sources to net CO 2 sink, in a manner similar to that already advocated for the coastal ocean (e.g., Mackenzie et al., 2004;Bauer et al., 2013;Regnier et al., 2013a).

Conclusions and outlook
A generic modeling approach was used to quantitatively explore the biogeochemical dynamics across three idealized alluvial estuaries (marine, mixed and riverine) and to relate their behavior to their main hydro-geometrical characteristics.To this end, a comprehensive literature review of biogeochemical parameter values (e.g., rate and half-saturation constants) used in estuarine model applications was performed.This allowed deriving a generic set of parameter values for our simulations including reasonable ranges to perform sensitivity tests.This large literature survey, to our knowledge the first of its kind, highlighted that modeling studies are largely biased towards temperate latitudes.Out of a total of 51 modeling applications, 49 are from temperate regions and two are from the tropics (see map provided in the Supplement).Results from our simulations performed for typical temperate conditions (30-60 • in either hemisphere) are in line with recent literature syntheses (Borges and Abril, 2011;Laruelle et al., 2013) showing that the vast majority of estuaries are net heterotrophic and act as net sources for at- mospheric CO 2 .The average C filtering capacities for baseline conditions are 40, 30 and 22 % for the marine, mixed and riverine estuary, respectively.Extrapolating these filtration rates to the carbon loads of all temperate tidal estuaries worldwide, as defined in Dürr et al. (2011), results in a regional outgassing flux between 0.01 and 0.02 PgC yr −1 , assuming that the total amount of carbon delivered by rivers to these systems is 0.06 PgC yr −1 (calculated from Hartmann et al., 2009, andMayorga et al., 2010).Moreover, a global CO 2 outgassing comprised between 0.04 and 0.07 PgC yr −1 can be calculated by extending these calculations to all tidal estuaries worldwide and assuming that the global amount of carbon delivered by rivers is then 0.17 PgC yr −1 (calculated from Hartmann et al., 2009, andMayorga et al., 2010).Since this estimate ignores that tropical and polar tidal estuaries may process riverine carbon differently than temperate systems, it only represents a first-order quantification of the CO 2 outgassing flux from tidal estuaries worldwide.Nonetheless, it is broadly in line with the recent global estimate calculated by Laruelle et al. (2013) (0.06 PgC yr −1 ).In addition, results for baseline conditions indicate that the three idealized estuaries can retain between 15 and 22 % of the total nitrogen input from the land.On the other hand, sensitivity analysis performed by varying the rate constants for aerobic degradation, denitrification and nitrification over the range of values reported in the literature significantly widens these ranges (Fig. 13).Although the bulk of published N and C retention rates are generally based on local studies (e.g., Soetaert and Herman, 1995;Regnier and Steefel, 1999;Arndt et al., 2009;Billen et al., 2009), while literature syntheses are limited (Nixon et al., 1996;Laruelle, 2009;Regnier at al., 2013b), our estimates confirm the importance of considering the estuarine filter in global and regional C and N budgets.
Overall, our simulations using three idealized estuaries allow identifying the main influence of a marine or riverine character of a system on its biogeochemical functioning (Fig. 13).For instance, simulation results suggest that marine estuaries with short width convergence length (banks that strongly converge in the landward direction) could be efficient filters for nitrogen but relatively limited sources for atmospheric CO 2 .In contrast, mixed and riverine estuaries, owing to their relatively large CO 2 outgassing, could contribute disproportionally to the global estuarine CO 2 budget.The presented approach could thus ultimately enhance our ability to transfer knowledge from well-known estuaries to poorly constrained systems and to predict their biogeochemical functioning even when available data are scarce.Helpful insights are also provided by the sensitivity analysis, which reveals that, for the parameter ranges tested, the uncertainty in aerobic degradation and denitrification rates may translate into large specific-system responses.Finally, prospective simulations for the year 2050 indicate that, while the riverine and mixed estuaries will remain heterotrophic and only marginally affected by river load changes and increase in atmospheric pCO 2 , the marine estuary is likely to become a significant CO 2 sink in its downstream section.Such a change in behavior might offset the current balance between the CO 2 emitted by estuaries and that taken up by continental shelf seas (Andersson et al., 2005;Laruelle et al., 2010Laruelle et al., , 2013Laruelle et al., , 2014;;Cai, 2011;Bauer et al., 2013;Regnier et al., 2013a).
Although this study provides advances on the qualitative and quantitative understanding of the estuarine biogeochemistry in tidal estuaries, our approach also relies on simplifications and abstractions.As a consequence, even if simulation results may be considered as representative of the main tidal, alluvial estuarine classes, they only represent a limited num-ber of estuarine cases that may not cover the extremely wide spectrum of hydro-geometrical and biogeochemical proprieties typically observed in estuaries.For instance, our sensitivity results may not be valid in autotrophic estuaries, where primary production rather than heterotrophic processes controls the estuarine biogeochemistry.In addition, the use of the generic modeling platform C-GEM, whose reaction network does not include an early diagenetic model at this stage, hampers the generalization of simulation results to estuaries that are typically subject to intense biogeochemical processing in sediments.For instance, the lack of a representation of denitrification in sediment might result in overestimations of the N export to the coastal ocean.As a consequence, the implementation of an early diagenetic model for carbon, nutrients and O 2 will be the next logical step to further improve our modeling strategy.
In the future, our generic approach could be applied to perform ensemble runs in order to reduce the predictive uncertainty by covering a wider range of hydro-geometrical characteristics and the uncertainty associated with biogeochemical parameters and conditions, as well as climatological regimes.This design will ultimately provide a better understanding of the estuarine dynamics.Moreover, we believe that a generic approach such as ours can be combined with continuously growing high-resolution environmental databases and, coupled to Earth system models, could be used to derive robust regional and/or global biogeochemical mass budgets by explicitly incorporating the estuarine filter into the estimates.

C. Volta et al.: Linking biogeochemistry to hydro-geometrical variability in tidal estuaries
Appendix A: Sensitivity tests Table A1.List of the parameter values for aerobic degradation (k ox ), denitrification (k denit ) and nitrification (k nit ) rate constants used for the two sensitivity analyses (SA1 and SA2).S1 to S10 refer to different simulations, each corresponding to a specific parameter combination, while BS indicates baseline parameter values (see Sect. 2.3.4,Table 5).

Sensitivity analysis 1 (SA1)
Sensitivity analysis 2 (SA2)          B2) for a reference T of 20 • C. Red lines indicate the T functions as used in the present study (Table 3).Note that only T functions associated with k ox and k denit values not recognized as outliers are displayed.

Figure 2 .
Figure 2. Conceptual scheme of the generic modeling approach (from Volta et al., 2014).Each estuarine type responds in a specific manner to the interdependence between geometry and hydrodynamics and to the first-order control of the hydrodynamics on transport and biogeochemistry.Small panels correspond to longitudinal distribution of (a) A (cross-section area, in m 2 ), B (width, in m), and H (depth, in m); (b) water flow velocity (in m s −1 ); (c) salinity; and (d) O 2 concentration (in µM O 2 ).

Figure 3 .
Figure 3. Variability along the estuarine axis of (a) width B (m), (b) cross-sectional area A (m 2 ) and (c) volume V (m 3) in the three idealized estuaries.Note that the estuarine length EL (km) varies across systems.A is calculated as the product of the tidally averaged depth H (H = 7 m along EL) and B and V is the product between A and x.Profiles are obtained by using geometrical parameters reported in Table1.

Figure 4 .
Figure 4. Conceptual scheme of the C-GEM reaction network as implemented in this study.The new variables and reactions compared to the version presented in Volta et al. (2014) are highlighted in bold.

Figure 6 .
Figure 6.Results from the outlier/distribution analysis for the inhibition term for denitrification (K in,O 2 in µM O 2 ; left panel) and for the nitrification rate constant (k nit in µM N s −1 ; right panel).Dots correspond to parameter values from our literature review.The straight black line represents the theoretical normal distribution.In the left panel, K in,O 2 reveals a normal distribution, while on the right panel k nit displays a distribution strongly skewed by the presence of outliers (blue dots).Data are plotted on similar scales to facilitate the comparison.In each panel, the inset represents the corresponding whisker plot, where boxes contain literature values included between the 25th and the 75th percentiles; the whiskers extend between the maximum and the minimum, beyond which a value is considered an outlier; and blue dots and red lines represent outlier values and the median of literature values, respectively.

Figure 7 .
Figure 7. Longitudinal distributions of tidal amplitude (a), computed as difference between simulated maximum/minimum and average water depths, salinity (b) and SPM (c), modeled in the three representative estuaries using parameters listed in Tables1 and 4and baseline boundary conditions reported in Table7.

Figure 8 .
Figure 8. Longitudinal distributions of chemical species involved in organic (a) and inorganic carbon (b) dynamics in the three representative estuaries.Simulations are forced with parameters reported in Tables 1, 4 and 5 and with the baseline biogeochemical boundary conditions summarized in Table7.

Figure 9 .
Figure 9. Distributions along the three estuarine types of process rates affecting TOC (a), O 2 (b), NO 3 (c), TAlk (d), DIC (e) and H + (f).Net: net rate; R: aerobic degradation; NPP NO 3 : net primary production consuming NO 3 ; NPP: total net primary production (NPP NH 4 + NPP NO 3 ); N: nitrification; D: denitrification; FO 2 :O 2 exchange with the atmosphere; M: phytoplankton mortality; FCO 2 : CO 2 exchange with the atmosphere.Negative gas exchange rates indicate gas evasion from water towards the atmosphere, while positive rates represent gas fluxes from atmosphere to water column.Process contributions are calculated as inRegnier et al. (1997Regnier et al. ( , 2013b)).

Figure 10 .
Figure10.Sensitivity analysis of the biogeochemical indicators (NEM, FCO 2 , FC TC , FC TN ) to changes in aerobic degradation, denitrification and nitrification rate constants over 2 orders of magnitude.For the NEM, values normalized by specific-system organic carbon loads are also provided (FC org ).SA1 and SA2 (black and grey circles, respectively) refer to the two sensitivity analyses performed, while BS (black dots) indicates baseline simulations.The x axes do not report the values of applied parameters but instead the factor by which the parameters are modified compared to their baseline values (BS on the x axes).Values of the reaction rate constants used are summarized in TableA1.All simulations are forced with the baseline biogeochemical scenario for the year 2000 (Table7).

Figure 11 .
Figure 11.Longitudinal profiles of (a) total heterotrophic degradation (TH, in mmol C m −1 day −1 ), (b) CO 2 outgassing (FCO 2 , in mmol C m −1 day −1 ) and (c) denitrification (D, in mmol C m −1 day −1 ) rates in the three idealized estuaries for each of the simulations of the first sensitivity analysis (SA1; see TableA1).S1 to S10 refer to different simulations, each corresponding to a specific parameter combination while BS indicates baseline parameter values (Table5).Black vertical lines indicate the position of the maximal rates simulated by using the BS biogeochemical parameter sets, while arrows displays the longitudinal shift of the biogeochemical fronts by varying reaction rate constants.

Figure 12 .
Figure 12.Longitudinal distributions of pH (a) and depth-as well as width-integrated FCO 2 (b) simulated using a biogeochemical scenario for years 2000 and 2050 in the three idealized estuaries.Note that results from simulations carried out until year 2100 are also reported.The bar graphs in the lower panels represent the volume-integrated FCO 2 for the three simulated years (data are plotted on a different y axis scale).Positive FCO 2 values represent CO 2 emissions, while negative values correspond to CO 2 flux from the atmosphere towards the water.

Figure 13 .
Figure 13.Main results of the sensitivity analysis and scenario simulations summarized as functions of the key hydro-geometrical parameters: freshwater discharge (Q in m 3 s −1 ) and width convergence length (b in m).Results of the sensitivity analysis report ranges corresponding to the maximum and minimum values obtained for any parameter combination of the sensitivity analysis SA1. Results of the scenario simulations are expressed in percentages of relative change between years 2000 and 2050.

k
Figure B1.Temperature functions for aerobic degradation (k ox , left panel) and denitrification (k denit , right panel) rate constants as reported in the reviewed modeling applications (Table B2) for a reference T of 20 • C. Red lines indicate the T functions as used in the present study (Table3).Note that only T functions associated with k ox and k denit values not recognized as outliers are displayed.

Table 2 .
State variables and processes explicitly implemented in C-GEM.

Table 4 .
Volta et al. (2014) used in C-GEM simulations.*indicatesthat a linear variation is applied in the tidal river zone.In such cases, reported parameter values correspond to those imposed at the estuarine upper limit.Bold parameters refer to generic assumptions valid in alluvial estuaries, while other parameters are adapted fromVolta et al. (2014).

Table 5 .
Values of the biogeochemical parameters used in this study.a indicates temperature-dependent parameters.b and c refer to generic parameters estimated from the average and the median of literature values, respectively.Redfield ratios are from Redfield et al.

Table 6 .
Values used for the climate forcings as representative of temperate regions.

Table 7 .
Boundary conditions used for the different scenario simulations.
(Seitzinger et al., 2010)6.0;Mossetal., 2010; IPCC Report, 2013), as well as a rising socioeconomic development and a reactive approach to environmental problems (global orchestration scenario;Seitzinger et al., 2010).The global orchestration scenario is used to constrain future riverine inputs of DIN, PO 4 , TOC, DSi and SPM by year 2050(Seitzinger et al., 2010).It predicts an increase of about 29 and 57 % in DIN and PO 4 , respectively, essentially induced by increased inputs of sewage, fertilizer and manure, and a decrease in TOC, DSi and SPM of about 6, 5 and 17 %, respectively, owing to the influence of damming.On the other hand, no generic predictions are available for riverine DIC and TAlk concentrations and long time-series analyses indicate diverging trends that are inconclusive (e.g., Jones

Table B2 .
Values of temperature-dependent biogeochemical parameters reported from the literature for which a parameter analysis was performed.T corresponds to the temperature in • C. T abs is the absolute temperature.All parameter values are normalized at 20 • C. P B max,L , k main,L , k mort,L , k ox,L , k denit,L and k nit,L represent the parameter values reported in the corresponding literature source (refer to references).indicates that the parameter value is calculated from the mean of values reported for reactive and refractory organic carbon.* * indicates that the parameter value refers to the mean of values associated with different salinities.Outlier values are indicated in italic. *