The global water resources and use model WaterGAP v2.2d: model description and evaluation

Abstract

WaterGAP is a global hydrological model that quantifies human use of groundwater and surface water as well as water flows and water storage and thus water resources on all land areas of the Earth. Since 1996, it has served to assess water resources and water stress both historically and in the future, in particular under climate change. It has improved our understanding of continental water storage variations, with a focus on overexploitation and depletion of water resources. In this paper, we describe the most recent model version WaterGAP 2.2d, including the water use models, the linking model that computes net abstractions from groundwater and surface water and the WaterGAP Global Hydrology Model (WGHM). Standard model output variables that are freely available at a data repository are explained. In addition, the most requested model outputs, total water storage anomalies, streamflow and water use, are evaluated against observation data. Finally, we show examples of assessments of the global freshwater system that can be achieved with WaterGAP 2.2d model output.

1 Introduction

A globalized world is characterized by large flows of virtual water among river basins and by international responsibilities for the sustainable development of the Earth system and its inhabitants. The foundation of a sustainable management of water, and more broadly the Earth system, are quantitative estimates of water flows and storages as well as of water demand by humans and freshwater biota on all continents of the Earth . During the last three decades, global hydrological models (GHMs) have been developed and continually improved to provide this information. They enable the determination of the spatial distribution and temporal development of water resources and water stress for both humans and other biota under the impact of global change (including climate change). In addition, global-scale knowledge about water flows and storages on land is necessary to understand the Earth system, including interactions with the ocean and the atmosphere as well as gravity distribution and crustal deformation (affecting GPS).

Such models are frequently used in large-scale assessments, such as the assessment of virtual water flows for products within the framework of the Intergovernmental Panel on Climate Change and the assessment of impacts based on scenarios for a sustainable future (such as the Sustainable Development Goals). Furthermore, global-scale modeling of water use and water availability is frequently used to evaluate large-scale water issues, for example water scarcity and droughts .

Some of these models are contributing to the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) where the focus is on both the model evaluation/improvement and the impact assessment of anthropogenic changes such as human water use or climate change. A series of evaluation exercises shows that high-performing simulation is challenging due to uncertain process representation at the given resolution, input data uncertainty and unequal data availability in terms of spatial and temporal distribution, e.g., river discharge observations . In this context, a proper model description is of great value for a better understanding of the process representation and parameterization of such models, and a related work is in progress .

A continuous improvement of process representations in GHMs is required to reduce uncertainty in assessments of water resources over historical periods and thus increase confidence in future projection assessments. In the recent past, some of the GHM approaches consider new processes such as the CO${}_{\mathrm{2}}$ fertilization effect or gradient-based groundwater models . Improved methods for the estimations of agricultural and other water use have been developed, and total water storage data from satellite observations are being increasingly employed either for evaluation or calibration/assimilation of models . Ultimately, there are attempts to achieve a finer spatial resolution than the typically used 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ grid cell .

Water – Global Assessment and Prognosis (WaterGAP), which has been developed since 1996, is one of the pioneers in this field. WaterGAP as described here operates with a spatial resolution of 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ and is part of the model family WaterGAP 2. Key model versions are WaterGAP 2.1d , 2.1e , 2.1f , 2.1g , 2.1h , 2.2 , 2.2a , 2.2(ISIMIP2a) , 2.2b , 2.2c () and 2.2d (this paper). In addition, a model family with 5${}^{\prime }$ $\times$ 5${}^{\prime }$ is named WaterGAP 3 . While the model family 3 has similar algorithms to the model family 2, this paper only refers to the recent model version WaterGAP 2.2d.

The major model purpose was to quantify global-scale water resources with a specific focus on anthropogenic inventions due to human water use and man-made reservoirs, to assess water stress. Furthermore, a lot of effort have been assigned to specific water storages like groundwater, lakes and wetlands. In the previously mentioned evaluation studies, WaterGAP has been qualified as a robust and qualitatively good-performing model in those key issues and for most climate zones worldwide.

Since the last complete model description of WaterGAP 2.2 , a number of modifications and improvements have been achieved. To be able to follow these changes and to transparently understand the process representation, a new model description can guide model output data users, especially in the case of discrepant model outputs from a GHM ensemble approach, and the GHM developing community in general. Hence, the aim of this paper is to provide an overview of the newest model version WaterGAP 2.2d by

• comprehensively describing the full model including all developments since WaterGAP 2.2 ,

• showing and discussing standard model output,

• providing insights into model evaluation, and

• giving guidance for the users of model output.

The framework of WaterGAP 2.2d is presented in Sect. , followed by the in-depth description of the water use models (Sect. ) and the global hydrological model (Sect. ). The description of standard model outputs is given in Sect.  including caveats of using the model outputs. In Sect. , model output is compared against multiple observation-based datasets, followed by typical model applications in Sect.  and the conclusions and outlook (Sect. ). The Supplement contains a table of symbols used in the equations (Table S1) and abbreviations, highlights the current fields of scientific use of WaterGAP, and shows additional figures (Figs. S1–S12).

2 WaterGAP 2 framework

WaterGAP 2 consists of three major components, the global water use models, the linking model Groundwater-Surface Water Use (GWSWUSE) and the WaterGAP Global Hydrology Model (WGHM) (Fig. ). Five global water use models for the sectors irrigation , livestock, domestic, manufacturing and cooling of thermal power plants compute consumptive water use and, in the case of the latter three sectors, also withdrawal water uses. Consumptive water use refers to the part of the withdrawn ($=$ abstracted) water that evapotranspirates during use. Whereas the output of the Global Irrigation Model (GIM) is available at monthly resolution, annual time series are calculated by all non-irrigation water use models (Sects. , ). The linking model GWSWUSE serves to distinguish water use from groundwater and from surface water bodies (Sect. ). It computes withdrawal water uses from and return flows to the two alternative water sources to generate monthly time series of net abstractions from surface water (NA${}_{\mathrm{pot},\mathrm{s}}$) and from groundwater (NA${}_{\mathrm{pot},\mathrm{g}}$) . These time series are input to the WGHM, affecting the daily water flows and storages computed by it (Sect. ).

Figure 1

The WaterGAP 2.2d framework with its water use models and the linking module GWSWUSE that provides potential net water abstraction from groundwater and surface water as input to the WaterGAP Global Hydrology Model (WGHM). Figure adapted from .

[Figure omitted. See PDF]

2.1 Spatial coverage and climate forcings

The WaterGAP 2 framework operates on the so-called CRU land–sea mask , which covers the global continental area (including small islands and Greenland but excluding Antarctica) with 67 420 grid cells in total, each 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ in size, which represents approx. 55 $\mathrm{km}$ $\times$ 55 $\mathrm{km}$ at the Equator. WaterGAP uses the continental area of the grid cell, which is defined as the cell area (calculated with equal area cylindrical projection) minus the ocean area with the borders according to the ESRI worldmask shapefile . The continental area comprises land area and surface water body area (lakes, reservoirs and wetlands only; river area is not considered). Since WaterGAP 2.2a, surface water body areas, and consequently land area, are dynamic and are updated in each time step.

Both GIM and WGHM use meteorological input data that consist of air temperature, precipitation, downward shortwave radiation and downward longwave radiation, all with daily temporal resolution. Various global meteorological datasets (hereafter referred to as climate forcings) were developed by the meteorological community at the 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ spatial resolution, such as WFD , WFDEI , GSWP3 , the Princeton meteorological forcing , and recently ERA5 and WFDE5 . Alternative climate forcings may lead to significantly different WaterGAP outputs .

2.2 Modifications of WaterGAP since version 2.2

The general framework of WaterGAP 2.2d does not differ from model version 2.2 described in . Improvements of water use modeling since WaterGAP 2.2 include, among others, deficit irrigation in regions with groundwater depletion (Sect. ) as well as integration of the Historical Irrigation Dataset (HID), which provides the historical cell-specific development of the area equipped for irrigation . Major improvements in WGHM include (1) a consistent river-storage-based method to compute river flow velocity; (2) simulation of land area dynamics in response to varying areas of lakes, reservoirs and wetlands; (3) groundwater recharge from these surface water bodies in (semi)arid grid cells; (4) if daily precipitation is below a threshold value, the potential groundwater recharge remains in the soil and does not (as in WaterGAP 2.2) become surface runoff; (5) return flows to groundwater from surface water use are corrected (by adjusting NA${}_{\mathrm{g}}$) by the amount of NA${}_{\mathrm{pot},\mathrm{s}}$ that cannot be satisfied; and (6) the integration of reservoirs by taking into account their commissioning year (and not assuming anymore that they have existed during the whole study period). Other changes concern model calibration or consist of the inclusion of new datasets and software improvements. A complete list of modifications of WaterGAP 2.2d compared to WaterGAP 2.2 is provided in Appendix .

3 WaterGAP water use models

3.1 Global Irrigation Model

Irrigation accounts for 60 %–70 % of global withdrawal water uses and 80 %–90 % of global consumptive water uses, and for even larger shares in almost all regions with severe water stress and groundwater depletion . Therefore, a reliable simulation of irrigation water use is decisive for the quality of WaterGAP simulations of streamflow and water storage in groundwater and surface water bodies as well as for the reliability of computed water stress indicators. Based on information on irrigated area and climate for each grid cell, GIM computes first cell-specific cropping patterns and growing periods and then irrigation consumptive water use (ICU), distinguishing only rice and non-rice crops . ICU can be regarded as the net irrigation requirement that would lead to optimal crop growth.

3.1.1 Computation of cropping patterns and growing periods of rice and non-rice crops

The cropping pattern for each cell with irrigated cropland describes whether only rice, non-rice crops or both are irrigated during either one or two growing seasons. The growing period for both crop types is assumed to be 150 d. A total of 17 cropping patterns are possible including simple variants (e.g., one cropping season with non-rice on the total irrigated area) and complex variants (non-rice after rice on one part of the total irrigated area and non-rice after non-rice on the other). The following data are used to model the cropping pattern: total irrigated area, long-term average temperature and soil suitability for paddy rice in each cell, harvested area of irrigated rice in each country, and cropping intensity in each of 19 world regions. In a second step, the optimal start date of each growing season is computed for each crop. To this end, each 150 d period within a year is ranked based on criteria of long-term average temperature, precipitation and potential evapotranspiration provided in . The most highly ranked 150 d period(s) is (are) defined as growing season(s).

3.1.2 Computation of consumptive water use due to irrigation

GIM implements the Food and Agriculture Organization of the United Nations (FAO) CROPWAT approach of to compute crop-specific ICU per unit irrigated area ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$) during the growing season as the difference between crop-specific optimal evapotranspiration $E_{{\mathrm{pot}}_{\mathrm{c}}}$ and effective precipitation $P_{\mathrm{irri},\mathrm{eff}}$ if the latter is smaller than the former, with

