Coupling global models for hydrology and nutrient loading to simulate nitrogen and phosphorus retention in surface water – description of IMAGE–GNM and analysis of performance

Abstract

The Integrated Model to Assess the Global Environment–Global Nutrient Model (IMAGE–GNM) is a global distributed, spatially explicit model using hydrology as the basis for describing nitrogen (N) and phosphorus (P) delivery to surface water, transport and in-stream retention in rivers, lakes, wetlands and reservoirs. It is part of the integrated assessment model IMAGE, which studies the interaction between society and the environment over prolonged time periods. In the IMAGE–GNM model, grid cells receive water with dissolved and suspended N and P from upstream grid cells; inside grid cells, N and P are delivered to water bodies via diffuse sources (surface runoff, shallow and deep groundwater, riparian zones; litterfall in floodplains; atmospheric deposition) and point sources (wastewater); N and P retention in a water body is calculated on the basis of the residence time of the water and nutrient uptake velocity; subsequently, water and nutrients are transported to downstream grid cells. Differences between model results and observed concentrations for a range of global rivers are acceptable given the global scale of the uncalibrated model. Sensitivity analysis with data for the year 2000 showed that runoff is a major factor for N and P delivery, retention and river export. For both N and P, uptake velocity and all factors used to compute the subgrid in-stream retention are important for total in-stream retention and river export. Soil N budgets, wastewater and all factors determining litterfall in floodplains are important for N delivery to surface water. For P the factors that determine the P content of the soil (soil P content and bulk density) are important factors for delivery and river export.

Introduction

Eutrophication, induced by a surge in anthropogenic nutrient loads to the global freshwater domain (e.g., rivers, lakes and estuaries), has an increasingly negative impact on aquatic ecosystems. In order to ameliorate and reverse this trend, ecological principles must be integrated into environmental management and restoration practices. These actions require a thorough understanding of the interactions between various human-induced disturbances (e.g., climate change, land use change, nutrient loadings and hydrology regulation) and their effects on freshwater systems (Stanley et al., 2010). To fully grasp the human impact on biogeochemical cycles, studies must collectively consider the biogeochemical turnover and exchange among the atmosphere, and the aquatic and terrestrial ecosystems.

Numerical models can assess the interaction between multiple processes in various river basin environments. They can furthermore improve predictions for the regional to global nutrient flux from the land to the ocean. Integrated assessment models (IAM) have established themselves as powerful tools to study future development of complex, large-scale environmental and sustainable development issues. There are at least two key reasons for this: (i) many of these issues are strongly interlinked and integrated models can capture important consequences of these linkages; and (ii) substantial inertia is an inherent property of these problems, which can only be captured using long-term scenarios.

The Integrated Model to Assess the Global Environment (IMAGE) (Stehfest et al., 2014) is an IAM. IMAGE is structured around key global sustainability problems (Fig. 1). Similar to other IAMs, it contains two main subcomponents: i.e., (i) the human system, describing the long-term development of human activities relevant for sustainable development issues, and (ii) the Earth system, describing changes in the natural environment. The two systems are coupled via the impact of human activities on the environment, and via the impacts of environmental change back on the human system.

This paper describes the IMAGE–Global Nutrient Model (GNM), which simulates the fate of nitrogen (N) and phosphorus (P) in surface water arising from concentrated point sources (wastewater from urban and rural populations, and industrial wastewater), and from dispersed (non-point) sources such as agricultural production systems with its fertilizer application and manure management, and natural ecosystems. This global-scale model focuses on prolonged historical periods for testing output results, and future scenarios to analyze consequences of future global change. IMAGE–GNM uses the grid-based global hydrological model PCRaster Global Water Balance (PCR-GLOBWB) (Van Beek et al., 2011) to quantify water stores and fluxes, volume, surface area, and thus depth of water bodies, and water travel time. IMAGE–GNM takes spatially explicit input from the IMAGE land model, including land cover and the annual surface N balance from inputs such as biological N fixation, atmospheric N deposition and the usage of synthetic N fertilizer and animal manure. The IMAGE–GNM model comprises processes such as N removal due to crop harvesting, hay and grass cutting and grazing (Fig. 1). Starting from the soil nutrient budgets, IMAGE–GNM simulates the outflow of nutrients from the soil in combination with emissions from point sources and direct atmospheric deposition to determine the nutrient load to surface water and its fate during transport via surface runoff. It furthermore tracks nutrient transport in groundwater, riparian zones, lakes and reservoirs and in-stream biogeochemical retention processes. Earlier versions of parts of this model, particularly for the nutrient flows towards surface water, have been described previously for N (Van Drecht et al., 2003; Bouwman et al., 2013a), where the retention of N in streams, rivers, lakes and reservoirs was represented by a single, global coefficient. A first step to improve these approaches was the coupling of IMAGE with a hydrological model at the global scale to analyze N retention as pioneered by Wollheim et al. (2008a). Following Wollheim et al. (2008a), the version of IMAGE–GNM presented here uses the nutrient spiraling approach (Newbold et al., 1981) to describe in-stream retention of both total N and total P with a yearly time step.

Scheme of the Integrated Model to Assess the Global Environment (IMAGE) modified from Stehfest et al. (2014).

[Figure omitted. See PDF]

Various other model approaches exist (Bouwman et al., 2013c). The widely used regression models lump the combined effects of nutrient transformations in the continental system into a set of parameters and equations, which can ultimately predict the drainage basin discharge of various geochemical species (e.g., dissolved inorganic and organic, and particulate N, P, C; Seitzinger et al., 2005; Mayorga et al., 2010; Seitzinger et al., 2010). For our purposes, these lumped regression models have limited value, because they both ignore spatial variability of sources and sinks within river basins, and amalgamate all processes in the river continuum. They thus cannot elucidate the nonlinear behavior that results from the interplay between nutrient sources and biogeochemical processes. The SPARROW (SPAtially Referenced Regression On Watershed attributes; Smith et al., 1997; Alexander et al., 2008) model and similar hybrid approaches correlate measured stream nutrient fluxes with spatial data on nutrient sources and landscape characteristics. However, the disadvantage of such an approach is that only a limited time period is covered, while many scientific questions regarding the anthropogenic pressures on the nutrient cycles require prolonged time periods. On the other extreme, there is a range of continuous or event-based distributed watershed-scale models available, which simulate all the components of a landscape, with the hydrology as the basis of calculations. An inventory of such mechanistic models was presented by Borah and Bera (2003). These models usually focus on N while ignoring P and tend to require extensive data that may be difficult to obtain at the spatiotemporal scales of human–climate interactions, and thus are less appropriate to implement in IMAGE–GNM.

In summary, IMAGE–GNM is a global, spatially explicit model, which uses hydrology as the basis for describing N and P delivery to surface water and in-stream transport and retention. It is part of the IAM IMAGE, and used to study the impact of multiple environmental changes over time frames, which capture the mutual feedbacks between humanity and the Earth system. In this manuscript, we compare the model behavior against observations for a number of rivers, and test its sensitivity to a range of model parameter variations to analyze the impact of changing nutrient loading, climate and hydrology.

Model description

General aspects

The IMAGE model utilizes historical data for testing the model behavior, and projections to describe direct and indirect drivers of future global environmental change. Most of these drivers (such as technology and lifestyle assumptions) are used as input in various subcomponents of IMAGE such as GNM (Fig. 1). Clearly, the exogenous assumptions made on these factors need to be consistent. To ensure this, so-called storylines are used, brief descriptions about how the future may unfold, that can be used to derive internally consistent assumptions for the main driving forces of each IMAGE module. Important categories of scenario drivers include demographic factors, economic development, lifestyle and technology change. Among these, population and economic development form a special category as they can be dealt with in a quantitative sense as exogenous model drivers.

The geographical resolution of IMAGE 3.0 is 26 socio-economic world regions (Stehfest et al., 2014). These regions are selected given their relevance for global environmental problems and a relatively high degree of internal coherence. In the Earth system, the key geographic scale is a 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ grid for plant growth, land cover, carbon, nutrient and water cycles. In terms of temporal scale, both systems are run at an annual time step, focusing on long-term trends to capture important inertia aspects of global environmental problems such as simultaneously changing climate and various human activities. Within the Earth system, much shorter time steps are used for water, crop and vegetation modeling. For many applications the IMAGE model deliberately runs over the historical period of 1970 until present-day in order to test model dynamics against key historical trends and then up to 2050, depending on the focus of the analysis. IMAGE–GNM is integrated in the IMAGE model framework, as it has to account for all the drivers that determine the nutrient emissions from point and diffuse sources and their transport. IMAGE–GNM is therefore a distributed model with temporal resolution of 1 year, and a spatially explicit resolution of 0.5 by 0.5 degrees.

Scheme of the model framework with PCR-GLOBWB and IMAGE and the data flows between the models.

[Figure omitted. See PDF]

IMAGE provides land cover and soil budgets for N and P and IMAGE–GNM outputs the nutrient delivery to surface water via surface and subsurface runoff (see Sect. 2.4.2 and 2.4.3) (Fig. 2). IMAGE distinguishes grid cells with natural vegetation or agriculture. Within each agricultural grid cell IMAGE computes distributions of seven crop groups that are aggregated in IMAGE–GNM to larger groups (pastoral grassland, grassland in mixed systems, wetland rice, legumes and upland crops). The soil N budget (N${}_{\mathrm{budget}}$) is calculated for each of these groups and then aggregated to the level of 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ grid cells for individual years as follows: $${\mathrm{N}}_{\mathrm{budget}}={\mathrm{N}}_{\mathrm{fix}}+{\mathrm{N}}_{\mathrm{dep}}+{\mathrm{N}}_{\mathrm{fert}}+{\mathrm{N}}_{\mathrm{man}}-{\mathrm{N}}_{\mathrm{withdr}}-{\mathrm{N}}_{\mathrm{vol}},$$ where N${}_{\mathrm{fix}}$ is biological N fixation (kg), N${}_{\mathrm{dep}}$ is atmospheric N deposition (kg), N${}_{\mathrm{fert}}$ is application of synthetic N fertilizer (kg), N${}_{\mathrm{man}}$ is animal manure (kg), N${}_{\mathrm{withdr}}$ is N removal from via crop harvesting, hay and grass cutting, and grass grazed by animals (kg), and N${}_{\mathrm{vol}}$ is ammonia (NH${}_{\mathrm{3}}$) volatilization (kg). The N budget is prone to erosion, leaching or denitrification, or can accumulate in the soil. Following the approach of Bouwman et al. (2013d), the P budget is assumed to depend on erosion, and soil accumulation. P inputs for the soil budget are fertilizer and animal manure, and outputs are crop and grass withdrawal.

Scheme of the flows of water and nutrients, and retention processes within a grid cell.

[Figure omitted. See PDF]

The data exchange between PCR-GLOBWB and IMAGE–GNM is presented in Fig. 2. Spatial land cover distributions from IMAGE and global climate data from ERA-40 reanalysis (Uppala et al., 2005) are used in PCR-GLOBWB for computing the water balance, runoff and discharge for each year. For each grid cell, IMAGE–GNM provides the delivery of N and P to water bodies via diffuse sources (surface runoff, shallow and deep groundwater, riparian zones) and point sources (wastewater) (Figs. 3 and 4). Grid cells receive water with dissolved and suspended N and P from upstream grid cells, and from diffuse and point sources within the grid cell. In each grid cell, N and P retention in a water body is calculated on the basis of the residence time of the water and nutrient uptake velocity, and subsequently, water and nutrients are transported to downstream grid cells. Discharge is routed to obtain the accumulated water and nutrient flux in each grid cell, through streams, rivers, lakes, wetlands and reservoirs (Fig. 4).

The various submodels for hydrology, spatially explicit nutrient delivery patterns and in-stream retention (Fig. 3), used within IMAGE–GNM are parameterized independently. Furthermore, these parameters are not calibrated in order to better understand the model behavior, identify the lacunae in the data used, and discern the influence of the various processes considered in the model. Instead, the sensitivity of different model outputs to changes in values of input data and model parameters is analyzed in order to explore our model and data.

