Human impact parameterizations in global hydrological models improve estimates of monthly discharges and hydrological extremes: a multi-model validation study

We used a modeled monthly discharge (0.5° × 0.5° spatial resolution) for the period 1971–2010 from five GHMs: H08 (Hanasaki et al 2008a, 2008b), LPJmL (Bondeau et al 2007, Rost et al 2008, Schaphoff et al 2013), MATSIRO (Pokhrel et al 2012, 2015, Takata et al 2003), PCR-GLOBWB (van Beek et al 2011, Wada et al 2011, 2014b) and WaterGAP2 (Müller Schmied et al 2016). All simulations were carried out under the modeling framework of phase 2a of the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP2a: www.isimip.org/protocol/#isimip2a). For each GHM, we used two simulations: (1) HIP: a model run including time-varying land use and land cover change, historical dam construction and operation, irrigation and upstream consumptive water abstractions; and (2) NOHIP: a 'naturalized' model run without HIP.

An overview of the model characteristics of each of the GHMs, and the methods used to parameterize hydrological processes and human impact, can be found in supplementary table 1, and details on each GHM can be found in the individual model references provided therein. In the following subsections, we briefly outline the most important characteristics of the hydrological and human impact parameterizations.

2.1.1. Parameterizations of hydrological processes

2.1.2. Parameterization of human impact

Observed monthly river discharge data was taken from the Global Runoff Data Centre (GRDC, 56068 Koblenz, Germany). From the 9051 gauging stations in the GRDC database, we selected stations that meet the following criteria: (1) a minimum of 5 years' coverage (not necessarily consecutive) during the period 1971–2010 with a completeness of observations of ≥95%; and (2) a minimum catchment area of 9000 km², to omit catchments whose hydrological processes cannot be adequately represented by models operating at 0.5° × 0.5° (Hunger and Döll 2008). We discarded the stations for which the difference in the catchment area in the GRDC database and as estimated by using the DDM30 river routing network (Döll and Lehner 2002) is >25%.

Zoom In Zoom Out Reset image size

We then made a distinction between near-natural and managed catchments. Following Beck et al (2016), a catchment is classified as near-natural if the share of land-area subject to irrigation is <2% and the total reservoir capacity is <10% of its long-term mean annual discharge. If these conditions were not met the catchment was classified as managed. The classification was based on the HYDE 3/MIRCA land cover dataset (Fader et al 2010, Klein Goldewijk and Van Drecht 2006, Portmann et al 2010, Ramankutty et al 2008) together with the Global Reservoir and Dam database (Lehner et al 2011). Two stations shifted from near-natural to human-impacted conditions between 1971 and 2010 and were discarded from further analysis.

The aforementioned steps resulted in 471 stations with a total catchment area covering 19.8% of the global land (figure 2), of which 92 are located at the outlet of a catchment area. The mean length of observations is 32.8 years for all stations. Of all the stations, 226 are located in managed and 245 in near-natural catchments. Of the stations located at the outlet of a catchment, 45 are managed (4.8% of the global land area) and 47 are near-natural (15.1% of the global land area).

Figure 2 shows that the majority of selected stations (blue) are located in northern and Latin America, Europe, southern Africa and Australia. The number of stations in northern and central Africa and Asia is relatively small. We selected 12 stations in river basins located in different geographic regions (green circles in figure 2: Amazonas, Amur, Colorado, Congo, Guadiana, Mackenzie, Murray, Ob, Rhine, Tocantins, Volga and the Zambezi) for which a detailed analysis is provided in the supplementary results section (supplementary).

To evaluate the GHM simulation of monthly discharge and hydrological extremes under HIP and NOHIP conditions, we compared the modeled results with observed river discharge data using several evaluation metrics described below. To ensure a consistent comparison between the modeled and observed data, we only used modeled data for the same years for which observations were available. We also corrected modeled discharges for potential over- and underestimations caused by the difference in catchment size between the model and the GRDC. To do this, we used a multiplier that represented the difference in the upstream area as reported by the GRDC and as estimated from the DDM30 network.