1$$\mathrm{ICU}=\left\{\begin{array}{ll}{E}_{{\mathrm{pot}}_{\mathrm{c}}}-P_{\mathrm{irri},\mathrm{eff}}& E_{{\mathrm{pot}}_{\mathrm{c}}}>P_{\mathrm{irri},\mathrm{eff}}\\ \mathrm{0}& \text{otherwise}\end{array}\right.,$$ where $E_{{\mathrm{pot}}_{\mathrm{c}}}$ is the product of potential evapotranspiration $E_{\mathrm{pot}}$ and the dimensionless crop coefficient $k_{\mathrm{c}}$ which depends on the crop and the crop development stage . As a standard, $E_{\mathrm{pot}}$ is calculated according to Eq. (). $P_{\mathrm{irri},\mathrm{eff}}$ is the fraction of the total precipitation $P$ (including rainfall and snowmelt) that is available to plants and is computed as a simple empirical function of precipitation. Equation () is implemented with a daily time step, but to take into account the storage capacity of the soil and to remain consistent with the CROPWAT approach, daily precipitation values are averaged over 10 d, except for rice-growing areas in Asia, where the averaging period is only 3 d to represent the limited soil water storage capacity in the case of paddy rice .

3.1.3 Irrigated area

In the standard version of WaterGAP 2.2d, irrigated area per grid cell used in GIM is based on the HID , which provides area equipped for irrigation (AEI) in 5 arcmin grid cells for 14 time slices between 1900 and 2005. HID data are aggregated to 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ and temporally interpolated to obtain an annual time series of AEI. Cropping patterns and growing periods are generated for every year, with an individual combination of year-specific AEI and harvested area of rice and the respective 30-year climate averages, which are then used to calculate ICU for every day of the same year (Sect. ). Harvested area of rice per country from the MIRCA2000 dataset, representative for the year 2000 , is scaled according to annual AEI country totals, ensuring consistency to AEI.

To take into account that not the whole AEI is actually used for irrigation in any year, country-specific values of the ratio of area actually irrigated (AAI) to AEI are used to estimate AAI in each grid cell. AAI is then applied for calculating the consumptive irrigation water use in volume per time. AAI $/$ AEI ratios were derived from the Global Map of Irrigation Area (GMIA) for 2005 . To set AAI from 2006 to 2016, we found country-specific AAI for 2006–2008 from the AQUASTAT database of the FAO, other international organizations, and national statistical services (e.g., EUROSTAT and USDA) for 61 countries. For these countries, the AAI values for 2009–2016 were set to the 2008 values, while for the rest of the countries, AAI was set to the 2005 values for the whole period 2006–2016.

Alternatively, as in previous WaterGAP versions, GIM in WaterGAP 2.2d can be executed based on a temporally constant dataset of AEI per grid cell, e.g., the GMIA for 2005 . Cropping patterns and growing periods are then computed for AEI and harvested area of rice in a reference year and the pertaining 30-year average climate. For more details and application examples, we refer to and .

3.2 Non-irrigation water uses

Although irrigation water use is the dominant water use sector globally, non-irrigation water uses, particularly in terms of withdrawal water uses, play a major role in Europe and America . Competition between agricultural and non-agricultural water uses are not uncommon , and the estimation of water demands becomes even more crucial when water resources are scarce. Statistical information on withdrawal water uses and consumptive water uses for domestic, industrial and livestock purposes are difficult to obtain on a country basis since no comprehensive global database does exist. However, the FAO collects relevant water-related data from national statistics and reports to provide a comprehensive view on the state of sectoral water uses. Unfortunately, the database lacks data in space and time, and hence modeling is of importance to fill these gaps .

3.2.1 Livestock

Withdrawal water uses for livestock are computed annually by multiplying the number of animals per grid cell by the livestock-specific water use intensity . The number of livestock are taken from . It is assumed that the withdrawal water uses for livestock are equal to their consumptive water use.

3.2.2 Domestic

Domestic water use comprises withdrawal water uses and consumptive water uses of households and small businesses and is estimated on a national level. The main concept is to first compute the domestic water use intensity (${\mathrm{m}}^{\mathrm{3}}$ per capita per year) and then to multiply this by the population of water users in a country. The domestic water use intensity is expressed by a sigmoid curve which indicates how water use intensity (per capita water use) changes with income (gross domestic product per capita) and is derived from historical data on a national or regional level . Besides changes driven by income and population, technological changes are considered to reflect improvement in water-use efficiency. Continuous improvements in technology make appliances more water efficient and, hence, contribute to reductions in water use. Detailed data on domestic consumptive water uses do not exist from statistics, but a simple balancing equation has been used in WaterGAP since the year 2000 to simulate consumptive water uses as the difference between withdrawal water use and wastewater volume (i.e., return flow) as the latter information is available from statistics. The calculation of consumptive water use before the year 2000 is based on the application of consumptive water use coefficients that accounts for the proportion of the withdrawal water use that is consumed. In order to allow for a spatially explicit analysis, country values of domestic water uses are allocated to grid cells (0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$) within the country based on the geo-referenced historical population density maps from HYDE version 3.1 . Additionally, population numbers beyond 2005 as well as information on the ratio of rural to urban population of each grid cell come from .

3.2.3 Manufacturing

The manufacturing sector is rather diverse in terms of water use and varies between countries and subsectors, for example highly water-intensive production processes in the chemical industry compared to the processes in the glass industry that use less water. In WaterGAP, the manufacturing water use model simulates the annual withdrawal water use and consumptive water use of water that is used for production and cooling processes, whereas the water used for power generation is modeled separately. A manufacturing structural water intensity that describes the ratio of water abstracted over the manufacturing gross value added (GVA) is derived per country for the base year 2005 (in ${\mathrm{m}}^{\mathrm{3}}\mathrm{USD}(\mathrm{constant}\mathrm{for}\mathrm{the}\mathrm{year}\mathrm{2000}{)}^{-\mathrm{1}}$) based on national statistics . GVA is found to be positively correlated with the sector's withdrawal water uses and is used as the driving force to reflect the time variant system. In addition, technological improvements are considered through a technological change factor.

The consumptive water use for this sector is obtained by using the same approach as described for the domestic sector, i.e., the calculation of the difference between the withdrawal water use and the return flows (starting in the year 2000) and the application of a consumption factor before the year 2000. Contrary to the domestic sector, return flows from the manufacturing sector are further subdivided into cooling water and wastewater. For countries where no data are available, the fraction of consumptive water use is derived from neighboring or economically comparable countries. Less information is available on the location of manufacturing industries; therefore country-level manufacturing water use is downscaled to grid cells proportional to its urban population .

3.2.4 Thermal power

Water is abstracted and consumed for the production of thermal electricity, particularly for cooling purposes where water is used to condense steam from the turbine exhaust. The volume of cooling withdrawal water use and consumptive water use is modeled on a grid-cell level based on input data on the location, type and size of power stations from the World Electric Power Plant Database . Here, the annual cooling water requirements in each grid cell are calculated by multiplying the annual thermal electricity production with the respective water-use intensity of each power station . A key driver is the annual thermal electricity production ($\mathrm{MWh}{\mathrm{yr}}^{-\mathrm{1}}$) on a country basis, which is downscaled to the level of thermal power plants according to their capacities. Time series on thermal electricity production per country until 2010 are available online from the Energy Information Administration . Cooling water intensities in terms of withdrawal water use and consumptive water use vary between plant types and cooling systems. Therefore, the model distinguishes between four plant types (biomass and waste, nuclear, natural gas and oil, coal, and petroleum) and three cooling systems (tower cooling, once-through cooling, ponds) . The approach is complemented by considering technological change leading to reduced intensities.

In general, water abstractions of once-through flow systems are considerably higher compared to the withdrawal intensities of pond cooling or tower cooling systems. In contrast, consumptive water use of tower cooling systems is much higher than water consumed by once-through cooling systems. In ordering plant-type-specific water intensities, i.e., water abstraction per unit electricity production, it becomes obvious that intensities are highest for nuclear power plants, followed by fossil, biomass, and waste-fuelled steam plants, while natural gas and oil combined-cycle plants have the lowest intensities, respectively. The model has been validated for the year 2005 by comparing modeled values with published thermoelectric withdrawal water uses .

3.3 GWSWUSE

The linking model GWSWUSE computes the fractions of all five sectoral water abstractions, or withdrawal water use, WU and consumptive water use CU in each grid cell that stem from either groundwater or surface water bodies (lakes, reservoirs and river). Time series for WU and CU from the sectoral water use models are an input to GWSWUSE except for WU for irrigation. The latter is computed within GWSWUSE as water use efficiencies CU/WU for irrigation are assumed to vary between surface water and groundwater. Country-specific efficiency values are used for surface water irrigation, while in the case of groundwater irrigation, water use efficiency is set to a relatively high value of 0.7 worldwide . In GWSWUSE, CU due to irrigation is decreased to 70 % of optimal CU in groundwater depletion areas; these areas were defined as grid cells with a groundwater depletion rate for 1980–2009 of more than 5 mm yr${}^{-\mathrm{1}}$ and a ratio of WU for irrigation over WU for all sectors of more than 5$\mathrm{\%}$ as computed for optimal irrigation in .

Sectoral groundwater fractions were derived individually for each grid cell in the case of irrigation and for each country in the case of the other four water use sectors . They are assumed to be temporally constant. Water for livestock and the cooling of thermal power plants is assumed to be extracted exclusively from surface water bodies.

Finally, GWSWUSE computes monthly time series of net abstraction from surface water NA${}_{\mathrm{pot},\mathrm{s}}$ and from groundwater NA${}_{\mathrm{pot},\mathrm{g}}$ which are used as input to WGHM. Net abstraction is the difference between total water abstraction from one of the two sources and the return flow to the respective source according to Eqs. (1), (3) and (4) in . In all sectors except irrigation, return flows are only directed to surface water bodies. The fraction of return flow to groundwater in the case of irrigation water use is estimated as a function of degree of artificial drainage in the grid cell (Sect. 2.1.3 in ). Positive net abstraction values refer to the situation where storage is reduced due to human water use, and negative values indicate an increase in storage. In the case of groundwater, the latter only occurs if there is irrigation with surface water in the grid cell. The approach of direct net abstractions implicitly assumes instantaneous return flows. The sum of NA${}_{\mathrm{pot},\mathrm{g}}$ and NA${}_{\mathrm{pot},\mathrm{s}}$ is equivalent to (potential) consumptive water use. NA${}_{\mathrm{pot},\mathrm{s}}$ and NA${}_{\mathrm{pot},\mathrm{g}}$ as computed by GWSWUSE are potential net abstractions that may be adjusted depending on the availability of surface water (Sect. ).

4 WaterGAP Global Hydrology Model (WGHM)

The WGHM simulates daily water flows and water storage in 10 compartments (Fig. ). The vertical water balance (dashed box in Fig. ) encompasses the canopy (Sect. ), snow (Sect. ) and soil (Sect. ) components. Water storage in glaciers is not simulated by WaterGAP 2.2d. The lateral water balance includes groundwater (Sect. ), lakes, man-made reservoirs, wetlands (Sect. ) and rivers (Sect. ). Different to the vertical water balances, where the water balance is calculated based on water height units ($\mathrm{mm}$), the lateral water balance is calculated in volumetric units (${\mathrm{m}}^{\mathrm{3}}$). Water height units are converted to volumetric units by considering the land area (for flows) or continental area (for storages) of the grid cell, respectively. Local surface water bodies are defined to be recharged only by runoff generated in the cell itself, while global ones additionally receive streamflow from upstream cells (Fig. ). Upstream–downstream relations among the grid cells are defined by the drainage direction map DDM30 . Each cell can drain only into one of the eight neighboring cells as streamflow. There is no groundwater flow between grid cells.

The amount of water reaching the soil is regulated by the canopy and snow water balance. Total runoff from the land fraction of the cell $R_{\mathrm{l}}$ is calculated from the soil water balance. $R_{\mathrm{l}}$ is then partitioned into fast surface and subsurface runoff $R_{\mathrm{s}}$ and diffuse groundwater recharge $R_{\mathrm{g}}$. Lateral routing of water through the storage compartments is based on the so-called fractional routing scheme and differs between (semi)arid and humid grid cells (red and green arrows in Fig. ). The definition of (semi)arid and humid cells is given in Appendix . To avoid that all runoff generated in the grid cell is added to local lake or wetland storage, only the fraction $f_{\mathrm{swb}}$ times $R_{\mathrm{s}}$ flows into surface water bodies, and the remainder discharges into the river. The factor $f_{\mathrm{swb}}$ is calculated as the relative area of wetlands and local lakes in a grid cell multiplied by 20 (representing the drainage area of surface water bodies), with its maximum value limited to the cell fraction of continental area. In humid cells, groundwater discharge $Q_{\mathrm{g}}$ is partitioned using $f_{\mathrm{swb}}$ into discharge to surface water bodies and discharge to the river segment. In (semi)arid cells, surface water bodies (excluding rivers) are assumed to recharge the groundwater to mimic point recharge. To avoid a short circuit between groundwater and surface water bodies, the whole amount of $Q_{\mathrm{g}}$ flows into the river. Loosing conditions, where river water recharges the groundwater, are not modeled in WGHM.

In WaterGAP, human water use is assumed to affect only the water storages in the lateral water balance. Increases in soil water storage in irrigated areas are not taken into account as the WaterGAP approach of direct net abstractions implicitly assumes instantaneous return flows. To consider anthropogenic consumptive water use in the output variable of actual evapotranspiration $E_{\mathrm{a}}$ (Table ), we sum up all evapo(transpi)ration components and actual consumptive water use WC${}_{\mathrm{a}}$ (see note 5 in Table ). NA${}_{\mathrm{s}}$ is abstracted from the different surface water bodies except wetlands with the priorities shown as numbers in Fig. .

Outflow from the final water storage compartment in each cell, the river compartment, is streamflow ($Q_{\mathrm{r},\mathrm{out}}$), which becomes inflow into the next downstream cell.

The ordinary differential equations describing the water balances of the 10 storage compartments simulated in WGHM are solved sequentially for each daily time step in the following order: canopy, snow, soil, groundwater, local lakes, local wetlands, global lakes, global reservoirs/regulated lakes, river (Fig. ). An explicit Eulerian method is used to numerically solve all differential equations except those for global lakes and rivers, where an analytical solution is applied to compute storage change during one daily time step, which allows daily time steps instead of smaller time steps that would have been required in the case of an explicit Eulerian method. As the water balances of global lakes, global reservoirs/regulated lakes and river of a grid cell are not independent from those of the upstream grid cells, the sequence of grid cell computations starts at the most upstream grid cells and continues downstream according to the drainage direction map DDM30 .

Figure 2

Schematic of WGHM in WaterGAP 2.2d. Boxes represent water storage compartments, and arrows represent water flows. Green (red) color indicates processes that occur only in grid cells with humid ((semi)arid) climate. For details the reader is referred to Sect.  to , in which the water balance equations of all 10 water storage compartments are presented.

[Figure omitted. See PDF]

4.1 General model variants of human water use and reservoirs

The standard model setup of WGHM in WaterGAP 2.2d simulates the effects of both human water use and man-made reservoirs (including their commissioning years) on flows and storages and is referred to as “ant” simulation (anthropogenic). These stressors can be turned off in alternative model setups to simulate a world without these two types of human activities and to quantify the direct impact of human water use and reservoirs.

• “Nat” simulations compute naturalized flows and storages that would occur if there where neither human water use nor global man-made reservoirs/regulated lakes.

• “Use only” simulations include human water use but exclude global man-made reservoirs/regulated lakes.

• “Reservoirs only” simulations exclude human water use but include global man-made reservoirs/regulated lakes.

The following sections generally refer to ant simulations.

4.2 Canopy

Canopy refers to the leaves and branches of terrestrial vegetation that intercept precipitation. Modeling of the canopy processes does not differentiate between rain and snow.

4.2.1 Water balance

The canopy storage $S_{\mathrm{c}}$ ($\mathrm{mm}$) is calculated as

2$${\displaystyle \frac{\mathrm{d}{S}_{\mathrm{c}}}{\mathrm{d}t}}=P-P_{\mathrm{t}}-E_{\mathrm{c}},$$ where $P$ is precipitation ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$); $P_{\mathrm{t}}$ is throughfall, the fraction of $P$ that reaches the soil ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$); and $E_{\mathrm{c}}$ is evaporation from the canopy ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$).

4.2.2 Inflows

Daily precipitation $P$ is read in from the selected climate forcing (see Sect. ).

4.2.3 Outflows

Throughfall $P_{\mathrm{t}}$ is calculated as

3$$P_{\mathrm{t}}=\left\{\begin{array}{ll}\mathrm{0}& P<(S_{\mathrm{c},\max }-S_{\mathrm{c}})\\ P-(S_{\mathrm{c},\max }-S_{\mathrm{c}})& \text{otherwise}\end{array}\right.,$$ where $S_{\mathrm{c},\max }$ is maximum canopy storage calculated as 4$$S_{\mathrm{c},\max }=m_{\mathrm{c}}\cdot L$$ where $m_{\mathrm{c}}$ is 0.3 $\mathrm{mm}$ , and $L$ $(-)$ is the one-side leaf area index. $L$ is a function of daily temperature and $P$ and limited to minimum or maximum values. Maximum $L$ values per land cover class (Table ) are based on and , whereas minimum $L$ values are calculated as 5$$L_{\min }=\mathrm{0.1}{f}_{\mathrm{d},\mathrm{lc}}+(\mathrm{1}-f_{\mathrm{d},\mathrm{lc}})c_{\mathrm{e},\mathrm{lc}}{L}_{\max },$$ where $f_{\mathrm{d},\mathrm{lc}}$ is the fraction of deciduous plants and $c_{\mathrm{e},\mathrm{lc}}$ is the reduction factor for evergreen plants per land cover type (Table ).

The growing season starts when daily temperature is above 8 ${}^{\circ }\mathrm{C}$ for a land-cover-specific number of $\mathrm{days}$ (Table ) and cumulative precipitation from the day where growing season starts reaches at least 40 $\mathrm{mm}$. In the beginning of the growing season, $L$ increases linearly for 30 $\mathrm{d}$ until it reaches $L_{max}$. For (semi)arid cells, at least 0.5 $\mathrm{mm}$ of daily $P$ is required to keep the growing season on-going. When growing season conditions are not fulfilled anymore, a senescence phase is initiated and $L$ linearly decreases to $L_{\min }$ within the next 30 $\mathrm{d}$ . It is noteworthy that in WaterGAP $L$ only affects the calculation of the canopy water balance. $L$ is not taken into account in computing consumptive water use for irrigated crops (Sect. ) and evapotranspiration from land (Sect. ).

Following , $E_{\mathrm{c}}$ is calculated as 6$$E_{\mathrm{c}}=E_{\mathrm{pot}}{\left({\displaystyle \frac{S_{\mathrm{c}}}{S_{\mathrm{c},\max }}}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}},$$ where $E_{\mathrm{pot}}$ is the potential evapotranspiration ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$) calculated with the Priestley–Taylor equation according to as 7$$E_{\mathrm{pot}}=\mathit{\alpha }\left({\displaystyle \frac{s_{\mathrm{a}}R}{s_{\mathrm{a}}+g}}\right),$$ where, following , $\mathit{\alpha }$ is set to 1.26 in humid and to 1.74 in (semi)arid cells (Appendix ). $R$ is net radiation ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$) that depends on land cover (Table ) (for details on the calculation of net radiation, the reader is referred to ), and $s_{\mathrm{a}}$ is the slope of the saturation vapor pressure–temperature relationship ($\mathrm{kPa}{}^{\circ }{\mathrm{C}}^{-\mathrm{1}}$) defined as 8$$s_{\mathrm{a}}={\displaystyle \frac{\mathrm{4098}\left(\mathrm{0.6108}{e}^{\frac{\mathrm{17.27}T}{T+\mathrm{237.3}}}\right)}{(T+\mathrm{237.3}{)}^{\mathrm{2}}}},$$ where $T$ (${}^{\circ }\mathrm{C}$) is the daily air temperature and $g$ is the psychrometric constant ($\mathrm{k}\mathrm{Pa}{}^{\circ }{\mathrm{C}}^{-\mathrm{1}}$). The latter is defined as 9$$g={\displaystyle \frac{\mathrm{0.0016286}{p}_{\mathrm{a}}}{l_{\mathrm{h}}}},$$ where $p_{\mathrm{a}}$ is atmospheric pressure of the standard atmosphere (101.3 $\mathrm{kPa}$), and $l_{\mathrm{h}}$ is latent heat ($\mathrm{MJ}{\mathrm{kg}}^{-\mathrm{1}}$). Latent heat is calculated as 10$$l_{\mathrm{h}}=\left\{\begin{array}{ll}\mathrm{2.501}-\mathrm{0.002361}T& \text{if}T>\mathrm{0}\\ \mathrm{2.501}+\mathrm{0.334}& \text{otherwise}\end{array}\right..$$

4.3 Snow

To simulate snow dynamics, each 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ grid cell is spatially disaggregated into 100 non-localized subcells that are assigned different land surface elevations according to GTOPO30 . Daily temperature at each subcell is calculated from daily temperature at the 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ cell by applying an adiabatic lapse rate of 0.6 ${}^{\circ }\mathrm{C}$ per 100 m . The daily snow water balance is computed for each of the subcells such that within a 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ cell there may be subcells with and without snow cover or snowfall. For model output, subcell values are aggregated to 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ cell values.

4.3.1 Water balance

Snow storage accumulates below snow freeze temperature and decreases by snow melt and sublimation. Snow storage $S_{\mathrm{sn}}$ ($\mathrm{mm}$) is calculated as

11$${\displaystyle \frac{\mathrm{d}{S}_{\mathrm{sn}}}{\mathrm{d}t}}=P_{\mathrm{sn}}-M-E_{\mathrm{sn}},$$ where $P_{\mathrm{sn}}$ is the part of $P_{\mathrm{t}}$ that falls as snow ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$), $M$ is snowmelt ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$) and $E_{\mathrm{sn}}$ is sublimation ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$).

4.3.2 Inflows

Snowfall $P_{\mathrm{sn}}$ ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$) is calculated as

12$$P_{\mathrm{sn}}=\left\{\begin{array}{ll}{P}_{\mathrm{t}}& T<T_{\mathrm{f}}\\ \mathrm{0}& \text{otherwise}\end{array}\right.,$$ where $T$ is daily air temperature (${}^{\circ }\mathrm{C}$), and $T_{\mathrm{f}}$ is snow freeze temperature, set to 0 ${}^{\circ }\mathrm{C}$. In order to prevent excessive snow accumulation, when snow storage $S_{\mathrm{sn}}$ reaches 1000 mm in a subcell, the temperature in this subcell is increased to the temperature in the highest subcell with a temperature above $T_{\mathrm{f}}$ .

4.3.3 Outflows

Snow melt $M$ is calculated with a land-cover-specific degree-day factor $D_{\mathrm{F}}$ ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}{}^{\circ }\mathrm{C}$) (Table ) when the temperature $T$ in a subgrid surpasses melting temperature $T_{\mathrm{m}}$ (${}^{\circ }\mathrm{C}$), set to 0 ${}^{\circ }\mathrm{C}$, as

13$$M=\left\{\begin{array}{ll}{D}_{\mathrm{F}}(T-T_{\mathrm{m}})& T>T_{\mathrm{m}},S_{\mathrm{sn}}>\mathrm{0}\\ \mathrm{0}& \text{otherwise}\end{array}\right..$$

Sublimation $E_{\mathrm{sn}}$ is calculated as the fraction of $E_{\mathrm{pot}}$ that remains available after $E_{\mathrm{c}}$. For calculating $E_{\mathrm{pot}}$ according to Eq. (), land-cover-specific albedo values are used if $S_{\mathrm{sn}}$ surpasses 3 mm in the 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ cell (Table ). 14$$E_{\mathrm{sn}}=\left\{\begin{array}{ll}{E}_{\mathrm{pot}}-E_{\mathrm{c}}& E_{\mathrm{pot}}-E_{\mathrm{c}}>E_{\mathrm{sn}}\\ S_{\mathrm{sn}}& \text{otherwise}\end{array}\right.$$

4.4 Soil

WaterGAP represents soil as a one-layer soil water storage compartment characterized by a land-cover- and soil-specific maximum storage capacity as well as soil texture. The simulated water storage represents soil moisture in the effective root zone.

4.4.1 Water balance

The change of soil water storage $S_{\mathrm{s}}$ ($\mathrm{mm}$) over time ($\mathrm{d}$) is calculated as

15$${\displaystyle \frac{\mathrm{d}{S}_{\mathrm{s}}}{\mathrm{d}t}}=P_{\mathrm{eff}}-R_{\mathrm{l}}-E_{\mathrm{s}},$$ where $P_{\mathrm{eff}}$ is effective precipitation ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$), $R_{\mathrm{l}}$ is runoff from land ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$) and $E_{\mathrm{s}}$ is actual evapotranspiration from the soil ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$). Once the water balance is computed, $R_{\mathrm{l}}$ is partitioned into (1) fast surface and subsurface runoff $R_{\mathrm{s}}$, representing direct surface runoff and interflow, and (2) groundwater recharge $R_{\mathrm{g}}$ (Fig. ) according to a heuristic scheme .