Although part of the IMAGE framework, GNM can also be used as a stand-alone version, provided that all the input data are in the correct format. For example, land cover data and soil N budgets from various modeling groups could be used (Van Drecht et al., 2005; Fekete et al., 2011). Here we use an update of the nutrient data covering the period 1900–2000 presented by Bouwman et al. (2013d). Also, output from different hydrological models (e.g., Alcamo et al., 2003; Fekete et al., 2011) could be compared.

IMAGE–GNM is written in Python 2.7 code. The complete code is available in the Supplement.

Scheme of the routing of water (with N and P) in a landscape with streams, rivers, lakes, wetlands and reservoirs; each type of water body within a grid cell is defined by an inflow or discharge, depth and area. Floodplains may be temporarily or permanently flooded.

[Figure omitted. See PDF]

Hydrology

Water balance

The land surface in PCR-GLOBWB is represented by a topsoil (0.3 m thick or less) and a subsoil (1.2 m thick or less). Precipitation falls as rain if air temperature exceeds 0 ${}^{\circ }$C, and as snow otherwise. Snow accumulates on the surface, and melt is temperature controlled. Potential evapotranspiration is broken down into canopy transpiration and bare-soil evaporation, which are reduced to an actual evapotranspiration rate based on soil moisture content. Vertical transport in the soil column arises from percolation or capillary rise, depending on the vertical hydraulic gradient present between these layers.

Precipitation and temperature are from New et al. (2000) and downscaled to daily values using the ERA-40 reanalysis (Uppala et al., 2005). Precipitation and temperature were fed directly into the model whereas secondary variables (vapor pressure, wind speed, cloud cover) were used to compute reference potential evapotranspiration using the Penman–Monteith equation according to guidelines of the Food and Agriculture Organization of the United Nations (FAO) (Allen et al., 1998). For the overlapping period 1960–2001, the actual sequence of ERA-40 years was used.

Water drains from the soil column and is delivered as specific runoff to the drainage network, consisting of direct runoff, interflow and base flow. PCR-GLOBWB simulates runoff and converts it to regulated discharge (i.e., including reservoirs; water extraction is ignored), which is used to simulate waterborne nutrient transport. First, total runoff $q_{\mathrm{tot}}$ (m yr${}^{-\mathrm{1}}$) is split into surface runoff ($q_{\mathrm{sro}}$, m yr${}^{-\mathrm{1}}$) and excess water flow ($q_{\mathrm{eff}}$, m yr${}^{-\mathrm{1}}$): $$q_{\mathrm{tot}}=q_{\mathrm{sro}}+q_{\mathrm{eff}}=f_{\mathrm{qsro}}{q}_{\mathrm{tot}}+q_{\mathrm{eff}},$$

where $f_{\mathrm{qsro}}$ is the fraction of surface runoff with respect to total runoff. Surface runoff represents a large proportion of total runoff in locations where drainage into soils is restricted (e.g., urban areas with sealed surfaces, areas covered with impermeable topsoil, and locations with a steep topography) and is represented as $$f_{\mathrm{qsro}}=f_{\mathrm{qsro}}\left(\mathrm{slope}\right)f_{\mathrm{qsro}}\left(\mathrm{texture}\right)f_{\mathrm{qsro}}\left(\mathrm{landuse}\right).$$ Surface runoff is assumed to not be limited ($f_{\mathrm{qsro}}$(texture) $=$ 1.0) in soils with very fine topsoil texture, whereas for loam and sandy loam, and for coarse sand and peat the value $f_{\mathrm{qsro}}$(texture) is adjusted to 0.75 and 0.25, respectively.

The slope-runoff classification for unconsolidated sediments is implemented following Bogena et al. (2005): $$f_{\mathrm{qsro}}\left(\mathrm{slope}\right)=\mathrm{1}-e^{-0.00617\mathrm{MAX}[\mathrm{1},S]},$$ where $S$ is the slope in m km${}^{-\mathrm{1}}$. Since this function is nonlinear, $f_{\mathrm{qsro}}$(slope) is the median value of all 90 m $\times$ 90 m cells within each 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ grid cell. Land use and soil texture can also influence the surface runoff, and these are implemented via the dimensionless factors $f_{\mathrm{qsro}}$(texture) and $f_{\mathrm{qsro}}$(land use), respectively (Velthof et al., 2007, 2009). The soil map used shows dominant soil texture, and has no bare rock class. In areas with bare rock such as in mountainous regions, slopes are generally steep, and Eq. (4) yields high values for $f_{\mathrm{qsro}}$(slope) and thus for $f_{\mathrm{qsro}}$.

Water stagnation may occur in flat land (slope < 20 m km${}^{-\mathrm{1}}$) where soils are saturated based on the Improved Arno Scheme (Todini, 1996; Hageman and Gates, 2003). Soils that are (semi-) permanently saturated are identified as poorly drained areas and are associated with the occurrence of bogs and peat lands. Also, where percolation at the interface between soil and the groundwater reservoir is impeded (e.g., in the case of permafrost), water can stagnate and drain as topographically driven saturated interflow.

When water infiltrates, it can either flow laterally to ditches and streams or vertically to groundwater. IMAGE–GNM implements two groundwater compartments, following Van Drecht et al. (2003), De Wit and Pebesma (2001) and De Wit (2001) (Fig. 3). The shallow groundwater system comprises the top 5 m of the saturated zone where water is retained over short residence times and can either enter the local surface water at short distances (< 1 m) or infiltrate into the deep groundwater system. A 50 m thick deep groundwater layer (Meinardi, 1994), is located below the shallow groundwater system and significantly contributes to the runoff. The water residence time in the deep groundwater system is much higher than that of the shallow groundwater system, as it flows more slowly at greater depths and drains into the fluvial system at greater distances (> 1 km). IMAGE–GNM assumes no deep groundwater presence (i) in areas with non-permeable, consolidated rocks; (ii) in sediments underlying surface waters (rivers, lakes, wetlands, reservoirs); and (iii) in coastal lowlands (< 5 m above sea level) where (artificial) drainage or a high groundwater level persists (Bouwman et al., 2013a).

The excess water flow $q_{\mathrm{eff}}$ (Eq. 5) splits into interflow through the shallow groundwater system ($q_{\mathrm{int}}$, m yr${}^{-\mathrm{1}}$) and deep groundwater runoff ($q_{\mathrm{gwb}}$, m yr${}^{-\mathrm{1}}$) as follows: $$q_{\mathrm{eff}}=(\mathrm{1}-f_{\mathrm{qsro}})q_{\mathrm{tot}}=q_{\mathrm{int}}+q_{\mathrm{gwb}}.$$ The partitioning $f_{\mathrm{qgwb}}\left(p\right)$ of the excess water flow $q_{\mathrm{eff}}$ between these two systems (Fig. 3) is based on the effective porosity ($p$) of the parent material (Table 1). The deep layer (if present) is assumed to have the same characteristics as the surface layer.

Porosity ($p$), the fraction of excess water $Q_{\mathrm{eff}}$ flowing to deep groundwater ($f_{\mathrm{qgwb}}\left(p\right)$), half-life of nitrate in groundwater (dt50${}_{\mathrm{den}}$), activation energy ($E_{\mathrm{a},\mathrm{w}}$) and background P concentration ($C_{\mathrm{PWeath}}$) for various lithological classes.

Lithological class${}^{\mathrm{a}}$	Porosity ($p$)${}^{\mathrm{b}}$	$f_{\mathrm{qgwb}}(p{)}^{\mathrm{c}}$	dt50${}_{\mathrm{den}}$	$E_{\mathrm{a},\mathrm{w}}$	$C_{\mathrm{PWeath}}^{\mathrm{d}}$
	m${}^{\mathrm{3}}$ m${}^{-\mathrm{3}}$	(–)	Year	kJ mol${}^{-\mathrm{1}}$	g m${}^{-\mathrm{3}}$
1. Alluvial deposits	0.15	0.50	2	50	0.0516
2. Loess	0.20	0.67	5	50	0.0256
3. Dunes and shifting sands	0.30	1.00	5	50	0.0790
4. Semi- to unconsolidated sedimentary	0.30	1.00	5	60	0.0248
5. Evaporites	0.20	0.67	5	0	0.0000
6. Carbonated consolidated sedimentary	0.10	0.33	5	0	0.0708
7. Mixed consolidated sedimentary	0.10	0.33	5	60	0.1032
8. Siliciclastic consolidated sediment${}^{\mathrm{e}}$	0.10	0.33	1	60	0.0568
9. Volcanic basic	0.05	0.17	5	50	0.0896
10. Plutonic basic	0.05	0.17	5	50	0.0896
11. Volcanic acid	0.05	0.17	5	60	0.0116
12. Complex lithology	0.02	0.07	5	60	0.0645
13. Plutonic acid	0.02	0.07	5	60	0.0224
14. Metamorphic rock	0.02	0.07	5	60	0.0336
15. Precambrian basement	0.02	0.07	5	60	0.0224

${}^{\mathrm{a}}$ Lithological classes as defined by Dürr et al. (2005). ${}^{\mathrm{b}}$ Porosity values from de Wit (1999). ${}^{\mathrm{c}}$ $f_{\mathrm{qgwb}}\left(p\right)=p/0.3$, 0.3 being maximum porosity. ${}^{\mathrm{d}}$ Background P concentrations ($C_{\mathrm{PWeath}}$) were calculated on the basis of Hartmann et al. (2014). ${}^{\mathrm{e}}$ Weathered shales containing pyrite.

IMAGE-GMN assumes that shallow groundwater interflow moves to the fluvial system via riparian zones (Fig. 3), except in (fractions of) grid cells with wetlands, lakes or large streams, where riparian zones are bypassed. Although riparian zones may only account for a small percentage of the drainage basin, they are critical control points for groundwater and N fluxes within many basins (Vidon and Hill, 2006). Riparian zones along small streams have long ecotone lengths within drainage networks, and may process groundwater N at faster rates than larger nearby water bodies (Bouwman et al., 2013a).

Vegetation and land cover

Vegetation effects are taken into account by partitioning the land surface by fraction into different types. Similarly, spatial variations in soil properties can be accounted for by considering effective values for each of these vegetation types. Soil characteristics are assumed to be constant under changing land cover, except for soil total available water capacity (tawc); the relative distribution of tawc varies with changing root depth distributions based on Canadell et al. (1996). All other soil parameters are from the FAO Digital Soil Map of the World (FAO, 1991) and the World Inventory of Soil property Estimates (WISE) data from the International Soil Reference and Information Center (ISRIC) World Soil Information (Batjes, 1997, 2002). Lithological properties (such as hydraulic conductivity) are derived from a global lithological map (Dürr et al., 2005).

Similar to earlier implementations of PCR-GLOBWB, vegetation parameters are taken from the Olson classification of the global land cover characterization (GLCC) data set with a resolution of 30 arcsec and values assigned using the parameter data set of Hagemann et al. (1999). The parameterization is adjusted to the reconstruction of agricultural land cover for 1900–2000 with 5-year time steps derived from the IMAGE model (Bouwman et al., 2013d) based on historical data (Klein Goldewijk et al., 2010, 2011) in order to achieve consistency between the simulated hydrology and imposed land use.

The land cover reconstruction for the 20th century specifies the fractions of arable land and grassland within each 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ grid cell. To combine this information with the Olson classification, three separate maps at the original resolution of 30 arcsec were created, including (i) Olson classes that were assumed to represent semi-natural vegetation and that were spatially extrapolated per Holdridge life zone (Holdridge, 1967); (ii) Olson classes representing cropland; and (iii) Olson classes representing grassland.

