Regionalization of catchment hydrological model parameters for global water resources simulations

Water is one of the most important natural resources that significantly influences the social and economic development of a region or globe (Parajka et al. 2007; Petelet-Giraud et al. 2018; Grill et al. 2019; Behboudian et al. 2021; Naghdi et al. 2021). With the development of society and economy, the conflict of water resources supply and demand has been intensified (Müller Schmied et al. 2016; Vargas Godoy et al. 2021). In addition, the changing climate and increasing water demand due to the growing population lead to a series of environmental and social problems, e.g. frequent floods and droughts, food crisis, and soil salinization (Zhang et al. 2020; Li et al. 2021; Roushangar et al. 2021; Zhou et al. 2022). Efficient water resource allocation could help to address water shortage problems around the world, which should be based on a full understanding of the spatial and temporal variation of global water resources (Shrestha et al. 2017; Zhang et al. 2020). However, due to the lack of observational data around the world, especially in poorly instrumented regions (e.g. the Middle East and West Africa), it is difficult to directly obtain global water information. In addition, problems caused by water conflicts in multinational catchments and virtual water trade both reflect the requirements of continental and global-scale hydrological simulations (Döll et al. 2003; Oki & Kanae 2006; Widén-Nilsson et al. 2007, 2009).

Global water resources are commonly approximated by using runoff coefficient (Liang & Greene 2020), gauge-based observations (Harris et al. 2014; Ghiggi et al. 2019), precipitation minus evaporation (Syed et al. 2010; Chandanpurkar et al. 2017) and water balance modelling, among which water balance modelling is one of the most popular solutions for global water resource estimations (Döll et al. 2003; Oki & Kanae 2006; Widén-Nilsson et al. 2007, 2009; Beck et al. 2016, 2020). Models developed to simulate continental or global water resources can be roughly classified into dynamic global vegetation models (DGVMs), land surface models (LSMs), and global hydrological models (GHMs). Most DGVMs do not include lateral water flows or surface water bodies and can therefore only be used to assess runoff generation but not streamflow discharge (Döll et al. 2015). The LSM is commonly used as a component of climate models in simulating the energy and water balance at soil, atmosphere, and vegetation interfaces (Haddeland et al. 2011; Bierkens 2015). However, global climate models are considerably biased in global runoff simulations (Sellers et al. 1986; Xu et al. 2005; Sood & Smakhtin 2015). Hence, GHMs focusing on the simulation of water resources have been developed to simulate (sub-) surface water fluxes and storages. Some of the widely used GHMs include Variable Infiltration Capacity model (VIC, Liang et al. 1994), Water Balance Model–Water Transport Model (WBM-WTM, Vörösmarty et al. 1989), PCRaster GLOBal Water Balance model (PCR-GLOBWB, Van Beek & Bierkens 2008; Van Beek et al. 2012), and Water And Snow balance MODeling system (WASMOD-M, Widén-Nilsson et al. 2007, 2009).

Although great progress has been made in global hydrological modelling during the past decades, several studies highlighted the considerable discrepancies among these simulation results, which is partly due to the difficulty of parameter estimation for hydrological models (Widén-Nilsson et al. 2007; Schellekens et al. 2016). The lack of observational data in many regions of the world (especially in poorly instrumented regions) and the low quality or ‘disinformative’ data make it difficult to directly calibrate global hydrological model parameters (Kauffeldt et al. 2013; Schellekens et al. 2016). Therefore, the majority of macroscale models (e.g. GHMs) were inclined to get attached to prior parameterizations, which may result in insufficient streamflow simulations (Beck et al. 2016). For example, the parameter values of the WBM-WTM were tuned by an adjustment factor, rather than optimization (Vörösmarty et al. 1989). Thus, the parameter values of WGHM were globally uniform or related to land cover and its associated properties, except the runoff coefficient which was tuned against the time series of measured annual discharges (Döll et al. 2003). The parameter values of Macro-PDM were set based on literature review or previous model applications, and 6 out of 13 parameters were globally uniform (Arnell 1999, 2003).

Considering the restriction of the observational data, parameter regionalization, the most common method to solve the Predictions in Ungauged Basins (PUB) problem, has been used in global hydrological modelling. For example, the WASMOD-M transferred the calibrated parameter sets from catchment scale to grid cells by searching for the most commonly occurring parameter set within a rectangular window and found that regionalized parameters produced better streamflow estimates than spatially uniform parameters (Widén-Nilsson et al. 2007). Beck et al. (2016) transferred the calibrated parameter sets of the HBV model from the selected donor catchments to 0.5° grid cell with the most similar climatic and physiographic characteristics and found that HBV with regionalized parameters outperformed nine state-of-the-art macroscale models. Liang & Greene (2020) employed the regression method to transfer the runoff coefficient and estimated global runoffs. The above studies concluded that regionalization methods have a great potential to be used for global hydrological modelling.