4.4.2 Inflows

$P_{\mathrm{eff}}$ is computed as

16$$P_{\mathrm{eff}}=P_{\mathrm{t}}-P_{\mathrm{sn}}+M,$$ where $P_{\mathrm{t}}$ is throughfall ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$; see Eq. ), $P_{\mathrm{sn}}$ is snowfall ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$; see Eq. ) and $M$ is snowmelt ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$; see Eq. ).

4.4.3 Outflows

$E_{\mathrm{s}}$ is calculated as

17$$E_{\mathrm{s}}=min\left((E_{\mathrm{pot}}-E_{\mathrm{c}}),(E_{\mathrm{pot},\max }-E_{\mathrm{c}}){\displaystyle \frac{S_{\mathrm{s}}}{S_{\mathrm{s},\max }}}\right),$$ where $E_{\mathrm{pot}}$ is potential evapotranspiration ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$), $E_{\mathrm{c}}$ is canopy evaporation ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$; Eq. ) and $S_{\mathrm{s},\max }$ is the maximum soil water content ($\mathrm{mm}$) derived as a product of total available water capacity in the upper meter of the soil and land-cover-specific rooting depth (Table ) . $E_{\mathrm{pot},\max }$ is set to 15 $\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$ globally. Following , runoff from land $R_{\mathrm{l}}$ is calculated as 18$$R_{\mathrm{l}}=P_{\mathrm{eff}}{\left({\displaystyle \frac{S_{\mathrm{s}}}{S_{\mathrm{s},\max }}}\right)}^{\mathit{\gamma }},$$ where $\mathit{\gamma }$ is the runoff coefficient (–). This parameter, which varies between 0.1 and 5.0, is used for calibration (Sect. ). Together with soil saturation, it determines the fraction of $P_{\mathrm{eff}}$ that becomes $R_{\mathrm{l}}$ (Fig. ). If the sum of $P_{\mathrm{eff}}$ and $S_{\mathrm{s}}$ of the previous day exceed $S_{\mathrm{s},\max }$, the exceeding fraction of $P_{\mathrm{eff}}$ is added to $R_{\mathrm{l}}$. In urban areas (defined from MODIS data, Sect. ), 50 % of $P_{\mathrm{eff}}$ is directly turned into $R_{\mathrm{l}}$.

Figure 3

Relation between runoff from land $R_{\mathrm{l}}$ as a fraction of effective precipitation $P_{\mathrm{eff}}$ and soil saturation $S_{\mathrm{s}}/S_{\mathrm{s},\max }$ for different values of the runoff coefficient $\mathit{\gamma }$ in WaterGAP.

[Figure omitted. See PDF]

$R_{\mathrm{l}}$ is partitioned into fast surface and subsurface runoff $R_{\mathrm{s}}$ and diffuse groundwater recharge $R_{\mathrm{g}}$ calculated as 19$$R_{\mathrm{g}}=min(R_{{\mathrm{g}}_{\max }},f_{\mathrm{g}}{R}_{\mathrm{l}}),$$ where $R_{{\mathrm{g}}_{\max }}$ is soil-texture-specific maximum groundwater recharge with values of 7, 4.5 and 2.5 $\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$ for sandy, loamy and clayey soils, respectively, and $f_{\mathrm{g}}$ is the groundwater recharge factor ranging between 0 and 1. $f_{\mathrm{g}}$ is determined based on relief, soil texture, aquifer type, and the existence of permafrost or glaciers . If a grid cell is defined as (semi)arid and has coarse (sandy) soil, groundwater recharge will only occur if precipitation exceeds a critical value of 12.5 $\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$, otherwise the water remains in the soil. The fraction of $R_{\mathrm{l}}$ that does not recharge the groundwater becomes $R_{\mathrm{s}}$, which recharges surface water bodies and the river compartment.

4.5 Groundwater

As there is no knowledge about the depth below the land surface where groundwater no longer occurs due to the lack of pore space, groundwater storage can only be computed in relative terms but is assumed to be unlimited. The groundwater storage $S_{\mathrm{g}}$ is always positive unless net abstractions from groundwater NA${}_{\mathrm{g}}$ are high and groundwater depletion occurs. Groundwater discharge is assumed to be proportional to (positive) $S_{\mathrm{g}}$ and to stop in the case of negative $S_{\mathrm{g}}$.

4.5.1 Water balance

The temporal development of groundwater storage $S_{\mathrm{g}}$ (${\mathrm{m}}^{\mathrm{3}}$) is calculated as

20$${\displaystyle \frac{\mathrm{d}{S}_{\mathrm{g}}}{\mathrm{d}t}}=R_{\mathrm{g}}+R_{{\mathrm{g}}_{\mathrm{l},\mathrm{res},\mathrm{w}}}-Q_{\mathrm{g}}-{\mathrm{NA}}_{\mathrm{g}},$$ where $R_{\mathrm{g}}$ is diffuse groundwater recharge from soil (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$, Eq. ), $R_{{\mathrm{g}}_{\mathrm{l},\mathrm{res},\mathrm{w}}}$ is point groundwater recharge from surface water bodies (lakes, reservoirs and wetlands) in (semi)arid areas (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$, Eq. ), $Q_{\mathrm{g}}$ is groundwater discharge (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) and NA${}_{\mathrm{g}}$ is net abstraction from groundwater (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$).

4.5.2 Inflows

$R_{\mathrm{g}}$ is the main inflow in most grid cells, except in (semi)arid grid cells with significant surface water bodies where $R_{{\mathrm{g}}_{\mathrm{l},\mathrm{res},\mathrm{w}}}$ may be dominant. $R_{{\mathrm{g}}_{\mathrm{l},\mathrm{res},\mathrm{w}}}$ varies temporally with the area of the surface water body, which depends on the respective water storage (Sect. ). In many cells with significant irrigation with surface water, NA${}_{\mathrm{g}}$ is negative, and irrigation causes a net inflow into the groundwater due to high return flows (Sect. ).

4.5.3 Outflows

$Q_{\mathrm{g}}$ quantifies the discharge from groundwater storage to surface water storage, with

21$$Q_{\mathrm{g}}=k_{\mathrm{g}}{S}_{\mathrm{g}},$$ where $k_{\mathrm{g}}=\mathrm{0.01}$ ${\mathrm{d}}^{-\mathrm{1}}$ is the globally constant groundwater discharge coefficient . The second outflow component NA${}_{\mathrm{g}}$ is described in Sect. .

4.6 Lakes, man-made reservoirs and wetlands

Where lakes, man-made reservoirs and wetlands (LResWs) of significant size exist, their water balances strongly affect the overall water balance of the grid cell due to their high evaporation and water retention capacity . WGHM uses the Global Lakes and Wetland Database (GLWD) and a preliminary but updated version of the Global Reservoir and Dam (GRanD) database to define location, area and other attributes of LResWs. It is assumed that surface areas given in the databases represent the maximum extent. Appendix describes how the information from these databases is integrated into WGHM. Two categories of LResWs are defined for WGHM, so-called “local” water bodies that receive inflow only from the runoff generated within the grid cell and so-called “global” water bodies that additionally receive the streamflow from the upstream grid cells (Fig. ). Six different LResW types are distinguished in WaterGAP.

• Local wetlands (wl) and global wetlands (wg). These cover a maximum area of 3.743 million ${\mathrm{km}}^{\mathrm{2}}$ and 3.752 million ${\mathrm{km}}^{\mathrm{2}}$, respectively, an area that is at its maximum at least 3 times larger than the combined maximum area of lakes and reservoirs (Appendix ). However, 0.3 million ${\mathrm{km}}^{\mathrm{2}}$ of floodplains along large rivers is included as global wetlands, and their dynamics are not simulated suitably by WGHM. They are assumed to receive the total streamflow as inflow while in reality only the part of the streamflow that does not fit in the river channel flows into the floodplain . All local (global) wetlands within a 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ grid cell are simulated as one local (global) wetland that covers a specified fraction of the cell.

• Local lakes (ll). These include about 250 000 small lakes and more than 5000 man-made reservoirs and are defined to have a surface area of less than 100 ${\mathrm{km}}^{\mathrm{2}}$ or a maximum storage capacity of less than 0.5 ${\mathrm{km}}^{\mathrm{3}}$. Like wetlands, all local lakes in a grid cell are aggregated and simulated as one storage compartment taking up a fraction of the grid cell area. Small reservoirs are simulated like lakes as (1) the required lumping of all local reservoirs within a grid cell into one local reservoir per cell necessarily leads to a “blurring” of the specific reservoir characteristics, and (2) small reservoirs are likely not on the main river simulated in the grid cell but on a tributary. Therefore, a reservoir algorithm is not expected to simulate water storage and flows better than the lake algorithm.

• 1355 global lakes (lg). These consist of lakes with an area of more than 100 ${\mathrm{km}}^{\mathrm{2}}$, are simulated in WaterGAP. Since a global lake may spread over more than one grid cell, the water balance of the whole lake is computed at the outflow cell (for consequences, see Sect. ). Only the maximum area of natural lakes is known, not the maximum water storage capacity.

• Global man-made reservoirs (res) and global regulated lakes. Global man-made reservoirs have a maximum storage capacity of at least 0.5 ${\mathrm{km}}^{\mathrm{3}}$, and global regulated lakes (lakes where outflow is controlled by a dam or weir) have a maximum storage capacity of at least 0.5 ${\mathrm{km}}^{\mathrm{3}}$ or an area of more than 100 ${\mathrm{km}}^{\mathrm{2}}$. Both are simulated by the same water balance equation. There can be only one global reservoir/regulated lake compartment per grid cell. Outflow from reservoirs/regulated lakes is simulated by a modified version of the algorithm, distinguishing reservoirs/regulated lakes with the main purpose of irrigation from others . Like in the case of global lakes, water balance of global reservoirs/regulated lakes is computed at the outflow cell (for consequences; see Sect. ). Different from lakes, information on maximum water storage capacity is available from the GRanD database, in addition to the main use and the commissioning year. In WGHM, reservoirs start filling at the beginning of the commissioning year, and regulated lakes then turn from global lakes into global regulated lakes (Appendix ). A total of 1082 global reservoirs and 85 regulated lakes are taken into account, but as those that have the same outflow cell are aggregated to one water storage compartment by adding maximum storages and areas, only 1109 global reservoirs/regulated lakes compartments are simulated in WGHM (Appendix ). Under naturalized conditions (Sect. ), there are no global man-made reservoirs, and regulated lakes are simulated as global lakes; however, local reservoirs remain in the model.

In each grid cell, there can be a maximum of one local wetland storage compartment, one global wetland compartment, one local lake compartment, one global lake compartment and one global reservoir/regulated lake compartment. The lateral water flow within the cell follows the sequence shown in Fig. . For example, if there is a local lake compartment in a grid cell, it is this compartment that receives, under a humid climate, a fraction of the outflow from the groundwater compartment and of the fast surface and subsurface outflow, and the outflow from the local lake becomes inflow to the local wetland if it exists (Fig. ). If there is no local wetland but a global lake, the outflow from the local lake becomes part of the inflow of the global lake. In the case of having a global lake and a global reservoir/regulated lake in one cell, water is routed first through the global lake.

4.6.1 Water balance

The water balance for the five types of LResW compartments is calculated as

22$${\displaystyle \frac{\mathrm{d}{S}_{\mathrm{l},\mathrm{res},\mathrm{w}}}{\mathrm{d}t}}=Q_{\mathrm{in}}+A(P-E_{\mathrm{pot}})-R_{{\mathrm{g}}_{\mathrm{l},\mathrm{res},\mathrm{w}}}-{\mathrm{NA}}_{\mathrm{l},\mathrm{res}}-Q_{\mathrm{out}},$$ where $S_{\mathrm{l},\mathrm{res},\mathrm{w}}$ is volume of water stored in the water body (${\mathrm{m}}^{\mathrm{3}}$), $Q_{\mathrm{in}}$ is inflow into the water body from upstream (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$), $A$ is global (or local) water body surface area (${\mathrm{m}}^{\mathrm{2}}$) in the grid cell at time $t$, $P$ is precipitation (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$), $E_{\mathrm{pot}}$ is potential evapotranspiration (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$, Eq. ), $R_{{\mathrm{g}}_{\mathrm{l},\mathrm{res},\mathrm{w}}}$ is groundwater recharge from the water body (only in arid/semiarid regions) (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$, Eq. ), NA${}_{\mathrm{l},\mathrm{res}}$ is the net abstraction from the lakes and reservoirs (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) (Fig.  and Sect. ), and $Q_{\mathrm{out}}$ is outflow from the water body to other surface water bodies including river storage (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) (Fig. ).

The temporally varying surface area $A$ of the water body is computed in each daily time step using the following equation: 23$$A=r\cdot A_{max},$$ where $r$ is reduction factor (–), and $A_{max}$ is maximum extent of the water body (${\mathrm{m}}^{\mathrm{2}}$) from GRanD or GLWD databases. In the case of local and global lakes 24$$r=\mathrm{1}-{\left({\displaystyle \frac{|S_{\mathrm{l}}-S_{\mathrm{l},\max }|}{\mathrm{2}{S}_{\mathrm{l},\max }}}\right)}^p,\mathrm{0}\le r\le \mathrm{1},$$ where $S_{\mathrm{l}}$ is the volume of the water (${\mathrm{m}}^{\mathrm{3}}$) stored in the lake at time $t$ ($\mathrm{d}$), $S_{\mathrm{l},\max }$ is the maximum storage of the lake (${\mathrm{m}}^{\mathrm{3}}$), $S_{\mathrm{l},\max }$ is computed based on $A_{\max }$ and a maximum storage depth of 5 $\mathrm{m}$, and $p$ is the reduction exponent (–), set to 3.32. According to the above equation, the area is reduced by 1 % if $S_{\mathrm{l}}=\mathrm{50}$ % of $S_{\mathrm{l},\max }$, by 10 % if $S_{\mathrm{l}}=\mathrm{0}$ and by 100 % if $S_{\mathrm{l}}=-S_{\mathrm{l},\max }$ . In the case of global reservoirs/regulated lakes and local and global wetlands 25$$r=\mathrm{1}-{\left({\displaystyle \frac{|S_{\mathrm{res},\mathrm{w}}-S_{\mathrm{res},\mathrm{w},\max }|}{S_{\mathrm{res},\mathrm{w},\max }}}\right)}^p,\mathrm{0}\le r\le \mathrm{1},$$ where $S_{\mathrm{res},\mathrm{w}}$ is the volume of the water (${\mathrm{m}}^{\mathrm{3}}$) stored in the reservoir/regulated lake or wetland, and $p$ is 2.814 and 3.32 for reservoirs/regulated lakes and wetlands, respectively. In the case of wetlands, $S_{\mathrm{res},\mathrm{w},\max }$ (${\mathrm{m}}^{\mathrm{3}}$) is computed based on $A_{\max }$ and a maximum storage depth of 2 $\mathrm{m}$. Wetland area is reduced by 10 % if $S_{\mathrm{w}}=\mathrm{50}$ % of $S_{\mathrm{res},\mathrm{w},\max }$ and by 70 % if $S_{\mathrm{w}}$ is only 10 % of $S_{\mathrm{res},\mathrm{w},\max }$. In the case of reservoirs/regulated lakes, storage capacity $S_{\mathrm{res},\mathrm{w},\max }$ is taken from the database. Reservoir area is reduced by 15 % if $S_{\mathrm{res}}$ is 50 % of $S_{\mathrm{res},\mathrm{w},\max }$ and by 75 % if $S_{\mathrm{res}}$ is only 10 % of $S_{\mathrm{res},\mathrm{w},\max }$. For regulated lakes without available maximum storage capacity, $S_{\mathrm{res},\mathrm{w},\max }$ is computed as in the case of global lakes.

While storage in reservoirs/regulated lakes and wetlands cannot drop below zero due to high outflows, high evaporation or NA${}_{\mathrm{s}}$, storage in lakes can become negative. This represents the situation where there is no more outflow from the lake to a downstream water body ($Q_{\mathrm{out}}=\mathrm{0}$). There, like groundwater storage, storage of local and global lakes is a relative and not an absolute water storage. Reservoir/regulated lake storage is not allowed to fall below 10 % of storage capacity.

With changing $A$ of the surface water compartments local wetland, global wetlands and local lakes, the land area fraction is adjusted accordingly. However, in the case of global lakes and reservoirs/regulated lakes, which may cover more than one 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ cell, such an adjustment is not made as it is not known in which grid cells the area reduction occurs. Therefore, land area fraction is not adjusted with changing $r$ and precipitation is assumed to fall on a surface water body with an area of $A_{\max }$ instead of $A$.

4.6.2 Inflows

Calculation of $Q_{\mathrm{in}}$ differs between local and global water bodies. In the case of local lakes and local wetlands, they are recharged only by local runoff generated within the same grid cell. A fraction $f_{\mathrm{swb}}$ of the fast surface and subsurface runoff generated within the grid cell $R_{\mathrm{s}}$ (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) and, only in the case of humid grid cells, a fraction $f_{\mathrm{swb}}$ of the base flow from groundwater $Q_{\mathrm{g}}$ (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) become inflow to local water bodies (Fig. , Sect. , ). In the case where one grid cell contains both local lake and wetland, then the outflow of the local lake will be the inflow to the local wetland according to Fig. . Global lakes, global wetlands, and global reservoirs/regulated lakes receive, in addition to local runoff, inflow from streamflow of the upstream grid cells as river inflow (Fig. ). In many cells with significant groundwater abstraction, NA${}_{\mathrm{s}}$ is negative, and return flow leads to a net inflow into surface water bodies (Sect. ).

4.6.3 Outflows