For the reconstructed land cover under the two agriculturally managed conditions, i.e., crops and pasture, all 30 arcsec cells within a 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ cell are ranked in order of decreasing suitability from 0 to 1. This is achieved by first delineating their current extent in the GLCC and ranking on the basis of slope, computed from the Hydro1k database (Verdin and Greenlee, 1996). Next, the adjoining cells are ranked on the basis of the slope parallel distance starting from the delineated areas. These rank orders are then normalized, values near zero indicating the most suitable locations, one indicating the poorest locations, and used to match the IMAGE derived fractions for each 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ cell. In this procedure, cropland has priority, followed by grassland. Any remaining areas are subsequently filled with semi-natural vegetation types. On the basis of the resulting patched land cover, the land cover parameterization for PCR-GLOBWB was then derived.

Drainage network

Drainage density is computed from the Hydro1k data set (Verdin and Greenlee, 1996). The drainage network is based on the DDM30 flow direction map of Döll and Lehner (2002) and the lake characteristics taken from the Global Lakes and Wetlands Database version 1 (GLWD1) product (Lehner and Döll, 2004). Reservoirs are from the Global Reservoir and Dam (GRaND) database (Lehner et al., 2011) and introduced dynamically on the basis of the reported construction year.

The water level in lakes is constant, as the through flow will increase with increasing discharge. The water travel time is determined by the discharge and the volume of the water body. Assuming that flooding occurs once a year and that all river discharge follows the main channel, the travel time in a river with floodplains is determined as follows: $$\mathit{\tau }={\displaystyle \frac{V}{Q-Q_{\mathrm{f}}}},$$ where $\mathit{\tau }$ is the travel time (year), $V$ is the volume of the water body (including river bed) (m${}^{\mathrm{3}}$), $Q$ is the discharge (m${}^{\mathrm{3}}$ yr${}^{-\mathrm{1}}$) and $Q_{\mathrm{f}}$ is the discharge into the flooded area (m${}^{\mathrm{3}}$ yr${}^{-\mathrm{1}}$). While the simulated discharge includes the regulating effect of reservoirs, consumptive water use has not been included as it is difficult to identify its source (groundwater, surface water) and to quantify its spatial distribution with certainty.

Water bodies such as lakes and reservoirs can extend over several 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ grid cells and are included if their volume exceeds that of the channel within a cell. Where more than one reservoir is located within the same grid cell, they are merged and the combined storage and volume assigned to the dominant reservoir. At the start of the simulation, in 1901, 107 out of a total of 132 reservoirs of the GRaND data set are included as 88 spatially individual water bodies, corresponding to 78 % of the reported total volume of 16.4 km${}^{\mathrm{3}}$. For 2000, 5595 out of a total of 6369 reservoirs are included as 3507 spatially individual water bodies, corresponding to 98 % of the reported total volume of 5848.4 km${}^{\mathrm{3}}$. No demand is imposed on the reservoirs and by default they are assigned the purpose of hydropower generation. In absence of pricing generation at the global scale (Haddeland et al., 2006; Adam and Lettenmaier, 2008), this results in an operation that maximizes the available potential energy. In this case, this conforms with 75 % of the maximum storage capacity in absence of detailed global data. The remaining 25 % are reserved to buffer inflow for flood control purposes. Reservoir release is linearly scaled to storage when reservoir storage falls below 30 % of the available capacity. This reduced outflow also results in a realistic, gradual filling of reservoirs after completion of dam construction.

Nutrient delivery to surface water

Surface and subsurface runoff are calculated from the soil N and P budgets on the basis of the hydrological flows provided by PCR-GLOBWB. Other nutrient sources that are directly delivered to surface water included in IMAGE–GNM are wastewater from urban areas, aquaculture, allochthonous organic matter, weathering and atmospheric deposition.

Nutrients directly delivered to surface water

N and P inputs from wastewater for the 20th century are from Morée et al. (2013), and those from freshwater aquaculture are calculated using the country-scale model estimates of Bouwman et al. (2013b) for finfish and Bouwman et al. (2011) for shellfish using data for the period 1950–2000 from FAO (2013); data indicate that prior to 1950 aquaculture production was negligible. N and P emissions from aquaculture are allocated within countries using three weighing factors, i.e., population density, presence of surface water bodies, and mean annual air temperature. For population density, all grid cells with no inhabitants and those with more than 10 000 inhabitants km${}^{-\mathrm{2}}$ are excluded; around an optimum density of 1000 inhabitants km${}^{-\mathrm{2}}$, a steep parabolic function on the left and less steep on the right are used to calculate the weighing. Lakes, reservoirs, rivers and wetlands have the maximum weight for water bodies, and floodplains and intermittent lakes only half of that; all other types have a weight of zero. Grid cells with mean annual air temperature <0 ${}^{\circ }$C are excluded for aquaculture. The three weighing factors are combined by multiplication to obtain the overall weight (range $=$ [0,1]). Then all grid cells with overall probability < 10 % are excluded for aquaculture, yielding the map for allocation for all years. Subsequently, the country production for shellfish and finfish are allocated separately. Grid cells with fish production less than a threshold are excluded for that particular year, and the remaining grid cells are used to allocate the N and P emissions from shellfish and finfish based on the weighing map.

Allochthonous organic matter input to surface water is an important flux in the global C cycle (Cole et al., 2007). This could be an important source of nutrients, but so far its magnitude has not been investigated. Here, estimates of NPP from IMAGE for wetlands and floodplains are used. Part of annual NPP is assumed to be deposited in the water during flooding, and where flooding is temporary, the litter from preceding periods is assumed to be available for transport in the flood water. The mass ratio of litter to belowground inputs of organic matter ranges from 30 : 70 to 70 : 30 (Vogt et al., 1986; Trumbore et al., 1995); 50 % of total NPP is assumed to end in the surface water. N and P inputs to the water are estimated based on a C : N ratio of 100 and a C : P ratio of 1200 (Vitousek, 1984; Vitousek et al., 1988).

The calculation of P release from weathering is based on a recent study (Hartmann et al., 2014), which uses the lithological classes distinguished by Dürr et al. (2005). The lithological classes are available on a 5 by 5 min resolution; hence, the weighted average P concentration within each 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ grid cell is calculated, and the P${}_{\mathrm{RivLoadWeath}}$ (kg P yr${}^{-\mathrm{1}}$) is computed as follows: $$\begin{array}{cc}{\mathrm{P}}_{\mathrm{RivLoadWeath}}& ={10}^{-\mathrm{3}}{C}_{\mathrm{PWeath}}{q}_{\mathrm{tot}}{A}_{\mathrm{gridcell}}{\mathrm{SS}}_{\mathrm{corr}}\\ & & \exp \left(-,{\displaystyle \frac{-E_{\mathrm{a},\mathrm{w}}}{R}},\left({\displaystyle \frac{\mathrm{1}}{K}},-,{\displaystyle \frac{\mathrm{1}}{284}}\right)\right),\end{array}$$ where $C_{\mathrm{PWeath}}$ (g m${}^{-\mathrm{3}}$) is the background concentration specified for each lithological class (Table 1) and derived from river runoff data, $q_{\mathrm{tot}}$ is the total runoff (m yr${}^{-\mathrm{1}}$), $A_{\mathrm{gridcell}}$ is the land area (m${}^{\mathrm{2}}$) in the grid cell considered, SS${}_{\mathrm{corr}}$ is a correction factor for soil shielding, $E_{\mathrm{a},\mathrm{w}}$ is the activation energy (J mol${}^{-\mathrm{1}}$) (Table 1), $K$ the local mean annual air temperature (Kelvin) and $R$ the molar gas constant (8.3144 J mol${}^{-\mathrm{1}}$ K${}^{-\mathrm{1}}$). The soil shielding correction SS${}_{\mathrm{corr}}$ is a correction factor of 0.1 leading to a 90 % reduction for FAO soil units (FAO/Unesco, 1988) Ferralsols, Acrisols, Nitosols, Lixisols, Gleysols (soils with hydromorphic properties) and Histosols (organic soils). For all other soils SS${}_{\mathrm{corr}}=\mathrm{1}$ (no reduction). With this approach, regions with the same lithology but with more precipitation have higher P-weathering losses than regions in dry climates.

Atmospheric N deposition to water bodies is from the ensemble of reactive-transport models for the year 2000 (Dentener et al., 2006), and the years before that were made by scaling the deposition with grid-based emissions of ammonia (Bouwman et al., 2013d). The deposition in floodplains, wetlands and river channels is ignored, because it is already part of the soil N budget, and does not need to be accounted for in periods of flooding.

Surface runoff

IMAGE–GNM distinguishes two surface runoff mobilization pathways for nutrients, i.e., losses from recent nutrient applications in the form of fertilizer, manure or organic matter (N${}_{\mathrm{sro},\mathrm{rec}}$, P${}_{\mathrm{sro},\mathrm{rec}}$) (Hart et al., 2004), and a memory effect (N${}_{\mathrm{sro},\mathrm{mem}}$, P${}_{\mathrm{sro},\mathrm{mem}}$) related to long-term historical changes in soil nutrient inventories (McDowell and Sharpley, 2001; Tarkalson and Mikkelsen, 2004): $${\mathrm{N}}_{\mathrm{sro}}={\mathrm{N}}_{\mathrm{sro},\mathrm{rec}}+{\mathrm{N}}_{\mathrm{sro},\mathrm{mem}}.$$ Estimates of soil loss by rainfall erosion from Cerdan et al. (2010) based on a large database of measurements were used as a basis for calculating P${}_{\mathrm{sro},\mathrm{mem}}$ and N${}_{\mathrm{sro},\mathrm{mem}}$. The approach presented by Cerdan et al. (2010) based on slope, soil texture and land cover type were used to estimate country aggregated soil-loss rates for arable land, grassland and natural vegetation. Soil loss from peat soils was assumed to be low (equal to fine texture). These estimates were then adjusted to obtain the mean erosion loss estimates for Europe (360 t of soil per km${}^{\mathrm{2}}$ for arable fields, 40 t km${}^{-\mathrm{2}}$ for grassland and 15 t km${}^{-\mathrm{2}}$ for natural vegetation). The model was then applied to all grid cells of the world. For global grasslands this yields an erosion rate of 60 t of soil per km${}^{\mathrm{2}}$, which exceeds the European rate by 50 % due to larger erosivity of grasslands in especially tropical and (semi-)arid climates.

As the model keeps track of all inputs and outputs in the soil P budget, the actual P content can be calculated. The initial P stock for the year 1900 in the top 30 cm is taken from Yang et al. (2013). All inputs and outputs of the soil balance are assumed to occur in the top 30 cm; the model replaces P enriched or depleted soil material lost at the surface by erosion with fresh soil material (with the initial soil P content) at the bottom. For N the soil organic C content, which is assumed to be constant over time, is used as a basis to calculate N in eroded soil material using land-use-specific C : N ratios (soil C : N for arable land is 12, for grassland 14 and for soils under natural vegetation 14) (based on Brady, 1990; Batjes, 1996; Guo and Gifford, 2002; McLauchlan, 2006). Hence, with changing land use, the N content in soil erosion loss will also change.

P${}_{\mathrm{sro},\mathrm{rec}}$ and N${}_{\mathrm{sro},\mathrm{rec}}$ are calculated from the N and P input terms (Eq. 1) on the basis of $f_{\mathrm{qsro}}$ (Eq. 4). For N the equation is $${\mathrm{N}}_{\mathrm{sro},\mathrm{rec}}=f_{\mathrm{cal}}{f}_{\mathrm{qsro}}{\mathrm{N}}_{\mathrm{inp}},\left(\mathrm{landuse}\right)$$ where $f_{\mathrm{cal}}$ is a correction coefficient of 0.3 to match the N runoff results of the Miterra model (Velthof et al., 2007, 2009).

Subsurface nitrogen removal and delivery