First, we applied the modified Kling–Gupta efficiency index (KGE) with its sub-components: the linear correlation coefficient (rKGE), the bias ratio (βKGE) and the variability ratio (γKGE) (Gupta et al 2009, Kling et al 2012). The KGE is a widely applied indicator for the validation of hydrological performance in modeling studies at the global and regional scale and provides a good representation of the 'closeness' of simulated discharges to observations (Huang et al 2017, Kuentz et al 2013, Nicolle et al 2014, Revilla-Romero et al 2015, Thiemig et al 2013, 2015, Thirel et al 2015, Wöhling et al 2013). Moreover, use of its three sub-components enables the identification of reasons for sub-optimal model performance (Gupta et al 2009, Kling et al 2012, Thiemig et al 2013). This was achieved by estimating for each sub-parameter its distance to optimal performance, and by subsequently comparing these distances across the different sub-parameters. The statistical significance of the change in KGE outcomes due to the inclusion of HIP was tested by means of regular bootstrapping (n = 1000, p ≤ 0.05 (two-tailed)), following the method of Livezey and Chen (1982) and Wilks (2006).

Second, we applied the Nash–Sutcliffe efficiency test (NSE, Nash and Sutcliffe 1970) to evaluate the representation of Q₁ (high-flow) and Q₉₉ (low-flow) conditions (e.g. Beck et al 2017a, Blösch et al 2013, Hejazi and Moglen 2008, Mohamoud 2008), obtained under fixed threshold level settings (van Loon 2015). By means of a two-sample Kolmogorov–Smirnov (KS) test (Massey 1951, p ≤ 0.05) we tested how often HIP leads to significant changes in the fit of the full modeled exceedance probability curve for hydrological extremes compared to the full observed exceedance probability curve.

Including the parameterization of human impacts in the GHMs leads to a large improvement in overall model performance. The hydrological performance under the HIP simulations shows a significant improvement compared to the NOHIP simulations for between 40.8% and 72.3% of the land area studied, depending on the GHM (figure 3(a)). For most GHMs, the positive effects of including HIP in the simulations outweigh the negative effects. This is the case for both near-natural and managed catchments, although the positive effects are more pronounced for the managed catchments (figures 3(a)–(d)). Near-natural catchments are only indirectly impacted by HIP, for example by receiving improved or altered water simulations from upstream managed catchments. The KGE sub-components show significant improvement in performance in large shares of the land area studied, especially for the bias and variability ratio. The bias ratio improves significantly for 36.1%–73.0% of the total land area for all catchments, compared to 64.8%–90.6% and 24.3%–70.4% in managed and near-natural catchments respectively (figure 3(b)). For the variability ratio, improvements were found for 31.4%–74.4% of the land area for all catchments (48.9%–92.6% for the managed, 23.0%–73.2% for the near-natural) (figure 3(c)). The lowest improvements are found for the correlation coefficient, with improvements for 15.9%–58.1% of the total land area for all catchments (22.1%–75.1% for the managed, 13.9%–61.4% for the near-natural) (figure 3(d)).

Zoom In Zoom Out Reset image size

The results are shown for each station in figure 4 for the overall model performance (KGE), and in supplementary figure 1 for the KGE sub-parameters. The results show particularly strong improvements in overall performance in Latin America, southern Africa and the northwest US. There are only a limited number of stations for which the inclusion of HIP leads to a significant decrease in overall hydrological performance for the majority of GHMs or where no to limited changes occur, for example in near-natural areas (e.g. the Amazonas).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When considering overall hydrological performance for each GHM under HIP conditions (figure 3(e)), WaterGAP2 and MATSIRO show the best performance globally. Even though the simulations with HIP include human impact parameterizations by definition, all GHMs still show better performance in near-natural catchments than in managed catchments (figures 3(e)–(h)). The KGE bias ratio values >1 indicate that all models systematically overestimate long-term mean monthly discharge (figure 3(f)), up to five-fold for LPJmL in managed catchments. For the variability ratio (figure 3(g)), WaterGAP2 is the only GHM that tends to slightly underestimate variability (variability ratio <1) in the monthly discharge, in both the managed and near-natural catchments. All other GHMs show overestimations, up to 1.55 fold for LPJmL for near-natural catchments. All GHMs show a reasonable correlation with the observed monthly discharge estimates (figure 3(h)), with values ranging between 0.49–0.69 in the managed catchments and 0.50–0.79 in the near-natural catchments. The highest correlation coefficients including HIP are found for WaterGAP2, with a global mean value across all catchments of 0.76 (0.69 for the managed catchments and 0.78 for the near-natural catchments).