LResWs lose water by evaporation $E_{\mathrm{pot}}$, which is assumed to be equal to the potential evapotranspiration computed using the Priestley–Taylor equation with an albedo of 0.08 according to Eq. (). In semiarid and arid grid cells (Appendix ), LResWs are assumed to recharge the groundwater with a focused groundwater recharge, $R_{{\mathrm{g}}_{\mathrm{l},\mathrm{res},\mathrm{w}}}$ with

26$$R_{{\mathrm{g}}_{\mathrm{l},\mathrm{res},\mathrm{w}}}=K_{{\mathrm{gw}}_{\mathrm{l},\mathrm{res},\mathrm{w}}}\cdot r\cdot A_{\max },$$ where $K_{{\mathrm{gw}}_{\mathrm{l},\mathrm{res},\mathrm{w}}}$ is the groundwater recharge constant below LResWs ($=\mathrm{0}$.01 $\mathrm{m}{\mathrm{d}}^{-\mathrm{1}}$). This process is applied only in the arid and semiarid grid cells, as in humid areas groundwater mostly recharges the surface water bodies as explained in Sect.  .

It is assumed that water can be abstracted from lakes and reservoirs but not from wetlands. An amount of NA${}_{\mathrm{l},\mathrm{res}}$ (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) is the net abstractions from lakes and reservoirs, which depends on the total unsatisfied water use Rem${}_{\mathrm{use}}$ and the water storage in the surface water compartment. In the case of a global lake and a reservoir within the same cell, NA${}_{\mathrm{l},\mathrm{res}}$ is distributed equally. In a reservoir, abstraction is only allowed until water storage reaches 10 % of storage capacity (after fulfilling $E$ and $R_{{\mathrm{g}}_{\mathrm{l},\mathrm{res}}}$ ). Outflow from LResWs to downstream water bodies including river storage (Fig. ) is calculated as a function of LResW water storage. The principal effect of a lake or wetland is to reduce the variability of streamflow, which can be simulated by computing outflow $Q_{\mathrm{out}}$ as 27$$Q_{\mathrm{out}}=k\cdot S_{\mathrm{ll},\mathrm{wl}}\cdot {\left({\displaystyle \frac{S_{\mathrm{ll},\mathrm{wl}}}{S_{\mathrm{ll},\mathrm{wl},\max }}}\right)}^a,$$ where $S_{\mathrm{ll},\mathrm{wl}}$ is the local lake or local wetland storage (${\mathrm{m}}^{\mathrm{3}}$), and $k$ is the surface water outflow coefficient ($=\mathrm{0}$.01 ${\mathrm{d}}^{-\mathrm{1}}$). $S_{\mathrm{ll},\mathrm{wl},\max }$ (${\mathrm{m}}^{\mathrm{3}}$) is computed based on $A_{\max }$ and a maximum storage depth of 2 m for local lakes and 5 m for local wetlands. The exponent $a$ is set to 1.5 in the case of local lakes, based on the theoretical value of outflow over a rectangular weir, while the exponent of 2.5 used for local wetlands leads to a slower outflow . The outflow of global lakes and global wetlands is computed as 28$$Q_{\mathrm{out}}=k\cdot S_{\lg ,\mathrm{wg}}.$$

Different from the commissioning year of a reservoir, which is the year the dam was finalized (Appendix ), the operational year of each reservoir is the 12-month period for which reservoir management is defined. It starts with the first month with a naturalized mean monthly streamflow that is lower than the annual mean. To compute daily outflow, e.g., release, from global reservoirs/regulated lakes, the total annual outflow during the reservoir-specific operational year is determined first as a function of reservoir storage at the beginning of the operational year. Total annual outflow during the operational year is assumed to be equal to the product of mean annual outflow and a reservoir release factor $k_{\mathrm{rele}}$ that is computed each year on the first day of the operational year as 29$$k_{\mathrm{rele}}={\displaystyle \frac{S_{\mathrm{res}}}{S_{\mathrm{res},\max }\cdot \mathrm{0.85}}},$$ where $S_{\mathrm{res}}$ is the reservoir/regulated lake storage (${\mathrm{m}}^{\mathrm{3}}$), and $S_{\mathrm{res},\max }$ is the storage capacity (${\mathrm{m}}^{\mathrm{3}}$). Thus, total release in an operational year with low reservoir storage at the beginning of the operational year will be smaller than in a year with high reservoir storage.

During the first filling phase of a reservoir after dam construction, $k_{\mathrm{rele}}$ = 0.1 until $S_{\mathrm{res}}$ exceeds 10 % of $S_{\mathrm{res},\max }$. If the storage capacity to mean total annual outflow ratio is larger than 0.5, then the outflow from the reservoir is independent of the actual inflow and temporally constant in the case of a non-irrigation reservoir. In the case of an irrigation reservoir, outflow is driven by monthly NA${}_{\mathrm{s}}$ in the next five downstream cells or down to the next reservoir . For reservoirs with a smaller ratio, the release additionally depends on daily inflow and is higher on days with high inflow . If reservoir storage drops below 10 % of $S_{\mathrm{res},\max }$, release is reduced to 10 % of the normal release to satisfy a minimum environmental flow requirement for ecosystems. Daily outflow may also include overflow, which occurs if reservoir storage capacity is exceeded due to high inflow into the reservoir.

4.7 Rivers

The water balance of the river compartment is computed to quantify streamflow, one of the most important output variables of hydrological models.

4.7.1 Water balance

The dynamic water balance of the river water storage in a cell is computed as

30$${\displaystyle \frac{\mathrm{d}{S}_{\mathrm{r}}}{\mathrm{d}t}}=Q_{\mathrm{r},\mathrm{in}}-Q_{\mathrm{r},\mathrm{out}}-{\mathrm{NA}}_{\mathrm{s},\mathrm{r}},$$ where $S_{\mathrm{r}}$ is the volume of water stored in the river (${\mathrm{m}}^{\mathrm{3}}$), $Q_{\mathrm{r},\mathrm{in}}$ is inflow into the river compartment (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$), $Q_{\mathrm{r},\mathrm{out}}$ is the streamflow (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) and NA${}_{\mathrm{s},\mathrm{r}}$ is the net abstraction of surface water from the river (${\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$).

4.7.2 Inflows

If there are no surface water bodies in a grid cell, $Q_{\mathrm{r},\mathrm{in}}$ is the sum of $R_{\mathrm{s}}$, $Q_{\mathrm{g}}$ and streamflow from existing upstream cell(s). Otherwise, part of $R_{\mathrm{s}}$, and in the case of humid cells also part of $Q_{\mathrm{g}}$, is routed through the surface water bodies (Fig. ). The outflow from the surface water body preceding the river compartment then becomes part of $Q_{\mathrm{r},\mathrm{in}}$. In addition, negative NA${}_{\mathrm{s}}$ values due to high return flows from irrigation with groundwater lead to a net increase in storage. Thus, if no surface water bodies exist in the cell, negative NA${}_{\mathrm{s}}$ is added to $Q_{\mathrm{r},\mathrm{in}}$ (Sect.  and Fig. ).

4.7.3 Outflows

$Q_{\mathrm{r},\mathrm{out}}$ is defined as the streamflow that leaves the cell and is transferred to the downstream cell.

It is calculated as

31$$Q_{\mathrm{r},\mathrm{out}}={\displaystyle \frac{v}{l}}\cdot S_{\mathrm{r}},$$ where $v$ ($\mathrm{m}{\mathrm{d}}^{-\mathrm{1}}$) is river flow velocity, and $l$ is the river length ($\mathrm{m}$). $l$ is calculated as the product of the cell's river segment length, derived from the HydroSHEDS drainage direction map , and a meandering ratio specific to that cell (method described in ). $v$ is calculated according to the Manning–Strickler equation as 32$$v=n^{-\mathrm{1}}\cdot R_{\mathrm{h}}^{\frac{\mathrm{2}}{\mathrm{3}}}\cdot s^{\frac{\mathrm{1}}{\mathrm{2}}},$$ where $n$ is river bed roughness (–), $R_{\mathrm{h}}$ is the hydraulic radius of the river channel ($\mathrm{m}$) and $s$ is river bed slope ($\mathrm{m}{\mathrm{m}}^{-\mathrm{1}}$). Calculation of $s$ is based on high-resolution elevation data (SRTM30), the HydroSHEDS drainage direction map and an individual meandering ratio. The predefined minimum $s$ is 0.0001 $\mathrm{m}{\mathrm{m}}^{-\mathrm{1}}$.

To compute the daily varying $R_{\mathrm{h}}$, a trapezoidal river cross section with a slope of 0.5 is assumed such that it can be calculated as a function of daily varying river depth $D_{\mathrm{r}}$ and temporally constant bottom width $W_{\mathrm{r},\mathrm{bottom}}$ . empirically derived equations relating river depth, river top width and streamflow for bankfull conditions. In former model versions, these equations were also applied at each time step, even if streamflow was not bankfull, to determine river width and depth required to compute $R_{\mathrm{h}}$ and thus $v$. As usage of these functions for any streamflow below bankfull is not backed by the data and method of , WaterGAP 2.2d implements a consistent method for determining daily width and depth as a function of river water storage.

As bankfull conditions are assumed to occur at the initial time step, the initial volume of water stored in the river is computed as 33$$S_{\mathrm{r},\max }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot l\cdot D_{\mathrm{r},\mathrm{bf}}\cdot (W_{\mathrm{r},\mathrm{bottom}}+W_{\mathrm{r},\mathrm{bf}}),$$ where $S_{\mathrm{r},\max }$ is the maximum volume of water that can be stored in the river at bankfull depth (${\mathrm{m}}^{\mathrm{3}}$), $D_{\mathrm{r},\mathrm{bf}}$ ($\mathrm{m}$) and $W_{\mathrm{r},\mathrm{bf}}$ ($\mathrm{m}$) are river depth and top width at bankfull conditions, respectively, and $W_{\mathrm{r},\mathrm{bottom}}$ is river bottom width ($\mathrm{m}$). River water depth $D_{\mathrm{r}}$ ($\mathrm{m}$) is simulated to change at each time step with actual $S_{\mathrm{r}}$ as 34$$D_{\mathrm{r}}=-{\displaystyle \frac{W_{\mathrm{r},\mathrm{bottom}}}{\mathrm{4}}}+\sqrt{W_{\mathrm{r},\mathrm{bottom}}\cdot {\displaystyle \frac{W_{\mathrm{r},\mathrm{bottom}}}{\mathrm{16}}}+\mathrm{0.5}\cdot {\displaystyle \frac{S_{\mathrm{r}}}{l}}}.$$ Using the equation for a trapezoid with a slope of 0.5, $R_{\mathrm{h}}$ is then calculated from $W_{\mathrm{r},\mathrm{bottom}}$ and $D_{\mathrm{r}}$. Bankfull flow is assumed to correspond to the maximum annual daily flow with a return period of 1.5 years and is derived from daily streamflow time series.

The roughness coefficient $n$ of each grid cell is calculated according to , who modeled $n$ as a function of various spatial characteristics (e.g., urban or rural area, vegetation in river bed, obstructions) and a river sinuosity factor to achieve an optimal fit to streamflow observations. Because of the implementation of a new algorithm to calculate $D_{\mathrm{r}}$, we had to adjust their gridded $n$ values to avoid excessively high river velocities . By trial and error, we determined optimal $n$-multipliers at the scale of 13 large river basins that lead to a good fit to monthly streamflow time series at the most downstream stations and basin-average total water storage anomalies from GRACE. We found that in 9 out of 13 basins, multiplying $n$ by 3 resulted in the best fit between observed and modeled data. We therefore set the multiplier to 3 globally, except for the remaining four basins, where other values proved to be more adequate; this concerns the Lena basin, where $n$ is multiplied by 2; the Amazon basin, where $n$ is multiplied by 10; and the Huang He and Yangtze basins, where $n$ is kept at its original value (Fig. S1).

Net cell runoff $R_{\mathrm{nc}}$ ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$), the part of the cell precipitation that has neither been evapotranspirated nor stored with a time step, is calculated as 35$$R_{\mathrm{nc}}={\displaystyle \frac{(Q_{\mathrm{r},\mathrm{out}}-Q_{\mathrm{r},\mathrm{in}})}{A_{\mathrm{cont}}}}\times {\mathrm{10}}^{\mathrm{9}},$$ where $A_{\mathrm{cont}}$ is the continental area (0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ grid cell area minus ocean area) of the grid cell (${\mathrm{m}}^{\mathrm{2}}$). Renewable water resources are calculated as long-term mean annual $R_{\mathrm{nc}}$ computed under naturalized conditions (Sect. ). Renewable water resources can be negative if evapotranspiration in a grid cell is higher than precipitation due to evapotranspiration from global lakes, reservoirs or wetlands that receive water from upstream cells.

4.8 Abstraction of human water use in WaterGAP Global Hydrological Model

The global water use models (Sect. ) together with GWSWUSE (Sect. ) calculate potential NA${}_{\mathrm{pot},\mathrm{g}}$ and NA${}_{\mathrm{pot},\mathrm{s}}$, which are independent of actual water availability. Potential NA${}_{\mathrm{pot},\mathrm{g}}$ is always satisfied in WGHM due to the assumed unlimited groundwater storage that can be depleted (with the exception described in last paragraph of this section).

Satisfaction of potential NA${}_{\mathrm{pot},\mathrm{s}}$ depends on the availability of water in surface water bodies including the river compartment, considering the abstraction priorities shown in Fig. . If the surface water in a grid cell cannot satisfy potential NA${}_{\mathrm{pot},\mathrm{s}}$ of the grid cell on a certain day, the unsatisfied NA${}_{\mathrm{s}}$ of the demand cell is distributed spatially and temporally to potentially increase the amount of satisfied NA${}_{\mathrm{s}}$. If the demand cell is a riparian cell of a global lake or reservoir, NA${}_{\mathrm{s}}$ can be satisfied from the lake/reservoir storage. Unsatisfied surface water demand of all other cells can be taken from the neighboring cell with the largest river and lake/reservoir storage (“second cell”). In both cases, negative values of consumptive use (sum of NA${}_{\mathrm{s}}$ and NA${}_{\mathrm{g}}$) can occur in the demand cells in case of irrigation with surface water. Here, a negative value of NA${}_{\mathrm{g}}$ in the demand cell may occur in the case of return flows from irrigation, while the positive value of NA${}_{\mathrm{s}}$ is allocated to a neighboring cell. Temporal distribution of unsatisfied NA${}_{\mathrm{s}}$ is achieved by adding it to NA${}_{\mathrm{s}}$ of the next day, but no longer than until the end of the calendar year (“delayed use”). If NA${}_{\text{pot,s}}$ still cannot be fulfilled, actual NA${}_{\mathrm{s}}$ becomes smaller than potential NA${}_{\mathrm{s}}$.

Delayed satisfaction aims at compensating for the fact that WaterGAP likely underestimates the storage of water, e.g., by small tanks and dams, and because of the generic reservoir operation scheme. Without delayed satisfaction, less than 50 % of potential NA${}_{\mathrm{pot},\mathrm{s}}$ could be satisfied in many semiarid regions (Fig. S2). The delayed satisfaction scheme may overestimate satisfaction of surface water demand in particular in highly seasonal flow regimes. However, this effect is hardly visible in the hydrograph of the monsoonal Yangtze River (Fig. S3) but more visible in semiarid regions (Figs. S4, S5). With delayed satisfaction of potential NA${}_{\mathrm{s}}$, 92.5 % of global potential NA${}_{\mathrm{s}}$ during 1981–2010 is satisfied, but only 82.2 % in the case of the alternative option that surface water demand needs to be satisfied by available surface water on the same day.

In the case of irrigation by surface water, it is assumed that any decrease in NA${}_{\mathrm{s}}$ is due to a decrease in withdrawal water uses for irrigation. This also reduces return flow to groundwater. Therefore, in WaterGAP 2.2d, NA${}_{\mathrm{g}}$ is increased in each time step in the water demand cell in accordance with the unfulfilled potential NA${}_{\mathrm{pot},\mathrm{s}}$ in the cell (after steps 1 and 2).

4.9 Calibration and regionalization

4.9.1 Calibration approach

The main purpose of WaterGAP is to quantify water resources and water stress for both historical time periods and scenarios of the future. Not only due to very uncertain global climate input data, uncalibrated global hydrological models may compute very biased runoff and streamflow values (e.g., ). To reduce the bias and simulate at least mean streamflow and thus renewable water resources with a reasonable reliability, WGHM has been calibrated to match observed long-term average annual streamflow at gauging stations on all continents . Calibration is required due to uncertain model parameters, input data (e.g., deviations of precipitation from meteorological forcings to observation networks; ) and model structure including the spatial resolution. The rationale behind the approach can be summed up by the phrase “if the model is not able to properly capture the average observed hydrological conditions, how well founded are future projections?” (see also the discussion in ). In order to minimize the problem of equifinality, WGHM is calibrated in a very simple basin-specific manner to match long-term mean annual observed streamflow ($Q_{\mathrm{obs}}$) at the outlet of 1319 drainage basins that cover $\sim$ 54 % of the global drainage area (except Antarctica and Greenland) (Fig. ). The runoff coefficient $\mathit{\gamma }$ (Eq. ) and up to two additional correction factors (the areal correction factor, CFA, and the station correction factor, CFS; for a brief description the reader is referred to the calibration status CS3 and CS4 below or to ), if needed, are adjusted homogeneously for all grid cells within the drainage basin. Calibration starts in upstream basins and proceeds to downstream basins, with the streamflow from the already calibrated upstream basin as inflow.

While the calibration approach in WaterGAP 2.2d is generally the same as in previous model versions , it was modified Appendix A3 to allow for a $\pm$10 % gauging station observation uncertainty following instead of $\pm$1 % in previous model versions. It is noteworthy that the discharge uncertainty (approximated here with $\pm$10 %) is unlikely to be stationary in space and time , but there are no further data available to better constrain the specific uncertainty of each gauging station. The source of streamflow data and selection criteria for stations is the same as in (their Appendix B2), but the 30-year period was shifted (if available) from 1971–2000 to 1980–2009 to capture a more recent time period.

