Modeling Process‐Based Biogeochemical Dynamics in Surface Fresh Waters of Large Watersheds With the IMAGE‐DGNM Framework

Over the last centuries, human activities have exerted increasing pressures on the environment, leading to drastic alterations in the functioning of freshwater bodies (e.g., eutrophication). Global biogeochemical models have proven crucial to investigate interactions between humans, hydrology, and water quality of surface fresh waters. However, most do not account for high‐resolution spatial and temporal variability within watersheds, and they typically lack any representation of benthic dynamics that can drive pollution legacy effects. We present here the Integrated Model to Assess the Global Environment‐Dynamic Global Nutrient Model (IMAGE‐DGNM), which couples global, spatially explicit hydrology and integrated assessment models with process‐based biogeochemistry in surface fresh waters. The new Dynamic In‐Stream Chemistry (DISC) module calculates advective transport from headwaters to estuaries, processes in the water column and in bed sediments, as well as the exchanges between these two compartments. As application examples of IMAGE‐DGNM, we simulate sediment dynamics and nitrogen cycling in two large river basins. We assess in‐stream concentration time series at specific locations, and identify governing processes in transfers along the aquatic continuum. Results highlight the importance of benthic dynamics in watersheds highly perturbed by damming. The implementation of such dynamics within IMAGE‐DGNM allows for including the temporal effect of pollution legacies in large scale water quality studies and shifts in pollutant speciation along river continua. This new framework therefore incorporates new features for large basin to global scale studies that are crucial to better predict the effects on receiving ecosystems and evaluate future environmental management pathways.