For each catchment (and therefore its associated land area), it is possible to distinguish which of the KGE sub-parameters contributes most to sub-optimal performance. These results are summarized in figure 5. The results show that under HIP conditions, the bias ratio contributes most to sub-optimal performance in managed catchments for most GHMs, except WaterGAP2 (for which the correlation coefficient contributes most). For near-natural catchments, sub-optimal performance is most often caused by the variability ratio for H08, LPJmL and WaterGAP2, by the bias ratio for MATSIRO and by the correlation coefficient for PCR-GLOBWB.

Spatially explicit results vary per GHM and are shown in supplementary figure 3. The distribution of dominant contributors to the sub-optimal overall hydrological performance is similar for H08, LPJmL, and PCR-GLOBWB. For these GHMs, we find dominant contributions from the bias ratio in southern Africa, Australia and inland US, and dominant contributions of the variability ratio and the correlation coefficient in Latin America as well as at higher latitude and altitude regions. For Europe, the dominant contributions for H08, LPJmL and PCR-GLOBWB are the variability ratio, the correlation coefficient, and the bias ratio respectively. The dominant contributors that cause sub-optimal overall hydrological performance for MATSIRO and WaterGAP2 are more equally distributed across the globe. While all sub-components contribute to sub-optimal overall hydrological model performance for MATSIRO, it is predominantly the correlation coefficient and the variability ratio that determines the sub-optimal performance in WaterGAP2.

The inclusion of HIP in the simulations affects the ability of GHMs to estimate hydrological extremes correctly in the majority of the land area studied (figure 6). The inclusion of HIP leads to better model performance for all GHMs, across a substantial share of the land area studied (figures 6(a)–(b)). For high-flows, HIP improves model performance significantly across 34.6%–77.0% of the land area for all catchments (36.4%–94.7% for managed, 24.1%–79.2% for near-natural). For low-flows, HIP improves model performance significantly across 39.4%–80.4% of the land area for all catchments (29.3%–81.8% for managed, 42.7%–90.3% for near-natural). The KS-test results (supplementary figure 4) show that HIP only leads to significant changes in the representation of the exceedance probability curve in a limited number of cases for H08 and LPJmL (up to 14.1% of the land area studied), predominantly in managed catchments.

Table 1. The performance metrics used in this study and their calculation procedure. Here, s_i and o_i are the simulated and observed monthly discharge at station i; μ_s and μ_o are the simulated and observed mean monthly discharge at station i; σ_s and σ_o are the standard deviation of the simulated and observed discharge at station i, respectively; Q_s and Q_o are the simulated and observed hydrological extremes.