Calibration follows a four-step scheme with specific calibration status (CS):

• CS1. Adjust the basin-wide uniform parameter $\mathit{\gamma }$ (Eq. ) in the range of [0.1–5.0] to match $Q_{\mathrm{obs}}$ within $\pm$1 %.

• CS2. Adjust $\mathit{\gamma }$ as for CS1, but within 10 % uncertainty range (90 %–110 % of observations).

• CS3. As CS2 but apply the areal correction factor CFA (adjusts runoff and, to conserve the mass balance, actual evapotranspiration as counterpart of each grid cell within the range of [0.5–1.5]) to match $Q_{\mathrm{obs}}$ with 10 % uncertainty.

• CS4. As CS3 but apply the station correction factor CFS (multiplies streamflow in the cell where the gauging station is located by an unconstrained factor) to match $Q_{\mathrm{obs}}$ with 10 % uncertainty to avoid error propagation to the downstream basin. Note that with CFS, actual evapotranspiration of this grid cell is not adapted accordingly to avoid unphysical values. Hence, mass is not conserved in the case of CS4 for the grid cell where CFS is applied in the upstream basin. For global water balance assessment, the mass balance is kept by adjusting the actual evapotranspiration component by the amount CFS modified streamflow.

For each basin, calibration steps 2–4 are only performed if the previous step was not successful.

4.9.2 Regionalization approach

The calibrated $\mathit{\gamma }$ values are regionalized to river basins without sufficient streamflow observations using a multiple linear regression approach that relates the natural logarithm of $\mathit{\gamma }$ to basin descriptors (mean annual temperature, mean available soil water capacity, fraction of local and global lakes and wetlands, mean basin land surface slope, fraction of permanent snow and ice, aquifer-related groundwater recharge factor). Just like the calibrated $\mathit{\gamma }$ values, the regionalized values are limited between 0.1 and 5.0; CFA and CFS are set to 1.0 in uncalibrated basins. A manual modification of the regionalized $\mathit{\gamma }$ value to 0.1 was done (from values of 3–5) for basins covering the North China Plain in northeastern China as groundwater depletion was overestimated by a factor of 4 in this region ; a lower $\mathit{\gamma }$ allows higher runoff generation that translates into higher groundwater recharge and thus a weaker overestimation.

4.9.3 Calibration and regionalization results

Calibration of WaterGAP 2.2d driven by the standard climate forcing (Sect. ) results in 485 basins with calibration status CS1, 185 basins with calibration status CS2, 277 basins with calibration status CS3 and 372 basins with calibration status CS4. This means that in 72 % of the calibration basins, the usage of the station correction factor CFS is not required to match the simulated long-term annual streamflow to observations. The spatial distribution of the calibration parameters and status is shown in Fig. .

Figure 4

Results of the calibration of WaterGAP 2.2d to the standard climate forcing with (a) the calibration status (see Sect. ) of each calibration basin, (b) calibration parameter $\mathit{\gamma }$, (c) areal correction factor CFA and (d) station correction factor CFS. Grey areas in (d) indicate regions with regionalized calibration parameter $\mathit{\gamma }$ and for (a)–(d) dark green outlines indicate the boundaries of the calibration basins.

[Figure omitted. See PDF]

5 Standard model output

5.1 Data provided at PANGAEA repository

A set of standard model outputs is provided via the data publisher and repository PANGAEA hosted by Alfred Wegener Institute, Helmholtz Center for Polar and Marine Research (AWI), Center for Marine Environmental Sciences and University of Bremen (MARUM), under the Creative Commons Attribution-NonCommercial 4.0 International license (CC-BY-NC-4.0). The data are stored using the network Common Data Form (netCDF) format developed by UCAR/Unidata and are available at https://doi.pangaea.de/10.1594/PANGAEA.918447.

The available storages and flows are listed in Table  and Table , respectively. To convert between equivalent water heights (e.w.h.) and volumetric units, the cell-specific continental area used in WaterGAP 2.2d is also provided. The assumed water density is 1 $\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$. The following additional static data used to produce the storages and flows are available: flow direction , land cover (Appendix ), location of outflow cells of global lakes and reservoirs/regulated lakes (Sect. ), rooting depth (Sect. ), maximum soil water storage ($S_{\mathrm{s},\max }$), and reservoir commissioning year (Sect. ). Additionally, the calibration factors $\mathit{\gamma }$, CFA, CFS and the calibration status CS (Sect. ) are provided. The netCDF files contain metadata with detailed information regarding characteristics of the data (e.g., whether a storage type contains anomaly or absolute values) and a legend where applicable.

Table 1

Standard WaterGAP output variables: water storages. Units are $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ ($\mathrm{mm}\mathrm{e}.\mathrm{w}.\mathrm{h}.$). Temporal resolution is monthly.

Storage type	PANGEA file	Symbol
Total water storage${}^{\mathrm{1},\mathrm{2}}$	tws	$S_{\mathrm{tws}}$
Canopy water storage	canopystor	$S_{\mathrm{c}}$
Snow water storage	swe	$S_{\mathrm{sn}}$
Soil water storage	soilmoist	$S_{\mathrm{s}}$
Groundwater storage${}^{\mathrm{2}}$	groundwstor	$S_{\mathrm{g}}$
Local lake storage${}^{\mathrm{2}}$	loclakestor	$S_{\mathrm{ll}}$
Global lake storage${}^{\mathrm{2}}$	glolakestor	$S_{\lg }$
Local wetland storage	locwetlandstor	$S_{\mathrm{wl}}$
Global wetland storage	glowetlandstor	$S_{\mathrm{wg}}$
Reservoir storage	reservoirstor	$S_{\mathrm{res}}$
River storage	riverstor	$S_{\mathrm{r}}$

${}^{\mathrm{1}}$ Sum of all compartments below. ${}^{\mathrm{2}}$ Relative water storages, only anomalies with respect to a reference period can be evaluated.

Table 2

Standard WaterGAP output variables: flows. Units are $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ ($\mathrm{mm}\mathrm{e}.\mathrm{w}.\mathrm{h}.{\mathrm{s}}^{-\mathrm{1}}$), except for $Q_{\mathrm{r},\mathrm{out}}$ and $Q_{\mathrm{r},\mathrm{out},\mathrm{nat}}$, which are in ${\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$. Temporal resolution is monthly.

Flow type	PANGEA file	Symbol
Monthly precipitation	precmon	$P$
Fast surface and fast subsurface runoff${}^{\mathrm{1}}$	qs	$R_{\mathrm{s}}$
Diffuse groundwater recharge	qrdif	$R_{\mathrm{g}}$
Groundwater recharge from surface water bodies	qrswb	$R_{{\mathrm{g}}_{\mathrm{l},\mathrm{res},\mathrm{w}}}$
Total groundwater recharge${}^{\mathrm{2}}$	qr	$R_{{\mathrm{g}}_{\mathrm{tot}}}$
Runoff from land${}^{\mathrm{3}}$	ql	$R_{\mathrm{l}}$
Groundwater discharge${}^{\mathrm{4}}$	qg	$Q_{\mathrm{g}}$
Actual evapotranspiration ${}^{\mathrm{5}}$	evap	$E_{\mathrm{a}}$
Potential evapotranspiration	potevap	$E_{\mathrm{p}}$
Net cell runoff	ncrun	$R_{\mathrm{nc}}$
Naturalized net cell runoff${}^{\mathrm{6}}$	natncrun	$R_{\mathrm{nc},\mathrm{nat}}$
Streamflow${}^{\mathrm{7}}$	dis	$Q_{\mathrm{r},\mathrm{out}}$
Naturalized streamflow${}^{\mathrm{7}}$	natdis	$Q_{\mathrm{r},\mathrm{out},\mathrm{nat}}$
Actual net abstraction from surface water	anas	NA${}_{\mathrm{s}}$
Actual net abstraction from groundwater	anag	NA${}_{\mathrm{g}}$
Actual consumptive water use ${}^{\mathrm{8}}$	atotuse	WC${}_{\mathrm{a}}$

${}^{\mathrm{1}}$ Fraction of total runoff from land that does not recharge the groundwater. ${}^{\mathrm{2}}$ Sum of qrdif and qrswb. ${}^{\mathrm{3}}$ Sum of qs and qrdif. ${}^{\mathrm{4}}$ Groundwater runoff. ${}^{\mathrm{5}}$ Sum of soil evapotranspiration $E_{\mathrm{s}}$, sublimation $E_{\mathrm{sn}}$, evaporation from canopy $E_{\mathrm{c}}$, evaporation from water bodies and actual consumptive water use WC${}_{\mathrm{a}}$. ${}^{\mathrm{6}}$ Equals renewable water resources if averaged over, for example, 30-year time period. ${}^{\mathrm{7}}$ River discharge. ${}^{\mathrm{8}}$ Sum of anas and anag.

5.2 Caveats in usage of WaterGAP model output

Based on feedback from data users and our own experience, here we describe caveats regarding analysis of specific WaterGAP 2.2d model output with the aim of guiding output users.

• WaterGAP does not consider leap years. This implies that model output (typically provided in netCDF file format) corresponding to leap years contains the “fill value” instead of a data value at the position of 29 February.

• The water balance of large lakes and reservoirs is calculated in the outflow cell only. Hence, large numerical values can occur for storages and flows, especially in the case of very large water bodies.

• In the case that the station correction factor CFS (Sect. ) is applied in the grid cell corresponding to the calibration station, multiplication of streamflow by CFS destroys the water balance for this particular grid cell. Hence, the calculation of water balance at various spatial units requires that the amount of reduced/increased streamflow is taken into account in order to close the water balance. A direct inclusion of modified streamflow in, for example, evapotranspiration is not done to avoid physically implausible values for this variable. Water balance is preserved in the case that CFA is used.

• Gridded model output always relates to the continental area (grid cell area minus ocean area within cell). If flows like runoff from land or diffuse groundwater recharge are simulated to occur only on the land area, i.e., the fraction of the continental area that is not covered by surface water bodies, these flow variables can be small in cells with large water bodies, e.g., groundwater recharge along the Amazon river with riparian wetlands (Fig. c).

• Groundwater recharge below surface water bodies (Eq. ) can lead to very high values in the case of large surface water bodies and especially in inland sinks that contain large lakes. Temporal changes of this variable can be implausibly high ($>$ 10${}^{\mathrm{3}}$ $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$).

• Renewable water resources (Fig. a) are defined as the amount of precipitation that is not evapotranspired in the long term (30 years) under naturalized conditions (no water use, no reservoirs). Data users should keep in mind that this variable can only be calculated from naturalized runs and the long-term average of the variable “net cell runoff” $R_{\mathrm{nc},\mathrm{nat}}$ (Table ). A calculation of renewable water resources using other model setups is not meaningful.

• Actual consumptive water use can become negative in those cases where water demand is satisfied by spatially distributed grid cells and in the case of irrigation with surface water (see Sect. ).

6 Model evaluation

This section comprises an evaluation of WaterGAP 2.2d using independent data of withdrawal water uses, streamflow and total water storage anomalies (TWSAs) as well as a comparison to the previous model version 2.2 .

6.1 Model setup and simulation experiments

In order to compare WaterGAP 2.2d with model version 2.2 (Sect. ), both versions were calibrated and run with the same climate forcing. However, version 2.2 was calibrated using the calibration routine of . The differences between model versions 2.2 and 2.2d are listed in Appendix .

A homogenized combination of WATCH Forcing Data based on ERA40 (for 1901–1978) and WATCH Forcing Data methodology applied to ERA-Interim reanalysis (for 1979–2016), with precipitation adjusted to monthly precipitation sums from GPCC , was used. The homogenization method is described in . The calibrated models have been run for the time period 1901–2016, with a spin-up of 5 years in which the model input for 1901 was used.

6.2 Evaluation datasets

6.2.1 AQUASTAT withdrawal water use data

AQUASTAT is the Food and Agriculture Organization of the United Nations Global Information System on Water and Agriculture . It contains information on country-level withdrawal water uses for different sectors. These data represent estimates mainly provided by the individual countries. In particular irrigation withdrawal water uses are, for most countries, not based on observations. Six different withdrawal water use variables (Table ) were available for comparison to WaterGAP 2.2d. For the evaluation, all database entries available on were used; hence it contains yearly values per country as data units. The evaluation metrics (Sect. ) are calculated using each single data point of AQUASTAT without any temporal aggregation by country.

Table 3

AQUASTAT variables used for evaluating WaterGAP 2.2d potential withdrawal water use WU, including variable ID reference of AQUASTAT.

No.	WU variable	Description	AQUASTAT equivalent (variable ID)
1	Total WU	Total WU from all sectors	Total freshwater withdrawal water use (4263)
2	Groundwater WU	As 1 but from groundwater resources only	Fresh groundwater withdrawal water use (4262)
3	Surface water WU	As 1 but from all surface water resources only	Fresh surface withdrawal water use (4261)
4	Irrigation WU	WU for irrigation	Irrigation withdrawal water use (4475)
5	Industrial WU	WU for manufacturing and cooling of thermal power plants	Industrial withdrawal water use (4252)
6	Domestic WU	WU for domestic sector	Municipal withdrawal water use (4251)

6.2.2 GRDC streamflow data

Monthly streamflow time series from 1319 calibration stations from the Global Runoff Data Centre (GRDC) were used for evaluating the performance of WaterGAP 2.2d and 2.2. As the GRDC archive has certain gaps in some regions and times and the calibration objective is to benefit from a maximum of observation data, the typical split-sampling calibration/validation is not appropriate. Even though the same observation data are used for calibration and validation, the validation against monthly time series is meaningful as only long-term mean annual streamflow values have been used for calibration.

6.2.3 GRACE total water storage anomalies

Three mascon solutions of monthly time series of TWSAs from the Gravity Recovery And Climate Experiment (GRACE) satellite mission are considered. The Jet Propulsion Laboratory (JPL) mascon dataset from the GRACE Tellus Website is based on the Level-1 product processed at JPL. A geocenter correction is applied to the degree-1 coefficients following the method from , the c${}_{\mathrm{20}}$ coefficient is replaced with the solutions from satellite laser ranging (SLR; ) and a glacial isostatic adjustment (GIA) correction is applied based on the ICE6G-D model published in . The Center of Space Research (CSR) RL05 GRACE mascon solution from the University of Texas website performs the same degree-1 and c${}_{\mathrm{20}}$ replacements (but following ) and removes the GIA signal based on the model from . Last, the Goddard Space Flight Center (GSFC) GRACE mascon solutions from the Geodesy and Geophysics Science Research Portal applies trend corrections for the $c_{\mathrm{21}}$ and $s_{\mathrm{21}}$ coefficients following in addition to the degree-1, $c_{\mathrm{20}}$ and GIA corrections described for CSR.

Monthly TWSA values are provided on 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ grid cells for JPL and CSR, while GSFC provides equal area grids with a spatial resolution of around 1${}^{\circ }\times \mathrm{1}{}^{\circ }$ at the Equator. In this study, the grid values are spatially averaged over 143 river basins with a total area of more than 200 000 ${\mathrm{km}}^{\mathrm{2}}$ each, out of the 1319 basins used for calibration. The considered time span for this study is 2003–2015 (full years of data), limited by available monthly solutions from GSFC between January 2003 and July 2016.

6.3 Evaluation metrics

6.3.1 Nash–Sutcliffe efficiency

The Nash–Sutcliffe efficiency metric (NSE) (–) is a traditional metric in hydrological modeling. It provides an integrated measure of modeling performance with respect to mean values and variability and is calculated as follows:

36$$\mathrm{NSE}=\mathrm{1}-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(O_i-S_i{)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n(O_i-\stackrel{\mathrm{‾}}{O}{)}^{\mathrm{2}}}},$$ where $O_i$ is the observed value (e.g., monthly streamflow), $S_i$ is the simulated value and $\stackrel{\mathrm{‾}}{O}$ is the mean observed value. The optimal value of NSE is 1. Values below 0 indicate that the mean value of observations is better than the simulation . For assessing the performance of low values of water abstraction (Sect. ), a logarithmic NSE was calculated in addition by applying logarithmic transformation before calculation of the performance indicator.

6.3.2 Kling–Gupta efficiency

The Kling–Gupta efficiency metric (KGE) transparently combines the evaluation of bias, variability and timing and is calculated (in its 2012 version) as follows:

37$$\begin{array}{rl}& \mathrm{KGE}=\mathrm{1}-\\ & \sqrt{({\mathrm{KGE}}_r-\mathrm{1}{)}^{\mathrm{2}}+({\mathrm{KGE}}_b-\mathrm{1}{)}^{\mathrm{2}}+({\mathrm{KGE}}_g-\mathrm{1}{)}^{\mathrm{2}}},\end{array}$$ where KGE${}_r$ is the correlation coefficient between simulated and observed values (–), an indicator for the timing; KGE${}_b$ is the ratio of mean values (Eq. ) (–), an indicator of biases regarding mean values; and KGE${}_g$ is the ratio of variability (Eq. ) (–), an indicator for the variability of simulated ($S$) and observed ($O$) values. 38$${\mathrm{KGE}}_b={\displaystyle \frac{{\mathit{\mu }}_S}{{\mathit{\mu }}_O}},$$ 39$${\mathrm{KGE}}_g={\displaystyle \frac{{\mathrm{CV}}_S}{{\mathrm{CV}}_O}}={\displaystyle \frac{{\mathit{\sigma }}_S/{\mathit{\mu }}_S}{{\mathit{\sigma }}_O/{\mathit{\mu }}_O}},$$ where $\mathit{\mu }$ is mean value, $\mathit{\sigma }$ is standard deviation and CV is coefficient of variation. The optimal value of KGE is 1.

6.3.3 TWSA-related metrics