Subsurface transport of P is neglected, as P is easily absorbed by soil minerals; leaching of P may occur only in P-saturated soils with long histories of heavy over-fertilization; below the saturated soil layer, P will be absorbed into the minerals occurring there, which are low in P. All the positive values of the soil N budget (Eq. 1) are subjected to leaching. Leaching from the top 1 m of soil (or less for thinner soils) is a fraction of the soil N budget excluding the N lost by surface runoff ($f_{\mathrm{leach},\mathrm{soil}}$; Van Drecht et al., 2003): $${\mathrm{N}}_{\mathrm{leach},\mathrm{soil}}=f_{\mathrm{leach},\mathrm{soil}}({\mathrm{N}}_{\mathrm{budget}}-{\mathrm{N}}_{\mathrm{sro}}),$$ where $f_{\mathrm{leach},\mathrm{soil}}$ is $$\begin{array}{cc}{f}_{\mathrm{leach},\mathrm{soil}}& =[\mathrm{1}-\mathrm{MIN}[(f_{\mathrm{climate}}+f_{\mathrm{text}}+f_{\mathrm{drain}}\\ & & +f_{\mathrm{soc}}),\mathrm{1}]]f_{\mathrm{landuse}}.\end{array}$$ The fraction of N lost by denitrification ($f_{\mathrm{den},\mathrm{soil}}$) complements $f_{\mathrm{leach},\mathrm{soil}}(f_{\mathrm{den},\mathrm{soil}}=\mathrm{1}-f_{\mathrm{leach},\mathrm{soil}})$. $f_{\mathrm{text}}$, $f_{\mathrm{drain}}$ and $f_{\mathrm{soc}}$ represent factors that address the soil texture, aeration and soil organic carbon (C) content, respectively (Table 2). Fine-textured soils are more susceptible to reach and maintain anoxia, which favors denitrification, as they are characterized by higher capillary pressures and hold water more tightly than sandy soils. Denitrification rates tend to be higher in poorly drained than in well-drained soils (Bouwman et al., 1993). The soil organic C content is used as a proxy for the C supply, which can have a direct impact on the soil oxygen concentrations. $f_{\mathrm{landuse}}$ is the land use effect on leaching, where arable land has a value of 1, and grassland and natural vegetation a value of 0.36 (Keuskamp et al., 2012).

Denitrification fractions for soil texture, soil organic carbon and soil drainage.

Soil texture	$f_{\mathrm{text}}$ (–)	Soil drainage	$f_{\mathrm{drain}}$ (–)	Soil organic	$f_{\mathrm{SOC}}$ (–)
class				carbon content
Coarse	0.0	Excessively well drained	0.0	< 1 %	0
Medium	0.1	Moderate well drained	0.1	1–3 %	0.1
Fine	0.2	Imperfectly drained	0.2	3–6 %	0.2
Very fine	0.3	Poorly drained	0.3	6–50 %	0.3
Organic	0.0	Very poorly drained	0.4	Organic	0.3

Source: Van Drecht et al. (2003).

The factor $f_{\mathrm{climate}}$ (–) combines the effects of temperature, water residence time, and NO${}_{\mathrm{3}}^{-}$ in the root zone on denitrification rates. $f_{\mathrm{climate}}$ is the product of the temperature effects on denitrification ($f_K$, –) and the mean annual residence time of water and NO${}_{\mathrm{3}}^{-}$ in the root zone ($T_{\mathrm{r},\mathrm{so}}$, yr): $$f_{\mathrm{climate}}=f_K{T}_{\mathrm{r},\mathrm{so}}.$$ The temperature effect $f_K$ follows the Arrhenius equation (Firestone, 1982; Kragt et al., 1990; Shaffer et al., 1991): $$f_K=7.94\times {10}^{12}\exp \left({\displaystyle \frac{-E_{\mathrm{a},\mathrm{d}}}{RK}}\right),$$ where $E_{\mathrm{a},\mathrm{d}}$ is the activation energy (74830 J mol${}^{-\mathrm{1}}$), $K$ the mean annual temperature (Kelvin) and $R$ is the molar gas constant (8.3144 J mol${}^{-\mathrm{1}}$ K${}^{-\mathrm{1}}$). $T_{\mathrm{r},\mathrm{so}}$ is calculated via: $$T_{\mathrm{r},\mathrm{so}}={\displaystyle \frac{\mathrm{tawc}}{q_{\mathrm{eff}}}},$$ where tawc (m) is the total available water capacity for the top 1 m (or less if thinner) of soil and $q_{\mathrm{eff}}$ is described in Eq. 5. Based on the negligible retardation of NO${}_{\mathrm{3}}^{-}$, the water and NO${}_{\mathrm{3}}^{-}$ residence times are assumed to be the same. Soils used for agricultural crops in dry regions with $T_{\mathrm{r},\mathrm{so}}$ < 1 receive a $T_{\mathrm{r},\mathrm{so}}$ value of 1.0 assuming that irrigation is required to grow crops in these locations.

Arid regions under grassland or natural vegetation have long residence times according to Eq. (14), and results in values of $f_{\mathrm{climate}}$ and $f_{\mathrm{den},\mathrm{soil}}$ equal 1, implying that denitrification removes all the N. This representation is not realistic, since N can accumulate in the vadose zone below the root zone as nitrate (Walvoord et al., 2003), and can escape via surface runoff, ammonia-N volatilization and denitrification (Peterjohn and Schlesinger, 1990). It is not possible to quantify the relative contribution of each process (Peterjohn and Schlesinger, 1990), but it is clear that only a negligible part of N surpluses in arid climates is lost by denitrification. Denitrification was thus neglected from the fate of N surplus in soils receiving an annual precipitation of < 3 mm and overlain with grasslands and natural vegetation. For the year 2000, N surplus in the 3100 Mha of global arid lands was 20 Tg.

The N concentration $C_{\mathrm{N}}$ in the excess water leaching from the root zone (depth $z=\mathrm{0}$) is represented by the ratio of leached N over $q_{\mathrm{eff}}$ (Eq. 5): $$C_{\mathrm{N}}(z=\mathrm{0})={\displaystyle \frac{{\mathrm{N}}_{\mathrm{leach}}}{q_{\mathrm{eff}}}}.$$ The groundwater N concentration varies according to the historical year of infiltration into the saturated zone and the denitrification (including anammox) during groundwater advection (Böhlke et al., 2002; Van Drecht et al., 2003). The time available for denitrification is represented by the mean travel time $T_{\mathrm{r},\mathrm{aq}}$, which is the ratio of the specific groundwater volume and the water recharge: $$T_{\mathrm{r},\mathrm{aq}}\left(t\right)=\mathrm{MIN}\left[{\displaystyle \frac{pD}{q_{\mathrm{inflow}}\left(t\right)}},,,,1000\right],$$ where $D$ is aquifer thickness (m) and can either be for shallow groundwater ($D_{\mathrm{sgrw}}=$ 5 m) or for deep groundwater ($D_{\mathrm{dgrw}}=$ 50 m) following Meinardi (1994). $q_{\mathrm{inflow}}$ is either the shallow groundwater recharge ($q_{\mathrm{int}}$, m yr${}^{-\mathrm{1}}$) or deep groundwater recharge, ($q_{\mathrm{gwb}}$, m yr${}^{-\mathrm{1}}$). The vertical drainage of the shallow groundwater feeds the deep groundwater (Fig. 3). The vertical flow distribution for the shallow system is uniform; therefore, the travel time can be equated to the mean travel time. In contrast, travel times for lateral flows to the fluvial system vary considerably. The travel time distribution for lateral flow in a vertical cross section is represented by Meinardi (1994): $$g_{\mathrm{age}}\left(z\right)=-T_{\mathrm{r},\mathrm{aq}}\ln (\mathrm{1}-(z/D\left)\right),$$ where $g_{\mathrm{age}}$ (yr) is the age of groundwater at a specific depth, and $z$ (m) is the depth in the aquifer (i.e., $z=\mathrm{0}$ at the top of the aquifer and $z=D$ at the bottom of the aquifer).

Denitrification takes place during transport in the shallow system along the various flow paths in a homogeneous and isotropic aquifer, drained by parallel rivers or streams. IMAGE–GNM simulates the effects of denitrification in N concentrations at time t and depth z ($C_{\mathrm{N}}(t,z)$) through a first-order degradation reaction, leading to an exponential decay Eq. for the nitrogen concentration: $$C_{\mathrm{N}}(t,z)=C_{\mathrm{N}}\left(t,-,g_{\mathrm{age}},\left(,z,\right),,,\mathrm{0}\right)e^{-{\mathrm{kg}}_{\mathrm{age}}\left(z\right)},$$ where $t$ is time and the decay rate $k$ is obtained via the half-life of nitrate ($dt{50}_{\mathrm{den}}$) due to denitrification: $$k={\displaystyle \frac{\ln \left(\mathrm{2}\right)}{dt{50}_{\mathrm{den}}}}.$$ Lithology can have a direct effect on denitrification, and thus $dt{50}_{\mathrm{den}}$ (Dürr et al., 2005). Siliciclastic material exhibits low $dt{50}_{\mathrm{den}}$ values of 1 yr${}^{-\mathrm{1}}$ , whereas alluvial material has $dt{50}_{\mathrm{den}}$ values of 2 yr${}^{-\mathrm{1}}$ and all other lithology classes have a $dt{50}_{\mathrm{den}}$ value of 5 yr${}^{-\mathrm{1}}$ (Table 1). The N concentration in water percolating to deep groundwater represents the outflow from shallow groundwater. IMAGE-GMN assumes that denitrification is absent in deep groundwater. Although denitrification could occur in organic matter- and/or pyrite-rich deep aquifers, denitrification measurements in the literature have a bias toward high rates (Green et al., 2008), which makes their global assessment difficult.

Following Beusen et al. (2013), nitrogen transported through submarine groundwater discharge (SGD) is excluded from the delivery to rivers and other water bodies. This assumption is justified, since, only 10 % of the gridded map could contribute to SGD. The remaining aquifer discharge in the grid box goes towards streams and rivers.

While urban areas can have an effect in the N loss to the environment (e.g., Foppen, 2002; Wakida and Lerner, 2005; Van den Brink et al., 2007; Nyenje et al., 2010), the total urbanized land represents 0.3 % of the total land area (Angel et al., 2005), and thus it is neglected from the model. The median NH${}_{\mathrm{4}}$ concentration in groundwater of 25 European aquifers is 0.15 mg L${}^{-\mathrm{1}}$ (Shand and Edmunds, 2008), which represents a small part (0.7–1.2 %) of the nitrogen concentration (EEA, 2012), and thus NH${}_{\mathrm{4}}$ in groundwater is also neglected.

N transport and removal in riparian zones

Modeling geochemical processes in riparian zones require a detailed hydrological and geographical information at very high spatial scales, since, even at 0.1 km resolution, the topography of the riparian area cannot be adequately assessed (Vidon and Hill, 2006). IMAGE–GNM therefore uses a conceptual approach.

In riparian zones, denitrification rates depend highly on the local pH (Knowles, 1982; Simek and Cooper, 2002), temperature, water saturation, NO${}_{\mathrm{3}}^{-}$ availability and soil organic carbon availability. Previous laboratory studies of pure cultures have shown that denitrification is maximized at a pH of 6.5 to 7.5, and decreases at both low (below 4) and high (above 10) pH values (Van Cleemput, 1998; Van den Heuvel et al., 2011).