Various regionalization methods were proposed and numerous studies have compared these methods over different regions during the past decades (Abdulla & Lettenmaier 1997; Hundecha & Bárdossy 2004; Oudin et al. 2008; Pokhrel et al. 2008; Jin et al. 2009; Samaniego et al. 2010; He et al. 2011; Razavi & Coulibaly 2013; Li & Zhang 2017; Yang et al. 2019; Yang et al. 2020). However, there is still no clear conclusion on the best-performed regionalization method for use in ungauged catchments. In particular, there is a lack of detailed information regarding the performance of regionalization methods on the global scale and the guidance for the use of regionalization methods for global hydrological modelling. Therefore, the objectives of this study are to (1) evaluate the performance of the most widely used regionalization methods over thousands of catchments in the world to determine an appropriate global-scale regionalization scheme (GSRS) for global hydrological modelling and (2) propose a robust regionalization-based global hydrological simulation scheme (RGHS) to estimate global runoff.

The observed daily streamflow data used in this study were obtained from three sources, i.e. the Global Runoff Data Centre (GRDC; http://www.bafg.de/GRDC/), the Canadian model parameter experiment (CANOPEX) database (Arsenault et al. 2016, 2020), and some catchments of China (Gu et al. 2018). The GRDC dataset comprises river discharge data for more than 9,500 stations from 161 countries. Basin boundaries and flow paths were taken from HYDRO1 K (Gong et al. 2009). Continent boundaries were taken from STN-30p (Vörösmarty et al. 2000).

There are some uncertainties that may impact the results of global hydrological modelling. For example, small catchments might result in less reliable simulations since the model resolution is low (i.e. 0.5°). Therefore, the following three criteria were used to choose catchments for our analysis:

• (1)

  The streamflow record length was not shorter than 5 years (not necessarily consecutive) during the 1982–2015 period.

• (2)

  The catchment size is over 2,500 km² to ensure that each catchment covers at least one 0.5° grid cell.

• (3)

  The upper limit of catchment size was set as 50,000 km² to minimize the effects of river regulation and channel routing effects (Gericke & Smithers 2014; Beck et al. 2016, 2020).

Based on the above criteria, 2,277 catchments were selected in total for the comparison of regionalization methods and GSRS. Figure 1 shows the distribution of these catchments and climate classification extracted from the World Map of Köppen-Geiger climate classification (Kottek et al. 2006). The majority of those catchments are located in North America, Europe, Southeast Asia, and central South America. However, catchments located in the Middle East, North Africa, central Australia, and the Russian Far East and Siberia are much fewer.

The meteorological data used in this study are daily precipitation, air temperature, and potential evaporation. The potential evaporation data at the global scale were obtained from the Global Land Evaporation Amsterdam Model (GLEAM V3) potential evaporation dataset (1980–2015) at the 0.5° resolution, which is calculated by using the Priestley and Taylor equation based on observations of surface net radiation and near-surface air temperature (Miralles et al. 2010; Martens et al. 2017). The daily temperature data were obtained from the European Centre for Medium range Weather Forecasts (ECMWF)–Interim Reanalysis (ERA-Interim) at the 0.5° resolution for the 1979–2019 period (Dee et al. 2011). The precipitation data were obtained from the Global Precipitation Climatology Centre (GPCC) V.2018 precipitation dataset at the 0.5° resolution for the 1982–2016 period (Fuchs et al. 2009), which was interpolated using gauged precipitation data provided by national meteorological and hydrological services, regional and global data collections, as well as the World Meteorological Organization (WMO) GTS-data (GPCC, http://gpcc.dwd.de).

In addition, comprehensive reviews of regionalization methods were made by He et al. (2011) and Razavi & Coulibaly (2013), in which the number of times that catchment descriptors used in other studies were counted after reviewing the regionalization methods. Based on this information, 13 catchment descriptors, classified as climate index, terrain characteristics, land use, and soil characteristics, were selected for regionalization (Table 1). Those datasets with spatial resolutions <0.5° were resampled to 0.5° using bilinear averaging.

Four conceptual hydrological models, i.e. Génie Rural à 4 paramètres Journalier model (GR4 J), simple lumped conceptual daily rainfall-runoff model (SIMHYD), Xinanjiang model (XAJ), and Hydrological Model of École de technologie supérieure (HMETS), were used to simulate runoff at the daily time step. They were chosen because of their proven effectiveness around the world and successful application in regionalization studies (Perrin et al. 2003; Oudin et al. 2008; Zhang & Chiew 2009; Li et al. 2014; Zhang et al. 2014, 2016; Li & Zhang 2017; Chen et al. 2018; Shen et al. 2018; Yang et al. 2020; Qi et al. 2021b). In addition, the four hydrological models used in this study have different numbers of parameters and vary in the generalization of the natural hydrological processes and physical mechanisms (Zhao et al. 1980, Zhao 1992; Chiew et al. 2002; Perrin et al. 2003; Chiew 2010; Boumenni et al. 2017; Martel et al. 2017; Shen et al. 2018). For example, the model structure for representing the hydrological process is more complicated and detailed in SIMHYD, XAJ, and HMETS than the most parsimonious GR4 J. For the simulation of evaporation, SIMHYD and HMETS use a one-layer evaporation model, while a three-layer model is used in XAJ. As for the streamflow generation, XAJ and HMETS use the saturation excess runoff generation mechanism based on the soil moisture content of the aeration zone reaching its field capacity, while SIMHYD considers both infiltration-excess and saturation-excess runoffs by using an interception store and a soil moisture store. In addition, GR4 J and HMETS take into account the groundwater exchange based on surface water–groundwater interaction functions, while this is not the case for SIMHYD and XAJ. Since there is no snow module in GR4 J, SIMHYD, and XAJ, the snow module-CemaNeige (Valéry 2010) was incorporated into the original models in this study. CemaNeige allows estimating the snowmelt and simulating the snowpack evolution by using two parameters, and the coupling of these models and CemaNeige has been tested in some other studies (e.g. Valéry 2010; Coron et al. 2014; Hublart et al. 2015; Guo et al. 2020; Wang et al. 2020; Yang et al. 2020). Table 2 shows the main structures of the four hydrological models. More detailed information about model parameters can be found in Supplementary Appendix Table A1.

The regionalization methods used in this study are the most commonly used step-wise regionalization methods, i.e. the global mean method, the regression-based methods, and the distance/attribute-based methods (Jin et al. 2009; He et al. 2011; Razavi & Coulibaly 2013; Yang et al. 2018, 2019, 2020). Moreover, the great performance of these methods has been observed for the ungauged runoff prediction in many studies (Merz & Blöschl 2004; Oudin et al. 2008; Yang et al. 2018, 2019).

The global mean method is a simple regionalization method. Generally, the arithmetic mean of the parameters of all gauged catchments in the region is directly applied to the ungauged catchments.

The regression-based methods assume that a well-behaved relationship exists in the observable catchment characteristics and model parameters (Burn & Boorman 1993). Parameters for an ungauged catchment are derived by using the relationship between catchment descriptors and model parameters for the donor catchments (Xu 1999a, 1999b, 2003; Muller-Wohlfeil et al. 2003; Arsenault & Brissette 2014; Yang et al. 2018, 2019). The multiple linear regression method was used in this study, in which the relationships among model parameters and the selected 13 catchment descriptors were established using multiple linear regressions and these functions were then used to estimate model parameters for ungauged catchments.

Considering the limitation of the two distance/attribute-based methods mentioned earlier, some studies (e.g. Samuel et al. 2011; Viviroli & Seibert 2015; Yang et al. 2018) integrated spatial proximity with physical similarity to improve the regionalization ability. In the present study, the physical similarity method considering distance in which the distance was considered as one of the catchment descriptors was used and then the similarity index was calculated.

For distance/attribute-based methods, there are two different averaging options to transfer the model parameter sets from donor catchments: (1) parameter average option, which transfers the averaged model parameters from donor catchments to ungauged catchments; and (2) output average option, which averages runoff simulations calculated by using individual parameter sets from donor catchments to ungauged catchments (Oudin et al. 2008; Yang et al. 2018). In addition, there are two different weighting approaches used to combine the model parameters or model outputs: (1) Arithmetic Mean (AM) method and (2) Inverse Distance Weighted (IDW) method (Parajka et al. 2007; Yang et al. 2018). All regionalization methods used in this study are summarized in Table 3.

To find suitable donor catchments and the optimal regionalization scheme (i.e. combination of different optional regionalization methods), the performance of regionalization methods under different thresholds of model efficiency in terms of KGE was tested. Threshold values of model efficiency (all, ≥0, ≥0.5, ≥0.6, ≥0.7, ≥0.8, or ≥0.9) were determined for the calibration period. In other words, when the efficiencies of the catchments were below a threshold, these catchments were not used as donor catchments to simulate runoffs for ungauged catchments. The threshold named ‘all’ means that all catchments are used as donor catchments regardless of their performance in the calibration stage. What is more, all catchments, whether poorly or well modelled, were all considered as pseudo-ungauged under different efficiency thresholds for cross-validation.

The commonly used leave-one-out cross-validation strategy was used to evaluate the performance of the regionalization methods. Specifically, the performance of five parameter regionalization methods with two different average options, two weighting approaches (14 regionalization schemes, Table 3), and seven efficiency thresholds was compared over thousands of catchments based on four conceptual hydrological models. The best-performed regionalization method (with the highest KGE value) was then selected as the GSRS to regionalize parameters of four hydrological models to the spatial resolution of 0.5°×0.5° grid cell all over the world except for Antarctica. This procedure is based on an assumption that the parameters at the catchment scale can be transferred to the grid scale. In other words, the 0.5°×0.5° grid cell was treated as a catchment (Widén-Nilsson et al. 2007; Beck et al. 2016).

The values of these parameters for calibration are [4,5,6,7,8,9,10] for and [0.2,0.3,0.4,0.5,0.6] for b, respectively. They were chosen based on computer capability limitations and the physical meaning of each parameter. Therefore, there are 35 routing parameter sets, composed of five b values and seven values. In this study, the enumeration method was used to calibrate these two routing parameters. The catchment streamflow was generated by converging grid runoff from all grids in the catchment based on each routing parameter set. Since two routing parameters are calibrated rather than regionalized, it is hard to obtain streamflow time series for real ungauged catchments. To overcome this problem, the commonly used aggregation was used to obtain the simulation performance of ungauged catchments without calibration of routing parameters (Widén-Nilsson et al. 2007; Beck et al. 2020).

The performance of the RGHSs was evaluated in four ways (Figure 2).

• (1)

  The effectiveness of RGHSs was first evaluated by comparing the evaluation metrics of the 2,277 catchments obtained by RGHSs (RGHS-NRF and RGHS-aggregation) with the results obtained by calibration and regionalization at the catchment scale.

• (2)

  The second evaluation was done by comparing the evaluation metrics of the 2,277 catchments obtained by RGHSs with the results obtained from nine state-of-the-art macroscale models, coming from the EartH2Observe project (Schellekens et al. 2016), which were driven by WATCH Forcing Data ERA-Interim (WFDEI) for the period 1979–2012 at the 0.5° resolution. The catchments’ streamflow time series of these nine models were calculated by using the NRF method in this study. For the comparison of RGHSs versus nine models from the EartH2Observe project, RGHSs were driven by WFDEI (RGHS-WFDEI).

• (3)

  In addition, 420 independent catchments with surface area ranging between 50,000 and 100,000 km² and with streamflow records being longer than 5 years were used as validation catchments. The KGE, NSE, and AVE values of the validation catchments obtained by RGHSs and the results obtained by nine models were further compared to assess the value of RGHSs.

• (4)

  The last way to evaluate the effectiveness of RGHSs was by comparing the long-term averaged annual total runoff (1983–2015, 1982 is considered as the warm-up period) from global land and each continent obtained by RGHSs with the results in literatures.

Figure 3 shows the cumulative density function curves of the percentage of catchments with KGE, NSE, and AVE values exceeding the given value (the value of y-axis) for four hydrological models over the calibration and validation periods. Generally, the performances of the four hydrological models are close to each other. For all hydrological models, the KGE value is higher than 0.7 for more than 60% catchments at the calibration period and 0.5 for the validation period. Specifically, SIMHYD shows the best performance, followed by XAJ, GR4 J, and HMETS. It is the same when using NSE and AVE as evaluation metrics. Figure 4 presents the spatial distribution of the model efficiency for the calibration period in terms of KGE values. The results show that the KGE values are above 0.8 for most catchments in eastern Canada, the USA, southern China, and along the Atlantic Coast of Europe. The catchments in the American tropics, the Andes (South America), and northwest China show low KGE values, which might be because of the complex topography and insufficiency of precipitation observations for the calibration of hydrological models (Demirel et al. 2015; Vetter et al. 2015; Beck et al. 2016). Among the five different climate regions, the highest median KGE values are obtained for the cold regions, and the lowest for the arid regions (Figure 5). The good performance in the cold regions may be because the hydrological processes are not that sensitive to precipitation variability. The large different performance among four hydrological models in the equatorial and arid regions is likely due to their different skills in simulating the highly nonlinear process from precipitation to runoff.

The median KGE values of these hydrological models are 0.748 for GR4 J, 0.774 for SIMHYD, 0.766 for XAJ, and 0.750 for HMETS in the calibration period and 0.632 for GR4 J, 0.644 for SIMHYD, 0.641 for XAJ, and 0.617 for HMETS in the validation period. The performance of hydrological models used in this study is comparable to or higher than those from previous studies. For example, Beck et al. (2020) calibrated the HBV model for 4,229 catchments globally and obtained median daily KGE values of 0.770 and 0.690 for the calibration and validation periods, respectively. Alfieri et al. (2020) calibrated the LISFLOOD model for 1,126 catchments and obtained median daily KGE values of 0.670 and 0.610 for the calibration and validation periods, respectively.

For illustrative purposes, the model performance (GR4 J) of distance/attribute-based methods using different numbers of donor catchments with different weighting and averaging options under each KGE threshold is shown in Figure 6 (similar results can be found for the other three models, which are shown in Supplementary Appendix Figure A1). Results show that the IDW approach performs better than the AM approach for all distance/attribute-based regionalization methods and all hydrological models. This is consistent with the results from other studies (e.g. Samuel et al. 2011; Arsenault & Brissette 2014; Li et al. 2014). The worse performance of the AM approach is probably caused by the large difference of distance or similarity among our studied catchments. However, the weighting scheme of IDW minimizes the negative impact caused by the farthest distance or the least similar donors. Furthermore, the optimal number of donor catchments for the parameter averaging option always ranges between 1 and 5. For the output averaging option, it is between 4 and 6 (see Figure 6). In general, the output averaging option outperforms the parameter averaging option globally in this study, which is also consistent with previous regional studies (e.g. Samuel et al. 2011; Arsenault & Brissette 2014; Li et al. 2014; Yang et al. 2018, 2020) at the regional scale. Therefore, to balance the effect and the amount of computation, five donor catchments are suggested to be used for the output averaging method.

Poorly modelled catchments may yield higher uncertain model parameter values, but they may add some diversities for modelling ungauged catchments as well (Oudin et al. 2008). Therefore, it is worth evaluating the necessity to consider poorly modelled catchments in regionalization. In this study, we compared the performance of regionalization methods under different KGE thresholds for all regionalization methods (including regression method, not shown) and found that the KGE threshold of 0.5 is the best (see Figure 6). Including catchments with KGE values below 0.5 as donors has little effect on the performance of hydrological models for ungauged catchments. However, when using the efficiency threshold of 0.9, the hydrological model performance remarkably drops for ungauged catchments. This indicated that there is a loss of information between those thresholds, which is attributable to the loss of donor catchments. When using the threshold of 0.9, only 170 (for GR4 J), 211 (for SIMHYD), 236 (for XAJ), and 193 (for HMETS) donor catchments out of 2,277 were available and most of them are located in the USA, southern China, and along the Atlantic Coast of Europe. This result is consistent with those in Neri et al. (2020), as they concluded that the performance of regionalization methods strongly depends on the informative content of the dataset of available donor catchments. Therefore, the threshold of 0.5 is taken as the best performance threshold, and the total number of donor catchments for the global regionalization scheme ranges from 1902 (for GR4 J) to 1985 (for SIMHYD) for different hydrological models.

The performance differences among regionalization methods are observed in Figure 7. Generally, the SPI-OUT method outperforms the others, followed by PSDI-OUT and PSI-OUT, and the GM method performs the worst. Overall, the performances of different regionalization methods are consistent in each hydrological model. This indicates that the ranking of regionalization methods is independent of model structures. As for hydrological models, the two parsimonious hydrological models (i.e. GR4 J and SIMHYD) slightly outperform the other two more complex models (i.e. XAJ and HMETS) in most situations. Figure 8 summarizes the median KGE values of all regionalization methods for five climate regions. The results show that SPI-OUT outperforms other regionalization methods for all climate regions, followed by PSDI-OUT. GM and RE perform worse than other regionalization methods for all climate regions. In general, the performance of SPI-OUT is relatively stable (except for arid regions) for different climate regions. However, all regionalization methods perform much worse for arid regions than for the other four climate regions.

Based on the comparison of regionalization methods, it can be concluded that the SPI-OUT outperforms others. However, according to the principle of this method (i.e. the nearby catchments should have similar behaviours), the performance of the SPI-OUT method may deteriorate with the increase of the mean distance between the donor and ungauged catchments. Therefore, it is important to explore how the mean distance impacts the performance of regionalization methods. Figure 9 shows the percentage of the best-performed regionalization method under different mean distances (between the donor and ungauged catchments) under the threshold of 0.5 (the SPI-OUT, PSI-OUT, and PSDI-OUT were selected for their outstanding performance). For all hydrological models, when the mean distance between donors and ungauged catchments is smaller than 700 km, the average proportion of catchments that the SPI-OUT method outperformed others is the largest. It is hard to decide which method outperforms others when the mean distance is larger than 700 km for GR4 J and HMETS, 900 km for SIMHYD, and 1,000 km for XAJ. The advantage of the SPI-OUT method is reduced with the increased mean distance between donors and ungauged catchments. Overall, the greatest gains of the SPI-OUT method in performance are achieved for ungauged catchments with mean distances no more than 700 km from the donors for most hydrological models.

The GSRS was selected based on the best-performed catchment scale regionalization methods (Section 3.2). That is, only catchments with KGE value being greater than 0.5 at the calibration period were used as donor catchments. For grid cells with a mean distance less than 700 km to donors, the calibrated parameter sets of the five nearest donor catchments were transferred by using the SPI-OUT method. For grid cells with a mean distance larger than 700 km, the parameters were extracted from the PSDI-OUT method. Figure 10 shows the mean distance to the five nearest donor catchments. The mean distances are generally ≤700 km for most parts of the world. However, there were still some regions where the mean distance is >700 km, like Greenland, central Siberia, southern South America, Southeast Asia, and Northeast Africa.

By using the GSRS, the runoff was first calculated for all global grid cells except for Antarctica. The NRF approach was then used to converge gird runoff to catchment streamflows. The 2,277 catchments were used to evaluate the performance of the RGHSs, since all of the grids were regionalized (Table 4). The median KGE value of SIMHYD-RGHS is the largest (0.385) and that of HMETS-RGHS is the smallest (0.374). When using NSE, HMETS-RGHS has the largest value (0.283) and GR4 J-RGHS has the smallest value (0.231). In terms of AVE, SIMHYD-RGHS performs the best (0.686) and XAJ-RGHS performs the worst (0.672). All above-mentioned differences are considered not large.

Generally, RGHSs perform better than the catchment scale regression and global mean methods, but worse than other regionalization methods. The differences between the median KGE values of RGHSs and that using the best-performed catchment scale regionalization method (i.e. SPI-OUT) are 0.166 (for GR4 J), 0.173 (for SIMHYD), 0.166 (for XAJ), and 0.155 (for HMETS). The differences between the median KGE value of RGHSs and that using the calibrated parameters are 0.372 (for GR4 J), 0.389 (for SIMHYD), 0.386 (for XAJ), and 0.376 (for HMETS). The performance of RGHSs is about half of that obtained using calibrated parameters in terms of the median value of KGE. Similar results are also observed when using NSE. However, the differences between RGHSs and those using the calibrated parameters are small when using AVE. In addition, the median KGE and NSE values of 2,277 catchments of RGHS-aggregation are smaller than that of RGHS-NRF. However, when using AVE, the difference between these two methods becomes small. This result indicates that the aggregation method is suitable to be used in the water balance analysis. However, effective estimation of streamflow time series and extreme flows in ungauged catchments still needs to be further studied.

Figure 11 shows the spatial distribution of KGE for RGHSs. For most of Europe and the east coast of North America, KGE is generally above 0.5, and even above 0.7 over many catchments. In contrast, most catchments in southwestern Africa and northwestern Australia perform much worse. To further evaluate the effectiveness of the RGHSs in various climate regimes, Table 5 summarizes the median KGE values over 2,277 catchments for five Köppen-Geiger climate types. Generally, all RGHSs perform much worse for arid climate regions than for the other four climate regions. Previous studies (e.g. Oki et al. 2001; Döll et al. 2003; Beck et al. 2016) indicated that most of the existing GHMs overestimate runoff in arid basins. This is because the high evaporative losses in arid regions, the highly nonlinear response of runoff to rainfall, and the flashy nature of the streamflow time series make it difficult to simulate streamflow time series (Pilgrim et al. 1988; Widén-Nilsson et al. 2007; Beck et al. 2016).

To further evaluate the effectiveness of the RGHSs, the KGE, NSE, and AVE values obtained from four RGHSs are compared to those obtained from nine macroscale hydrologic models from the EartH2Observe project. Figures 12 and 13 present evaluation metrics over 2,277 catchments that participated in global-scale regionalizations and 402 independent validation catchments, respectively. The overall performance of RGHS-GPCC is better than that of the nine models, with varying degrees of improvement in the median, upper and lower quartiles, and whiskers of each performance metric (except for the HTESSEL model for the lower quartiles and whiskers). When comparing the models in terms of the median of each performance metric, the median KGE of RGHS-GPCC ranges between 0.374 (HMETS-RGHS) and 0.385 (SIMHYD-RGHS), while that of nine models ranges between −0.061 (PCR-GLOBWB) and 0.124 (HBV-SIMREG) over 2,277 catchments. As for the independent validation catchments, the median KGE of RGHSs is between 0.213 (HMETS-RGHS) and 0.231 (XAJ-RGHS), while that of nine models is between −0.154 (ORCHIDEE) and 0.082 (PCR-GLOBWB). Similar results are observed when using NSE and AVE as evaluation metrics.

Compared to RGHS-GPCC, the performance of RGHSs driven by WFDEI (RGHS-WFDEI) somewhat deteriorates. For example, the median KGE of RGHS-WFDEI is between 0.209 (GR4 J-RGHS) and 0.260 (HMETS-RGHS) over 2,277 catchments, and 0.056 (GR4 J-RGHS) and 0.143 (SIMHYD-RGHS) for the independent validation catchments, which is lower than that of RGHS-GPCC. However, RGHS-WFDEI still outperforms the other nine models from the EartH2Observe project in most cases. For example, over 2,277 catchments, the 25% quantile and median evaluation metrics of RGHS-WFDEI are always bigger than the other nine models and only the 75% quantile evaluation metrics of RGHS-WFDEI are lower than that of the HTESSEL. As for the independent validation catchments, the 25% quantile and median KGE and NSE of RGHS-WFDEI (GR4 J-RGHS) are only lower than that of the PCR-GLOBWB. The median AVE does not vary greatly among the models and the median AVE of RGHS-WFDEI is only lower than that of the HTESSEL. However, the 75% quantile NSE and AVE of RGHS-WFDEI have no advantages.

The spatial variability of long-term (1983–2015) averaged annual runoff is presented in Figure 14, which shows a very high spatial variability around the world. However, the spatial distributions of the long-term annual runoff from four RGHSs in this study are consistent. High long-term annual runoffs are obtained for the Amazon, maritime Southeast Asia, central Africa, and the southern section of the Andes in Patagonia. The long-term averaged annual runoffs in central Australia, North Africa, the Middle East, and the northwest of North America are low. The global long-term annual runoffs estimated by RGHSs are 46,810 km³/yr for GR4 J-RGHS, 42,733 km³/yr for SIMHYD-RGHS, 42,592 km³/yr for XAJ-RGHS, and 45,100 km³/yr for HMETS-RGHS. Overall, the simulated continental runoffs from the GR4 J-RGHS and HMETS-RGHS are larger than those from the SIMHYD-RGHS and XAJ-RGHS, especially for Asia and South America.

The global and continental water resources are calculated using RGHSs and compared to those calculated using the other five GHMs in literatures (Table 6). The six individual simulations (including the results of this study) of global long-term average water resources encompass a range between 29,485 and 46,810 km³/yr. The difference between these individual global water resource simulations is 17,325 km³/yr, which is even much higher than the global consumptive water use estimated using Global Water Use Model of WaterGAP 2 (1,250 km³/yr in 1995). A comparative study of seven models has shown that the range of predicted global runoff from different models was about 45% of the mean simulated runoff, which means that only limited agreement exists between GHMs (Haddeland et al. 2011). In fact, it is difficult to compare different GHMs, for there are different periods, data quality, spatial and temporal resolution, and so on of the data used in the building of the GHMs.

Therefore, we further compared the inter-annual variation of global runoff from RGHSs and the EartH2Observe project for the same period (1983–2012, see Figure 15). The results show that inter-annual variability of global runoff is relatively consistent. However, there is a large difference in magnitudes. In particular, large differences exist for nine global models from the EartH2Observe project, even if they were all driven by the same precipitation dataset (i.e. WFDEI). However, the annual global runoff from RGHSs driven by GPCC and WFDEI are similar and in the range of nine global models, probably due to the fact that the same GSRS were used for global runoff simulations.

The regionalization methods used in this study are called step-wise regionalization, namely, the model parameters are calibrated, and then transferred to the ungauged regions (Samaniego et al. 2017). The equifinality problem of hydrological models may be one of the problems when using this method (Bárdossy 2007; Götzinger & Bárdossy 2007), because only one parameter set is used for the regionalization, even though the hydrological model may have multiple sets of parameters that lead to equally acceptable model performance (Beven & Freer 2001; Samaniego et al. 2010). To illustrate the influence of equifinality on the regionalization performance, another nine calibration sets were generated during model calibration with different initial random seeds. The difference of median KGE values between these sets is no more than 0.05. For each of 2,277 catchments, one of the nine calibrated parameter sets was then randomly selected for the catchment scale regionalization. Figure 16 shows the performance of the randomly selected parameter set and the original parameter set under the threshold of 0.5. The results show that the regionalization performance of the randomly selected parameter set is consistent with the results of the original parameter set. This indicates that the equifinality does not contribute significantly to the overall uncertainty in the applications of regionalization methods. A similar conclusion was also drawn by Arsenault & Brissette (2014).

The comparison of different models showed that great uncertainties exist in global water resources simulation. The uncertainties are usually induced by meteorological input, hydrological model structures and parameters (Nijssen et al. 2001; Tuo et al. 2016; Sawunyama & Hughes 2018). This study preliminarily investigated the impacts of precipitation inputs on global water resources by running RGHSs using two different precipitation datasets (i.e. GPCC and WFDEI). The global long-term annual runoffs estimated by RGHS-GPCC range from 42,592 to 46,810 km³/yr and those estimated by RGHS-WFDEI range from 48,132 to 53,365 km³/yr. This indicates that the precipitation input is one of the largest uncertainty sources for global water resources estimation.

In addition, the difference of the global long-term average annual runoff is 4,218 km³/yr between RGHS-GR4 J and RGHS-XAJ (driven by GPCC), indicating that the instability of hydrological models could also result in uncertainty in global water resources estimation. Similar conclusions were also drawn in Schellekens et al. (2016). They simulated the global water budget by using ten global models driven by the same precipitation data and found significant differences in global runoff among these models. In fact, due to the complexity of the water cycle and the lack of precise understanding of natural hydrological rules, each hydrological model emphasizes different aspects and has its own conceptualization of hydrological processes, which leads to uncertainty in water resources estimation. In addition, multiple objective functions and various optimization algorithms in hydrological models would result in different parameter sets and may lead to a certain level of uncertainty in water resource modelling (Kunnath-Poovakka et al. 2021). Therefore, the reduction of uncertainty for global water resources estimation is an avenue for future studies.

Moreover, the observed discharge rather than natural discharge was used to calculate the parameter maps and global water resources. Almost all large rivers are regulated (Vörösmarty et al. 2004; Nilsson et al. 2005), which may result in modified high-flow and low-flow and affect the performance of nature catchment streamflow simulation and the performance of regionalization of hydrological models (Widén-Nilsson et al. 2007; Lehner et al. 2011). Therefore, catchments larger than 50,000 km² were discarded in this study, as regulation by dams exists in larger catchments with high probability (Vörösmarty et al. 2004; Nilsson et al. 2005). In general, considering the impact caused by regulated catchments is a challenge in global hydrological modelling due to the lack of global information to an adequate accuracy. Many GHMs, therefore, do not include regulation effects. For example, Widén-Nilsson et al. (2009) assumed that regulation has less influence on average flow volumes and used long-term average runoff instead of time series to minimize the effect of the regulation problem on model calibration. Beck et al. (2016) used the Global Reservoir and Dam (GRanD) database (v1.1) (Lehner et al. 2011) to exclude catchments influenced by reservoirs. Arheimer et al. (2020) concluded the potential and difficulty to improve the global model by individual modelling river regulation. Therefore, further studies need to be carried out to consider the influence of human activities (e.g. regulation) to improve the performance of regionalization and GHMs.

This study first evaluated the performance of multiple regionalization methods over 2,277 catchments distributed around the world based on four catchment hydrological models. A GSRS was then selected to regionalize hydrological parameters at 0.5° grid cell, and four RGHSs (GR4 J-RGHS, SIMHYD-RGHS, XAJ-RGHS, and HMETS-RGHS) were presented for global water resource estimation. Finally, nine global models from the EartH2Observe project were selected to further validate the performance of RGHSs in global water resources estimation. The following conclusions can be drawn.

• (1)

  The SPI-OUT method offered the best results with the highest KGE value when using the 0.5 efficiency threshold for selecting donor catchments, and the optimal number of donor catchments ranges between 3 and 6 for the output averaging option.

• (2)

  Compared to the models from the Earth2Observe project, the presented RGHSs provide better simulation results. For example, the median KGE values of RGHSs were from 0.374 (HMETS-RGHS) to 0.385 (SIMHYD-RGHS). However, the median KGE values of EartH2Observe were −0.06 (PCR-GLOBWB) to 0.124 (HBV-SIMREG). Similar results were also observed when using NSE and AVE as evaluation metrics. On the whole, RGHSs offered great performance in the cold and temperate regions, while the poor performance of the RGHSs was observed in the arid climate regions.

• (3)

  The RGHSs showed reasonable estimations of global water resources. The global long-term average annual runoff estimated by RGHSs are 46,810 (for GR4 J-RGHS), 42,733 (for SIMHYD-RGHS), 42,592 (for XAJ-RGHS), and 45,100 (for HMETS-RGHS) km³/yr, which is in the good range of those estimated by other GHMs in the literature.

This work was partially supported by the National Key Research and Development Program of China (No. 2017YFA0603704), the National Natural Science Foundation of China (Grant Nos. 51779176, 51539009), the Overseas Expertise Introduction Project for Discipline Innovation (111 Project) funded by the Ministry of Education of China and State Administration of Foreign Experts Affairs P.R. China (Grant No. B18037), the Research Council of Norway (FRINATEK Project 274310), and the Thousand Youth Talents Plan from the Organization Department of CCP Central Committee (Wuhan University, China).

Data cannot be made publicly available; readers should contact the corresponding author for details.