For the evaluation of total water storage anomaly performance, the following metrics were used: $R^{\mathrm{2}}$ (coefficient of determination) as the strength of the linear relationship between simulated and observed variables, and the amplitude ratio as the indicator for variability and trends of both GRACE and WaterGAP data. Amplitude and trends were determined by a linear regression for estimating the most dominant temporal components of the GRACE time series. The time series of monthly TWSAs was approximated by a constant $a$, a linear trend $b$, and an annual and a semiannual sinusoidal curve as follows:

40$$\begin{array}{rl}y\left(t\right)=& a+b\cdot t+c\cdot \sin (\mathrm{2}\cdot \mathit{\pi }\cdot t)+d\cdot \cos (\mathrm{2}\cdot \mathit{\pi }\cdot t)\\ & +e\cdot \sin (\mathrm{4}\cdot \mathit{\pi }\cdot t)+f\cdot \cos (\mathrm{4}\cdot \mathit{\pi }\cdot t)+r,\end{array}$$ where $r$ denotes the residuals. The parameters $a$ to $f$ were estimated via least-squares adjustment. The annual amplitude can be computed by $A=\mathrm{sqrt}(c^{\mathrm{2}}+d^{\mathrm{2}})$, and thus the annual ratio was calculated by $A_{\mathrm{WGHM}}/A_{\mathrm{GRACE}}$.

6.4 Evaluation results

6.4.1 Water withdrawals

The performance of WaterGAP potential withdrawal water uses is generally of reasonable quality (Fig. , for a non-logarithmic graph see Fig. S6). The highest agreement in terms of performance indicator is shown for the total withdrawal water uses with both efficiency metrics close to the optimum value. Slightly less agreement is visible for the separation into groundwater withdrawals (underestimation by WaterGAP) and surface water withdrawals (overestimation by WaterGAP). The domestic sectoral withdrawal water uses are best simulated with WaterGAP, followed by the industrial sector. Here, large differences between NSE and logarithmic NSE are visible, indicating that WaterGAP has specific problems in representing the small values and tending to a general overestimation of industrial withdrawal water uses. Comparing simulated industrial water uses from WaterGAP with data of the FAO AQUASTAT database reveals inconsistencies due to overestimation (i.e., for values $>$ 200 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$) as well as underestimation (i.e., for small values) (Fig.  and Fig. S6). In terms of overestimated values, values for India and Germany dominate the differences in the time intervals 2008–2012 and 2013–2016, respectively. Water withdrawals of 56 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ for the industry sector (including thermoelectric) were assessed by India's National Commission on Integrated Water Resources Development for 2010 . Here AQUASTAT reports 17 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ and WaterGAP simulates 72 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$. In the case of Germany, AQUASTATs reports only the water use of the manufacturing sector but omits the water abstractions of cooling water for thermal electricity production that is included in the WaterGAP results. The underestimation of industrial water uses $>$ 200 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ (Fig. S6) is particularly biased by the reported numbers from the US statistics. While AQUASTAT data include both freshwater and saline water abstractions from manufacturing, thermoelectric abstractions and mining, WaterGAP only accounts for the freshwater part of the manufacturing and thermoelectric abstractions.

WaterGAP performs reasonably well in the irrigation sector with a slightly better logarithmic NSE metric but with the overall lowest sectoral performance in terms of NSE (no visible direction in under- or overestimation).

Figure 5

Comparison of potential withdrawal water uses from WaterGAP 2.2d with AQUASTAT . Each data point represents one yearly value (if present in the database) per country for the time span 1962–2016.

[Figure omitted. See PDF]

6.4.2 Streamflow

The performance of WaterGAP 2.2d in terms of monthly streamflow time series at 1319 gauging stations (Fig. ) reaches a median NSE (KGE) of 0.52 (0.61). However, NSE values below 0 for 259 stations show that WaterGAP 2.2d cannot reproduce monthly and annual streamflow dynamics in one-fifth of the evaluated basins, although the simulated mean annual streamflow fits to the observations due to the calibration. The median for KGE${}_r$ of 0.79 indicates a relatively satisfactory simulation of the timing of monthly streamflows both seasonally and interannually. As the model is calibrated to match long-term annual river discharge (Sect. ), the median of the bias measure KGE${}_b$ is, with a value of 1.01, close to the optimum value. In rare cases, values outside the range of 0.9–1.1 occur as for calibration the individual basins were run for the calibration time period (plus 5 initialization years) while the evaluation run was a global run from 1901 to 2016. In the normal global runs, water demand can be fulfilled from neighboring grid cells while this is not possible in the calibration runs. This partially explains the larger biases also seen in Fig. . Streamflow variability is mostly underestimated by WaterGAP 2.2d, and median KGE${}_g$ is 0.85 (Fig. ).

When analyzing the spatial distribution of streamflow performance indicators, note that a highly seasonal streamflow regime tends to lead to high NSE and KGE${}_g$ not due to the quality of the evaluated hydrological model but due to the highly seasonal precipitation input. The global distribution of NSE classes shows a diverse pattern (Fig. ). Whereas large parts of central Europe, Asia and southern America are simulated reasonably well, the performance in northern America and large parts of Africa is in many cases below a value of 0.5. Based on NSE alone it remains unclear why WaterGAP consistently fails to satisfactorily simulate large parts of the well observed northern American region. Further insights can be gained by assessing the spatial distribution of KGE and its components (Fig. ). The broad picture of overall KGE (Fig. a) is similar to the NSE spatial distribution (Fig. ). In a large fraction of river basins with low NSE and KGE, the timing is off, with KGE${}_r<\mathrm{0}$.5. One reason could be the inappropriate modeling of the dynamics of lakes and wetland (mainly in Canada) and of reservoir regulations. As most snow-dominated basins in Alaska, Europe and Asia show a reasonably high KGE${}_r$ of $>\mathrm{0}$.8, it is not likely that snow dynamics are the dominant cause for low correlations between observed and simulated streamflow. For many other regions (e.g., central Asia and the Nile Basin), streamflow regulations due to reservoirs as well as the timing of water abstractions are most likely to cause low performance in timing. The indicator of variability KGE${}_g$ shows a medium to strong underestimation of streamflow variability in most of the northern snow-dominated basins. Underestimation in the Amazon basin is caused by the inability of WaterGAP to simulate wetland dynamics there. There are also many gauging stations for which WaterGAP overestimates seasonality, even by more than 50 %. Further research and development is needed for improving the GHMs in this respect .

Figure 6

Efficiency metrics for monthly streamflow of WaterGAP 2.2d at the 1319 GRDC stations with NSE, KGE and its components. Outliers (outside 1.5$\times$ inter-quartile range) are excluded but the number of stations that are defined as outliers are indicated after the metric.

[Figure omitted. See PDF]

Figure 7

Classified NSE efficiency metric for the 1319 river basins in WaterGAP 2.2d.

[Figure omitted. See PDF]

Figure 8

Classified KGE efficiency metric and its components for the 1319 river basins in WaterGAP 2.2d.

[Figure omitted. See PDF]

6.4.3 TWSA

WaterGAP 2.2d underestimates the mean annual TWSA amplitude in 54 % of the 143 investigated river basins by more than 10 % (Fig. ). Most of these basins are located in Africa, in the northern and monsoon regions of Asia, in Brazil, and in western North America. In contrast, the mean annual amplitude is overestimated in western Russia as well as in eastern and central North America. The correlation coefficient exceeds 0.7 in almost 75 % of the river basins and 0.9 in 22 %. Only 8 % of the basins show a correlation coefficient below 0.5.

The comparison of the TWSA trends shows that GRACE and WaterGAP 2.2d agree in the sign of the trend for 63 % of the 143 basins, for example most European basins; nearly the entire South American continent; and several basins in North America, Asia and Australia, but trends are often underestimated, e.g., in the Amazon and western Russia. Basins with different signs of the trend are scattered around the globe. GRACE suggests strong decreases in water storage in Alaskan basins, which is likely due to glacier mass loss, while WaterGAP determines a small mass increase, likely because WaterGAP does not simulate glaciers. Comparing the spatial pattern of Figs.  and , no obvious interrelation can be derived between the performances of streamflow and TWSAs.

Figure 9

Comparison of basin-average TWSAs of WaterGAP 2.2d and the average values of three GRACE mascon products for 143 basins larger than 200 000 ${\mathrm{km}}^{\mathrm{2}}$, with (a) ratio of amplitude (reddish colors indicate underestimated amplitude of WaterGAP, vice versa for bluish), (b) correlation coefficient, (c) trend of GRACE and (d) trend of WaterGAP 2.2d. All values based on the time series January 2003–December 2015.

[Figure omitted. See PDF]

6.5 Performance comparison between WaterGAP 2.2d and WaterGAP 2.2

Performance differences are expected due to modifications in model algorithms and the calibration routine (for details on modifications see Appendix ). When comparing the NSE of monthly streamflow (Figs.  and S7), the broad picture is similar. WaterGAP 2.2d shows some improvements in northern South America (especially the Amazon) but at the same time gets worse in southern South America. Slight decreases in performance for WaterGAP 2.2d are observed in southern Africa. No major changes are visible in North America, Europe and Asia, with small bidirectional changes. KGE patterns are also relatively similar for both versions (Figs.  and S8) and generally follow the differences in NSE. However, there are more regions in Europe and Asia where WaterGAP 2.2d performs better in overall KGE, resulting mainly from an improvement of KGE${}_r$. This is also visible in the number of basins per Köppen climate zone, where especially in the tropical A and dry B climates WaterGAP 2.2d has higher performance in KGE${}_r$ (Table ). The differences of KGE${}_b$ are negligible.

KGE${}_g$ shows significant differences between both model versions, in both directions, but performance of WaterGAP 2.2d is significantly better. Summarizing the basin statistics per Köppen climate zone, 272 instead of only 241 basins are within $\pm$10 % of observed variability in WaterGAP 2.2d in all climate zones except E (Table ). Fewer river basins (56 % compared to 61 % in 2.2) are subject to an underestimation of streamflow variability. However, the number of basins with overestimation increases slightly from 21 % for WaterGAP 2.2 to 23 % for WaterGAP 2.2d.

The performance of streamflow of the 1319 basins (Fig. S9) is similar for most indicators. The higher variation in KGE${}_b$ stems from modifications in the calibration routine, where up to $\pm$10 % uncertainty of observed streamflow is allowed. Similarly, the performance statistics of both streamflow and TWSAs (for the 143 basins $>$ 200 000 ${\mathrm{km}}^{\mathrm{2}}$) are very similar for both model versions (Fig. S10).

Table 4

Model performance with respect to streamflow timing: number of calibration basins per KGE${}_r$ category and Köppen–Geiger climate zone.

Model	Class	KGE${}_r$	A	B	C	D	E	Sum
2.2d	1	$>\mathrm{0}$.8	159	35	173	251	16	634
	2	0.5–0.8	109	47	77	200	17	450
	3	$<\mathrm{0}$.5	17	45	18	146	9	235
2.2	1	$>\mathrm{0}$.8	160	28	169	250	16	623
	2	0.5–0.8	104	46	80	202	18	450
	3	$<\mathrm{0}$.5	21	53	19	145	8	246

Table 5

Model performance with respect to streamflow variability: number of calibration basins per KGE${}_g$ category and Köppen–Geiger climate zone.

Model	Class	KGE${}_g$	A	B	C	D	E	Sum
2.2d	1	$>\mathrm{1}$.5	37	15	22	29	4	107
	2	1.1–1.5	46	22	71	58	5	202
	3	0.9–1.1	59	26	78	99	10	272
	4	0.5–0.9	124	51	88	281	10	554
	5	$<\mathrm{0}$.5	19	13	9	130	13	184
2.2	1	$>\mathrm{1}$.5	29	16	19	27	3	94
	2	1.1–1.5	46	18	57	54	6	181
	3	0.9–1.1	48	21	74	88	10	241
	4	0.5–0.9	141	49	109	277	10	586
	5	$<\mathrm{0}$.5	21	20	12	151	13	217

A comparison of simulated seasonality of streamflow and TWSAs in 12 selected large river basins across climate zones shows that performance with respect to both variables are improved in WaterGAP 2.2d for the Lena, Amazon and Yangtze basins (Fig. ). Simulations for the Congo, Mekong, Mackenzie and Murray basins do not differ. In some basins (Orange, Volga) the simulation of streamflow is improved in WaterGAP 2.2d whereas TWSA seasonality remains similar. In other basins (Rio Parana) seasonality agreement of TWSAs remains the same for WaterGAP 2.2d but streamflow seasonality agreement decreases.

Figure 10

Seasonality of streamflow and TWSAs of selected large river basins: model results of WaterGAP 2.2d and WaterGAP 2.2 as well as streamflow and TWSA observations.

[Figure omitted. See PDF]

7 Examples of model application

This section provides some examples of the WaterGAP 2.2d model applications for characterizing historical freshwater conditions at the global scale.

7.1 Model setup

The model setup is similar to those for the evaluation (Sect. ). For the purpose of model examples, the model was run in both the naturalized (nat) and the anthropogenic (ant) variant (Sect. ).

7.2 Spatial patterns of the global freshwater system

7.2.1 Renewable water resources

The quantification of (total) renewable water resources is one of the key elements of WaterGAP model application. They are defined as the long-term annual difference between precipitation and actual evapotranspiration of a spatial unit, or long-term annual net cell runoff. As runoff and evapotranspiration are influenced by human interference, renewable water resources are calculated based on the naturalized model variant, by averaging $R_{\mathrm{nc}}$ (Sect. ) over e.g., a 30-yr time period, resulting in $R_{\mathrm{nc},\mathrm{lta},\mathrm{nat}}$. On around 42.6 % of the global land area (excluding Greenland and Antarctica), total water resources are calculated to be $<$ 100 mm yr${}^{-\mathrm{1}}$ during the period 1981–2010, whereas on 19.8 % values are $>$ 500 mm yr${}^{-\mathrm{1}}$ (Fig. a). Globally averaged renewable water resources are computed to be 307 $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$ or 40 678 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$. The global map of inter-annual variability of runoff production (Fig. b), here defined as the ratio of runoff in a 1-in-10 dry year to total renewable water resources, shows regions with relatively constant and relatively variable annual runoff generation, in bluish and reddish colors, respectively. High variability is linked with low renewable water resources.

Total renewable water resources include renewable groundwater resources which are the sum of long-term average diffuse groundwater recharge $R_{\mathrm{g}}$ (Fig. c) and long-term average point (or focused) groundwater recharge from surface water bodies $R_{{\mathrm{g}}_{\mathrm{l},\mathrm{res},\mathrm{w}}}$ (Fig. d). While focused recharge is the major type of groundwater recharge in some (semi)arid grid cells, its quantification is highly uncertain, and diffuse groundwater recharge dominates in most cells. For 1981–2010, global mean diffuse groundwater recharge is calculated as 111.0 $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$, and global mean focused recharge as 12.8 $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$. Note that as $R_{\mathrm{g}}$ is calculated on (time-variable) land area (continental area minus fraction of lakes, reservoirs, wetlands) but is related to continental area in the standard output (Sect. ), grid cells with large gaining surface water bodies, e.g., wetlands along the Amazon river, show significant lower $R_{\mathrm{g}}$ values than surrounding grid cells.

The sum of diffuse and focused renewable groundwater resources amounts to 40 % of total renewable water resources, highlighting the important contribution of groundwater resources. There have been a number of studies on the potential impact of climate change on renewable groundwater resources (either including or excluding focused recharge), in which WaterGAP was applied as the impact model .

Figure 11

Water resources assessment 1981-2010 using WaterGAP 2.2d, with (a) total renewable water resources defined as long-term annual net cell runoff $R_{\mathrm{nc},\mathrm{lta},\mathrm{nat}}$ $\left[\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}\right]$, (b) 1-in-10 dry-year runoff generation in percent of total renewable water resources $\left[\mathrm{\%}\right]$, (c) long-term annual diffuse groundwater recharge $R_{\mathrm{g}}$ $\left[\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}\right]$, (d) long-term annual focused groundwater recharge $R_{{\mathrm{g}}_{\mathrm{l},\mathrm{res},\mathrm{w}}}$ $\left[\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}\right]$. Results are based on naturalized model runs. In (a) note that negative values for total water resources are possible (Sect.  ). In (b) areas where the denominator is $<{\mathrm{10}}^{-\mathrm{5}}$ are labeled as not defined.

[Figure omitted. See PDF]

Figure 12