As with soil denitrification, riparian zone denitrification is calculated using dimensionless reduction factors and is based on the characteristics of the groundwater flow, soil and climate. Heterotrophic denitrification is assumed to be highest at pH > 7 (Van den Heuvel et al., 2010). A pH reduction factor $f_{\mathrm{denpH},\mathrm{rip}}$ is then used to reduce the value with decreasing pH, such that $f_{\mathrm{denpH},\mathrm{rip}}=\mathrm{1}$ at pH > 7 and 0 at pH < 3 (Fig. 5). $${\mathrm{N}}_{\mathrm{den},\mathrm{rip}}=f_{\mathrm{den},\mathrm{rip}}{\mathrm{N}}_{\mathrm{in}},$$ where N${}_{\mathrm{in}}$ is the nitrogen that enters the riparian zone from the shallow groundwater. $$\begin{array}{cc}{f}_{\mathrm{den},\mathrm{rip}}& =\mathrm{MIN}\left[\left(f_{\mathrm{climate}},+,f_{\mathrm{text}},+,f_{\mathrm{drain}}\right.\right.\\ & & \left..\left..+,f_{\mathrm{soc}}\right),,,\mathrm{1}\right]f_{\mathrm{denpH},\mathrm{rip}},\end{array}$$ where $f_{\mathrm{climate}}$ is the product of $f_K$ (Eq. 13) and the water (and NO${}_{\mathrm{3}}^{-}$) travel time through the riparian zone ($T_{\mathrm{r},\mathrm{rip}}$). $T_{\mathrm{r},\mathrm{rip}}$ depends on the thickness of the riparian zone ($D_{\mathrm{rip}}\le$ 0.3 m, depending on the soil thickness), on the available water capacity for the top 1m of the riparian zone (tawc), and on the flow of water entering the riparian zone from the shallow groundwater ($q_{\mathrm{int}}$) : $$T_{\mathrm{r},\mathrm{rip}}={\displaystyle \frac{D_{\mathrm{rip}}\mathrm{tawc}}{q_{\mathrm{int}}}}.$$

Reduction fraction ($f_{\mathrm{denpH},\mathrm{rip}}$) of riparian denitrification as a function of soil pH modified from Bouwman et al. (2013a).

[Figure omitted. See PDF]

In-stream nutrient retention

Three processes contribute to N retention, i.e., denitrification, sedimentation and uptake by aquatic plants. Denitrification is generally the major component of N retention (Saunders and Kalff, 2001). P is removed by sedimentation and sorption by sediment (Reddy et al., 1999). Retention in a grid cell is calculated as a first-order approximation according to $$R=\mathrm{1}-\exp \left({\displaystyle \frac{v_{f,E}}{H_{\mathrm{L}}}}\right),$$

where $R$ is the fraction of the nutrient load that is removed (–), $v_f$ is the net uptake velocity (m yr${}^{-\mathrm{1}}$), $E$ is the nutrient considered (N or P), and $H_{\mathrm{L}}$ is the hydraulic load (m yr${}^{-\mathrm{1}}$) obtained from $$H_{\mathrm{L}}={\displaystyle \frac{D}{\mathit{\tau }}},$$ where $D$ is the depth of the water body (m), $\mathit{\tau }$ is the residence time (yr) and $\mathit{\tau }$is calculated from the volume $V$ (m${}^{\mathrm{3}}$) of the water body and the discharge $Q$ (m${}^{\mathrm{3}}$ yr${}^{-\mathrm{1}}$): $$\mathit{\tau }={\displaystyle \frac{V}{Q}}$$ for all water bodies except for river channels and floodplains where the discharge $Q$ is reduced by the water volume in the floodplains (see Eq. 6). In this approach hydrological (defined by $H_{\mathrm{L}}$) and biological and chemical factors (defined by $v_f$) controlling retention are isolated, assuming first-order kinetics is applicable (i.e., areal uptake changes linearly with concentration).

Net uptake velocity is different for each element $E$ (N or P). For N, the basic value for all water body types of 35 m yr${}^{-\mathrm{1}}$ taken from (Wollheim et al., 2006, 2008a) is modified based on temperature and N concentration: $$v_{f,\mathrm{N}}=35f\left(t\right)f\left(C_{\mathrm{N}}\right),$$ where $t$ is annual mean temperature (${}^{\circ }$C) and $C_{\mathrm{N}}$ is the N concentration in the water. $f(C_{\mathrm{N}}$) describes the effect of concentration on denitrification as a result of electron donor limitation in the case of high N loads; the results of Mulholland et al. (2008) were mimicked by assuming a decrease of $f(C_{\mathrm{N}}$) from a value of 7.2 at $C_{\mathrm{N}}=0.0001$ mg L${}^{-\mathrm{1}}$ to 1 for $C_{\mathrm{N}}=\mathrm{1}$ mg L${}^{-\mathrm{1}}$, a further decrease to 0.37 for $C_{\mathrm{N}}=$ 100 mg L${}^{-\mathrm{1}}$ and constant at higher concentrations.

The temperature effect $f\left(t\right)$ is calculated as $$f\left(t\right)={\propto }^{t-20}(t-20),$$ where $\mathit{\alpha }=1.0717$ for N (following Wollheim et al., 2008a and references therein) and $\mathit{\alpha }=1.06$ for P (following Marcé and Armengol, 2009).

For P, the basic value for $v_f$ of 44.5 m yr${}^{-\mathrm{1}}$ taken from Marcé and Armengol (2009) is used for all water body types, with a modification based on temperature: $$v_{f,\mathrm{P}}=44.5f\left(t\right).$$ The drainage network of PCR-GLOBWB represents streams and rivers of Strahler order (Strahler, 1957) 6 and higher. The parameterization of lower-order streams follows the approach presented by Wollheim et al. (2008b). A globally uniform subgrid river network is included for all grid cells without lakes or reservoirs. It is assumed that PCR-GLOBWB has one river of order 6 in each grid cell, and all lower-order rivers are lacking. The river network is then defined on the basis of stream length and basin area of the first-order river. The mean length ratio $R_{\mathrm{L}}$ (–) is used to calculate the stream length of the next higher order the river according to $$L_n=L_{\mathrm{1}}{R}_{\mathrm{L}}^{(n-\mathrm{1})},$$ with $L_n$ being the stream length of order $n$ (km); $L_{\mathrm{1}}=1.6$ km. The drainage area ratio $R_{\mathrm{a}}$ (–) is used to calculate the basin area for higher-order stream as follows: $$A_n=A_{\mathrm{1}}{R}_{\mathrm{a}}^{(n-\mathrm{1})},$$ where $A_n$ is basin area of order $n$ in km${}^{\mathrm{2}}$; $A_{\mathrm{1}}=2.6$ km${}^{\mathrm{2}}$. With the stream number ratio $R_b$ (–) the number of lower-order streams is calculated as $$R_n=R_b^{(\mathrm{6}-n)},$$ with $R_n$ being the number of streams of order $n$ in this grid cell; $R_b=4.5$. The discharge for each stream is calculated with the runoff ($q$): $$Q_n=q{A}_n{C}_Q,$$ with the discharge of stream order $n$ ($Q_n$) in m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$, runoff in mm yr${}^{-\mathrm{1}}$ and $C_Q$ the unit conversion ($C_Q=1000/(3600\times 24\times 365)$). The midpoint discharge of a stream length of order n is calculated as $$Q_{\mathrm{mid},n}=Q_n+0.5Q_{n-\mathrm{1}}.$$ The width of the stream of order n is calculated as $$W_n=A(Q_{\mathrm{mid},n}{)}^B,$$

where $W_n=$ width (m), $A$ is a constant ($A=8.3$ m) and coefficient $B=0.52$. It is now possible to calculate the hydrologic load ($H_{\mathrm{L}}$) and thus the retention of the stream according to $$H_{\mathrm{L}}={\displaystyle \frac{C_{Q\mathrm{1}}{Q}_{\mathrm{mid},n}}{L_n{W}_n{C}_{Q\mathrm{2}}}},$$ with $C_{Q\mathrm{1}}$ being the conversion from seconds to years ($C_{Q\mathrm{1}}=3600\times 24\times 365$), $C_{Q\mathrm{2}}$ the conversion from km to m (1000) and $H_{\mathrm{L}}$ in m yr${}^{-\mathrm{1}}$. The local diffuse load in a grid cell is spatially uniformly distributed over the streams. Here, the fraction of the total stream length per order is used to calculate the distribution of the load. The direct load is allocated to stream-order n as follows: $$F_{\mathrm{d},n}={\displaystyle \frac{R_n{L}_n}{{\sum }_{i=\mathrm{1}}^{\mathrm{6}}{R}_i{L}_i}},$$ where $F_{\mathrm{d},n}$ is the fraction of the total load, which is direct input for streams of order $n$. The pathway of the outflow of the streams is determined according to a matrix ${\mathbf{T}}_{i,j}$ representing the fraction of the outflow of stream-order $i$ to stream-order $j$, whereby ${\mathbf{T}}_{i,j}=0.0$ for $i\ge j$. For $i$ < $j$, ${\mathbf{T}}_{i,j}$ is calculated as follows: $${\mathbf{T}}_{i,j}={\displaystyle \frac{R_j{L}_j}{{\sum }_{k=i+\mathrm{1}}^{\mathrm{6}}{R}_k{L}_k}}.$$ The calculation of the retention is performed for each stream order, starting with order $n=\mathrm{1}$, and is identical to the calculation of the PCR-GLOBWB schematization. The load of a stream is the sum of the direct load and the sum of the outflow from lower-order streams.

Model parameters included in the sensitivity analysis, their symbol and description, for which nutrient it is used, and the standard, minimum, mode and maximum value considered for the sampling procedure. Parameters are listed in alphabetical order of their symbol.