Introduction
Over the past centuries, humans have subjected the environment to increasing pressures, by transforming the land, controlling water and sediment flows, and accelerating nutrient and pollutant transfers via inputs to agriculture and emissions from industries and households (Vitousek et al., 1997). Given the extent of global ecosystems' perturbation, there is growing scientific consensus that the Earth has entered a new, human-dominated geological era, the Anthropocene (Lewis & Maslin, 2015). Among these alterations, river damming, chemical contamination, and acceleration of natural nutrient cycles, leading to widespread eutrophication, constitute major global features of changes in river networks (Meybeck, 2003). From an Earth system point of view, surface fresh waters, even though they represent a very small proportion of the water over the globe (<0.1‰), constitute the link between the atmosphere, the pedosphere, the biosphere and oceans (Meybeck, 2003). From an anthropogenic point of view, they constitute a vital natural resource, notably for drinking water and irrigation to support global food production. They are used in various other sectors as well, for example, for navigation, wastewater dilution, hydroelectricity production, cooling water, or leisure activities. In the context of global climate change, growing populations and increasing food and energy needs, modern societies must define integrated strategies that sustainably coordinate water use without compromising the health of ecosystems (Baron et al., 2002).
Achieving this requires a thorough understanding of large-scale, complex interactions between changes in population, land use, climate, and fluxes of water, particles, and dissolved (biogeo)chemical compounds in the Earth system. Global models have proved useful to investigate the main controls of river water quality. Using a lumped regression approach, the Global NEWS model suite constituted a first tool to assess the export of different forms of nitrogen (N), phosphorus (P), and carbon (C) from large river basins to coastal seas (Seitzinger et al., 2005(Seitzinger et al., , 2010. However, Bouwman et al. (2013) pointed out the necessity of including mechanistic process descriptions in large-scale biogeochemical models to identify spatial and temporal variability of water quality within river basins, and to be able to project responses to disturbances in hydrology or nutrient inputs. In this context, Beusen et al. (2015) developed a coupled hydrology and nutrient model, the Integrated Model to Assess the Global Environment-Global Nutrient Model (IMAGE-GNM). IMAGE-GNM mechanistically simulates total nitrogen (TN) and total phosphorus (TP) delivery to and transport in global river networks at a half-degree spatial resolution and at a yearly time step, using a spiraling approach (Newbold et al., 1981;Stream Solute Workshop, 1990;Wollheim et al., 2008) to account for retention along the river continuum. Its embedding within an integrated assessment framework (IMAGE, Stehfest et al., 2014), simulating interactions and feedbacks between the Earth and human systems, allows for a better understanding and prediction of large-scale interactions between land use/human activities and major nutrient transfers in river systems. Nevertheless, the spiraling approach is a lumped, steady-state concept that represents total nutrient cycling (e.g., assimilation from the water column into benthic biomass, temporary retention, and mineralization back into the water column), using a single uptake velocity parameter. This approach therefore does not allow for quantifying the relative amounts of specific biogeochemical compounds (e.g., various forms of nutrients, which have different effects on ecosystem functioning), and the interactions between these compounds (e.g., transformations from one form to another, sorption interactions, substrate-limited processes). Another key feature lacking in the spiraling approach and current large-scale biogeochemical models is the potential deposition and remobilization of accumulated biogeochemical compounds within river networks. This could induce spatial heterogeneity in dominating processes (e.g., onset of denitrification in anoxic sediments at locations of high particle accumulation), and time lags in transfers due to erosion or diffusive release from the benthic layer. This is referred to as in-stream legacies, or internal loading, and can contribute significantly to exports from aquatic systems (James & Barko, 2004;Sharpley et al., 2013). Existing biogeochemical models, such as RIVE (Billen et al., 1994) or QUAL2K (Pelletier et al., 2006), that account for in-stream transformations and benthic dynamics, can simulate such processes at small scales. A major challenge, restraining large-scale applications so far, is to balance detail in process descriptions with availability of knowledge and input/validation data.
With the goal to bridge these knowledge gaps, we present a new modeling framework, the IMAGE-Dynamic Global Nutrient Model (IMAGE-DGNM), designed to simulate biogeochemical cycling of multiple elements and forms in global watersheds, by using the speciation of different nutrient sources from Vilmin et al. (2018). It replaces IMAGE-GNM's spiraling approach with a process-based biogeochemical module: the Dynamic In-Stream Chemistry (DISC) module. DISC computes transfer and transformation processes along the river continuum. Simulated biogeochemical species and transformations are user defined. DISC implements species in the water column and in the benthic layer, in addition to the interactions within and between both of these compartments. We present two applications of the IMAGE-DGNM framework, of differing complexity: (1) the simulation of total suspended solids (TSS), and (2) the simulation of in-stream nitrogen (N) cycling. Results are presented and discussed for the Mississippi and Rhine, two major river basins of different size and land use, where extensive data are available for comparison. These examples illustrate how the framework can be used to investigate the controlling factors of biogeochemical functioning of systems at different spatial and temporal scales, how they can affect nutrient forms in watersheds, and, as a consequence, influence eutrophication risks. By coupling spatially resolved Human and Earth system models with process-based biogeochemistry in global watersheds, this new framework is useful to understand and predict the effects of large-scale socioeconomic forcings and management decisions on the environment.

The IMAGE-DGNM Framework
The IMAGE-DGNM framework results from the coupling of the following ( Figure 1): • IMAGE-GNM (Beusen et al., 2015) that simulates the long-term interactions between human activities and transfers of biogeochemical species in the hydrosphere, that is, delivery processes of nutrients from the land to surface freshwater bodies; • a global, mechanistic, hydrological model (PCR-GLOBWB, Sutanudjaja et al., 2018;van Beek et al., 2011); and • the newly developed process-based DISC module, simulating in-stream biogeochemical cycling. IMAGE-GNM assesses global changes in land use, and natural and anthropogenic agricultural, domestic, or industrial emissions at a yearly time step over the twentieth century. Among other results, it provides estimates of TN and TP delivery from different sources to global river networks. These total nutrient loads are split into different forms, based on a literature review on source speciation, as described by Vilmin et al. (2018). Within IMAGE-GNM, total sediments reaching surface fresh waters are assumed to mainly originate from soil loss. Annual soil loss by rainfall erosion is calculated based on the results from Cerdan et al. (2010), as described in Beusen et al. (2015). Rates of soil loss from arable lands, grasslands, and lands covered with natural vegetation are estimated based on slope, soil texture, and land cover type. These relationships are adjusted to fit the mean erosion loss estimates for Europe and are applied globally. Vilmin et al. (2018) presented a method to add monthly patterns to these inputs, based on the identified controlling hydrological drivers. We further assume that diffuse sources are distributed to the entire river network (based on the relative lengths of the different stream orders in a cell, see Text S1 and Tables S1 and S2 in the supporting information), while point sources are only discharged into streams with Strahler orders ≥6 (herein called "high-order streams"), lakes, or reservoirs.
PCR-GLOBWB simulates the water cycle in the Earth system (van Beek et al., 2011;Sutanudjaja et al., 2018). It dynamically simulates the volumes, surface areas, and discharges of the different waterbodies of river networks, including lakes, reservoirs, and high-order streams. The discharge and characteristics of smaller streams are estimated for each half-degree continental grid cell, based on runoff and the properties of high-order streams, with the parameterization from Wollheim et al. (2008) (see Text S1 and Tables S1 and S2).
With the description of subannual variations and speciation of nutrient delivery to river networks, and the integration of the new DISC module, the IMAGE-DGNM framework now allows for the simulation of transfers of multiple nutrient forms from land to estuaries, and with a higher time resolution, than IMAGE-GNM alone (using the spiraling approach). A major novelty is that it can furthermore quantify the effect of remobilization processes. In the current version of the framework, monthly to yearly output time steps can be used, to be consistent with inputs from IMAGE-GNM and PCR-GLOBWB. The DISC module is here used at a half-degree latitude/longitude spatial resolution, matching that of IMAGE-GNM estimates of the delivery of biogeochemical compounds to surface fresh waters.

The New Globally Applicable Disc Module
The Dynamic In-Stream Chemistry (DISC) module is developed in Python 3. It simulates advection from upstream to downstream cells for high-order streams, based on a local drain direction map, and from headwaters to high-order streams within each grid cell (see Text S1 and Tables S1 and S2). In every cell, in-stream processes are calculated mechanistically for every stream order, or lake/reservoir present in the cell. In the module's current version, vertical stratification is not accounted for.
Consistent with mechanistic models commonly used for basin-scale applications (e.g., SWAT, Neitsch et al., 2011;INCA, Jackson-Blake et al., 2016;RIVERSTRAHLER, Billen et al., 1994), we assume that amounts/concentrations of the different simulated compounds are, together with natural and anthropogenic inputs, controlled by the same physical and biogeochemical processes in all waterbodies, from headwaters to estuaries.
At each computation time step, the amount of species i (X i in moles) for each waterbody type, and in each grid cell, is calculated dynamically by implicitly solving: with V the volume of the waterbody in m 3 ). L i (mol hr −1 ) is the load of i entering the waterbody from upstream grid cells, from lower stream orders within the grid cell, and delivery from external natural and anthropogenic sources within the cell. DISC distinguishes between species that are dissolved or suspended in the water column, transported with the flow by advection, and those in the benthic compartment, only subject to biogeochemical transformations and/or vertical exchanges at the sediment-water interface. Dispersion is assumed negligible, which may limit applications at small spatial scales.
The parameters v p,i (-) and r p (mol hr −1 ) represent the stoichiometric coefficient for i in the process p and the transformation/exchange rate of p. The implementation of these processes, their stoichiometry, and the parameters used to describe them are user defined. Examples of expressions for v p,i and corresponding r p to simulate sediment transfers and N cycling are presented in the following section.
DISC runs can be initialized using a spinup (i.e., permanent state at the first simulation time step) or using a startup file containing the amounts of the simulated species for each waterbody type in every grid cell.
DISC is applicable to individual, multiple, or all river basins over the globe. This new, grid-cell-based, biogeochemical module could be applied in stand-alone mode, using input data from spatially explicit hydrology models and input estimates for the selected studied compounds to river networks. The implicit scheme ensures that the time step used in the model is adapted to correctly assess for the variations in process rates. The output time step of the module is user-defined, and should be consistent with available hydrology and load input data.
The DISC module is generic and can deal with multiple customizable processes. Hereafter, we show applications within the IMAGE-DGNM framework, to explore its potential to quantify sediment dynamics and N cycling at the scale of large watersheds.
For these examples, a sensitivity analysis of major model outputs to the model inputs and parameters is performed using the Latin Hypercube Sampling method (Beusen et al., 2015;Saltelli et al., 2000). The sensitivity of the model output to each parameter is calculated using a standardized regression coefficient (SRC), which represents the relative change of the output due to the relative change of the parameter. For details in the methods and results, see Text S3 and Tables S4 and S5. The results are used in our discussion to support our model results' analysis and hypothesis on sources of discrepancies with observation data.

Modeling Examples
Here, we present two applications of the IMAGE-DGNM framework: the simulation of sediment dynamics and in-stream N cycling within the Mississippi and Rhine river basins. The goal is to present potential uses of the new framework, to better understand the functioning of large river networks. In the following examples, we use a level of complexity for the representation of biogeochemical processes adapted to large-scale applications. All parameters originate from previous studies; we do not perform any calibration.

Sediment Dynamics
Sediment dynamics play a significant role in river biogeochemical functioning. We present an example of the implementation of physical particle transfers within the IMAGE-DGNM framework. We do not include biogeochemical transformations that could affect the amounts of sediments in the water column or the benthic layer (e.g., breakdown of particulate organic matter). This description of sediment dynamics could further be used to better describe the transport of sediment-bound compounds, to assess TSS concentrations as a proxy for light availability for primary production, or to estimate particle deposition and accumulation in river beds.

Expression of Process Rates and Stoichiometry in DISC
Two species are represented: TSS in the water column, and total benthic sediments (TBS) that can be remobilized by the flow. We use a simplified implementation of the formalisms from Hairsine and Rose (1992), with a complexity adapted to global scale applications. Accordingly, we assume that sedimentation and reentrainment occurs simultaneously ( Figure 2 and Table 1). Over a certain threshold of accumulated sediments in the benthic layer, compaction occurs, and part of the bed sediments become "unavailable" for erosion ( Figure 2 and Table 1). Total sediments have no unique molar mass; calculations are performed in kg for this application (i.e., X i , C i , L i , and r p from Equation 1 are in kg, kg m −3 , kg hr −1 , and kg hr −1 , respectively).
Particle deposition is proportional to the amount of TSS in the water column and the settling velocity (u sed,TSS , in m hr −1 ): where H (m) is the depth of the water column.
The equation of the erosion rate describes the fact that a fraction (η, nondimensional) of the total hydraulic power dissipated by the

Journal of Advances in Modeling Earth Systems
flow is used to maintain particles in suspension in the water column. With the assumption, based on the transport capacity theory, that particles tend to be at equilibrium, Vilmin et al. (2015) expressed the total erosion rate as where ρ w (kg m −3 ) and ρ TBS (kg m −3 ) refer to the density of water and TBS. J (-) and U (m hr −1 ) are the flow energy slope and the water velocity, respectively. A (m 2 ) is the total bed area.
We simplify this expression by assuming that sediments have a density of ρ TBS ¼ 2 · ρ w , and that the energy slope of the flow is constant. To avoid numerical oscillations, we include a limitation by the amount of total benthic sediments (TBS) per bed area (with half-saturation constant of K sed ¼ 2 kg m −2 , which corresponds to 1 mm of bed sediments, assuming a porosity ϵ ¼ 0.9). In this way, the total erosion rate is calculated with a single erosion coefficient parameter (k ero , in kg m −3 ): Compaction of benthic sediments is expressed as in Billen et al. (2015): We use a settling velocity for TSS of u sed,TSS ¼1 m hr −1 (Meybeck et al., 1998). k ero is set to 0.02 kg m −3 , using η ¼ 0.01 (optimal value from Sander et al., 2002) and J ¼ 0.001. Compaction occurs over a minimum accumulated sediments per bed area of (TBS/A) lim ¼ 0.5 kg m −2 , at a rate of c max ¼0.0005 hr −1 , as suggested by Billen et al. (2015).

Application of IMAGE-DGNM to the Mississippi and Rhine River Basins
The model is applied to the Mississippi and Rhine river basins. Runs are performed over the whole twentieth century and compared to measurements after 1950. For initialization, we perform a model spinup for the water column processes in 1900. Pelagic and benthic processes are then run for the entire twentieth century.
Our results for the 1950-2000 period therefore account for the effect of 50 years of accumulated sediment legacies.
The Mississippi drains ∼40% of the continental United States, between the Rocky and Appalachian mountains, to the Gulf of Mexico. It is historically subject to high sediment loads, transporting, in pristine conditions ∼400 · 10 6 tonnes of sediments per year from the interior United States to the coasts (Meade & Moody, 2010), predominantly originating from the muddy parts of the Missouri basin (Meade & Moody, 2010;Turner & Rabalais, 2004). Water and sediment transfers have been strongly modified by humans over the past centuries, notably due to the construction of dams, with nowadays thousands of them dispersed throughout the basin (Lehner et al., 2011). TSS data are available at 11 United States Geological Survey monitoring stations between 1970 and 2000 (USGS, 2013). The station closest to the outlet is located in St Francisville (Louisiana).
The Rhine flows from the Northern Alps, and the Black Forest and Vosges mountains to its lowland area, from the German/Dutch border to the mouth, in Hoek van Holland (The Netherlands). The lowland part comprises extensive sedimentary areas (Asselman et al., 2003). Around 50% of the total basin area is used for agricultural production (Buiteveld et al., 2005). The high-order streams of the Rhine basin are less subjected to damming than those of the Mississippi (16 dams identified for 2000 by Lehner et al., 2011 and are represented along the mainstream in PCR-GLOBWB); flows in small streams are however controlled by numerous weirs in the Alpine part of the basin (Asselman et al., 2003). Weekly TSS measurements are available at Lobith, just downstream of the German/Dutch border for the 1950-2000 period (Rijkswaterstaat, 2018).

Journal of Advances in Modeling Earth Systems
VILMIN ET AL.
In both basins, we compare yearly model outputs with discharge-weighted yearly measured concentrations ( Figure 3 and supporting information Text S2 and Table S3).
The model provides good estimates of TSS concentrations throughout the Mississippi river basin ( Figure 3 and Table S3). On average, concentrations are slightly overestimated, by 27% for the 11 USGS monitoring stations. For reference, according to the "Model evaluation guidelines for systematic quantification of accuracy in watershed simulations" from Moriasi et al. (2007), who derived criteria based on the results of smaller-scale simulations for a monthly time step, a bias of 15-30% on sediment estimates corresponds to a "good" model performance. However, this trend is not systematic: Concentrations are underestimated for 6 of the 11 stations. These discrepancies do not systematically concur with underestimates in the main channel discharge estimates from PCR-GLOBWB (data not shown), nor they follow an upstream to downstream trend. This implies that they may arise from the following: (i) local errors in the estimated soil loss reaching the river network or (ii) errors in the estimates of TSS retention in small streams before reaching the mainstream stations. Soil loss delivery to the river network, resulting from the IMAGE-GNM model, constitutes an input to the DISC module and is not the scope of this paper. The second hypothesis is confirmed by the sensitivity analysis, which shows that TSS accumulation in the Mississippi river basin is strongly influenced (SRC ¼ 0.33) by the length of small streams (parameterized by the length ratio R Lsee supporting information Text S1 and Tables S1 and S2). The settling velocity (u sed, TSS ) influences TSS accumulation to a similar extent (SRC ¼ 0.30). Differentiating between several grain size classes for soil loss, with their own settling velocities could help better reproducing spatial patterns of TSS concentrations within the river basin. The sediment export at the river mouth is well reproduced for the 1970-2000 period, with a mean underestimation of TSS concentrations of 21% at the most downstream station (Figure 3).
Even though the order of magnitude of TSS concentrations in the Rhine River at Lobith is well captured, the model overestimates them by a factor of 5 for the 1950-2000 period ( Figure 3 and Table S3). A first possible cause for this discrepancy is an overestimation of the soil loss reaching the river network. In fact, Asselman et al. (2003) estimated that 11.7 · 10 6 tonnes of sediments are supplied to the channels of the Rhine river basin over 1 year, while our soil loss estimates reach 22 · 10 6 tonnes for the year 2000. This suggests that we overestimate sediment delivery to the river network by a factor of 2 in this watershed. Another potential source of error is the underestimation of TSS retention in low-order streams. IMAGE-DGNM estimates very little retention in the river basin, with 21 · 10 6 tonnes of TSS reaching the Lobith monitoring station in 2000. This is in line with previous modeling studies on particle retention in global basins: Beusen et al. (2005) report a trapping efficiency of zero for the Rhine river basin. However, Asselman et al. (2003) found that 27% of the total sediment inputs to the basin reaches the German-Dutch border. They relate this high sediment retention to trapping behind weirs in small headwater streams in the upper parts of the basin, and to floodplain sedimentation in the lower part of the basin (downstream Andernach, Germany, at the entrance of the Lower Rhine stretch), which are not included in the model. This shows that sediment dynamics in the Rhine basin are complex, and points out that the model could be improved by including the representation of slow flowing areas in small streams (connected lakes and small dams/weirs) and floodplain particle retention to more accurately represent the functioning of river basins such as the Rhine.

Particle Accumulation in the Mississippi River Basin Over the Twentieth Century
Since the framework provides good estimates of TSS concentrations throughout the Mississippi river basin, we assume that sediment inputs and process rates within the basin are well represented.
In addition to estimating concentrations at specific locations within a river basin, the model can further be used to assess stocks of biogeochemical compounds in the environment (e.g., sediment accumulation), or identify the major controlling processes for the transfers of these compounds. Here, we calculate the amounts of accumulating bed sediments over the twentieth century over the entire basin for different waterbody types (i.e., small streams, high-order streams, and lakes and reservoirs; Figure 4).
While sediment accumulation in headwaters and the mainstream of the Mississippi have remained stable over time, accumulation in lakes and reservoirs has escalated, especially from the 1940s to the 1970s, when dam construction was highest ( Figure 4). This has been identified as the cause to the decrease in sediment delivery to the coastal zone, and the shrinking of the delta area (Day Jr. et al., 2007).
Particle accumulation can lead to legacy effects within the network. Sediment-bound pollutants can be stored in the system within bed sediments, and released later on due to erosion events or transformations in the benthic layer, leading to time lags in pollutant delivery to downstream systems, and potentially masking the effect of aquatic pollution mitigation strategies (Sharpley et al., 2013). High particle accumulation could furthermore lead to a depletion of oxygen in the benthic layer, triggering changes in the biogeochemical cycling of elements (e.g., desorption, switch in redox potentials). These processes could be implemented in the model to investigate their contribution to pollutant retention and transfers within the basin. With its half-degree spatial resolution, the framework also allows for locating accumulation zones. For the present example, the upstream area of the Missouri basin (northwestern part of the Mississippi river basin) has had large bed sediment stocks over the whole twentieth century, while the Red River basin (southwest) has increasingly accumulated bed sediments, locally reaching similar amounts as in the former basin at the end of the century ( Figure 5).

In-Stream N Cycling
Nitrogen is a key element, controlling species' growth, composition, and diversity in the environment. Implementation of N cycling processes in the IMAGE-DGNM framework allows for the assessment of changes in N speciation and concentrations along the aquatic continuum, and thereby better appraising the potential effect on ecosystems (e.g., primary production, risk of hypoxic events). Recent work has shown that processes occurring both in the water column and in the benthic layer are important for N budgets, not only in the upper reaches (Flipo et al., 2007) but also in large lowland river systems (Raimonet et al., 2015). The simulations of these processes in IMAGE-DGNM further allows for the assessment of the effect of human infrastructure on N retention, for example, increased denitrification in accumulated sediments behind dams.

Expression of Process Rates and Stoichiometry in DISC
We represent four different N forms: , detrital organic N (ON DET ), and organic N (ON) constituting living primary producers (ON PP ). Dissolved oxygen (O 2 ) dynamics are simulated as well, since they drive oxidation (i.e., nitrification: NH þ 4 þ 2O 2 →H 2 O + 2H + +NO − 3 ) and reduction (i.e., denitrification at low O 2 levels: 4NO − 3 þ 5CH 2 O + 4H + →2N 2 + 5CO 2 +7H 2 O) processes. We use the Redfield molar C:N ratio of 106:16 for organic matter constitution and assume that the C:O 2 molar ratio for primary production and respiration is 1. We keep track of ON DET accumulated in the benthic layer (ON DET,b ) for the calculation of the uptake and release of dissolved N forms and O 2 by benthic biogeochemical processes. Gaseous N compounds (e.g., N 2 and N 2 O produced during denitrification) are not represented. Denitrification constitutes a net N loss to the system. The simulated processes are shown in Figure 6. Description of process-rate notations and formalisms and corresponding stoichiometry are provided in Tables 2 and 3. Details on notations and values used for parameters of the N cycle can be found in Table 4. The formalisms used are inspired by those from the RIVE biogeochemical model (www.fire.upmc.fr/rive/).
In the water column, dissolved inorganic forms (DIN, i.e., NH þ 4 and NO − x ) can be taken up for primary production to form ON PP . Primary production is directly related to the maximum photosynthesis rate (Phot, in hr −1 ) that is calculated for each day and for each latitude, based on the photosynthesis-irradiance relation from (Billen et al., 1994). It corresponds to an integration of the photosynthesis at time t and depth z (Phot(t, z), in hr −1 ) over the daylight hours and over the total depth of the water column. Hourly irradiance (I, in J cm −2 hr −1 ) is estimated based on latitude, as described in Billen et al. (1994). Light availability in the water column can be reduced by the presence of suspended sediments. An additional light extinction factor L ext (-) is included to represent this. The growth of primary producers is limited by DIN availability. NH þ 4 is released via respiration/excretion processes, and ON PP is transformed into ON DET upon mortality. N is furthermore mineralized (transformation from ON DET to NH þ 4 ) through the degradation of ON DET (Figure 6). Mineralization is represented by a first-order formalism with a temperature dependency. The proportion of O 2 and NO − x used as electron donors for this oxidation process depends on O 2 availability (Table 3). In oxic conditions, NH þ 4 is oxidized to NO − x via nitrification. O 2 is exchanged at the water-atmosphere interface through re-aeration, thus moving its in-stream concentration toward saturation. Particulate species (i.e., ON DET , ON PP ) can settle to and be incorporated in the benthic layer, and be resuspended by the flow. We consider that, once the primary producer biomass settles to the bottom of the waterbody, it cannot survive and then adds to the ON DET,b pool. Finally, benthic biogeochemical processes can lead to exchanges of DIN and O 2 at the sediment-water interface.
We assume that over a certain threshold of accumulated bed sediments, compaction occurs, and part of the benthic particulate species then become unavailable for erosion or biogeochemical transformations. The available part of ON DET,b can be removed from the benthic layer by organic matter degradation ( Figure 6). This mineralization process, together with benthic nitrification, lead to exchanges of NH þ 4 , NO − x , and O 2 at the sediment-water interface. To simulate these exchanges, we adopt a "vertically integrated dynamic" representation (Soetaert et al., 2000), assuming that dissolved species in the water column are directly produced/consumed during the degradation of benthic particulate organic matter. As a simplification, we consider that a fixed proportion (p nit ) of the NH − x produced by mineralization is instantly oxidized to NO − x via nitrification. We use p nit ¼ 0.5, consistent with the findings from Soetaert et al. (2000) for a large range of mineralization rates.

Application of IMAGE-DGNM to the Mississippi and Rhine River Basins
As with sediment dynamics, we apply the presented implementation of N cycling to the Mississippi and Rhine river basins for 1900-2000, using the same initialization method. The delivery of NH þ 4 , NO − x , and ON DET to global river networks was recently assessed by Vilmin et al. (2018). We  assume that all runoff water reaching the river network is saturated in O 2 and contains 1 μgChl a.L −1 , which we convert to 6.2 · 10 −3 mgN L −1 (using a C:Chl a ratio of 0.035 mgC μgChl a, as suggested by Garnier et al., 1998, and the Redfield C:N molar ratio of 6.625). In both basins, N loads are strongly affected by inputs from agricultural lands (Loos et al., 2009;Turner & Rabalais, 2003). The Rhine river basin, however, receives a greater proportion of DIN loads due to its higher population density, and thus higher contribution of sewage emissions to N delivery to the river network (see Figure S1). TN, ON, NH þ 4 , and NO − x concentration measurements are available at the 11 USGS stations of the Mississippi river basin for the 1970period (USGS, 2013. This allows for the evaluation of the capacity of the model to assess N speciation along the river continuum, in relation with N delivery to the river network and in-stream transformations. For the Rhine, the Netherlands Rijkswaterstaat provides concentrations measurements of TN, NH þ 4 , NO − x , and O 2 at the Lobith monitoring station for the 1950-2000 period (Rijkswaterstaat, 2018). We can thereby evaluate the capacity of the model to mimic oxic conditions as well. Yearly simulated concentrations are compared to measurements in these two basins, using the same methods as for TSS (Figure 7 and Table S3).

Journal of Advances in Modeling Earth Systems
The model captures well TN concentrations and speciation close to the outlets of the Mississippi and Rhine river basins (Figure 7). TN concentrations are overestimated at the outlets by 7% and 38%, respectively. For reference, Moriasi et al. (2007) define a bias <40% for nutrients as a "good" model performance.
Total N concentrations tend to be overestimated in the Mississippi river basin (by 75% for all 11 USGS monitoring stations, and by 38% near the basin outlet, see Figure 7 and Table S3). This is induced by a quasi-systematic overestimation of NO − x concentrations over the basin, the dynamics of other N forms being better simulated. NO − x overestimation is most likely due to an overestimation of its inputs to the river network. This is supported by the results from the sensitivity analysis that show that NO − x export from the Mississippi river basin is mainly controlled by NO − x diffuse inputs (SRC ¼ 0.82). NO − x mainly originates from groundwater input ( Figure S3). Sprague et al. (2011) showed that the flow-normalized N export from the Mississippi has been rising since 1980, increased N river concentrations at low flows being a strong indication of NO − x delivery by groundwater. Decline in groundwater N concentration during transfer through aquifers is incorporated in IMAGE-DGNM, and thus accounted for when estimating the load to the river network. Beusen et al. (2015) found that the TN outflow from groundwater and overall river TN export is extremely sensitive to the assumed aquifer thickness and porosity. A better representation of the physical characteristics of aquifers and their hydrogeochemistry will most likely lead to a more realistic simulation of NO − x reduction by denitrification in these systems. NO − x can be further removed during its transit through the hyporheic zone, which is not represented in the model. Kiel and Cardenas (2014) showed that one quarter of the lateral hyporheic zones (in relation to total river length) in the Mississippi river basin favor denitrification. Denitrification in the benthic layer is simulated here, but enhanced denitrification during groundwater exfiltration, and potential recirculation through the hyporheic zone during transfers along the aquatic continuum are difficult to include in large-scale models (Boulton et al., 2010), and these pathways for N removal are missing in the current implementation. Furthermore, the peak in NH þ 4 and ON concentrations that occurred in the 1980s is not well captured by the model. Since there is no such trend in the calculated loads from sources containing high proportions of these forms, that is, sewage and surface runoff, this could also be due to errors in these inputs. The sensitivity analysis shows that both NH þ 4 and ON exports are strongly influenced by diffuse ON inputs to the river network (SRC ¼ 0.62 and 0.49, respectively), supporting once more that the simulation of N cycling in the river basin would substantially gain from better understanding groundwater flows within the river basin. This is however not the focus of this work. Mineralization and nitrification might also be overestimated for this period, potentially due to an overestimation of dissolved O 2 concentrations.
At the Lobith station along the Rhine, the simulated TN concentrations match measurements (only available after 1970) very well, with a small bias of +7% (Figure 7 and Table S3). TN mainly consists of NO − x , which is also accurately modeled (bias of +2%). The sensitivity analysis shows that, similar to the Mississippi river basin, NO − x and ON export is mainly controlled by diffuse inputs to the river network.
However, in the Rhine basin, point sources have a significant effect as well, not only on NH þ 4 export but also on NO − x , which is produced via nitrification throughout the river network. The magnitude of NH þ 4 concentrations is well reproduced by the model. However, the succession of an increase and a decrease in NH þ 4 concentrations in the 1960s to the 1970s is not captured. This is most likely due to an overestimation of the nitrification rate in the model, presumably because of the overestimation of O 2 concentrations (Figure 7). NH þ 4 export by the Rhine river is positively influenced by T opt,nit , and negatively influenced by k nit (SRC ¼ 0.72 and −0.32, respectively). These nitrification parameters are therefore also worth further investigating to improve the simulation of NH þ 4 dynamics in basins with high point loads. The drop around 1970 in O 2 concentrations is indeed not mimicked by the model. This could be explained by an overestimation of the surface water re-aeration (e.g., k wind ), or an underestimation of the benthic O 2 demand. In the model, O 2 concentrations at the outlet are positively influenced by the number, drainage area, and length of small streams. This might indicate that processes consuming O 2 are underestimated in these water systems.
Finally, in both basins, N export is sensitive to the optimal temperature parameter for mineralization (T opt, min ). NH þ 4 exports are extremely sensitive to the optimal nitrification temperature (T opt,nit ) in the Rhine river basin. This highlights the importance of a correct representation of water temperature and of its influence of biogeochemical processes.

Processes Controlling N Speciation in the Different Waterbodies of the Mississippi and Rhine River Basins
In both the Mississippi and the Rhine river basins, TN speciation evolves along the river continuum, due to in-stream processes. Thereby, the relative proportions of NH þ 4 , NO − x , and ON in total inputs to the river networks and exports to the estuaries differ (Figure 8). To investigate the main controls of these changes, we  Garnier et al. (2000), using diatom parameters for processes related to primary producers. b Billen et al. (1994). c Sum of base respiration and excretion rates. d Calculated using a C:Chl a ratio of 0.035 mgC μgChl a −1 , and the Redfield C:N molar ratio of 106:16. e Assuming 40%, 40%, and 20% of highly biodegradable, moderately biodegradable, and refractory organic matter, respectively. f Seitzinger et al. (2006): upper limit of "suboxic" conditions, for which denitrification occurs. g Thomann and Mueller (1987); Cox (2003). h Helder and de Vries (1983).

Journal of Advances in Modeling Earth Systems
aggregate the rates of the different simulated processes of the N cycle over the entire Mississippi and Rhine river basins, for different waterbody types (small streams-orders <6, high-order streams-orders ≥6, and lakes and reservoirs), and over 10-year periods between 1950 and 2000. Retention of the species X is calculated as 100 · X in − X out X in , where X in and X out are the total inputs of X to the river basin and export to the sea, respectively. As defined, retention of an individual N form can either be caused by accumulation within the system, or transformation to another N form; it can be negative in case of net production of this specific form.
In the Mississippi river basin, N inputs are mainly organic (52-74% of the TN), for the whole study period. ON retention is very high within the basin (58-83%). Moreover, along its transit from the upper to the lower part of the river network, N speciation shifts to a dominance of NO − x (60-85% of TN export, see Figure 8). N cycling in the watershed is controlled by processing in lakes and reservoirs. Due to high water residence times, large amounts of ON are degraded in these systems, both in the water column and in the benthic layer, releasing NH þ 4 . NH þ 4 from external inputs and produced in-stream via mineralization is consumed by nitrification in the water column, leading to a net retention of NH þ 4 of 62-87%. In the model, NO − x production by nitrification is higher than its consumption in the benthic layer for organic matter degradation (i.e., denitrification), which leads to a net production of NO − x over the basin. The latter result might however be inaccurate, since the model overestimates NO − x throughout the Mississippi river basin. In the Rhine river basin, TN inputs are dominated by ON and NH þ 4 forms in 1950 (accounting for 41% and 42% of the TN, respectively, see Figure 8), due to the dominance of sewage as N source to the basin ( Figure S1). During the second half of the twentieth century, ON and NO − x became dominant (both constituting 41% of the TN inputs in 2000) as a result of N abatement in wastewater treatment plants and increasing loads from groundwater bodies draining agricultural lands. In contrast with the Mississippi, the Rhine river basin is not dominated by river damming. N cycling is thus mainly controlled by processes occurring in high-order streams. Since these waterbodies are subject to higher flow velocities and therefore less particle accumulation than lakes and reservoirs, benthic processes have less influence on N cycling within these systems. Over the study period, 13-22% of the ON is retained, mainly due to mineralization in the water column, thereby producing NH þ 4 . NH þ 4 from external inputs and in-stream generation is nitrified, resulting in NO − x dominating N export from the Rhine (48-62% of the TN at the outlet for 1950-2000).
In both basins, in-stream transformations lead to an increase in the proportion of DIN forms along the transfer from land to estuary. The mechanistic implementation of processes of the N cycle allows for the assessment of in-stream concentrations in different types of waterbodies throughout the simulated basins, and of the export loads of different forms to the coasts. These can further be used to better appraise the potential effect of anthropogenic inputs on ecosystems (e.g., increased primary production/eutrophication). The model also allows for the assessment of the effect of human infrastructure, that is, damming, on N transfers. Our results notably show the importance of accounting for benthic processes in reservoirs. For example, in the Mississippi river basin for the year 2000, 260 Gg of ON are mineralized in the benthic layer in these systems (which represents 26% of the overall ON retention in the basin). This benthic mineralization is furthermore associated with NO − x consumption by denitrification due to lower O 2 levels.

Future Framework Improvements and Outlooks
The comparison of simulated and measured concentrations of TSS and different N forms, and its connection with simulated process rates, allows for the identification of potential sources of error in the current state of the framework, which should be addressed in the future.
First, the simulation of retention in small streams has been identified as a potential source of uncertainty in the assessment of particle transfers (e.g., underestimation in the Rhine river basin). This retention notably depends on the water residence times in these systems, driven by their different characteristics (drainage area, stream geometry). The study from Wollheim et al. (2008), whose parameterization we use to describe small stream characteristics, included a validation for TN transfers in the Mississippi river basin. We can therefore assume that this parameterization is appropriate for applications to this specific basin. However, Beusen et al. (2015) showed that IMAGE-GNM was sensitive to characteristics of low-order streams to model TN and TP retention. Data on the characterization of different Strahler orders is emerging in the literature

10.1029/2019MS001796
Journal of Advances in Modeling Earth Systems (Schneider et al., 2017) and should be integrated to the framework. Modifications in the characterization of small streams will most likely have a significant effect in basins draining areas with extremely spatially variable topography. Moreover, in the case of the Rhine watershed, large amounts of TSS are trapped behind weirs in headwater streams (Asselman et al., 2003). In the current implementation of the framework, only lakes and reservoirs connected to high-order streams are represented. Therefore, integrating, when available, data on lakes and reservoirs connected to small streams could further improve the model.
The application of the framework to the simulation of sediment transfers in the Rhine river basin also showed that floodplain accumulation could be a missing sediment sink in the model. Floodplains could constitute accumulation zones for sediments and associated compounds that can stay immobile for long periods after drying, and can be remobilized during flooding events (Ciszewski & Grygar, 2016;Rudorff et al., 2018;Venterink et al., 2003). They have been shown to be active in nutrient retention (e.g., through denitrification or uptake by vegetation) or release (e.g., by litter decomposition), depending on the wetting-drying regime (Baldwin & Mitchell, 2000;Noe et al., 2013). Currently, floodplains are only considered as a source of organic matter due to the scouring of vegetation. Including physical and biogeochemical processes in floodplains would allow for a better simulation of hydrobiogeochemical dynamics in flood-prone river basins, such as the Amazon (Melack et al., 2009) or the Mekong (Hung, 2011). Sediment exchanges with floodplains could be inspired by the work from Cohen et al. (2014), who introduced the effect of floodplains in a global sediment model (regression model). Floodplains could also be introduced as an additional waterbody type, where processes are described with the same complexity as in streams, lakes, and reservoirs.
Finally, thanks to its coupled structure, IMAGE-DGNM could benefit from improvements in the conceptualization of species' inputs and hydrological dynamics. These could, for example, include advances on the assessment of soil loss (potentially overestimated by a factor of ∼2 in the Rhine river basin) and of inputs through groundwater exfiltration. This latter aspect could be refined by improving the physical and hydrogeochemical description of aquifers, and accounting for the effect of biogeochemical activity during transfers through the hyporheic zone. Both soil loss processes and biogeochemical transformations in the hyporheic zone, resulting from complex interactions, are highly variable over space and time (García-Ruiz et al., 2015;Harvey et al., 2013) and are therefore not fully understood at large scales.
The examples presented in this manuscript are used to illustrate potential applications of IMAGE-DGNM. As shown here, we recommend future studies using this framework to be carried out together with sensitivity analysis to identify how model outputs depend on input data and process parameterization (Refsgaard et al., 2007;Saltelli et al., 2004). This is indispensable to communicate results in a critical way and pinpoint potential sources of uncertainty.

Summary and Conclusion
The new DISC module mechanistically simulates transfer and transformation processes along the river network, both in the water column and in the benthic layer. Simulated species, biogeochemical transformations, and corresponding equations are user defined, which makes it a highly flexible module, adaptable to a wide range of applications, such as physical description of sediment transfers, nutrient cycling, coupled cycles, transfer of chemicals, or of pathogens.
Incorporating this new module, IMAGE-DGNM is a globally applicable framework to model the cycling of biogeochemical compounds in watersheds, during their transit from land to estuaries. The applications we presented here to the simulation of sediment dynamics and N cycling in different river basins shows how the framework can be used to • dynamically assess in-stream concentrations of simulated compounds at a half-degree spatial resolution, • dynamically estimate the export of these compounds to the coasts, • identify the processes controlling their transfer along the river continuum (from upstream to downstream areas, and from headwaters to large streams, lakes, and reservoirs), and • identify particle accumulation hotspots within the river network, potentially constituting future sources of sediment-bound biogeochemical species, or triggering anoxic processes.
the incorporation of state-of-the-art knowledge on large-scale hydrological processes and interactions between the Human and Earth systems in the PCR-GLOBWB and IMAGE models. By mechanistically incorporating the processes linking human activities to water quality deterioration risks in watersheds, this new framework constitutes a major step forward in the development of tools for future sustainable land use and water resources management studies.

Data Availability Statement
Data sets used for this research are included in Beusen et al. (2015) and Vilmin et al. (2018).