Streamflow indicators of WaterGAP 2.2d for 1981–2010 with (a) long-term average annual streamflow $Q_{\mathrm{r},\mathrm{out},\mathrm{lta}}$ (${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$); (b) indication of streamflow alteration due to human water use and man-made reservoirs, where a reddish color indicates less streamflow for ant conditions, blue the opposite; (c) statistical monthly low flow $Q_{\mathrm{r},\mathrm{out},\mathrm{90}}$ in percent of $Q_{\mathrm{r},\mathrm{out},\mathrm{lta}}$; (d) differences of long-term average statistical monthly low flows as indication of low flow alteration due to human water use and man-made reservoirs. Not defined are areas where the denominator is smaller than ${\mathrm{10}}^{-\mathrm{5}}$ ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$.

[Figure omitted. See PDF]

7.2.2 Streamflow

Streamflow (or river discharge) $Q_{\mathrm{r},\mathrm{out}}$ is the model output that integrates all model components and human intervention, routing runoff along the river network. The global map of long-term average annual streamflow under anthropogenic conditions distinctly shows the very high spatial variability of streamflow and very distinctly the large river systems of the Earth (Fig. a). Temporal variability of monthly streamflow is much higher in the (semi)arid areas than in humid areas, increasing the spatial discrepancy of streamflow; this can be seen in Fig. c, which presents the ratio of the statistical low flow $Q_{\mathrm{90}}$ (the streamflow that is exceeded in 9 out of 10 months) to long-term average annual streamflow. The regions with a ratio of less than 5 % of low flow contribution on average streamflow (the hydrologically highly variable regions) follow in general the definition of (semi)arid grid cells (Fig. ) with some exceptions such as northern Asia. Different from the spatial pattern of interannual variability of long-term average net cell runoff (Fig. b), the spatial pattern of streamflow is characterized by low temporal variability in cells with large rivers, due to the integration of runoff from diverse grid cells as well as large water storage capacities in lakes, reservoirs or wetlands.

The impact of human interventions (human water use and man-made reservoirs) on streamflow is assessed in Fig. b for long-term averages and Fig. d for the statistical low flow indicator $Q_{\mathrm{90}}$ (please note the different legend for both subfigures). In general, human interventions reduce long-term average streamflow by at least 10 % (50 %) in 11.3 % (1.8 %) of the global land area, mainly due to reduced groundwater discharge to lakes, reservoirs, wetlands and rivers as a consequence of groundwater abstractions, in particular groundwater depletion (compare the red pattern with net abstraction from groundwater in Fig. a). There is only a minor share (0.7 %) of global land area, where long-term annual streamflow has been increased by more than 10 % due to human interventions (mainly return flow from groundwater abstractions). The impact of human interventions on $Q_{\mathrm{90}}$ is more pronounced (Fig. d). Large reddish patterns (consistent to net abstraction from groundwater in Fig. a) indicate the reduction of low flows by at least 10 % (90 %) on 29.7 % (14.4 %) of the global land area. However, there are also bluish river systems visible which represent a global land area of 5.3 % with increase in low flows of more than 10 %. Those areas are located downstream of large reservoirs that due to their storage capacity attenuate the flow regime towards a temporally less variable streamflow. As WaterGAP 2.2d considers only the largest reservoirs with reservoir management algorithm and handles the remaining $\sim \mathrm{6000}$ reservoirs of GRanD as unmanaged water bodies, the impact of streamflow regulation is most likely underestimated.

7.2.3 Water stress

A major motivation for the initial WaterGAP development was to consistently assess water stress on all land areas of the globe . A common water stress indicator (WSI) is calculated as the ratio of long-term average annual withdrawal water uses (or water abstractions of withdrawal water use) (Sect. ) and total renewable water resources for different spatial units (e.g., river basins). Renewable water resources in a basin are equal to long-term average naturalized annual streamflow at the outlet of the basin. WSI of 0.2–0.4 is generally assumed to indicate mild water stress and WSI $>$ 0.4 severe water stress e.g.,, while WSI $>$ 1.0 represents a situation where withdrawal water uses are larger than renewable water resources, indicating extreme water scarcity e.g.,. For this example, zero-order river basins (basins that drain to the oceans or inland sinks) were chosen as spatial units (Fig. ). River basins covering 73.6 % of global land area have a WSI $<$ 0.2 and thus are calculated to have none to only minor water stress. Mild (severe, but below extreme) water stress is represented in river basins that cover 9.7 % (6.9 %) of global land area. Extreme water stress (WSI $>$ 1.0) is simulated in river basins that cover 9.9 % of global land area (red colors in Fig. ). The spatial pattern of river basins with water stress is similar to the pattern of modification of statistically low flow alteration due to human interventions (Fig. d).

Figure 13

Water stress in zero-order river basins for 1981–2010, computed as the ratio of the basin sum of long-term average annual potential total withdrawal water uses (Sect. ) to long-term average annual streamflow $Q_{\mathrm{r},\mathrm{out},\mathrm{lta},\mathrm{nat}}$ of the basin (i.e., at its outflow cell to the ocean or at its inland sink).

[Figure omitted. See PDF]

Output of global models is usually shown in the form of two-dimensional planar global maps, which are necessarily distorted. While the Robinson projection that we normally use when presenting WaterGAP results is pleasing to the eye, it does not preserve the actual area of the land surface, and areas closer to Equator are shown relatively smaller than the areas closer to the poles. Using an equal-area projection as in Fig. b, Africa is shown larger than in the traditional Robinson maps. For Africa, large blue areas indicate high total renewable water resources per capita. However, very few people live in these large areas. For representing water resources for people instead of on areas, cartograms with population numbers as a distorter can be used (Fig. b). In cartograms, map polygons representing spatial units on the Earth's surface are distorted in a way that the units' polygon areas on the map are proportional to a quantitative attribute of the spatial unit , here the population in 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ grid cells in 2010. The latter was derived by aggregating 2010 GPWv3 gridded population estimate for the year 2010 from its original resolution of 2.5 arcmin. Clearly, with a higher share of red areas, the cartogram indicates a world with less water availability than the “normal” map, and it leads the eye to regions where humans are affected by water scarcity.

Figure 14

Water availability indicator per capita renewable water resources $Q_{\mathrm{r},\mathrm{out},\mathrm{lta},\mathrm{nat}}$ (${\mathrm{m}}^{\mathrm{3}}{\mathrm{cap}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) for 1981–2010 visualized in (a) an equal area projection and (b) as a cartogram with population in 2010 as distorter. In the cartogram each half-degree grid cell is distorted such that its area is proportional to the population of the grid cell.

[Figure omitted. See PDF]

Figure 15

Long-term (1981–2010) annual net abstractions: potential net water abstractions from surface water bodies (a), potential net water abstractions from groundwater (b), ratio of actual net water abstractions from surface water bodies to its potential value (c) and ratio of actual net water abstractions from ground water to its potential value. In (a) and (b) negative values indicate a net recharge of surface water and groundwater, respectively, due to return flows caused by human water use, while positive values indicate a net removal of water from the sources. In (c) and (d), cells with potential net water abstractions smaller than $\left|\mathrm{1}\right|$ $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$ are greyed out. Furthermore, grid cells where the sign of water abstractions changes between potential and actual net abstractions are displayed in red.

[Figure omitted. See PDF]

7.2.4 Water abstractions

With human water use being essential for the estimation of water stress, quantification of sectoral water uses was a focus already in the initial stages of WaterGAP development . However, a distinction of the sources of water abstractions and the sinks of return flows (groundwater or surface water) was only implemented later, such that potential net abstractions from groundwater and from surface water could be computed . Model refinements (see Appendix ) have lead to a more consistent computation of actual net abstractions from both sources. The general patterns of potential net abstractions (Fig. a and b) are consistent with the earlier assessment of . Positive values of NA${}_{\mathrm{s}}$ and NA${}_{\mathrm{g}}$ indicate that human water use results in a net subtraction of water from surface water bodies and groundwater, while negative values indicate a man-made addition of water to these water storage compartments. As noted in Sect. , the actual net abstractions can differ from their potential values. The ratio of actual to potential net surface water abstractions NA${}_{\mathrm{s}}$ (Fig. c) shows a heterogeneous pattern, with adjacent grid cells with values below 0.9 and above 1.1. This is explained by the option to satisfy water demand from a neighboring grid cell. In the case of negative NA${}_{\mathrm{s}}$, potential and actual values are always the same, as it is assumed in the model that NA${}_{\mathrm{g}}$ can always be fulfilled so that return flows to surface water are not changed. There are only a few longer river stretches where actual NA${}_{\mathrm{s}}$ is smaller than the potential value.

Actual NA${}_{\mathrm{g}}$ is equal to potential NA${}_{\mathrm{pot},\mathrm{g}}$ except in a few grid cells where potential NA${}_{\mathrm{pot},\mathrm{s}}$ cannot be fulfilled and there is irrigation with surface water (Fig. d). In these cells, return flows to groundwater decrease and actual values of NA${}_{\mathrm{g}}$ increase compared to their potential values. For example, in the case of a positive (negative) potential NA${}_{\mathrm{g}}$, a ratio of 1.1 (0.9) means that the difference between actual and potential NA${}_{\mathrm{g}}$ is 10 % of the absolute value of potential NA${}_{\mathrm{g}}$. In most grid cells, actual NA${}_{\mathrm{g}}$ is equal to the potential value.

Table 6

Global-scale (excluding Antarctica and Greenland) water balance components for different time spans as simulated with WaterGAP 2.2d. All units in ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$. Long-term average volume balance error is calculated as the difference of component 1 and the sum of components 2, 3 and 7.

No.	Component	1961–1990	1971–2000	1981–2010	1991–2016	2001–2016
1	Precipitation	111 388	111 582	111 616	112 052	112 559
2	Actual evapotranspiration${}^{\mathrm{1}}$	70 734	71 604	71 979	72 225	72 328
3	Streamflow into oceans and inland sinks	40 659	40 09	39 678	39 930	40 357
4	Actual consumptive water use${}^{\mathrm{2}}$	906	1023	1146	1238	1302
5	Actual net abstraction from surface water	1002	1108	1220	1304	1353
6	Actual net abstraction from groundwater	$-$96	$-$85	$-$74	$-$66	$-$50
7	Change of total water storage	$-$6	$-$31	$-$40	$-$104	$-$125
8	Long-term average volume balance error	0.34	0.23	0.11	0.03	0.01

${}^{\mathrm{1}}$ Including actual consumptive water use. ${}^{\mathrm{2}}$ Sum of rows 5 and 6.

Table 7

Globally aggregated (excluding Antarctica and Greenland) water storage component changes during different time periods as simulated by WaterGAP 2.2d. All units in ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$.

No.	Component	1961–1990	1971–2000	1981–2010	1991–2016	2001–2016
1	Canopy	0.0	0.0	0.1	0.0	0.0
2	Snow	16.6	$-$6.3	3.7	$-$12.6	5.0
3	Soil	$-$9.4	$-$2.2	16	14.5	17.3
4	Groundwater	$-$62.9	$-$62.7	$-$90.8	$-$108.8	$-$138.2
5	Local lakes	$-$1.1	$-$0.8	2.8	$-$0.3	$-$1.9
6	Local wetlands	$-$1.4	$-$3.0	3.5	0.0	4.0
7	Global lakes	$-$4.3	$-$5.2	$-$0.4	4.0	9.9
8	Global wetlands	$-$5.8	2.4	0.2	0.1	$-$10.3
9	Reservoirs and regulated lakes	68.2	43.6	28.1	5.7	$-$3.6
10	River	$-$5.6	3.3	$-$3.2	$-$6.4	$-$7.7
11	Total water storage	$-$5.8	$-$31.0	$-$40.0	$-$103.9	$-$125.3

Table 8

Globally aggregated (excluding Antarctica and Greenland) sectoral potential withdrawal water use WU and consumptive water use CU (${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$) as well as use fractions from groundwater ($\mathrm{\%}$) as simulated by GWSWUSE of WaterGAP 2.2d for the time period 1991–2016. These values represent demands for water that cannot be completely satisfied in WGHM due to lack of surface water resources (row 5 in Table ).

Water use sector	WU	Percent of WU	CU	Percent of CU
		from groundwater		from groundwater
Irrigation	2363	25	1100	37
Thermal power plants	599	0	16	0
Domestic	348	36	56	35
Manufacturing	272	27	53	26
Livestock	29	0	29	0
Total	3610	22	1253	36

7.3 Globally aggregated components of the land water balance components

7.3.1 Major water balance components

Estimation of globally aggregated components of the land water balance components is an intrinsic application field of GHMs. Independent of the time span assessed in Table , streamflow into oceans and inland sinks, equivalent to global renewable water resources, amounts to around 40 000 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ (with a range of around 1000 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$). Actual evapotranspiration is estimated to be around 71 000 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ (with a range of 1200 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$). Renewable water resources estimates are in the range of the estimates of previous WaterGAP model versions and of other global assessments (compare , their Table 3). Temporal trends of precipitation, actual evapotranspiration and streamflow may not be reliable due to uncertainty of the climate forcing and WaterGAP 2.2d. With less than 10${}^{-\mathrm{1}}$ ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$, the water balance error is negligible (Table ), which is an improvement compared to earlier model versions (see , their Table 2).

7.3.2 Water storage components

Total actual consumptive water use has increased over time and reaches the maximum in the most recent time period 2001–2016. The negative value of actual net abstraction from groundwater in Table  indicates that, globally aggregated, the groundwater compartment is recharged by return flows from irrigation with surface water (addition of the positive and negative values of NA${}_{\mathrm{g}}$ in Fig. b). A globally averaged anthropogenic increase in groundwater recharge is consistent with a decrease in groundwater storage that is mainly caused by the net groundwater abstractions. The global groundwater storage, however, has decreased (Table ), mainly due to groundwater depletion in those grid cells where (positive) NA${}_{\mathrm{g}}$ is higher than groundwater recharge . The anthropogenic net recharge of groundwater in the grid cells with negative NA${}_{\mathrm{g}}$ in Fig. b does not lead to a substantial increase in groundwater storage but mainly increases groundwater discharge to surface water bodies. The decreasing trend of total water storage is dominated by increasing water storage losses that were balanced in earlier periods by increased water storage in newly constructed reservoirs while dam construction became less during the last three decades (Table , ). However, WaterGAP 2.2d underestimates water storage increases because only the largest reservoirs are simulated as reservoirs including their commissioning year and because the GRanD v1.1 database used in WaterGAP 2.2d does not include some of the major reservoirs that were put into operation after 2000 . Soil water storage also contributes significantly to total water storage changes, showing increases since 1981. Different from what may be expected due to global warming, simulated global snow storage does not decrease over time (Table ).

7.3.3 Water use components

For the time period 1991–2016, Table  presents global sums of annual sectoral potential withdrawal water uses and consumptive water uses as well as the respective fractions that are taken from groundwater (Sect. ). Potential net abstractions from surface water (groundwater) are calculated by GWSWUSE to be 1406 ($-$153) ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ (Sect. ). Actual net abstractions from surface water (groundwater) are computed by WGHM to be 1304 ($-$66) ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ due to restricted surface water availability and consequently less return flows to groundwater from irrigation with surface water. It is thus estimated that 98.8 $\mathrm{\%}$ of potential consumptive water use of 1253 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ could be fulfilled during 1991–2016, albeit causing groundwater depletion.

8 Conclusions and outlook

A globally consistent quantification of water flows and storages as well as of human water use is needed but challenging, not only due to a lack of observation data but also the difficulty of appropriate process representation in necessarily coarse grid cells . This study fully describes the state-of-the-art GHM WaterGAP in its newest version 2.2d. Evaluation of model performance using independent data or observations of the key output variables, namely withdrawal water uses, streamflow and total water storage, indicates a reasonable model performance and points to potential areas of model improvement. Model output has been widely used for studying diverse research problems but also for informing the public about the state of the global freshwater system (see Supplement). The description of model algorithms, model outputs and related caveats will allow for better usage of model outputs by other researchers, who can now access these data from the PANGAEA repository.

Ongoing WaterGAP development aims to fully integrate a gradient-based groundwater model , improve the floodplain dynamics of large river basins (e.g., the Amazon) as proposed by and integrate glacier mass data . In addition, an update of the data basis for water use computations is planned. To enhance cross-sectoral integration in the framework of ISIMIP, modeling of river water temperature according to and will be implemented.

Appendix A Description of changes between the model versions 2.2 and 2.2d

A1 Modifications of water use models compared to WaterGAP 2.2

The modifications of water use models compared to WaterGAP 2.2 were as follows:

• Deficit irrigation with 70 % of optimal (standard) consumptive irrigation water use was applied in grid cells, which were selected based on and have (1) groundwater depletion of $>\mathrm{5}$ $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$ over 1989–2009 and (2) a $>\mathrm{5}$ % fraction of mean annual irrigation withdrawal water uses in total withdrawal water uses over 1989–2009 (Sect. ). In WaterGAP 2.2, optimal irrigation allowing the plants to evapotranspirate at 100 % of PET was assumed to be done everywhere.

• The time series of the Historical Irrigation Dataset (HID) for 1900 to 2005 was integrated into the Global Irrigation Model (GIM) (Sect. ) . In WaterGAP 2.2, irrigated areas of the static Global Map of Irrigation Area (GMIA) were scaled by time series of irrigated area per country. In addition to that, the newly available country-specific area actually irrigated (AAI), which is available for 47 countries, was used to update computed ICU until 2010. Version 2.2d enables the cell-specific AAI $/$ AEI ratio to be considered (for details see ).

• Non-irrigation water uses (domestic, manufacturing) were corrected to plausible values for coastal cells with small continental areas to avoid unrealistically high total water storage values in those cells.

A2 Modifications of WGHM compared to WaterGAP 2.2

A2.1 General

The following general modifications were made:

• With the introduction of dynamic extents of surface water bodies, land area fractions became variable in time as well (Sect. ).

• A modified routing approach where water is routed through the storages depends upon the fraction of surface water bodies; otherwise water is routed directly into the river (Sect. ) .

• Since WaterGAP 2.2b, net cell runoff $R_{\mathrm{nc}}$ is the difference between the outflow of a cell and inflow from upstream cells at the end of a time step (Sect. ). In the versions before, cell runoff was defined as outflow minus inflow into the river storage.

• In a modified calibration routine, an uncertainty of 10 % of long-term average river discharge is allowed (following ), meaning that calibration runs in four steps as described in Sect. .

• Since WaterGAP 2.2b, all model parameters which are potentially used for the calibration/data assimilation integration (including also parameter multiplicators) are read from a text file in JavaScript object notation (JSON) format.

• The differentiation into semiarid/humid grid cells are defined with a new standard methodology (Appendix ).

• For WaterGAP 2.2d, the return flows from surface water resources are scaled according to actual NA${}_{\mathrm{s}}$ (see results in Sect.  and Fig. ). Return flows induced by irrigation from surface water resources were calculated in WaterGAP 2.2 under the assumption that NA${}_{\mathrm{s}}$ can be fully satisfied. However, this can lead to implausible negative total actual consumptive water use, if surface water availability leads to smaller actual NA${}_{\mathrm{s}}$ than the return flows.

• A new storage-based river velocity algorithm was implemented (Sect. ).

• The realization of naturalized runs was improved. In WaterGAP 2.2, reservoirs were treated like global lakes in naturalized runs, while now, global reservoirs are completely removed (but local reservoirs are still handled as local lakes) (Sect. ). Please note that the studies of and were performed with an even older model version, in which all reservoirs were removed in naturalized runs.

A2.2 Soil

The following modification was made with respect to soil:

• The total water capacity input was newly derived and is now based on (Sect. ) whereas in WaterGAP 2.2 it was based on .

A2.3 Groundwater

The following modifications were made with respect to groundwater:

• Groundwater recharge below surface water bodies (LResWs) is implemented in semiarid and arid regions of in WaterGAP 2.2d.

• Regional changes since WaterGAP 2.2b are based on : (1) for Mississippi Embayment Regional Aquifer, groundwater recharge was overestimated, and thus the fraction of runoff from land recharging groundwater was reduced from 80 %–90 % to 10 % in these cells by adapting the groundwater factor $f_{\mathrm{g}}$ (Fig. S11); (2) groundwater depletion in the North China Plain was overestimated by a factor of 4, and thus runoff coefficient $\mathit{\gamma }$ was reduced from 3–5 to 0.1 in this area (Fig. S12); (3) all wetlands in Bangladesh were removed since diffuse groundwater recharge was unrealistically low.

• In WaterGAP 2.2d and for semiarid/arid grid cells, in the case of less precipitation than 12.5 $\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$, groundwater recharge remains in the soil column and is not handled as runoff anymore as in the versions before (Sect. ).

A2.4 LResWs

The following modifications were made with respect to LResWs:

• Precipitation on surface water bodies is now also multiplied with the evaporation reduction factor (like evaporation) to keep the water balance consistent (Sect. ).

• Reservoir information was updated, including the year when reservoir began operation (commissioning year; Sect. ) .

• Reservoir commissioning years were implemented in the reservoir algorithm (Sect. ) ; before this year, the reservoir is not present, and in the case of a regulated lake it is simulated as a global lake. In the versions before 2.2d, reservoirs and regulated lakes are simulated to be always present.

• For global lakes and reservoirs (where the water balance is calculated in the outflow cell), water demand of all riparian cells is included in the water balance of the outflow cell and thus can be satisfied by global lake or reservoir storage (Sect. ).

• All water storage equations in horizontal water balance are solved analytically in WaterGAP 2.2d (except for local lakes). Those equations now include net abstractions from surface water or groundwater. As a consequence, the sequence of net abstractions has been changed to (1) global lakes, regulated lakes or reservoirs, (2) rivers, and (3) local lakes (Sect. ).

• The areal correction factor (CFA) is included in the water balance of lakes and wetlands in WaterGAP 2.2d (Sect. ).

• In WaterGAP 2.2d (as in versions before WaterGAP 2.2), local and global lake storage can drop to $-S_{\max }$ as described in . The area reduction factor (corresponding to the evaporation reduction factor in (their Eq. 1) has been changed accordingly (denominator: $\mathrm{2}\times S_{\max }$). If lake storage $S$ equals $S_{\max }$, the reduction factor is 1; if $S$ equals $-S_{\max }$, the reduction factor is 0 (Sect. ).

• Active reservoir storage is no longer assumed to be 85 % but rather 100 % of reported storage (based on comparisons with the literature) (Sect. ).

Appendix B Definition of arid and humid grid cells

The definition of semiarid and arid grid cells is the basis for fractional routing (Sect. ), a groundwater recharge scheme (Sect. , ) and a PET equation (Sect. ), for example. In the model versions before WaterGAP 2.2c as used in , we defined the input file for semiarid/arid or humid grid cells according to the climate forcing used. However, it turned out that this leads to problems when comparing model outputs from different model versions and climate forcings. For example, if well-known non-humid regions (e.g., the High Plains Aquifer and the North China Plain) are classified as humid to a large extent due to uncertain climate forcing (and the approach used), this is not representing reality and can lead to implausible calculation of hydrological processes in those regions. Therefore, a static definition of semiarid/arid and humid grid cells was developed (Fig. ).

Following , the Priestley–Taylor $\mathit{\alpha }$ is set to a value of 1.26 for humid regions and of 1.74 for semiarid/arid regions. WaterGAP 2.2c was run with EWEMBI for 1981–2010 with all grid cells defined as humid to avoid predefinition of areas with high or low PET due to the initial setup of the $\mathit{\alpha }$. Following , drylands were defined based on an aridity index $(\mathrm{AI}=P/\mathrm{PET})$ with $\mathrm{AI}<\mathrm{0.65}$ and non-drylands with $\mathrm{AI}\ge \mathrm{0.65}$. Due to the definition of $\mathit{\alpha }$ as a humid value globally, PET might be too low, especially for transitional zones between drylands and non-drylands. Therefore, and based on visual inspection, we defined all grid cells with $\mathrm{AI}<\mathrm{0.75}$ as semiarid/arid grid cells. Furthermore, we defined all grid cells north of 55${}^{\circ }$ N as humid grid cells.

Figure B1

Static definition of humid and semiarid/arid grid cells.

[Figure omitted. See PDF]

Appendix C Land cover input

WGHM is using a static land cover input map (Fig. ) which is derived from Moderate Resolution Imaging Spectroradiometer () data for the year 2004 . The primary land cover attribute at the original resolution of 500 m is used as a basis. In the case of 500 m MODIS primary land cover being defined as “urban area”, “permanent wetland” or “water body”, the secondary land cover was used instead as those land cover types are included as a separate input (for lakes/wetlands the GLWD dataset, Sect. ; urban areas are implemented as impervious areas, Sect. ). Finally, the dominant IGBP (International Geosphere-Biosphere Programme) land cover type (primary land cover) was selected for each 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ grid cell.

Figure C1

Land cover classification of WaterGAP 2.2d.

[Figure omitted. See PDF]

Table C1

Parameters of the leaf area index model from .

No.	Land cover		Fraction of	$L$ reduction factor for	Initial days to start/end
	type		deciduous plants $f_{\mathrm{d}}$	evergreen plants $C_{\mathrm{e}}$	with growing season (d)
1	Evergreen needleleaf forest	${\mathrm{4.02}}^{\mathrm{a}}$	0	1	1
2	Evergreen broadleaf forest	${\mathrm{4.78}}^{\mathrm{b}}$	0	0.8	1
3	Deciduous needleleaf forest	4.63	1	0.8	10
4	Deciduous broadleaf forest	${\mathrm{4.49}}^{\mathrm{c}}$	1	0.8	10
5	Mixed forest	${\mathrm{4.34}}^{\mathrm{d}}$	0.25	0.8	10
6	Closed shrubland	2.08	0.5	0.8	10
7	Open shrubland	1.88	0.5	0.8	10
8	Woody savanna	2.08	0.5	0.3	10
9	Savanna	1.71	0.5	0.5	10
10	Grassland	1.71	0	0.5	10
11	Cropland	3.62	0	0.1	10
12	Cropland/natural vegetation mosaic	3.62	0.5	0.5	10
13	Snow and ice	0	0	0	0
14	Bare ground	1.31	0	1	10

${}^{\mathrm{a}}$ $L_{\max }$ is assumed to be the mean value of TeENL and BoENL land cover classes of . ${}^{\mathrm{b}}$ Only value for TrEBL and not TeEBL from as in WaterGAP this class is mainly in the tropics. ${}^{\mathrm{c}}$ Mean value from TeDBL and TrDBL from . ${}^{\mathrm{d}}$ Mean value of all forest classes. Fraction of deciduous plants and $L$ reduction factor for evergreen plants based on IMAGE () initial days to start/end with growing season are estimated.

Table C2

Attributes for IGBP land cover classes used in WaterGAP 2.2d from . Water has an albedo of 0.08, snow 0.6.

No.	Land cover	Rooting depth${}^{\mathrm{a}}$	Albedo${}^{\mathrm{a}}$	Snow albedo	Emissivity${}^{\mathrm{b}}$	Degree-day factor
	type	($\mathrm{m}$)	(–)	(–)	(–)	$D_{\mathrm{F}}$${}^{\mathrm{c}}$ ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}{}^{\circ }{\mathrm{C}}^{-\mathrm{1}}$)
1	Evergreen needleleaf forest	2	0.11	0.278	0.9956	1.5
2	Evergreen broadleaf forest	4	0.07	0.3	0.9956	3
3	Deciduous needleleaf forest	2	0.13	0.406	0.99	1.5
4	Deciduous broadleaf forest	2	0.13	0.558	0.99	3
5	Mixed forest	2	0.12	0.406	0.9928	2
6	Closed shrubland	1	0.13	0.7	0.9837	3
7	Open shrubland	0.5	0.2	0.7	0.9541	4
8	Woody savanna	1.5	0.2	0.558	0.9932	4
9	Savanna	1.5	0.3	0.7	0.9932	4
10	Grassland	1	0.25	0.7	0.9932	5
11	Cropland	1	0.23	0.376	0.9813	4
12	Cropland/natural vegetation mosaic	1	0.18	0.3	0.983	4
13	Snow and ice	1	0.6	0.7	0.9999	6
14	Bare ground	0.1	0.35	0.7	0.9412	6

${}^{\mathrm{a}}$ Adapted from the IMAGE model . ${}^{\mathrm{b}}$ . ${}^{\mathrm{c}}$ , .

Appendix D Integration of GLWD and GRanD data of lakes, reservoirs and wetlands (LResWs) into WGHM

WGHM uses the Global Lakes and Wetland Database (GLWD) and a preliminary but updated version of the Global Reservoir and Dam (GRanD) database to define location, area and other attributes of LResWs. The GLWD database consists of three datasets. GLWD-1 contains shoreline polygons of 3067 large lakes (area is $\mathit{>=}\mathrm{50}$ ${\mathrm{km}}^{\mathrm{2}}$) and 645 large reservoirs (capacity $\mathit{>=}\mathrm{0.5}$ ${\mathrm{km}}^{\mathrm{3}}$); GLWD-2 contains shoreline polygons of approximately 2 500 000 smaller lakes, reservoirs and rivers; and GLWD-3 is a 30 arcsec raster dataset with lakes, reservoirs, rivers and wetland types, including both GLWD-1 and GLWD-2 water bodies. The GRanD v1.1 database includes 6824 reservoir polygons . Information from these databases was translated to the six categories of LResWs implemented in WaterGAP and assigned to the 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ grid cells (see Table ). Figure  shows the spatial distribution of the maximum extent of all LResWs (all six categories) in terms of fractional coverage.

Table D1

LResW representation in WGHM. The total continental area represented in WaterGAP is 136.782 million ${\mathrm{km}}^{\mathrm{2}}$ (Antarctica is not included in WaterGAP) and 134.396 million ${\mathrm{km}}^{\mathrm{2}}$ without Greenland. The minimum land area (without Greenland), i.e., continental area minus maximum LResW area, is 124.449 million ${\mathrm{km}}^{\mathrm{2}}$.

No.	Surface water	Data source	Area description	Maximum global area	Definition
	body type			[million ${\mathrm{km}}^{\mathrm{2}}$]
1	Local wetland	GLWD-3	% of cell area	3.743	Wetland types 10, 11, 12, part of wetland types 4, 5, 7 and 8 of GLWD-3 (see description in Appendix )${}^{*}$
2	Global wetland	GLWD-3	% of cell area	3.752	Part of wetland types 4, 5, 7 and 8${}^{*}$
3	Local lake	GLWD-1, GLWD-2	% of cell area	0.850	Lakes with area $<$ 100 ${\mathrm{km}}^{\mathrm{2}}$ and reservoirs where a maximum storage capacity $<$ 0.5 ${\mathrm{km}}^{\mathrm{3}}$
4	Global lake	GLWD 1	% of cell area, total area of water body	1.010	Lakes with area $\mathit{>=}\mathrm{100}$ ${\mathrm{km}}^{\mathrm{2}}$
5	Global reservoir	GRanD	% of cell area, total area of water body	0.404	Man-made reservoirs with a maximum storage capacity $\mathit{>=}\mathrm{0}$.5 ${\mathrm{km}}^{\mathrm{3}}$
6	Global regulated lake	GRanD	% of cell area, total area of water body	0.188	Global lakes that are regulated and simulated like global reservoirs. Maximum storage capacity provided by GRanD is only the additional storage due to dam construction

${}^{*}$ Wetland categories of GLWD-3: 4 – freshwater marsh, floodplain, 5 – swamp forest, flooded forest, 7 – pan, brackish/saline wetland, 8 – bog, fen, mire, 10 – 50 %–100 % wetland (using 75 % of area as local wetland), 11 – 25 %–50 % wetland (using 35 % of area as local wetland), 12 – wetland complex (0 %–25 % wetland) (using 15 % of area as local wetland).

Figure D1

Fraction of local lakes (a), local wetlands (b), global lakes (c), global wetlands (d), global reservoirs (e), regulated lakes (f), grid cell area covered by LResWs (represents the maximum extent of LResWs) and land fraction (represents minimum extent of LResWs).

[Figure omitted. See PDF]

• Implementation of wetlands. GLWD-3 provides approximately the temporal maximum of wetland extent as wetland outlines were mainly derived from maps and are used to determine $A_{\max }$. In the case of various input datasets, a wetland was assumed to be present if at least one of the datasets showed one. The wetland types “coastal wetland” (covering 660 000 ${\mathrm{km}}^{\mathrm{2}}$) and “intermittent wetland/lake” (690 000 ${\mathrm{km}}^{\mathrm{2}}$) which are in GLWD-3 are not included in WGHM. Inclusion of coastal wetlands would require the simulation of ocean–land interaction, while intermittent wetlands/lakes of GLWD-3 cover very large parts of the deserts (compare Fig. 5 in ) that cannot be assumed to be covered totally by water at any time but rather represent areas where very rarely and at different points in time some parts may be flooded. Rivers shown in GLWD-3 are considered to be (lotic) wetlands and included as wetlands in WGHM. It is assumed that only a river with adjacent wetlands (floodplain) is wide enough to appear as a polygon on the coarse-scale source maps . For the fractional wetland type “50 %–100 % wetland”, an arbitrary value of 75 % grid cell coverage with wetland is assumed, for “25 %–50 % wetland” a value of 35 % and for “wetland complex” a value of 15 %. The large floodplain wetland of the lower Ganges–Brahmaputra in GLWD-3, covering almost all of Bangladesh, is not simulated as a wetland in WGHM, as during most of the time, only a small part of Bangladesh is inundated.

  All wetlands subsumed in fractional classes are assumed to be local, i.e., locally fed. In the case of all other wetland types, global wetlands fed by the whole catchment were identified as follows. All wetland polygons with a direct connection to a major river (as defined by the big_river.shp file available from ESRI) are assumed to receive inflow from a large upstream area and are therefore categorized as global. However, if rivers in this file are categorized as intermittent, the adjacent wetlands are categorized as local in WGHM. All other wetlands are first buffered (to the inside, using a GIS) by a 10 $\mathrm{km}$ wide ring such that the outer 10 $\mathrm{km}$ of a wetland is considered to be local and the core wetland area inside this buffer ring is considered to be global.

• Implementation of lakes, man-made reservoirs and regulated lakes. The 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ outflow cell of each global lake is determined based on the GLWD lake polygon and the DDM30 drainage direction map. If more than one global lake has the same outflow cell, the lakes are treated as one lake by adding the lake areas. The same procedure is done in the case of reservoirs/regulated lakes. There are 43 grid cells with 2 reservoirs, 6 grid cells with 3 reservoirs, 2 grid cells with 1 regulated lake and 1 reservoir, 1 grid cell with 2 regulated lakes, and 1 grid cell with 1 global lake and 1 regulated lake. Each cell can be the outflow cell of both a global lake and a global reservoir/regulated lake but if there is a regulated lake and a reservoir in one outflow cell, then they are aggregated. The commissioning year and main purpose of the larger reservoir/regulated lake is used. The commissioning year of the resulting 1109 reservoirs/regulated lakes that are simulated as individual reservoirs/regulated lakes was obtained mainly from the GRanD database but also other sources. In the commissioning year, the reservoir area is increased to its full extent (thus land area fraction is adjusted), the reservoir starts filling and reservoir algorithm is enabled. The storage capacity of the reservoirs which are in operation in the model initialization year is set to the maximum value .

Code and data availability

WaterGAP 2.2d is on the way to open source but still in the process of clarifying licensing and copyright issues. Hence, source code cannot be made publicly available but has been available for referees and editors. The standard model output data is available at 10.1594/PANGAEA.918447 and described in Sect. . For latest papers published based on WaterGAP 2, we refer to http://www.watergap.de .

The supplement related to this article is available online at: https://doi.org/10.5194/gmd-14-1037-2021-supplement.

Author contributions

HMS and PD led the development of WaterGAP 2.2d. HMS led the software development, supported by DC, CH, CN, TAP, EP, FTP, RR, SS, TT, and PD. The paper was conceptualized by HMS and PD. HMS did the calibrations, simulations, data analysis, visualization and model validation, supported by MS regarding validation against GRACE TWS. CN prepared model output for the PANGAEA data repository. The original draft was written by HMS, with specific parts drafted and reviewed by all authors. All authors contributed to the final draft. The revised version was written by HMS, with specific contribution from PD, MF, FTP and DC.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

We thank Tim Schön for generating Fig.  and processing data for Fig.  and Hans-Peter Ruhlhof-Döll for processing and generating Fig. . We furthermore thank Florian Herz for polishing the reference list and for technical support during manuscript preparation. We are grateful for Edwin Sutanudjaja for providing insights into the withdrawal water use comparison of PCR-GLOBWB. We acknowledge the evaluation datasets from GRDC (The Global Runoff Data Centre, 56068 Koblenz, Germany), AQUASTAT and GRACE (CSR RL05 GRACE mascon solutions were downloaded from http://www2.csr.utexas.edu/grace (last access: 5 March 2020), and JPL GRACE mascon data are available from http://grace.jpl.nasa.gov (last access: 5 March 2020), supported by the NASA MEaSUREs Program). We are grateful for valuable comments and suggestions from one anonymous referee and Gemma Coxon which helped to streamline and improve the consistency of the paper.

Financial support

The publication of this article was funded by the Open Access Fund of the Leibniz Association.

Review statement

This paper was edited by Jeffrey Neal and reviewed by Gemma Coxon and one anonymous referee.