Symbol	Description	Nutrient	Distribution${}^{\ast }$	Standard	Min.	Max.
$A$	Width factor	N/P	U3	8.3	7.5	9.1
$A_{\mathrm{1}}$	Drainage area first-order stream	N/P	U3	2.6	2.3	2.9
$A_{\mathrm{flooding}}$	Area of flooding areas	N/P	U1	1.0	0.9	1.1
$B$	Width exponent	N/P	U3	0.52	0.47	0.57
$B_{\mathrm{soil}}$	Bulk density of the soil	N/P	U1	1.0	0.9	1.1
CN${}_{\mathrm{gnpp}}$	CN weight ratio of gnpp in flooding areas	N	U3	100	90	110
CN${}_{\mathrm{soil},\mathrm{crop}}$	CN weight ratio of soil loss under crops	N	U3	12	11	13
CN${}_{\mathrm{soil},\mathrm{grass}}$	CN weight ratio of soil loss under grassland	N	U3	14	12.5	15.5
CN${}_{\mathrm{soil},\mathrm{nat}}$	CN weight ratio of soil loss under natural ecosystems	N	U3	14	12.5	15.5
CP${}_{\mathrm{aomi}}$	CP weight ratio of gnpp in flooding areas	P	U3	1200	1080	1320
$C_{\mathrm{sro},\mathrm{N}}$	Correction coefficient for N in surface runoff	N	U3	0.3	0.27	0.33
$C_{\mathrm{sro},\mathrm{P}}$	Correction constant for P in surface runoff	P	U3	0.3	0.27	0.33
$D_{\mathrm{dgrw}}$	Thickness of deep groundwater system	N	U3	50.0	45	55
$D_{\mathrm{flooding}}$	Depth of flooding areas	N/P	U1	1.0	0.9	1.1
$D_{\mathrm{rip}}$	Thickness of riparian zone	N	U3	0.3	0.27	0.33
$D_{\mathrm{sgrw}}$	Thickness of shallow groundwater system	N	U3	5.0	4.5	5.5
dt50${}_{\mathrm{den},\mathrm{dgrw}}$	Half-life of nitrate in deep groundwater	N	U3	$\mathrm{\infty }$	50.0	100.0
dt50${}_{\mathrm{den},\mathrm{sgrw}}$	Half-life of nitrate in shallow groundwater	N	U1	1.0	0.9	1.1
$F_{\mathrm{aomi}}$	Reduction factor for litter load to surface water	N/P	U1	0.5	0.45	0.55
$F_{\mathrm{leach},\mathrm{crop}}$	Reduction fraction of N towards the shallow groundwater system	N	U3	1.0	0.9	1.0
$F_{\mathrm{leach},\mathrm{grass}}$	Reduction fraction of N towards the shallow groundwater system	N	U3	0.36	0.32	0.4
$F_{\mathrm{leach},\mathrm{nat}}$	Reduction fraction of N towards the shallow groundwater system	N	U3	0.36	0.32	0.4
$f_{\mathrm{qgwb}}$	Fraction of $q_{\mathrm{eff}}$ that flows towards the deep system	N	U1	1.0	0.9	1.1
$f_{\mathrm{qsro}}$	Overall runoff fraction	N/P	U1	1.0	0.9	1.1
$f_{\mathrm{qsro}}$(crops)	Land use effect on surface runoff for soils under crops	N/P	T2	1.0	0.75	1.0
$f_{\mathrm{qsro}}$(grass)	Land use effect on surface runoff for soils under grassland	N/P	T1	0.25	0.125	0.5
$f_{\mathrm{qsro}}$(nat)	Land use effect on surface runoff for soils in natural ecosystems	N/P	T3	0.125	0.1	0.3
AOMI	Litterfall in flooding areas	N/P	U1	1.0	0.9	1.1
$L_{\mathrm{1}}$	Mean length first-order stream	N/P	U3	1.6	1.4	1.8
N${}_{\mathrm{aqua}}$	N load from aquaculture	N	U1	1.0	0.9	1.1
N${}_{\mathrm{budget},\mathrm{crops}}$	N budgets in croplands	N	U1	1.0	0.9	1.1
N${}_{\mathrm{budget},\mathrm{grass}}$	N budget in grasslands	N	U1	1.0	0.9	1.1
N${}_{\mathrm{budget},\mathrm{nat}}$	N budget in natural ecosystems	N	U1	1.0	0.9	1.1
N${}_{\mathrm{conc},\mathrm{high}}$	Retention multiplier for retention at high N concentrations	N	U3	0.3	0.2	0.4
N${}_{\mathrm{conc},\mathrm{low}}$	Retention multiplier for retention at low N concentrations	N	U3	7	6	9
N${}_{\mathrm{depo}}$	N deposition on surface water	N	U1	1.0	0.9	1.1
N${}_{\mathrm{point}}$	N from point sources	N	U1	1.0	0.9	1.1
N${}_{\mathrm{uptake},\mathrm{crops}}$	N uptake in croplands	N	U1	1.0	0.9	1.1
N${}_{\mathrm{uptake},\mathrm{grass}}$	N uptake in grasslands	N	U1	1.0	0.9	1.1
P${}_{\mathrm{aqua}}$	P load from aquaculture	P	U1	1.0	0.9	1.1
P${}_{\mathrm{budget},\mathrm{crops}}$	P budgets in croplands	P	U1	1.0	0.9	1.1
P${}_{\mathrm{budget},\mathrm{grass}}$	P budget in grasslands	P	U1	1.0	0.9	1.1
P${}_{\mathrm{budget},\mathrm{nat}}$	P budget in natural ecosystems	P	U1	1.0	0.9	1.1
Poros	Porosity of aquifer material	N	U1	1.0	0.9	1.1
P${}_{\mathrm{point}}$	P from point sources	P	U1	1.0	0.9	1.1
P${}_{\mathrm{soil}}$	P content of the soil	P	U1	1.0	0.9	1.1
P${}_{\mathrm{uptake},\mathrm{crops}}$	P uptake in croplands	P	U1	1.0	0.9	1.1
P${}_{\mathrm{uptake},\mathrm{grass}}$	P uptake in grasslands	P	U1	1.0	0.9	1.1
P$v_{f,\mathrm{wetland}}$	Net uptake velocity for wetlands	P	U3	44.5	40	49
CP${}_{\mathrm{Weath}}$	P content of per lithology class	N	U1	1.0	0.9	1.1
${}_{\mathrm{qtot}}$	Runoff (total)	N/P	U1	1.0	0.9	1.1
$R_{\mathrm{a}}$	Drainage area ratio	N/P	U3	4.7	4.2	5.2
$R_b$	Stream number ratio	N/P	U3	4.5	4.05	4.95
$R_{\mathrm{L}}$	Mean length ratio	N/P	U3	2.3	2.0	2.6
Temp	Mean annual air temperature	N/P	U2	0.0	$-$1.0	1.0
$v_{f,\mathrm{lake}}$	Net uptake velocity for lakes	N	U3	35	32	38
$v_{f,\mathrm{lake}}$	Net uptake velocity for lakes	P	U3	44.5	40	49
$v_{f,\mathrm{reservoir}}$	Net uptake velocity for reservoirs	N	U3	35	32	38
$v_{f,\mathrm{reservoir}}$	Net uptake velocity for reservoirs	P	U3	44.5	40	49
$v_{f,\mathrm{river}}$	Net uptake velocity for rivers	N	U3	35	32	38
$v_{f,\mathrm{river}}$	Net uptake velocity for rivers	P	U3	44.5	40	49
$v_{f,\mathrm{wetland}}$	Net uptake velocity for wetlands	N	U3	35	32	38
$V_{\mathrm{water}}$	Water volume of all water bodies	N/P	U1	1.0	0.9	1.1

${}^{\ast }$ Samples values are applied to all grid cells. For sampling, either uniform of triangular distributions are used. A triangular distribution is a continuous probability distribution with lower limit $a$, upper limit $b$ and mode $c$, where $a\le c\le b$. The probability to sample a point depends on the skewness of the triangle. In the case of dt50${}_{\mathrm{den},\mathrm{dgrw}}$, $ac=bc$, and probability to sample a point on the left and right hand side of $c$ is the same. In other cases, for example, $f_{\mathrm{qsro}}$(crops) is a fraction (range $=$ [0,1]), with standard value of 1.0. To achieve a high probability to sample close to 1.0, the triangle is designed with $b=\mathrm{1}$ and $c$ is close to 1. For some of the above distributions the expected value is not equal to the standard. Since the calculated $R^{\mathrm{2}}$ for all output parameters exceeds 0.99, this approach for analyzing the sensitivity is still valid. The distributions used are U1. Uniform; values are multipliers for standard values on a grid cell basis. U2. Uniform; values are added to the standard values on a grid cell basis. U3. Uniform; values are used as such. T1. Triangular; values between 0.125 and 0.5 with an expected value of 0.25. T2. Triangular; values between 0.75 and 1.0 with an expected value of 0.995. T3. Triangular; values between 0.1 and 0.3 with an expected value of 0.125.

Data analysis

For the comparison of observations for individual monitoring stations or ad hoc measurements in rivers and simulated concentrations of river water, we use the root mean squared error (RMSE) expressed as a percentage. RMSE is calculated as follows: $$\mathrm{RMSE}={\displaystyle \frac{100}{\stackrel{\mathrm{‾}}{O}}}\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(O_i-M_i{)}^{\mathrm{2}}}{n}}},$$ where $\stackrel{\mathrm{‾}}{O}$ is the mean of the observations, $O_i$ is observation $i$, $M_i$ is the simulated value $i$ and $n$ is the number of data pairs. We consider values of 50 % acceptable in view of the global scale of the model.

The sensitivity of the modeled delivery, retention and river export for the year 2000 to variation of 48 model parameters for N and 34 for P is based on parameter-specific distributions between a minimum and maximum value around the standard parameter values (Table 3). The sensitivity analysis was performed using the Latin hypercube sampling (LHS) technique (Saltelli et al., 2000). LHS is a multi-parameter, stratified sample method based on subdividing the range of each of the $k$ parameters into disjunct equiprobable intervals or runs (Num). By sampling one value in each of the Num intervals according to the associated distribution in this interval, Num sampled values are obtained for each parameter. Num was 500 for P and 750 for N.

The sampled values for the first model parameter are randomly paired to the samples of the second parameter, and these pairs are subsequently randomly combined with the samples of the third source, and so forth. This results in an LHS consisting of Num combinations of $k$ parameters. The parameter space is thus representatively sampled with a limited number of samples.

The uncertainty contributions of each input parameter ($X_i$) can be further assessed by combining LHS with linear regressions with respect to the model outputs ($Y_i$) via Saltelli et al. (2000, 2004): $$Y={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}{X}_{\mathrm{1}}+{\mathit{\beta }}_{\mathrm{2}}{X}_{\mathrm{2}}\mathrm{\cdots }+{\mathit{\beta }}_n{X}_n+e,$$ where ${\mathit{\beta }}_i$ is the so-called ordinary regression coefficient and $e$ is the error of the approximation. The linear regression model can be evaluated using the coefficient of determination ($R^{\mathrm{2}}$), which represents the $Y$ variation as explained by $Y-e$. ${\mathit{\beta }}_i$ depends on the scale and dimension of $X_i$, the sensitivity study can be normalized by rescaling the regression equation using of the standard deviations for $Y$ and $X$ (${\mathit{\sigma }}_Y$ and ${\mathit{\sigma }}_{Xi}$, respectively) and calculating the standardized regression coefficient (SRC${}_i$): $${\mathrm{SRC}}_i={\mathit{\beta }}_i{\displaystyle \frac{{\mathit{\sigma }}_{X_i}}{{\mathit{\sigma }}_Y}}.$$ SRC${}_i$ can take values in the interval [$-$1,1]. SRC is the relative change $\mathrm{\Delta }Y/{\mathit{\sigma }}_Y$ of $Y$ due to the relative change $\mathrm{\Delta }{X}_i/{\mathit{\sigma }}_{X_i}$ of the parameter $X_i$ considered (both with respect to their standard deviation $\mathit{\sigma }$). Hence, SRC${}_i$ is independent of the units, scale and size of the parameters, and thus sensitivity analysis comes close to an uncertainty analysis. A positive SRC${}_i$ value indicates that increasing a parameter value will cause an increase in the calculated model output, while a negative value indicates a decrease in the output considered caused by a parameter increase.

The sum of squares of SRC${}_i$values of all parameters equals the coefficient of determination ($R^{\mathrm{2}}$), which for a perfect fit equals 1. Hence, SRC${}_i^{\mathrm{2}}/R^{\mathrm{2}}$ yields the contribution of parameter $X_i$ to $Y$. For example, a parameter $X_i$ with SRC${}_i=0.1$ adds 0.01 or 1 % to $Y$ in case $R^{\mathrm{2}}$ equals 1.

Comparison of modeled (black line) and measured (light blue, and aggregated yearly) discharge-weighed concentrations of total N in the rivers Mississippi (a–c), Rhine (d–f) and Meuse (g–i). Panels on the left are comparisons over time; panels in the center represent plots of simulations vs. observations with a 1 : 1 line, and panels on the right are the concentrations vs. the residues (observation minus simulation) with a regression line.

[Figure omitted. See PDF]

Analysis of the model results

Comparison with measurement data

We first compared the IMAGE–GNM model results with observed concentrations for two stations in the rivers Rhine and Meuse and at 11 stations in the Mississippi, USA (see Supplement). Stations near the river mouth (Lobith at the Rhine, Eysden at the Meuse, and St. Francisville, Louisiana, for the Mississippi) are shown first. The latter station was selected for comparison due to its widespread use in literature, for example by the US Geological Survey analysis of water quality (US Geological Survey, 2009). The measured concentrations were first aggregated to annual discharge-weighed concentrations, whereby for the US data years with < 6 observations were excluded. The model performance for the river Rhine for N concentrations (RMSE $=$ 15 %) is better than for the Meuse and Mississippi (Fig. 6a, b, d, e, g, h). IMAGE–GNM overestimates N concentrations in the river Meuse (RMSE $=$ 31 %) in almost all years; the model underestimates N concentrations in the early 1980s for the Mississippi, while its performance is better from the second half of the 1980s (RMSE for Mississippi $=$ 23 %). P concentrations in the Mississippi are consistently underestimated (RMSE $=$ 51 %) (Fig. 7a, b, d, e, g, h). P concentrations are overestimated in the Rhine in all years with data, although the declining trend is well captured (RMSE $=$ 28 %). The modeled P concentrations are close to observations in the Meuse, with deviations in both directions (RMSE $=$ 36 %).

The residues (observation minus simulation) for the observed vs. simulated concentrations of N and P (Figs. 6c and 7c) in the Mississippi show a very clear trend from overestimation at low concentrations to underestimation at high concentrations. The residues show a trend in the Rhine, with a slight increase along with increasing concentrations (Figs. 6f and 7f). The Meuse also shows such trends, although less clear. For P the residue increases with increasing concentration, and for N the opposite occurs (Figs. 6i and 7i).