Abbreviation	Name	Calculation procedure	Range and ideal value
KGE	Modified Kling–Gupta efficiency index	${\rm KGE} = 1 - \sqrt {(r{\rm KGE}^{*}-1)^{2} + (\beta {\rm KGE}^{*} - 1)^{2} + (\gamma {\rm KGE}^{*} - 1)^{2}}$	−∞ – 1 (ideal value: 1)
rKGE	KGE correlation coefficient (Pearson)	$r{\rm{KGE}} = \frac{{\mathop \sum \nolimits_{i = 1}^n \left( {{s_i} - {\mu _{s,i}}} \right)\left( {{o_i} - {\mu _{o,i}}} \right)}}{{\sqrt {\mathop \sum \nolimits_{i = 1}^n {{\left( {{s_i} - {\mu _{s,i}}} \right)}^2}} \sqrt {\mathop \sum \limits_{i = 1}^n {{\left( {{o_i} - {\mu _{o,i}}} \right)}^2}} }}$	−1 – 1 (ideal value: 1)
βKGE	KGE bias ratio	βKGE = μs,i/μo,i	0 – ∞ (ideal value: 1)
γKGE	KGE variability ratio	$\gamma {\rm{KGE}} = \frac{{{\raise0.7ex\hbox{${{\sigma _{s,i}}}$} \!\mathord{\left/ {\vphantom {{{\sigma _{s,i}}} {{\mu _{s,i}}}}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{${{\mu _{s,i}}}$}}}}{{{\raise0.7ex\hbox{${{\sigma _{o,i}}}$} \!\mathord{\left/ {\vphantom {{{\sigma _{o,i}}} {{\mu _{o,i}}}}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{${{\mu _{o,i}}}$}}}}$	0 – ∞ (ideal value: 1)
NSE	Nash–Sutcliffe model efficiency	${\rm{NSE}} = 1 - \frac{{\sum {{({Q_s} - {Q_o})}^2}}}{{\sum {{({Q_o} - \overline {{Q_o}} )}^2}}}$	−∞ – 1 (ideal value: 1)
Q1	High-flow indicator	Monthly discharge (m s−3) that is exceeded on average in 1 out of 100 months
Q99	Low-flow indicator	Monthly discharge (m s−3) that is exceeded on average in 99 out of 100 months
KS	Two-sample Kolmogorov–Smirnov test	[h, p] = kstest2 (cdf (Qs,), cdf (Qo), 'Alpha', 0.05)a	For p > 0.05, H0 (the two cdfs coming from the same distribution) is not rejected.

^aThe calculation procedure for the two-sample Kolmogorov–Smirnov test presented in the table is the Matlab function for the KS-test.

Zoom In Zoom Out Reset image size

Overall, hydrological extremes are represented reasonably well under HIP conditions, with globally weighted mean NSE values ranging between 0.8–0.98 for high-flows, and 0.84–0.98 for low-flows (figures 6(c)–(d)). However, there is a significant difference in the ability of the GHMs to represent hydrological extremes between managed and near-natural catchments.

Figure 7 indicates that for the majority of stations, the inclusion of HIP leads to an improvement in the representation of hydrological extremes, for most GHMs. A deterioration in the representation of hydrological extremes across the majority of GHMs as a result of the inclusion of HIP was only found in selected areas, for example at higher latitudes and along the east coast of the US. When comparing the results for Q₁ high-flows with Q₉₉ low-flows, no large differences in the spatial distribution of the number of GHMs are found with a significant improvement or deterioration.

Zoom In Zoom Out Reset image size

The effects of HIP on the magnitude of extreme discharge differ for low-flows and high-flows (supplementary figure 5). Whilst the magnitude of high-flows mostly decreases with the inclusion of HIP, the effects on the magnitude of low-flows are both positive and negative. The convergence of results towards higher observed discharges, in both high- and low-flow estimates (as identified for all models in supplementary figure 5), indicates that HIP becomes less important for the correct representation of hydrological extremes with increasing discharge volumes.

Whilst the inclusion of HIP predominantly leads to the largest improvements in simulated discharge in the managed catchments, simulated discharge is also improved in a large share of the near-natural catchments. Improvements in model performance associated with the inclusion of HIP can be attributed to improvements in the different KGE sub-components, and in turn to different model components parameterizing hydrological and human processes. In addition, insights into those factors bounding the optimal hydrological model performance under HIP conditions may help to identify priorities for further model improvement.

4.1.1. Representation of long-term mean discharges (bias ratio)

4.1.2. Representation of hydrological variability (variability ratio)

4.1.3. Representation of the goodness-of-fit (correlation coefficient)

4.1.4. Representation of hydrological extremes

As the GHMs have very different parameterizations of hydrological and human processes, the current study does not allow a systematic assessment of specific cause–effect relations between HIP and the observed improvements in performance (Döll et al 2016, Haddeland et al 2014, Hagemann et al 2013, Schewe et al 2014, Beck et al 2016). To do this, a substantial Monte Carlo analysis would be required, whereby individual parameters and combinations of parameters are systematically modified for all GHMs (Döll et al 2016). Undertaking such an analysis in parallel to the different GHMs incorporated is computationally expensive and requires a strict modeling protocol. It may, however, provide additional information on how to adapt and improve the individual GHMs and would be a valuable addition to the results presented in this study.

When interpreting the results of this study, one must take into account that we only evaluated the GHMs with respect to the monthly discharge. Whilst the monthly discharge may be sufficient for the assessment and management of low-flows, droughts and freshwater resource availability, flood risk assessment and management require information on the daily peak discharge. Further research should therefore attempt to validate GHMs using daily peak discharge and assess how it is affected by the inclusion of HIP.

The spatial resolution of the GHMs applied in this study is 0.5° × 0.5° (~50 km × 50 km at the equator), dictated by the resolution of the GSWP3 input dataset. At a 0.5° × 0.5° spatial resolution, hydrological processes are often represented by GHMs in a simplified or generalized form which are not fit for local applications (Bierkens 2015). To account for this, we applied a minimum catchment size of 9000 km², thereby omitting catchments that are too small to be adequately represented by GHMs (Hunger and Döll 2008). Newer versions of several of the GHMs now operate at higher resolutions; for example WaterGAP and PCR-GLOBWB have recently been published in 5 m/6 m versions respectively (Verzano et al 2012, Wada et al 2016b). Future research can investigate whether the inclusion of these high-resolution model-runs improves the representation of discharges and hydrological extremes in the selected catchments and whether these high-resolution runs also allow for the inclusion of smaller catchments.

In this study, a relatively simple distinction was made between managed and near-natural catchments using two parameters: irrigated agriculture and reservoirs. These parameters were chosen as they have been reported to be the most significant human parameters on river hydrology (Beck et al 2016, 2017a). However, to make a more detailed distinction between catchments that are impacted by human activity and those that are not, future studies could consider incorporating additional criteria, such as the share of sectoral water abstractions and return flows, and the share of built-up land area. Additional catchment descriptors (Eisner 2016), like climate conditions and the physiographic properties of the drainage area, could also be applied to further assess the important controls on modeled discharges.

When evaluating the impact of HIP on hydrological extremes, we only incorporated results for the Q₁ high-flow and Q₉₉ low-flow. In this study, we did not consider other ranges of the extreme value distribution explicitly. Although the inclusion of HIP influences these hydrological extremes substantially, we found very few instances for which this led to a significant change in the full exceedance probability curve. Future research should therefore also incorporate other ranges of the probability exceedance curve in order to do a full assessment of the influence of HIP on high- and low-flow extremes.

Next to the parameterizations and representation of hydrological processes and human impacts, other sources contribute to uncertainty in the modeling of discharges and hydrological extremes. These include the quality of the data, the uncertainties in the input data and observation datasets and the calibration/validation strategy (Döll et al 2016, Sood and Smakhtin 2015). The quality of the selected forcing data, for example, may limit the representation of monthly discharges and hydrological extremes significantly (Döll et al 2016, Beck et al 2016), although this not been evaluated explicitly here. However, climate forcing uncertainty is probably a dominant driver for model outputs (Müller Schmied et al 2014, 2016). Benchmarking of the GSWP3 dataset against historical observations of precipitation and temperature, or against other forcing datasets (e.g. similar to Beck et al 2017b, Sun et al 2018), may therefore be of added value.

Differences in the quality and trustworthiness of the historical discharge observations (e.g. due to sampling, measurement and interpretation errors), may potentially result in artificial biases in the validation results (Renard et al 2010). The spatial representativeness of our results is limited by the availability of consistent publicly available in situ observations of sufficient quality. Future research should therefore consider extending the GRDC data-points to the regional repositories of observed discharges, as recently attempted by Beck et al (2016), Do et al (2017) and Gudmundsson et al (2017). However, increasing spatial representation comes at the cost of consistency, and special attention should be paid to the harmonization of these different databases. The use of remotely sensed data could also provide a valuable way of carrying out calibration and validation in ungauged regions (Döll et al 2014a, 2014b, Scanlon et al 2018). Remotely sensed data can also be of added value in the assessment of water consumed by agricultural irrigation (Peña-Arancibia et al 2016), operational drought monitoring and early warning (Ahmadalipour et al 2017), as well as the estimation of terrestrial water budgets (Zhang et al 2017). Moreover, clear potential exists for assimilating remotely sensed data into the models (Eicker et al 2014).

Calibration and validation are essential to compensate for factors such as the impossibility of measuring all required model parameters at the applied scale, the lack of process understanding, the simplistic process representation in GHMs and errors in forcing data (Beck et al 2016, Bierkens 2015, Döll et al 2016, Liu et al 2017). Hence, calibration/validation is key for obtaining realistic model performance. It should be acknowledged, though, that the representation of hydrological and/or human processes is artificially altered by means of calibration/validation processes and that limited calibration may introduce uncertainties to the model output (Sood and Smakhtin 2015). Before using any calibrated/validated model data, one should therefore critically reflect on whether the calibration/validation procedure executed—together with their optimization objectives—are fit for the specific application in mind.