Since the deviations from observed concentrations can stem from errors in the hydrology, we compared the simulated vs. observed discharges (Fig. 8). Results for the Mississippi (Fig. 8a) show a good agreement but with overestimation in most years. While the RMSE is 19 % for the Mississippi, there is no consistent trend between residue and discharge, indicating no systematic error (Fig. 8b). The RMSE for the discharge of the Rhine is 14 %, with a consistent underestimation by the model (Fig. 8c), and the residues show a clear increase with observed discharge (Fig. 8d), indicating a systematic error in the model. For the Meuse, the RMSE for the discharge is 23 %, the discharge seems to be overestimated (Fig. 8e), and there is only a small trend between discharge and residue (Fig. 8f).

Overall, while discharge is overestimated in the Mississippi, N and P concentrations are underestimated in most years, indicating that part of the problem is in the hydrology. The hydrology model consistently underestimates discharge, while N concentrations are underestimated in most years, and P is underestimated in the first period up till about 1980, and after this year there is a slight overestimation. So apparently errors in the hydrology cannot explain those in the nutrient concentrations. The discharge of the Meuse is overestimated; simulated P concentrations are in good agreement with observations, while N concentrations are overestimated; hence, there is no clear connection between the model errors in discharge and nutrients.

Comparison of modeled (black line) and measured (light blue, and aggregated yearly) discharge-weighed concentrations of total P in the rivers Mississippi (a–c), Rhine (d–f) and Meuse (g–i). Panels on the left are comparisons over time; panels in the center represent plots of simulations vs. observations with a 1 : 1 line, and panels on the right are the concentrations vs. the residues (observation minus simulation) with a regression line.

[Figure omitted. See PDF]

We also investigated the model performance for 10 more stations in various states within the Mississippi River basin (Table 4). These stations along with the St. Francisville station form the monitoring network for nine subbasins in the Mississippi (US Geological Survey, 2007). The plotted data for all 11 stations in Mississippi River basin are available as separate graphs in the Supplement. The model performance is acceptable (RMSE < 50 %) for eight stations for N concentrations and five stations for P concentrations. There are some stations where the model poorly simulates the N concentrations such as Arkansas River and Red River (Table 4). Such high RMSE values do not occur for P. In general, simulated P concentrations are closer to observed values than N concentrations.

One of the reasons for poor agreement is the large fluctuation of discharge, load and concentration at some stations. Apparently, these peaks are associated with periods of high rainfall. We do not know if these peak values represent the full period of the measurement interval. For example, a peak value that represents 2 months (in the case there are six measurements per year) also yields a peak in the aggregated annual value. However, it is not known if this peak actually represents 1 day (with a much lower aggregated annual value) or 2 months. In contrast to St. Francisville, P concentrations (and N concentrations) at the other stations are not consistently underestimated or overestimated. Furthermore, at this level of comparison, the spatial data for land use and wastewater discharge locations in urban areas may not be realistic. For example, our wastewater discharge occurs in all grid cells with urban population, while in reality discharge takes place in discrete locations such as wastewater treatment plants.

Comparison of simulated and observed annual discharge (left-hand graphs with 1 : 1 lines) and residues (observation minus simulation) vs. observation (right-hand graphs with regression lines) for Mississippi (a, b), Rhine (c, d) and Meuse (e, f).

[Figure omitted. See PDF]

A further performance test involves a direct comparison between aggregated data and model results for a large number of European rivers (see Supplement) (European Environment Agency, 2013). This data set includes monitoring data at different stations for 125 rivers, 49 for N and 76 for P. River basins with less than four grid cells, of $\sim$ 2500 km${}^{\mathrm{2}}$ each, were removed because river basin areas of < 10 000 km${}^{\mathrm{2}}$ do not have adequate spatial data representation. This is an arbitrary choice, and probably many river basins with 4–10 grid cells also suffer the problem of poor spatial data. Measurements for some stations were removed from the data set as outliers (Table S1). Results for all measurements show a coefficient of determination of 0.59 and RMSE of 124 % for N ($n=709$) and 0.58 and RMSE of 184 % for P ($n=1010$) (Fig. 9a and b). Results show reasonable coefficients of determination ($r^{\mathrm{2}}$) of 0.79 and RMSE of 112 % for P and 0.55 and RMSE of 95 % for N (Fig. 9c and d). The average of all measurements for N and P is slightly lower than the simulated concentrations (0.16 vs. 0.25 mg P L${}^{-\mathrm{1}}$ and 1.25 vs. 1.78 mg N L${}^{-\mathrm{1}}$). The mean of observations and model values over the monitoring period for each station showed good agreement (Fig. 9e and f). There is also good agreement between model and data for the mean for all stations for each year with deviations never exceeding 1 mg N L${}^{-\mathrm{1}}$ and 0.2 mg P L${}^{-\mathrm{1}}$ (Fig. 9e and f). It is clear that the model has problems when modeling individual stations in small rivers in the database. The plotted data for all stations in the European rivers (available as separate graphs in the Supplement) show that the model results for the Danube, for example, are in good agreement with observations for two stations. Most simulated concentrations are within a factor of 2 of the observed concentrations in the EEA database (EEA, 2012).

Our model results also show a fair agreement with the validation data set for the early 1990s for total N collected by Van Drecht et al. (2003) (Fig. 10). Modeled total N concentrations for the Amazon for the early to mid-1980s (0.7–0.9 mg L${}^{-\mathrm{1}}$) are close to measured values (0.4–0.5), and results for total P (0.07 mg L${}^{-\mathrm{1}}$) are also close to observations (0.06 mg L${}^{-\mathrm{1}}$) (Forsberg et al., 1988; Meybeck and Ragu, 1995).

These comparisons of our model output with data at various aggregation levels show that IMAGE–GNM based on three calibrated submodels (hydrology, nutrient input and in-stream removal) performs very well without any tuning of the overall, integrated model. We have deliberately chosen to not further tune the model so that we can identify its shortcomings. Further improvement of model performance requires a sensitivity analysis.

Model sensitivity

The influence of a range of parameters on model sensitivity was investigated for modeled N and P delivery, retention and river export. Here we discuss only those parameters that are significant and have an SRC value > 0.2 or < $-$0.2 (parameters that add > 4 % to the delivery, retention or river export). Results presented in Tables 5 and 6 show that the sensitivity of N delivery, retention and river export for the year 2000 differs from that of P in many aspects.

RMSE for simulated vs. measured N concentrations, N load, discharge, P concentration and P load for 11 stations in the Mississippi River, Ohio River, Red River, Missouri River and Arkansas River. Measurement frequency ranges from 28 per year to 3. Years with less than 6 observations were excluded.

Station id	Name	RMSE (%)
		Discharge	N	N	P	P
			concentration	load	concentration	load
5420500	Mississippi River at Clinton, IA	60	36	72	23	66
3612500	Ohio River at dam 53 near Grand Chain, IL	32	19	44	48	53
5587550	Mississippi River below Alton, IL	56	48	47	53	71
7355500	Red River near Alexandria, LA	18	119	152	69	72
7022000	Mississippi River at Thebes, IL	67	49	34	64	52
5587455	Mississippi River below Grafton, IL	51	46	27	44	26
3303280	Ohio River at Cannelton dam, KY	56	10	59	58	89
6610000	Missouri River at Omaha, NE	35	74	76	88	78
6934500	Missouri River at Hermann, MO	19	53	56	73	82
7263620	Arkansas River at David D. Terry L&D BL Little Rock, AR	53	244	369	52	92
7373420	Mississippi River near St. Francisville, LA	19	23	26	51	44

Comparison of simulated total N and P concentration with the EEA data set for the period 1970–2000 (EEA, 2013). (a) N concentration for all stations, rivers and years; (b) P concentration for all stations, rivers and years; (c) mean N concentration of all years per station; (d) mean P concentration of all years per station; (e) mean N concentration of all rivers per year; (f) mean P concentration of all rivers per year. Please note that the European coverage is not constant and the trend is not representative of European rivers, because the number and location of stations has changed in time, causing changes in the trend. The 1 : 1 lines are also shown in panels (a–d). Comparison of modeled and observed concentrations for all individual EEA stations is in the Supplement.

[Figure omitted. See PDF]

Total runoff ($q_{\mathrm{tot}}$; Eq. 5) is significant for retention and river export of both N and P; runoff largely determines all transport pathways and flows of N (runoff, leaching, groundwater flow and also in-stream retention), and it determines P runoff, the major transport pathway for P. The soil N budget in natural ecosystems and arable land (N${}_{\mathrm{budget},\mathrm{crops}}$; N${}_{\mathrm{budget},\mathrm{nat}}$; Eq. 1) are important factors for the N delivery, but not for the retention and river export. For P the soil budgets are less important, because soil P content (P${}_{\mathrm{soil}}$) and bulk density ($B_{\mathrm{soil}}$) govern the runoff of P more than the budget; actually, soil P content is actually a result of the long-term soil P budget.

Our model results suggest that allochthonous organic matter input to stream is an important but uncertain nutrient source. The factors determining the allochthonous organic matter input of N to streams and rivers (flooded area, $A_{\mathrm{flooding}}$; litterfall, AOMI; its reduction factor for litterfall, $F_{\mathrm{AOMI}}$; and its C : N ratio) are similarly important for the delivery and river export of N. For P both the parameters determining allochthonous inputs and weathering ($C_{\mathrm{PWeath}}$; Eq. 7) are not significant nor important, as the biomass from litterfall contains only small amounts of P and because the anthropogenic sources are dominant.

For the modeling of river retention, the sensitivity analysis for a range of parameters shows that net uptake velocity ($V_{f,\mathrm{river},\mathrm{N}}$; $V_{f,\mathrm{river},\mathrm{P}}$; Eqs. 23, 26, 28) and mean length ratio ($R_{\mathrm{L}}$; Eq. 29) are important for retention and river export for both N and P, and logically not for nutrient delivery. The assumption that N retention depends on N concentrations (N${}_{\mathrm{conc},\mathrm{low}}$; Eq. 26) is significant in all years for the retention and river export. Temperature (Temp; Eq. 27) is important for retention of P, and for retention and river export of N.

Comparison of simulated total N concentrations for the year 1990 with the validation data set for the early 1990s for total N collected by Van Drecht et al. (2003) with a 1 : 1 line.

[Figure omitted. See PDF]

Results of the sensitivity analysis differ from previous studies (e.g., Bouwman et al., 2013a), primarily because the current model includes additional sources (allochtonous inputs) and changes in the model for surface runoff and leaching.

Future improvements

On the basis of the comparison with measurements and the model sensitivity, we can now analyze what parts of the model need improvement. Improvements are possible in both data and model components. Many components and data are ignored in this discussion, including all the data stemming from the IMAGE on soils, lithology, land use, vegetation distribution, nutrient cycles in agriculture and natural ecosystems and climate. We recognize that updates of the data used in this paper are now available. For example, soil data (http://www.isric.org/content/soilgrids), hydrographic information (http://hydrosheds.cr.usgs.gov/index.php) and lithology (Hartmann and Moosdorf, 2012) and associated porosity and permeability data (Gleeson et al., 2014). These updates have a finer resolution, allowing more specific calculation of surface characteristics (bare rock, more detailed soil texture classes, etc.). Hence, these updates and additional data sets will be considered for future improved versions of the model, and tested with new sensitivity analyses.

Standardized regression coefficient (SRC)${}^{\ast }$ representing the relative sensitivity of N delivery, N retention and river N export representing global model results (columns) for the year 2000 to variation in 48 parameters.

Parameter	N delivery	N retention	N export
${}_{\mathrm{qtot}}$	0.24	$\mathbf{-}\mathbf{0.23}$	0.28
$D_{\mathrm{rip}}$	$-$0.02	0.01	$-$0.02
N${}_{\mathrm{budget},\mathrm{crops}}$	0.26	$-$0.06	0.16
N${}_{\mathrm{budget},\mathrm{grass}}$	0.05		0.02
N${}_{\mathrm{budget},\mathrm{nat}}$	0.20	$-$0.02	0.10
N${}_{\mathrm{uptake},\mathrm{crops}}$	0.06		0.03
N${}_{\mathrm{uptake},\mathrm{grass}}$	0.03		0.01
$B_{\mathrm{soil}}$
CN${}_{\mathrm{soil},\mathrm{crop}}$	$-$0.13		$-$0.06
CN${}_{\mathrm{soil},\mathrm{grass}}$	$-$0.03		$-$0.01
CN${}_{\mathrm{soil},\mathrm{nat}}$	$-$0.04		$-$0.02
$C_{\mathrm{sro}}$	0.18	$-$0.01	0.09
$f_{\mathrm{qgwb}}$	$-$0.09	0.02	$-$0.06
$f_{\mathrm{qsro}}$	0.15	$-$0.01	0.07
$f_{\mathrm{qsro}}$(crops)	0.11	$-$0.01	0.06
$f_{\mathrm{qsro}}$(grass)	0.16		0.07
$f_{\mathrm{qsro}}$(nat)	0.07		0.03
$F_{\mathrm{leach},\mathrm{crop}}$	0.10	$-$0.02	0.06
$F_{\mathrm{leach},\mathrm{grass}}$	0.04	$-$0.01	0.03
$F_{\mathrm{leach},\mathrm{nat}}$	0.19	$-$0.02	0.10
$D_{\mathrm{dgrw}}$	$-$0.02	0.01	$-$0.02
$D_{\mathrm{sgrw}}$	$-$0.13	0.01	$-$0.07
dt50${}_{\mathrm{den},\mathrm{dgrw}}$	0.02
dt50${}_{\mathrm{den},\mathrm{sgrw}}$	0.14	$-$0.01	0.07
Poros	$-$0.15	0.01	$-$0.08
$A_{\mathrm{flooding}}$	0.34	$-$0.11	0.23
AOMI	0.35	$-$0.10	0.24
CN${}_{\mathrm{aomi}}$	$\mathbf{-}\mathbf{0.35}$	0.10	$\mathbf{-}\mathbf{0.24}$
$F_{\mathrm{aomi}}$	0.35	$-$0.10	0.24
$A$		0.16	$-$0.12
$A_{\mathrm{1}}$		$-$0.04	0.03
$B$		0.09	$-$0.07
$D_{\mathrm{flooding}}$		$-$0.01	0.01
$L_{\mathrm{1}}$		0.21	$-$0.16
N${}_{\mathrm{conc},\mathrm{high}}$		0.16	$-$0.12
N${}_{\mathrm{conc},\mathrm{low}}$	-0.01	0.40	$\mathbf{-}\mathbf{0.31}$
$R_{\mathrm{a}}$		$-$0.08	0.06
$R_b$		0.08	$-$0.06
$R_{\mathrm{L}}$		0.53	$\mathbf{-}\mathbf{0.41}$
Temp	$-$0.09	0.41	$\mathbf{-}\mathbf{0.36}$
$v_{f,\mathrm{lake},\mathrm{N}}$		0.06	$-$0.04
$v_{f,\mathrm{reservoir},\mathrm{N}}$		0.07	$-$0.05
$v_{f,\mathrm{river},\mathrm{N}}$		0.38	$\mathbf{-}\mathbf{0.30}$
$v_{f,\mathrm{wetland},\mathrm{N}}$
$V_{\mathrm{water}}$		0.01
N${}_{\mathrm{aqua}}$	0.03	$-$0.01	0.02
N${}_{\mathrm{depo}}$	0.03	0.01
N${}_{\mathrm{point}}$	0.22	$-$0.06	0.14

${}^{\ast }$ Cells with no values represent insignificant SRC values; all cells with values have significant SRC, numbers with normal font indicate values $-$0.2 < SRC < 0.2; numbers with bold and italic font indicate values exceeding $+$0.2 and $-$0.2, respectively. An SRC value of 0.2 indicates that the parameter concerned has an influence of 0.2${}^{\mathrm{2}}=$ 0.04 (4 %) on the model variable considered.

Standardized regression coefficient (SRC)${}^{\ast }$ representing the relative sensitivity of P delivery, P retention and river P export representing global model results (columns) for the year 2000 to variation in 34 parameters.

Parameter	P delivery	P retention	P export
${}_{\mathrm{qtot}}$	0.17	$\mathbf{-}\mathbf{0.47}$	0.48
P${}_{\mathrm{budget},\mathrm{crops}}$	0.07		0.05
P${}_{\mathrm{budget},\mathrm{grass}}$
P${}_{\mathrm{budget},\mathrm{nat}}$
P${}_{\mathrm{uptake},\mathrm{crops}}$	0.06		0.04
P${}_{\mathrm{uptake},\mathrm{grass}}$	0.02		0.01
$B_{\mathrm{soil}}$	$\mathbf{-}\mathbf{0.62}$	$-$0.13	$\mathbf{-}\mathbf{0.36}$
$C_{\mathrm{sro}}$	0.13		0.10
$f_{\mathrm{qsro}}$	0.13		0.10
P${}_{\mathrm{soil}}$	0.63	0.13	0.36
$F_{\mathrm{leach},\mathrm{crop}}$
$F_{\mathrm{leach},\mathrm{grass}}$
$F_{\mathrm{leach},\mathrm{nat}}$
P${}_{\mathrm{weathering}}$	0.17	$-$0.04	0.15
$A_{\mathrm{flooding}}$	0.13	$-$0.02	0.11
AOMI	0.14	$-$0.02	0.12
CP${}_{\mathrm{aomi}}$	$-$0.14	0.02	$-$0.11
$F_{\mathrm{aomi}}$	0.14	$-$0.02	0.12
$A$		0.22	$-$0.17
$A_{\mathrm{1}}$		$-$0.13	0.10
$B$			0.01
$D_{\mathrm{flooding}}$		$-$0.01
$L_{\mathrm{1}}$		0.28	$\mathbf{-}\mathbf{0.22}$
$R_{\mathrm{a}}$		$\mathbf{-}\mathbf{0.24}$	0.19
$R_b$		0.16	$-$0.12
$R_{\mathrm{L}}$		0.49	$\mathbf{-}\mathbf{0.38}$
Temp	0.12	0.27	$-$0.12
$v_{f,\mathrm{lake},\mathrm{P}}$		0.06	$-$0.04
$v_{f,\mathrm{reservoir},\mathrm{P}}$		0.10	$-$0.08
$v_{f,\mathrm{river},\mathrm{P}}$		0.40	$\mathbf{-}\mathbf{0.30}$
$v_{f,\mathrm{wetland},\mathrm{P}}$
$V_{\mathrm{water}}$		0.01
P${}_{\mathrm{aqua}}$	0.01		0.02
P${}_{\mathrm{point}}$	0.14	$-$0.06	0.15

${}^{\ast }$ See Table 5.

It is difficult to know from the available analyses what could be done to improve the model, because error may be the result of uncertainties in the input data (land use, climate, hydrology, wastewater flows, etc.), in surface and subsurface processes or in-stream processes. However, two parts of the model have a dominant importance for the model results, i.e., total runoff from the water balance model PCR-GLOBWB and the factors determining the in-stream biogeochemistry including the uptake velocity and factors used in the parameterization of sub-grid processes for streams and rivers of Strahler orders of 6 and less. Here we do not touch upon improvements of the hydrology model and focus on the nutrient-related processes, but see a clear need for improvement of the way the water flow in lakes and reservoirs is simulated, i.e., only the water that actually enters and leaves the lake is considered, with no role for the total water mass. Also, there is a need to improve the geohydrological information in order to better describe global aquifers, their thickness and their denitrification potential.

To improve the in-stream process description, the first short-term improvement is to add processes in sediments to allow for simulating P saturation of sediments and desorption in case of decreasing river P loads.

The current model version uses air temperature as a proxy for water temperature. A clear improvement would be to use water temperatures in the spiraling approach, since there may be large differences, especially in low-order streams. Other examples are large rivers influenced by cooling water from nuclear or other power plants. The river Meuse is such an example, and the overestimation of N concentrations may be caused by underestimation of the water temperature.

The importance of factors such as the P content of the soil call for attention to the description of the processes determining P (and N) transport to surface water via surface runoff. Our approach distinguishes an instant transport route, and the transport of soil material with the memory simulated by changing P content of the soil. The delay of the transport may be an important aspect to consider, but at present we have no data available to do so.

Longer-term improvements center on the incorporation of a mechanistic model for describing biogeochemical processes in the water column and sediment. This allows for further analysis of individual processes and their interplay (plant uptake, sedimentation, benthic processes, denitrification). This will involve a change to a temporal resolution that matches the requirements of the description of the biogeochemical processes (day, week, month). Mechanistic modeling of in-stream processes with shorter time steps requires a further refinement of the processes on the land, i.e., the temporal distribution of fertilizer application, manure spreading and grazing. This will also allow us to analyze the delay between rainfall events causing runoff and the discharge to the surface water. Also, such mechanistic models require a delivery and in-stream model that distinguishes different nutrient forms.

Mechanistic modeling also allows for the coupling of the processes of C with the nutrients N, P and Si, which may lead to better understanding of the C and nutrient fluxes to and from river basins. Regarding spatial scale, the current 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ resolution is large enough to assume that there are no interactions between grid cells. Future models at finer resolutions need to consider the fact that transport and processes may cross boundaries of grid cells.

Conclusions

The performance of our global nutrient model is similar to that of the more commonly used empirical approaches. The comparisons of our model output with data at various aggregation levels show that our model based on three submodels (hydrology, nutrient delivery and in-stream retention) performs very well without any calibration. We have deliberately chosen to not further tune the model so that we can identify its shortcomings.

IMAGE–GNM can simulate the present-day river nutrient export at the basin and global scales with acceptable deviations from observed values for large rivers, and generally within a factor of 2 for small European rivers. The model can also be used to explore changes in various processes and interactions between them during the 20th century. More specifically, the IMAGE–GNM model allows for attributing changes in nutrient transport, retention and export to changes in hydrology and nutrient delivery or their interactions (Beusen et al., 2015). It will therefore be a very valuable research tools to examine the effect of hydrological measures or climate-induced changes on nutrient processing and export and therefore on the functioning of downstream ecosystems.

Moreover, GNM is fully integrated into the integrated assessment model IMAGE and can thus provide nutrient transport and processing estimates fully consistent with scenarios based on, for example, the story lines of the shared socio-economic pathways currently in use by the global climate change community (Kriegler et al., 2014).

An interesting application of IMAGE–GNM is to study the impacts of increasing river export, i.e., eutrophication of coastal marine ecosystems leading to phenomena such as increased production and hypoxia. The changing nutrient stoichiometry in freshwater and coastal systems may lead to phenomena such as harmful algal blooms. Such analyses require coupling our model to coastal biogeochemistry models.

The Supplement related to this article is available online at doi:10.5194/gmd-8-4045-2015-supplement.

A. H. W. Beusen and A. F. Bouwman developed the model for the delivery of nutrients to surface water on the basis of the work presented in Van Drecht et al. (2003). A. H. W. Beusen, J. M. Mogollón and L. P. H. Van Beek integrated IMAGE–GNM and PCR-GLOBWB and further developed the routing, J. J. Middelburg, J. M. Mogollón, A. F. Bouwman and L. P. H. Van Beek wrote the text.

Acknowledgements

This paper was supported by the Water, Climate and Ecosystems project, part of the Sustainability strategic theme of Utrecht University (http://wce.uu.nl/), and contributes to the Netherlands Earth System Science Centre (NESSC, http://www.nessc.nl/). We gratefully acknowledge financial support from the Global Environment Facility (GEF), United Nations Environment Programme (UNEP), Intergovernmental Oceanographic Commission of the UNESCO (IOC/UNESCO) and other partners through the UNEP/GEF project Global Foundations for Reducing Nutrient Enrichment and Oxygen Depletion from Land-based Pollution in Support of Global Nutrient Cycle (GNC project). Additional funding was provided by the EU H2020 (MSCA award 661163 to J. M. Mogollón). Edited by: S. Arndt