Interactive comment on “ Global modeling of withdrawal , allocation and consumptive use of surface water and groundwater resources ”

Abstract. To sustain growing food demand and increasing standard of living, global water withdrawal and consumptive water use have been increasing rapidly. To analyze the human perturbation on water resources consistently over large scales, a number of macro-scale hydrological models (MHMs) have been developed in recent decades. However, few models consider the interaction between terrestrial water fluxes, and human activities and associated water use, and even fewer models distinguish water use from surface water and groundwater resources. Here, we couple a global water demand model with a global hydrological model and dynamically simulate daily water withdrawal and consumptive water use over the period 1979–2010, using two re-analysis products: ERA-Interim and MERRA. We explicitly take into account the mutual feedback between supply and demand, and implement a newly developed water allocation scheme to distinguish surface water and groundwater use. Moreover, we include a new irrigation scheme, which works dynamically with a daily surface and soil water balance, and incorporate the newly available extensive Global Reservoir and Dams data set (GRanD). Simulated surface water and groundwater withdrawals generally show good agreement with reported national and subnational statistics. The results show a consistent increase in both surface water and groundwater use worldwide, with a more rapid increase in groundwater use since the 1990s. Human impacts on terrestrial water storage (TWS) signals are evident, altering the seasonal and interannual variability. This alteration is particularly large over heavily regulated basins such as the Colorado and the Columbia, and over the major irrigated basins such as the Mississippi, the Indus, and the Ganges. Including human water use and associated reservoir operations generally improves the correlation of simulated TWS anomalies with those of the GRACE observations.


Introduction
In 1900, global population was less than 1.7 billion, but grew by more than 4 times during the 20th century, currently exceeding 7 billion.To sustain growing food demand and increasing standard of living, global water withdrawal increased by nearly 6 times from ∼ 500 km 3 yr −1 in 1900 to ∼ 3000 km 3 yr −1 in 2000, of which agriculture is the dominant water user (≈ 70 %) (Falkenmark et al., 1997;Shiklomanov, 2000a, b;Döll and Siebert, 2002;Vörösmarty et al., 2005;Haddeland et al., 2006;Bondeau et al., 2007;Wisser et al., 2010;Wada et al., 2013).Soaring water withdrawal worsens water scarcity conditions already prevalent in semiarid and arid regions (e.g., India, Pakistan, northeastern China, the Middle East and North Africa), where available surface water is limited due to lower precipitation, increasing uncertainty for sustainable food production and economic development (World Water Assessment Programme, 2003;Hanasaki et al., 2008b;Döll et al., 2009;Kummu et al., 2010;Vörösmarty et al., 2010;Wada et al., 2011b).In these regions, the water demand often exceeds the available surface water resources due to intense irrigation which requires large Published by Copernicus Publications on behalf of the European Geosciences Union.
To quantify the surface water balance, i.e., water in rivers, lakes, wetlands, and reservoirs, and to analyze the human perturbation on water resources consistently over a large scale, a number of macro-scale hydrological models (MHMs) have been developed in recent decades.Yates (1997) and Nijssen et al. (2001a, b) applied MHMs to calculate runoff and river discharge over river basin to continental scales at a relatively coarse spatial grid (1-2 • ).Arnell (1999Arnell ( , 2004) ) and Vörösmarty et al. (2000b) used respectively the Macro-PDM and WBM to simulate global surface water balance at a finer scale (0.5 • ).Oki et al. (2001) used the TRIP (0.5 • ) to route global local runoff simulated by land surface models (LSMs).These models, however, do not include the effect of water withdrawal on the surface water balance.Alcamo et al. (2003a, b) developed the Water-GAP model (0.5 • ), which simulates the global surface water balance and global water use, i.e., water withdrawal and consumptive water use, from agricultural, industrial, and domestic sectors.Döll et al. (2003Döll et al. ( , 2009) ) used the WGHM (0.5 • ) (Alcamo et al., 2007;Flörke et al., 2013;Portmann et al., 2013) to simulate globally the reduction of river discharge by human water consumption.Hanasaki et al. (2008aHanasaki et al. ( , b, 2010) ) and Pokhrel et al. (2012a, b) developed the H08 (0.5 • ) and MATSIRO (0.5 • ) respectively, both of which incorporate the anthropogenic effects (e.g., irrigation, reservoir regulation) into global surface water balance calculation.Wada et al. (2010Wada et al. ( , 2011a, b) , b) and Van Beek et al. (2011) developed the PCR-GLOBWB model (0.5 • ) to calculate the surface water balance and monthly sectoral water demand, and incorporated groundwater abstraction at the global scale.However, these models generally calculate water demand separately and independent of water availability, i.e., there is no feedback between human water use and terrestrial water fluxes, and equate water demand with either water withdrawals or consumptive water use (Döll and Siebert, 2002;Wisser et al., 2010;Wada et al., 2011b).In addition, water allocation or water use per source (surface water and groundwater) has rarely been dynamically incorporated in the models.
Here, we substantially improve the PCR-GLOBWB model (version 2.0) on the basis of the previous version of the model (version 1.0) presented in Wada et al. (2010Wada et al. ( , 2011a, b) , b) and Van Beek et al. (2011).We first couple the global water demand model developed by Wada et al. (2011a, b) with the global hydrological model PCR-GLOBWB (Wada et al., 2010;Van Beek et al., 2011).In the previous version of the model, water availability (water in rivers, lakes, and reservoirs) and water demand (agriculture, industry, and households) were calculated independently, and the simulation results were compared afterwards (as a post-process) to estimate, for example, water scarcity (Van Beek et al., 2011;Wada et al., 2011a, b).In the present version (version 2.0) of the model, water availability and water demand calculation is integrated to dynamically simulate water use at a daily time step and to account for the interactions between human water use and terrestrial water fluxes.The main goal of this integrated modeling framework is to estimate actual water use (i.e., withdrawal and consumption) rather than potential water demand (that is independent of available water).To enable this modeling framework, we implement a new irrigation scheme, in which irrigation water is supplied based on daily surface water and soil water balance and deficit.This will consider the mutual feedback from irrigation water supply to the soil and groundwater system, and the associated evapotranspiration over irrigated areas.Another improvement is that to satisfy the (potential) demands we consider water allocation and use from available surface water and groundwater resources at a daily time step.This allows us to distinguish the different response of human water use impacts on surface water (faster) and groundwater (slower) systems, and the mutual feedback between them due to, for instance, irrigation return flow.To improve the simulations of surface water availability, we also include the newly available extensive Global Reservoir and Dams data set (GRanD).Moreover, we update the climate forcing and use two newly available climate reanalysis data sets (ERA-Interim and MERRA) over the period 1979-2010, extending beyond most global analyses.
The overall objectives of this study are (1) to develop a coupled global hydrological and water demand model, (2) to evaluate the performance of the integrated modeling approach in terms of simulated water withdrawal and consumptive water use from surface water and groundwater resources, and (3) to quantify the impact of human perturbation (human water use and reservoir regulation) on terrestrial water resources consistently across large scales (e.g., basin).
Section 2 of this paper presents the integrated modeling framework which describes the coupling of the global hydrological model and the global water demand model at a daily temporal resolution.The section includes a brief introduction of the global hydrological model, but the other parts of the section are limited to our improved approaches modeling daily water demand or requirement for irrigation and other sectors, routing and surface water retention, and water allocation from surface water and groundwater resources and associated return flow.After an introduction of the simulation protocol in Sect.3, Sect. 4 presents the simulation results and evaluates their performance by comparing them to available statistics and satellite information.Section 5 discusses the advantages and the limitations of our modeling framework and the associated uncertainties, and provides conclusions from this study.

Water balance
The global hydrological model PCR-GLOBWB simulates for each grid cell (0.5 • × 0.5 • globally over the land) and for each time step (daily) the water storage in two vertically stacked soil layers and an underlying groundwater layer, as well as the water exchange between the layers (infiltration, percolation, and capillary rise) and between the top layer and the atmosphere (rainfall, evapotranspiration, and snowmelt).The model also calculates canopy interception and snow storage.Subgrid variability is taken into account by considering separately tall and short vegetation, open water (lakes, reservoirs, floodplains and wetlands), different soil types based on the FAO Digital Soil Map of the World (FAO, 2003), and the area fraction of saturated soil calculated by the improved Arno scheme (Todini, 1996;Hagemann and Gates, 2003) as well as the frequency distribution of groundwater depth based on the surface elevations of the HYDRO1k Elevation Derivative Database (HYDRO1k; US Geological Survey Center for Earth Resources Observation and Science; http://eros.usgs.gov/#/Find_Data/Products_and_Data_Available/HYDRO1K).The groundwater layer represents the deeper part of the soil that is exempt from any direct influence of vegetation and constitutes a groundwater reservoir fed by active recharge.The groundwater store is explicitly parameterized based on lithology and topography, and represented as a linear reservoir model (Kraaijenhoff van de Leur, 1958).Natural groundwater recharge fed by net precipitation and additional recharge from irrigation, i.e., return flow, fed by irrigation water (see Sect. 2.2) occurs as the net flux from the lowest soil layer to the groundwater layer, i.e., deep percolation minus capillary rise.Groundwater recharge interacts with groundwater storage as it can be balanced by capillary rise if the top of the groundwater level is within 5 m of the topographical surface (calculated as the height of the groundwater storage over the storage coefficient on top of the streambed elevation and the subgrid distribution of elevation).Groundwater storage is fed by groundwater recharge and drained by a reservoir coefficient that includes information on lithology and topography (e.g., hydraulic conductivity of the subsoil).The ensuing capillary rise is calculated as the upward moisture flux that can be sustained when an upward gradient exists and the moisture content of the soil is below field capacity.Also, it cannot exceed the available storage in the underlying groundwater reservoir.
The detailed description of the basic hydrologic model structure, and associated calculation and parameterization is given in Appendix A, and only newly developed parts of the model are described in the following sections.Figure 1 shows a schematic diagram of the integrated modeling framework that couples the hydrological model with human activities including water use and reservoir regulation.

Irrigation water requirement
A new irrigation scheme was implemented that separately parameterizes paddy and nonpaddy crops and that dynamically links with the daily surface and soil water balance considering the feedback between the application of irrigation water and the corresponding changes in surface and soil water balance.This in turn affects the amount of soil moisture and irrigation water requirement over the paddy and nonpaddy fields in following days.This enables to simulate more realistically the state of daily soil moisture condition, and associated evaporation and crop transpiration over irrigated areas.Previous studies used various methods simulating irrigation water requirement (IWR) as shown in Table 1.However, few models separately parameterize paddy and nonpaddy crops, and explicitly consider the feedback between irrigation water application, and associated change in surface and soil water balance.
The losses during water transport and irrigation application are included in the calculation of IWR (∼ gross irrigation water requirements).To account for such losses, other MHMs (e.g., H08, MATSIRO, WaterGAP, and WBM) use irrigation or project efficiency taken from available country statistics (Döll and Siebert, 2002;Rohwer et al., 2007;Rost et al., 2008), whereas we dynamically calculate the efficiency based on daily evaporative and percolation losses per unit crop area based on the surface and soil water balance (i.e., susceptible to the amount of soil moisture).
Crop-specific calendars and growing season lengths were obtained from the MIRCA2000 data set (Portmann et al., 2010), which accounts for various growing seasons of different crops and regional cropping practices under different climatic conditions, and distinguishes up to nine subcrops that represent multi-cropping systems in different seasons in different areas per grid cell.The corresponding crop coefficient per crop development stage and maximum crop rooting depth were additionally obtained from the Global Crop Water Model (Siebert and Döll, 2010).Although the MIRCA2000 data set considers 26 crop classes, we aggregated these to paddy and nonpaddy crop classes since distinct flooding irrigation is applied over most of paddy fields.The crop-specific parameters were aggregated by weighing the area of each crop class.
Daily (potential) crop evapotranspiration, ET c [m d −1 ], was calculated combining a crop coefficient, k c [dimensionless], that accounts for crop-specific transpiration and bare soil evaporation over the surface, with reference (potential) evapotranspiration, ET 0 [m d −1 ], computed with the Penman-Monteith equation according to FAO guidelines (Doorenbos and Pruitt, 1977;Allen et al., 1998): (1) Irrigation water [m d  Döll and Siebert (2002) CRU TS 1.0 (New et al., 2000) Priestley and Taylor   (Kalnay et al., 1996) FAO Penman-Monteith Siebert et al. (2005Siebert et al. ( , 2007) ) FAO Thenkabail et al. (2006) Siebert et al. (2007) 20 crops (You et al., 2006) FAO CROPWAT with some adjustments JRA-25 Reanalysis (Kim et al., 2009;Onogi et al., 2007) FAO Penman-Monteith Siebert et al. (2007) Freydank andSiebert (2008) 18 crops (Leff et al., 2004) SWIM model (Krysanova et al., 1998) Energy paddy fields, we maintained a 50 mm surface water depth, S max , (Wisser et al., 2008(Wisser et al., , 2010) ) until the late crop development stage (∼ 20 days) before the harvest.We opted for no irrigation approximately 20 days before the harvest based on irrigation practices that generally occur over paddy fields (Allen et al., 1998;Aslam, 1998).The duration of the nonirrigation period (∼ late crop development stage) varies depending on a region and local practices.For some regions (e.g., Africa), water is drained from the paddy field ∼ 10 days before the expected harvest date as draining hastens maturity and improves harvesting conditions.It is also common that irrigation is ceased a few weeks before harvest over the paddy fields to dry and for the rice to transfer maximum nutrients into the grains (e.g., Asia).Paddy irrigation water requirement and associated surface water balance are estimated as IWR paddy,t = max(0, S max − (S 0,t−1 + P net,t )), (2) where S 0,t is the surface water layer [m] over the paddy fields at a given time, t, and P net is the net liquid precipitation [m d −1 ], precipitation reduced by interception losses and snowfall.q i is the infiltration from the surface water layer, S 0 , to the first soil layer, S 1 , at a rate of saturated hydraulic conductivity of the first soil layer (k sat ) [m d −1 ].The saturated hydraulic conductivity was reduced by a factor ∼ 10 considering compacted soil preventing high percolation losses that is commonly practiced over paddy fields (Bhadoria, 1986).
EW is the open water evaporation from the surface water layer (S 0 ) [m d −1 ], assumed to occur at the potential rate over shallow water (Allen et al., 1998).t denotes time step [day].
We assumed that no direct runoff occurs over the paddy fields as farmers tend to irrigate much less before expected (heavy) rainy days or periods.However, this may underestimate direct runoff that occurs over flooded paddy fields during substantial rainfall particularly in humid regions (e.g., southern China, Indonesia, Bangladesh).
For the nonpaddy crop type, we estimated IWR nonpaddy by taking the difference between total (TAW) and readily available water (RAW) in the first and second soil layer with no surface water layer (Allen et al., 1998): where TAW is the total soil moisture available to irrigated crops in the soil column and RAW is for each time step the actual soil moisture available in the root zone (see Fig. 1).
effective degree of saturation at wilting point [all dimensionless].θ sat is the saturated (volumetric) water content, and θ res is the residual (volumetric) water content [all in m 3 m −3 ].SC is the storage capacity of the soil layer, and Z r is the rooting depth assuming an exponential growth to the maximum rooting depth over the growing season (Jackson et al., 1996) [all in m].S 1 and S 2 denote the first and second soil layer respectively.The parameter, p, is the soil water depletion fraction that is a function of daily crop evapotranspiration [m d −1 ], and p ref is the reference soil water depletion fraction per crop type (0.2 for paddy and 0.5 for nonpaddy).Although water in root zone is theoretically available until wilting point, crop water uptake is reduced well before wilting point is reached (Allen et al., 1998).When the soil is sufficiently wet, the soil supplies water fast enough to meet the atmospheric demand of the crop, and water uptake equals ET c (crop evapotranspiration; Eq. 5), however, as the soil water content decreases, water becomes more strongly bound to the soil matrix and is more difficult to extract (Allen et al., 1998).Thus, when the soil water content drops below a threshold value, soil water can no longer be transported quickly enough towards the roots to respond to the transpiration demand and the crop begins to experience stress.The soil water depletion fraction determines the fraction of TAW that a crop can extract from the root zone without suffering the water stress (∼ RAW; Eq. 4).
Historical growth of irrigated areas  was estimated using country-specific statistics of irrigated areas  for ∼ 230 countries (FAOSTAT; http://faostat.fao.org/) and by downscaling these to 0.5 • using the spatial distribution of the gridded irrigated areas from the MIRCA2000 data set (Portmann et al., 2010).This method is unable to reproduce changes in the distribution within countries, but it adequately reflects the large-scale dynamics of the expanding irrigated areas over the past decades (Wisser et al., 2010).

Other sectoral water demands
Other sectoral water demands include those from livestock, industry, and households and were estimated at a daily time step [all in m d −1 ] over the period 1979-2010, considering the past change in population, socioeconomic and technological development, and livestock densities.Livestock water demand was calculated by multiplying the number of livestock in a grid cell with its corresponding daily drinking water requirement, which is a function of daily air temperature (Wada et al., 2011b).The gridded global livestock densities of cattle, buffalo, sheep, goats, pigs and poultry in 2000, and their corresponding drinking water requirements were obtained from FAO (2007) and Steinfeld et al. (2006) respectively.For the other years , the numbers of each livestock type per country (FAOSTAT; http://faostat.fao.org/) were downscaled to a grid scale using the distribution of each gridded livestock density in 2000.
Gridded industrial water demand data for 2000 was obtained from Shiklomanov (1997), WRI (1998), andVörösmarty et al. (2005).Due to limited available data in order to identify the seasonal trends, daily industrial water demand was kept constant over the year similar to the study of Hanasaki et al. (2006Hanasaki et al. ( , 2008a, b) , b) and Wada et al. (2011b).However, in reality daily industrial water demand likely fluctuates over the year, although the seasonal amplitude may not be large.To calculate time series  of industrial water demand, we multiplied the gridded industrial water demand for 2000 with water use intensities calculated with an algorithm developed by Wada et al. (2011a).This algorithm calculates country-specific economic development based on four socioeconomic variables: gross domestic product (GDP), electricity production, energy consumption, and household consumption.Associated technological development per country was then approximated by energy consumption per unit electricity production, which accounts for industrial restructuring or improved water use efficiency.
Household water demand was estimated multiplying the number of persons in a grid cell with the country-specific per capita domestic water withdrawal.The daily course of household water demand was estimated using daily air temperature as a proxy (Wada et al., 2011a).The country per capita domestic water withdrawals in 2000 were taken from the FAO AQUASTAT database (http://www.fao.org/nr/water/aquastat/main/index.stm) and Gleick et al. (2009), which were multiplied with water use intensities to account for economic and technological development.Available gridded global population maps per decade (Klein Goldewijk and van Drecht, 2006) were used to downscale the yearly country population data (FAOSTAT) to produce gridded population maps for each year.

Routing and surface water retention
The simulated local direct runoff, interflow, and baseflow (see Appendix A) were routed along the river network based on the Simulated Topological Networks (STN30; Vörösmarty et al., 2000a).The routing is based on the characteristic distances, where volumes of water are transported over a distance, R cd , along the drainage network.R cd is given by where b and z are the channel width and channel depth respectively [m], G is the gradient derived from the elevation and the drainage network, and n is Manning's roughness coefficient.
Reservoirs are located on the drainage or river network based on the newly available and extensive GRanD data set (Lehner et al., 2011), which contains 6862 reservoirs with a total storage capacity of 6197 km 3 .The reservoirs were placed over the river network based on the years of their construction.If more than one reservoir fell into the same grid cell, we aggregated the storage capacities and modeled a single reservoir.In case no reported value was available, reservoir surface area [m 2 ], A, was calculated using the storage volume (V )-reservoir depth (h) relationship (Campos, 2010): where α is the reservoir specific shape factor [dimensionless], computed from the reported dam height and the reported storage capacity or S max .Similar to Hanasaki et al. (2006) and Van Beek et al. (2011), reservoir release was simulated to satisfy local and downstream water demands that could be reached within ∼ 600 km (∼ a week with an average discharge velocity of 1 m s −1 ) or a next downstream reservoir if present.In case of no water demand, the reservoir release, R r [m 3 day −1 ], was simulated as a function of minimum, S min (set to ∼ 10 % of storage capacity), maximum, S max (set to ∼ 100 % of storage capacity), and actual reservoir storage [all in m 3 ], S r , and mean average inflow, where I is the inflow to the reservoir, P local is the local precipitation over the reservoir surface, and EW r is the open water evaporation from the reservoir surface, assumed to occur at a rate of potential evapotranspiration [all in m 3 ].Reservoir spills occur when the reservoir storage exceeds the maximum reservoir storage.

Water allocation and return flow
Water demands for irrigation, livestock, industry, and households can be met from three water resources: (1) desalination, (2) groundwater, and/or (3) surface water.Around the globe, more than 10 000 desalination plants in 120 countries are in operation (World Water Assessment Programme, 2003).Although energy and economic costs to process sea water to produce purified water is still much higher than conventional water supply measures such as groundwater pumping, the amount of desalinated water use has been rising since the 1990s.Desalinated water use is generally limited to coastal areas and provides a stable amount of water supply over arid regions such as the Middle East and North Africa, where over 70 % of the global desalination capacity is installed and people receive ∼ 1 % of the global runoff.We used available country statistics of desalination water withdrawal for the period 1960-2010 from two data sources.
The country statistics were primarily obtained from the FAO AQUASTAT database, but were supplemented by the WRI EarthTrends (http://www.wri.org/project/earthtrends/;World Resources Institute, 1998) where applicable (global total ≈ 15 km 3 yr −1 ).The data are given in 5 yr intervals and we linearly interpolated these to estimate annual values.We then spatially downscaled the country values onto a global coastal ribbon of ∼ 40 km based on the gridded population intensities considering the fact that desalinated water is used in coastal areas (Wada et al., 2011b), and assumed constant withdrawals of desalination over the year.Allocation of surface water and groundwater to satisfy the remaining water demand (after subtracting desalinated water withdrawal) depends on available surface water including local and upstream reservoirs and readily extractable groundwater reserves.Since the absolute amount of available groundwater resources is not known at the global scale, we used the simulated daily (accumulated) baseflow, Q base [m 3 day −1 ], against the long-term average river discharge, Q avg [m 3 day −1 ], as a proxy to infer the readily available amount of renewable groundwater reserves, WA gw [m 3 day WA gw was then extracted daily from renewable groundwater storage, S 3 [m 3 day −1 ], to meet part of the water demand (Eq.13), WD tot [m 3 day −1 ].To avoid no local groundwater withdrawal over arid regions with negligible local baseflow, we used accumulated baseflow over a catchment, which allows regions with no local baseflow to extract local groundwater resources.The remaining water demand was then withdrawn from the simulated surface water.However, in case reservoirs are present at local or upstream grid cells over the river network, we first allocated surface water rather than groundwater (WA gw ) to meet the water demand, and the remaining water demand was met from available groundwater storage or S 3 .In case of no outstanding water demand, no groundwater is abstracted.
In case of lack of accumulated baseflow due to extremely dry conditions, surface water availability is also expected to be very small.The unmet water demand is then imposed on (nonrenewable) groundwater (e.g., groundwater withdrawal in excess of available groundwater storage, S 3 ).The available water is allocated proportionally to the amount of sectoral water demands.No priority is given to a specific sector, but a competition of water use among the sectors likely occurs over many water scarce regions, particularly for surface water resources.
Return flow from water that is withdrawn for the industrial and domestic sectors is assumed to occur to the river system on the same day (no retention due to waste water treatment).For the domestic sector, the return flow occurs only from the areas where urban and rural population have access to water (UNEP; http://www.unep.org/),whereas for the industry sector, the return flow occurs from all areas where water is withdrawn.For both sectors, the amount of return flow is determined by recycling ratios developed per country.
The country-specific water recycling was calculated according to the method developed by Wada et al. (2011a, b) who interpolated recycling ratios on the basis of GDP and the level of economic development, i.e., high income (80 %; 20 % of water is actually consumed.),middle income (65 %; 35 % of water is consumed.),and low income economies (40 %; 60 % of water is consumed.).The ratio was kept at 80 % if a country reached the high income economy, and the ratio of 40 % was assigned to countries with no GDP data.For the irrigation sector, return flow occurs to the soil layers as infiltration and to the groundwater layer as additional recharge (see Sect. 2.2).No return flow to the soil or river system occurs from the livestock sector.For completeness, we note that consumptive water use is equal to water withdrawal minus return flow.

Model simulation
To simulate global water use, i.e., water withdrawal and consumptive water use, we obtained daily climate drivers (e.g., precipitation and mean air temperature) over the period 1979-2010.We retrieved the data from the ERA-Interim reanalysis, where the precipitation was corrected with GPCP precipitation (GPCP: Global Precipitation Climatology Project; http://www.gewex.org/gpcp.html)(Dee et al., 2011).To account for climate uncertainty, we also retrieved the data from the MERRA reanalysis product (available at http://gmao.gsfc.nasa.gov/merra).Over the same period, we calculated reference evapotranspiration based on the Penman-Monteith equation according to FAO guidelines for a hypothetical grass surface with a specified height of 0.12 m, an albedo of 0.23, and a surface resistance of 70 s m −1 (Allen et al., 1998) with relevant climate fields (e.g., cloud cover, vapor pressure, wind speed) retrieved from the ERA-Interim and MERRA data sets.
For compatibility with our overall analysis, we biascorrected these data sets (precipitation, reference evapotranspiration, and temperature) on a grid-by-grid basis (0.5 degree grid) by scaling the long-term monthly means of these fields to those of the CRU TS 2.1 data set (Mitchell and Jones, 2005) over the overlapping period   (Wada et al., 2012b).For temperature, we calculated per month the long-term mean temperature  for each climate forcing (ERA-Interim and MERRA) and for the CRU data, and attributed the difference (additive) to the mean daily temperature from each climate forcing.For reference evapotranspiration, we corrected per month the amount for each climate forcing by attributing the ratio (multiplicative) of the long-term mean of the CRU data over that of each climate forcing.For precipitation, we first corrected the number of wet days for each climate forcing (ERA-Interim and  MERRA).We estimated per month the mean threshold precipitation by equalizing the number of wet days for each climate forcing to that for the CRU data over the 1979-2001 period (Wada et al., 2012b).The daily precipitation below these thresholds was removed.We then corrected per month the amount of precipitation for each climate forcing by attributing the ratio (multiplicative) of the long-term mean precipitation of the CRU data over that of each climate forcing.The resulting monthly additive (temperature), multiplicative bias-correction factors (reference evapotranspiration and precipitation), and the wet days correction (precipitation) were subsequently applied to the daily climate fields for the entire simulation period .We applied this method over regions wherever at least two CRU stations are present.Otherwise the original ERA-Interim and MERRA climate data were returned by default.

Results
To evaluate our modeling approach, we first compared our simulated water use to available reported national and subnational statistics.Since simulated river discharge, total water withdrawal and total consumptive water use have been extensively validated in earlier work (Van Beek et al., 2011;Wada et al., 2011aWada et al., , 2012a)), we, here, focus on validating simulated water withdrawal per source (surface water and groundwater), to assess our water allocation scheme.Reported statistics on consumptive water use per water source rarely exists even at a national or subnational level.After the validation, we provide a regional overview of water withdrawal and consumptive water use trends over the period 1979-2010.A limited validation exercise is also provided to assess the impact of human-induced change on simulated river discharge per river basin.We then compare our simulated terrestrial water storage (TWS) anomalies with those of the GRACE observations over the period 2003-2010 to assess the impacts of human water use and associated reservoir operations on TWS over the selected catchments.

Accuracy of simulated irrigation water requirement
Figure 2 compares our simulated IWR with reported country statistics obtained from the FAO AQUASTAT database.IWR was simulated with the CRU TS2.1, ERA-Interim and MERRA climate respectively.Table 2 shows the correlation between the simulated IWR and reported statistics per country, and Table 3 shows the reported and simulated IWR for major irrigated countries of the world.The results show generally good agreement with R 2 (the coefficient of determination) above 0.95 (p value < 0.001).Our estimates are also comparable to those of previous studies as shown in Table 1.With the CRU TS2.1 climate, our model tends to overestimate the IWR particularly in India, the USA, China, Pakistan, and Mexico.With the ERA-Interim and MERRA climate, we slightly overestimate IWR, but the magnitude is less compared to that of the CRU TS2.1 climate.With the ERA-Interim climate, IWR is generally overestimated over South and East Asia, e.g., India, Pakistan, China, Japan, and is underestimated over Europe, Africa, and South America, e.g., Spain, France, Germany, Egypt, South Africa, Brazil, and Argentina.With the MERRA climate, the overestimation is less obvious due to the wetter climate compared to the CRU TS2.1 and ERA-Interim climate, and our simulated IWR is rather underestimated over many regions, e.g., Europe, Africa, Asia except East Asia, and North America.When we use the average of the two or the three simulated IWRs, the correlation generally improves and the deviation between the simulated and reported values decreases.We thus used the average of the simulated results with the ERA-Interim and MERRA climate for the following analysis.

Accuracy of simulated surface water and groundwater withdrawal
Figure 3 and Table 4 show the comparison of our simulated water withdrawal per water source (surface water and groundwater), to reported country and state values for the year 2005 over the globe and for Europe, the USA, and Mexico.The comparison shows good agreement for both surface water and groundwater withdrawal over the globe (R 2 ≥ 0.96, p value < 0.001).However, our model tends to overestimate surface water withdrawal over South, Central, and East Asia (≈ +30 %), and tends to underestimate it over Southeast Asia and Africa (≈ −20 %).Simulated groundwater withdrawal shows good agreement with reported value over most of the regions of the world except Africa where the deviation is rather large (≈ ±30 %).Over Europe, the comparison shows reasonable agreement for surface water withdrawal and groundwater use with R 2 above 0.93 (p value < 0.001).However, our simulated surface water withdrawal is generally overestimated with α (the slope of regression line) being 0.85.Conversely, our simulated groundwater withdrawal is underestimated (α = 1.08).The overestimation of surface water withdrawal and the underestimation of groundwater withdrawal is large for the UK, and central and eastern Europe (> ±20 %) respectively.Over the conterminous USA and Mexico, the correlation is lower (R 2 < 0.9, p value < 0.001) compared to that over the global average and Europe, although regional variations of surface water and groundwater withdrawal are captured reasonably well.Our model generally overestimates both surface water and groundwater withdrawal for the central and eastern USA, whereas the deviation between the simulated and reported water use is smaller over the western USA.For Mexico, the comparison shows a contrasted trend compared to that of Europe in which surface water withdrawal is under-estimated, but groundwater withdrawal is overestimated over northern and southern Mexico.
In Fig. 4 we compare simulated and reported trends of groundwater withdrawal per country over the period 1980-2005 in 5 yr intervals when the reported statistics are available (the statistical data is not available before 1980).The comparison for 19 countries indicates that our approach is able to capture the decadal trends of groundwater withdrawal (R 2 > 0.95, p value < 0.001).The simulated trends of groundwater withdrawal match reasonably well with the observed trends not only for major groundwater users including the USA, China, and Mexico, but also for other countries including Poland, Greece, Spain, and Slovakia.However, the discrepancy between reported and observed trends tends to be larger for developed countries such as France, the UK, Austria, the Netherlands, and Finland.This suggests a limitation of our global application in which the partitioning between surface water and groundwater withdrawal represented by our approach needs further consideration or adjustment for these countries.

Regional trends of surface water and groundwater withdrawal and consumption
In Figs.  by growth in the agricultural sector (mostly irrigation), accounting for as much as ∼ 80 % of the total.Most of industrial and domestic water that is withdrawn from surface water and groundwater returns to river systems (40-80 %).Surface water and groundwater withdrawal increased respectively from ∼ 1350 and ∼ 650 km 3 yr −1 in 1979 to ∼ 2100 and ∼ 1200 km 3 yr −1 in 2010.During the period 1979-1990, groundwater withdrawal increased by ∼ 1 % per year, while surface water use rose by ∼ 2 % per year.However, during the recent period 1990-2010, the rate of groundwater withdrawal increased to ∼ 3 % per year, while that of surface water use decreased to ∼ 1 %.This is likely due to the fact that surface water has been extensively exploited in response to the consistent increase of global water demands, while the construction of new (large) reservoirs has been decreasing since the 1990s (Chao et al., 2008).The results suggest that the net increase in the demand has been mostly supplemented by groundwater withdrawal.These trends can also be seen  The regional trends of surface water and groundwater withdrawal and consumption exhibit very different trajectories over the period 1979-2010.Over Europe, groundwater withdrawal and consumption accounts for ∼ 30 % of the total and has not increased substantially over the past decades.However, over North and Central America, groundwater withdrawal and consumption account for ∼ 60 and ∼ 70 % of the total, and have increased by more than 40 % over the last 30 yr.Over western Asia, groundwater withdrawal has tripled and accounts close to ∼ 70 % of the total.Desalination water withdrawal accounts for 5 % of the total and is rapidly increasing over the region.Over North and Central America, and Asia, irrigation is the dominant water use sector and is predominantly relying on groundwater resources (∼ 70 %).Over South and East Asia, surface water and groundwater withdrawal nearly doubled from ∼ 600 and ∼ 360 km 3 yr −1 in 1979 to ∼ 1100 and ∼ 600 km 3 yr −1 in 2010, respectively.Total surface water and groundwater withdrawal over these regions accounts for more than half of the global surface water and groundwater withdrawal respectively.Over the other regions, e.g., Southeast Asia and South America, surface water withdrawal exceeds ∼ 80 % of the total except in North Africa where groundwater withdrawal is substantial (> 30 %).These trends are also visible from the development of consumptive use of surface water and groundwater (Fig. 6).

The impact of human-induced change on river discharge and terrestrial water storage change
Table 5 compares simulated river discharge under the pristine conditions (natural climate variability only) and under the human-induced change (human water use and reservoir operations) with observed river discharge taken from the selected GRDC stations (http://www.bafg.de/GRDC).For the comparisons, we selected major basins of the world that cover a wide range in climate and human impacts including reservoir regulation.Human-induced change is clearly observable for the rivers crossing major irrigated areas of the world, given the number of existing reservoirs, including the Nile, the Orange, the Murray, the Mekong, the Ganges, the Indus, the Yangtze, the Huang He, the Mississippi, the Columbia, and the Volga.For the other river basins, the human impact is less obvious, but still noticeable such as the Orinoco, the Parana, the Brahmaputra, the Danube, the Rhine, the Dnieper, and the Elbe.For the Amazon, the Congo, the Niger, the Zambezi, the Mckenzie, and the Lena, the river discharge is hardly affected because of small reservoir capacity and lower human water use.For those river basins where human impacts are large, the performance of the simulated river discharge under the pristine conditions tends to be lower compared to that of the simulated river discharge under the human-induced change, except for the Huang He where our overall model performance is low.Overall, the correlation between the simulated and the observed river discharge is high for most of the river basins, while the Nash-Sutcliffe model efficiency coefficient is high for some river basins but low for several basins including the Nile, the Niger, and the Orange, where the number of observation records are limited.Figure 7 compares the simulated monthly terrestrial water storage (TWS) anomalies with those of the GRACE observations (Liu et al., 2010) for a number of major river basins over the period 2003-2010.The selection of the basins is rather arbitrary, but is based on the fact that they are heavily affected by human activities, which enables to quantify the impact of human water use and reservoir operations on terrestrial water resources (e.g., surface water and groundwater).Simulated TWS was calculated from the sum of simulated snow, surface water, soil water, and groundwater storage.The TWS anomalies were computed over the overlapping period of 2003-2010 with the GRACE data.Here, we compared two simulation runs: one for pristine conditions (no human water use and no reservoirs) or natural climate variability only, and the other including human-induced change such as human water use (water withdrawal and consumptive water use) from surface water and groundwater storage, and reservoir operations.The comparison shows that human activities have a  noticeable impact on regional TWS signal and alter the seasonal and interannual TWS change.Over the Colorado and the Columbia basins, the seasonal TWS amplitude slightly decreased, which is explained by a combined effect of human water use and reservoir operations.The peak TWS signals are reduced due to human water extraction from surface and groundwater storage, while reservoirs release more water during the low flow period to satisfy the water demands downstream.The results indicate the large impact of human water use and regulation over these basins.Wang et al. (2011) also reported a large mass redistribution due to the presence of the Three Gorges Reservoir in the Yangtze Basin.Including human-induced change subsequently improves R 2 (between the simulated and observed TWS) from 0.75 to 0.80 (p value < 0.001) for the Columbia, but not for the Colorado where R 2 does not change substantially (∼ 0.65, p value < 0.001).Over the Mississippi and the Nile basins, human water use, primarily for irrigation purpose, decreases the peak TWS signals, which coincides with the crop growing season.Human water use extracts a large amount of water from groundwater and surface water storage, most of which evapotranspires over irrigated areas.This is less obvious for Table 5.Comparison of simulated to observed river discharge for selected major basins of the world.Simulated river discharge was derived under the pristine conditions (P; natural climate variability only) and under the human-induced change (H; human water use and reservoir operations).The results were obtained from the mean of the simulation with the ERA-Interim and MERRA climate.The observed river discharge was taken from the selected GRDC stations (http://www.bafg.de/GRDC)closest to outlets based on available records  for each basin.R 2 , α, and NSC denote the coefficient of determination, the slope (x coordinate: simulated discharge; y coordinate: observed discharge), and the Nash-Sutcliffe model efficiency coefficient (Nash and Sutcliffe, 1970)

The sensitivity of a water allocation scheme in simulated groundwater and surface water withdrawal
To evaluate our water allocation algorithm, we performed a sensitivity analysis for simulated groundwater and surface water withdrawal using three different scenarios.We recomputed the amount of groundwater and surface water withdrawal under each scenario, which specifies different wateruse behavior in groundwater and surface water resources including reservoirs.We calculated this amount at a spatial resolution of 0.5 • under each scenario, but the result of the global figure is presented in Fig. 8. Scenario 1 corresponds to our water allocation algorithm (Eq.13).To briefly recap our algorithm, the fraction of daily accumulated baseflow to the long-term average discharge was used to estimate the amount of water demand that is met from the renewable groundwater storage (S 3 ).The remaining water demand was then withdrawn from surface water availability including reservoirs.However, in case reservoirs are present at local or upstream grid cells over the river network, we first allocated surface water rather than groundwater to meet the water demand, and the remaining water demand was met from the renewable groundwater storage.Under Scenario 2, water demand is first met from the renewable groundwater storage (regardless of the presence of local or upstream reservoirs) and the remaining water demand is withdrawn from surface water availability.Under Scenario 3, water demand is first met from surface water availability and the remaining water demand is withdrawn from the renewable groundwater storage.In all scenarios, unmet water demand is imposed on (nonrenewable) groundwater, that is added to groundwater withdrawal.
The results show a clear difference in global trends of groundwater and surface water withdrawal.Under Scenario 2, the volume of groundwater withdrawal is more than 40 % larger than that simulated by our allocation algorithm (Scenario 1), while under Scenario 3, it is 16 % lower compared to the global groundwater withdrawal computed in our algorithm (Fig. 8a).The large difference in simulated groundwater withdrawals under these scenarios is also reflected in the difference in simulated surface water withdrawals (Fig. 8b).Note that most withdrawals occur over major irrigated regions that locate semiarid or arid regions where surface water availability is limited.Therefore, under all scenarios increase in surface water withdrawals is slowing down, but groundwater withdrawals are consistently increasing to supplement the unmet demand, particularly during the last decade.In terms of model performance, under Scenario 2 groundwater and www.earth-syst-dynam.net/5/15/2014/Earth Syst.Dynam., 5, 15-40, 2014 surface water withdrawals are likely very overestimated and underestimated respectively, while under Scenario 3 groundwater and surface water withdrawals are likely underestimated and overestimated respectively.

Discussion and conclusions
In this study, we coupled a global water demand model to a global hydrological model, and dynamically simulated daily water use, i.e., water withdrawal and consumptive water use, considering water allocation from surface water and groundwater resources.We implemented a new irrigation scheme, which interacts with the daily surface and soil water balance, and included a newly available extensive reservoir data set.
To simulate global water use, we used the newly available two climate reanalysis data sets (ERA-Interim and MERRA) over the period 1979-2010.The simulation period extended beyond most previous global analyses and the results provided new insights of the trends in global surface water and groundwater use over the recent decades.
To evaluate simulated water withdrawals, we compared our results with available reported statistics.Comparison of simulated IWR to reported statistics showed good agreement for most of the countries of the world.Although our model tends to overestimate IWR over some regions (e.g., Asia), the deviation is not substantial.The results showed substantial variability over country IWR depending on the climate input used (ERA-Interim and MERRA).As a result, we opted to use the average of the two simulated results for the subsequent analyses.However, compared to the ERA-Interim climate, the MERRA produces lower IWR due to the wetter climate over many regions, e.g., Europe, Africa, and North America.This subsequently lowered the mean of simulated water withdrawals and our water use numbers (Figs. 5,6).
Simulated water withdrawals per source (surface water and groundwater) were compared to reported statistics per country, and showed good agreement with R 2 above 0.93 (p value < 0.001).However, simulated surface water withdrawal was overestimated over Asia, and central and eastern Europe.Contrarily, groundwater withdrawal was underestimated over the same regions.To evaluate the spatial variability within a country, we then compared our estimates to reported subnational statistics.Results for the USA and Mexico showed that regional variations of surface water and groundwater withdrawal were captured reasonably well, although the correlation was lower compared to that for the country comparison.Comparison of simulated trends of groundwater withdrawal to reported trends also showed generally good agreement, but reported statistics were limited for only ∼ 20 countries of the world.Our simulated global groundwater withdrawal of ∼ 1000 km 3 yr −1 for 2000 is around the average when comparing to previous global estimates varying between ∼ 600 and ∼ 1700 km 3 yr −1 (Table 6).The large range is explained primarily by the use of different global hydrological models resulting in different runoff and water demand estimates that were used to calculate global groundwater abstraction (model based estimates).Validation of simulated consumptive water use (per source) remains difficult due to a lack of reliable information in many regions of the world.A recent study by Anderson et al. (2012) combined remotely sensed precipitation and satellite observations of evapotranspiration and groundwater depletion to estimate surface water consumption by irrigated agriculture in California's Central Valley.This approach may be promising and opens up new ways to measure surface water consumption, particularly over data poor regions.
A global and regional overview of water use showed a solid increase of surface water and groundwater use over the period 1979-2010.Global water withdrawal increased by more than ∼ 60 % from ∼ 2000 km 3 yr −1 in 1979 to ∼ 3300 km 3 yr −1 in 2010.The agricultural, mostly irrigation, sector accounts for as much as 80 % of the total.Surface water and groundwater withdrawals increased respectively from ∼ 1350 and ∼ 650 km 3 yr −1 in 1979 to ∼ 2100 and ∼ 1200 km 3 yr −1 in 2010, respectively.Although the decadal increase of water withdrawal decreased from ∼ 20 % during the 1990s to ∼ 14 % during the 2000s, water withdrawal has been consistently increasing over most of the regions of the world, e.g., Asia, and Central America, primarily due to a growing population and their water and food demand over the period 1979-2010.Our results also suggest that during the recent period 1990-2010 people have increasingly relied on groundwater, as surface water has been extensively exploited during the past periods.While readily accessible groundwater is an obvious choice to fill the gap between the increasing demand and limited surface water availability, this increasing dependence on groundwater likely worsens groundwater depletion already reported in various regions, e.g., northwestern India, northeastern Pakistan, northeastern China, the western and central USA, Mexico, and northern Iran (Konikow and Kendy, 2005;Rodell et al., 2009;Wada et al., 2010;Konikow, 2011;Famiglietti et al., 2011).
In our model, the allocation of surface water and groundwater resources to meet the demands is affected by the amount of simulated daily baseflow that is imposed over the long-term average discharge except for areas where local or upstream reservoirs are present.Water allocation is thus affected by the number of factors.For example, baseflow stems from the groundwater storage that is fed by daily groundwater recharge subject to seasonal and interannual climate variability.The number of upstream reservoirs affect the groundwater withdrawal downstream since surface water is allocated first to meet the demands in case reservoirs are present upstream.Moreover, groundwater withdrawal reduces the local groundwater storage and the associated baseflow, which in turn change the amount of groundwater withdrawal downstream due to the change in baseflow amount.Surface water withdrawal upstream also affects water availability downstream.These impacts of human water use and reservoirs are accumulated along the river network and become substantial over heavily affected basins (e.g., the Indus, the Colorado, and the Mississippi).Note that the temporal increase in simulated groundwater withdrawal is driven strongly by the increase in total water demand and the variability in surface water availability including reservoirs over the period 1979-2010.The analysis of simulated TWS anomalies revealed that human water use and associated reservoir operation alter the seasonal and interannual variability of TWS change.The alteration is particularly large over heavily regulated basins, e.g., the Colorado and Columbia basins, and over basins with major irrigated regions, e.g., the Mississippi, Indus, and Ganges basins.Including human water use generally improves the correlation of simulated TWS anomalies with those of the GRACE observations over basins (e.g., the Columbia, the Mississippi, and the Ganges).Note that the model performance is low over some of the simulation period for some basins, e.g., 2003 for the Indus, the Syr Darya, and the Euphrates, and 2006 for the Indus and the Euphrates.Nevertheless, the model reproduces TWS adequately for most of the basins.
To account for the climate uncertainty, we used two independent climate data sets: ERA-Interim and MERRA.However, model uncertainty can be large since model outputs can vary substantially among different global hydrological models (GHMs) with different model structure (Gosling et al., 2010(Gosling et al., , 2011;;Haddeland et al., 2011).Nevertheless, our simulated water use shows good agreement with reported statistics for most countries of the world.Moreover, our simulated TWS anomalies also show reasonable agreement with observed TWS data, but the comparison is limited to the selected major basins of the world where human impacts are substantial.
Our integrated modeling framework is capable of simulating human water use more realistically by including newly developed approaches.Our irrigation scheme is based on the surface and soil moisture deficit for paddy and nonpaddy fields respectively, and considers the feedback between daily irrigation and associated changes in surface and soil moisture condition.This, in turn, influences the amount of daily evapotranspiration on the same day and the soil moisture in following days over irrigated areas.Compared to earlier work (Wada et al., 2012a, b), the calculation of return flow from irrigation is more realistic considering the soil water balance and associated percolation losses underlying irrigated areas, which substantially contributes to groundwater recharge.Including the fine temporal dynamics of irrigation water requirement is critical as it concerns the daily water withdrawal to satisfy the demand, which affects the amount of water available downstream through river network.The model also calculates the other daily sectoral water demands including livestock, industry, and households, and considers the seasonal pattern of the demands (except the industry), while most of other GHMs calculate these demands on an annual basis.Our improved approach also includes a water allocation scheme that distinguishes surface water and groundwater withdrawals dynamically at a daily temporal resolution.Water is withdrawn and consumed from groundwater and available river discharge including the effect of reservoir operations that are parameterized using a newly available extensive global reservoir data set.Groundwater withdrawal also affects the amount of baseflow that is calculated with the available groundwater storage with the reservoir coefficient (see Appendix A), which in turn changes the amount of groundwater withdrawal downstream.These series of anthropogenic impacts on surface water and groundwater resources are reflected on the terrestrial water storage change, which Y. Wada et al.: Surface water and groundwater resources was well captured in our analysis using the GRACE observation.Thus, our new modeling framework enables one to comprehensively assess human-induced change in global water systems and to track those changes over time.The sensitivity analysis to assess our water allocation scheme showed that our scheme is preferred for the simulation of water use behavior compared to the other scenarios that prioritize either groundwater or surface water withdrawal first to meet the water demand.
However, there are certain limitations and major assumptions, that should be sufficiently addressed.First, allocation of surface water and groundwater to satisfy the sectoral water demands is currently simulated in a simplistic way using the faction of daily accumulated baseflow to the long-term average discharge and the number of available local and upstream reservoirs.Simulated (accumulated) baseflow is used to infer the readily extractable groundwater reserves.This assumption may be realistic in semiarid and arid regions where people largely rely on groundwater resources to satisfy the demands (e.g., northern India, Pakistan, northern China, Iran, and Mexico).However, there are likely discrepancies over regions where people predominantly rely on surface water resources despite the presence of the readily accessible groundwater reserves over shallow groundwater tables (e.g., humid regions).In addition, although the influence may not be large at the global scale, groundwater pumping is regulated in many developed countries (e.g., Japan, the Netherlands, Germany, and France), while water use management including pumping regulation is not accounted for in our model.Moreover, a realistic representation of the groundwater table considering lateral flow (Fan et al., 2013) is not incorporated, but it may substantially affect the simulation of groundwater storage and associated baseflow.Such process is important over aquifers with high transmissivity.Second, long-distance and cross-basin water diversions (e.g., aqueducts) provide additional surface water availability, which may substantially contribute to supply irrigation water requirement over regions such as the Indo-Gangetic plains, California's Central Valley, and the Colorado River, where extensive diversion works are present.Some information is available, e.g., the Periyar Project (maximum capacity: 40 m 3 s −1 ) and the Kurnool Cudappah Canal (maximum capacity: 85 m 3 s −1 ) in India, and the Irtysh-Karaganda Canal (maximum capacity: 75 m 3 s −1 in Central Asia (World Bank; http://www.worldbank.org/;UNDP; http://www.undp.org).However, artificial diversion networks and the actual amount of water transferred are difficult to be parameterized, and are not represented in our modeling framework due to limited data available worldwide.Third, we assumed that return flow from the industrial and domestic sectors to the river system occurs on the same day as water is withdrawn, but retention likely occurs due to waste water treatment, particularly in developed countries.Finally, although we used the recycling ratios developed on the basis of country GDP, the amount of water recycling to calculate consumptive water use is difficult to verify due to the lack of (sub-)national statistics over many regions of the world.
This study builds upon previous modeling efforts and contributes to improve a current modeling framework that quantifies the impact of anthropogenic impacts on global water resources.Despite its limitations, our modeling framework advances an important step beyond earlier work by attempting to account more realistically for the behavior of human water use and associated impacts on the terrestrial water system.It can be also used to assess future increase in water use per source due to population growth and economic development that will pose a serious threat to regions currently under substantial water scarcity and groundwater depletion, and to identify regions of looming water scarcity under future climate or under envisaged socioeconomic developments.

Appendix A
Here we present the essential and recently updated features of the model (Van Beek et al., 2011).The newly developed parts of the model are described in the main manuscript.The model is a grid-based model of global terrestrial hydrology (excluding Antarctica), which is essentially a leaky-bucket type of model, but a certain consideration is given representing the groundwater reservoir.

A1 Snow accumulation and melt
Over each grid cell, precipitation falls in the form of rain or snow.Any precipitation that falls on the soil surface can be intercepted by vegetation and in part or in whole evaporated.Snow accumulation and melt are temperature driven and modeled according to the snow module of the HBV model (Bergström, 1995).To represent the rain-snow transition over subgrid elevation dependent gradients of temperature, 10 elevation zones were made on each grid cell based on the HYDRO1k, and scaled the 0.5 • grid temperate fields with a lapse rate of 0.65 • C per 100 m.Over the 10 elevation zones, precipitation accumulates as snow if the temperature, T , is below the melt temperature (0 • C), T m .The snowmelt [m], SC melt , is then modeled using a degree day factor [m • C −1 day −1 ], f d : Above the melt temperature, precipitation and melt water are stored as liquid water in the available pore space in the snow cover.Melt water in the snow cover can refreeze depending on the water holding capacity of the snow (10 % of snow water equivalent) or evaporate.Snow is accumulated when the temperature is sufficiently low, otherwise it will melt and be added to the net liquid precipitation (P net ) that reaches the soil as rain or throughfall.Excess water from snowmelt and rainfall forms direct runoff or infiltrates into the first soil layer (S 1 ), which can further infiltrate into the second soil layer (S 2 ) and percolates into the third groundwater reservoir (S 3 ).

A2 Infiltration and direct runoff
Liquid water passed on from the snow cover to the soil surface will infiltrate if sufficient water storage is available, else it will drain over the surface as direct runoff.The partitioning into infiltration and direct runoff is dependent on the degree of saturation and the distribution of available water storage in the soil (Todini, 1996).
where W max is the total water storage capacity (SC S1 + SC S2 ) [m].W min is the minimum water storage capacity [m] according to the improved Arno scheme (IA) (Hagemann and Gates, 2003).W is the actual water storage (S S1 + S S2 ) [m].
To determine the water storages, the sum of the two upper soil layers and the average of each grid cell (0.5 • ) are considered.The parameter, b, is the dimensionless shape factor that defines the distribution of soil water storage within the cell.b is calculated based on the distribution of maximum rooting depths, which is derived from the 1 km by 1 km Global Land Cover Characteristics database version 2.0 (GLCC 2.0; http://edc2.usgs.gov/glcc/globe_int.php) and the land surface parameter data set (LSP2; Hagemann, 2002).To avoid an overestimation due to the zero (minimum) water storage capacity that always yields direct runoff when water falls onto the soil, the minimum water storage capacity was adjusted according to the IA (Hagemann and Gates, 2003).This allows a reasonable response time over which a minimum storage needs to be filled before any direct runoff will occur.Direct runoff, Q direct [m d −1 ], can occur for each time step [day] when P net [m d −1 ] falls over the surface and the sum of the actual water storage and the liquid precipitation exceeds W min .Moreover, any liquid precipitation is converted into direct runoff once pervious area is completely saturated.
where W max is the range between the maximum and the minimum water storage (W max -W min ) [m], and W is the range between the maximum and the actual water storage (W max -W ) [m].
where q i is the infiltration [m d −1 ] from P net to S 1 calculated from the difference between P net , and Q direct , augmented with any initial water storage to replenish the storage to W min .When the infiltration rate exceeds the saturated hydraulic conductivity (k sat ) of S 1 [m d −1 ], the infiltration excess is passed on to the direct runoff.

A3 Soil evaporation and transpiration
Soil moisture is liable to soil evaporation when the surface is bare and to transpiration when vegetated.Actual evapotranspiration is partitioned into soil evaporation and plant transpiration.Soil evaporation can occur from the top soil layer (S 1 ) over which the evaporation rate is limited by the saturated conductivity, k sat , in the saturated area, A sat (x), and by the unsaturated conductivity, k(θ E ), in the unsaturated area, A uns (1-x).The evaporation amount is considered separately for soil, E soil [m d −1 ], and melt water stored in snow cover, The actual T c is partitioned over the two soil layers (S 1 and S 2 ) according to the root fraction distribution, r f [dimensionless], in the two soil layers.
When the available soil moisture is limited to accommodate the overall fluxes, the overall fluxes are reduced proportionally to the amount of their original fluxes.

A4 Vertical water fluxes in soil and groundwater stores
Water moves through S 1 and S 2 and the third groundwater layer (S 3 Otherwise, the capillary flux is set to zero (θ E,S1 ≥ θ E,S2 ).
The unsaturated hydraulic conductivity of each soil layer is dependent on the effective degree of saturation (Clapp and Hornberger, 1978) and is calculated as where β is a dimensionless empirical exponent that varies on average between ∼ 4 and ∼ 11 over the range from sand to clay soil.It is based on a soil water retention curve parameter developed by Clapp and Hornberger (1978), who describe the relationship between the effective degree of saturation and the soil matric suction, ψ [m], as where ψ sat is the soil matric suction at saturation [m].Vertical water exchange, i.e., deep percolation, q p,S2→S3 [m d −1 ], and capillary rise, C r,S3→S2 [m d −1 ], between S 2 and S 3 is calculated in a similar way, except that the rate is given in the geometric mean of the unsaturated hydraulic conductivity of S 2 and the saturated hydraulic conductivity of S 3 : Capillary rise only occurs given the proximity of the water table and the resulting moisture content of S 2 cannot rise above field capacity, θ E FC,S2 (with ψ fc = 1.0 m).S S3 is the water storage in S 3 [m].f 5 m is the fraction of the grid cell with a groundwater depth within 5 m.Capillary rise is at maximum when the groundwater level is at the surface and it becomes zero when the groundwater level is below 5 m from the surface.The factor 0.5 is an estimate of the average capillary flux that occurs over the area fraction (f 5 m ) with a groundwater table within 5 m depth.The fraction (f 5 m ) is determined from the groundwater depth distribution.First all 1 km × 1 km grid cells are determined within the 0.5 • × 0.5 • grid that belongs to the perennial drainage network.Using the Perennial Inland Water Areas of the World (Vmap0; FAO, 1997), the average drainage density, D [m −1 ] is estimated (i.e., total length of perennial water courses divided by catchment area).For each of these grid cells, the upstream drainage area of the catchment is determined using the HYDRO1k.Next, taking the actual water levels (simulated by the model) of the perennial stream cells as a reference, the groundwater height [m], H = S S3 /f d (with f d drainable porosity or specific yield) is added to arrive at a local groundwater level and groundwater depth.From the groundwater depth distribution for each catchment, the area with a groundwater depth smaller than 5 m is determined.Adding these areas for all catchments and dividing by the total 0.5 • × 0.5 • grid cell area gives the estimates of the fraction (f 5 m ).

A5 Interflow or subsurface storm flow
In mountainous areas soils develop in regolith (unconsolidated solid material) covering the bedrock.The steep gradient of high (soil) to low (bedrock) hydraulic conductivity results in the occurrence of perched groundwater bodies during wet periods, which will cause a fast downslope flux of water through the soils down to the water courses.The model calculates this lateral drainage from S 2 over the height of the saturated wedge that may form over the contact with S 3 .This occurs when percolation from S 1 is high and percolation to S 3 is low, which results in the high gradient that drives lateral flow along the slope.This lateral drainage, known as interflow or subsurface storm flow, is modeled according to Sloan and Moore (1984): (q p,S1→S2,t + C r,S3→S2,t ) − (q p,S2→S3,t + C r,S2→S1,t ) , Earth Syst.Dynam., 5, 15-40, 2014 where L s is the average slope length [m] and t is the time step [day].TCL is the centroid-lag time or the characteristic response time [days].θ av is the available pore space based on the difference between the saturated (volumetric) moisture content (θ sat ) and the (volumetric) soil moisture content at field capacity (θ FC ) or θ av,S2 = θ sat,S2 − θ FC,S2 [all in m 3 m −3 ]. tan(α s ) is the gradient equal to the tangent of the slope angle [dimensionless].The average slope is determined from the average of calculated slopes from the 1 km × 1 km HYDRO1k within the 0.5 • × 0.5 • grid cell, excluding the lowest 10 % of the elevations, which are assumed to be part of the floodplain.
Interflow, Q int [m 3 d −1 ], is modeled along the slope (L s ) as a function of the lateral drainage (interflow) over the previous time step and the present recharge adding to or drawing from saturated wedge, taking into account the centroid-lag time.It is assumed that the saturated wedge responds to the recharge without delay and that it always will be draining with the available pore space.The recharge here is the net percolation, being the total of the gains to the second store due to percolation from the top soil layer and the capillary rise from the underlying groundwater layer as well as that of the losses due to percolation to the groundwater layer and the capillary rise to the top soil layer.Interflow is assumed to occur over areas with steep slopes and bedrock (i.e., mountainous areas).We calculated within the 0.5 • × 0.5 • grid cell the fraction of soils with a soil depth smaller than 1.5 m (maximum soil depth in the model) and used this as a proxy for the areas with the occurrence of interflow.Q int is calculated in terms of area per time step [m 2 d −1 ], but rather expressed as height per time step [m d −1 ] by dividing both terms of Eq. (A19) by L s [m].

A6 Baseflow
S 3 represents the deeper part of the soil that is exempt from any direct influence of vegetation and constitutes a groundwater reservoir fed by active recharge from S 2 .The active groundwater recharge consists of the net percolation that is calculated from the difference between the deep percolation (q p,S2→S3 ) and the capillary rise (C r,S3→S2 ).Drainage from the groundwater reservoir contributes as baseflow to the total river discharge.Groundwater discharge (baseflow) contributes an important part to streamflow in many parts of the world, particularly during low flow conditions.In the model, the groundwater storage and discharge are modeled by a first order linear reservoir approach.Baseflow, Q base [m d −1 ], is modeled by multiplying groundwater store, S S3 [m], and the reservoir or recession coefficient parameterized based on drainage theory (liner reservoir) developed by Kraaijenhoff van de Leur (1958), J [days].J is the spatial variable with a large regional variability, representing the average residence time of water in the groundwater store.S S3,t = S S3,t−1 − Q base,t−1 + (q p,S2→S3,t − C r.S3→S2,t ) (A21) Q base,t = S S3,t × J (A22) where k sat,S3 is the saturated hydraulic conductivity of the aquifer (S 3 ), f d is the drainable porosity, D c is the aquifer depth [m] and B c is the drainage length [m].The parameter B c is obtained from the drainage density analysis.The saturated hydraulic conductivity and drainable porosity have been related to a simplified version (7 classes) of the lithological map of the world (Dürr et al., 2005) and a literature search.Since there is no reliable information about aquifer thickness in relation to e.g., drainage distance and lithology, the aquifer thickness is arbitrarily assumed to be a constant of 50 m, this being the order of magnitude of the groundwater in contact with the surface water at the timescale of our simulations (several decades).By crossing the drainage length map with the lithological map and using the literature values, a global map of the global reservoir coefficient (groundwater residence time) can be estimated through Eq. (A23).This parameterization can be used as an initial estimate of global residence time, which can be further calibrated by comparing models results with low flows from discharge data and tuning k sat,S3 and f d for each lithological class.We refer to Sutanudjaja et al. (2011) for the sensitivity analysis of k sat,S3 and D c values of the PCR-GLOBWB model, and associated outcome.
Overall water balance for each soil layer can be written as S j = q p,j −1 + C r,j +1 − q p,j + C r,j + E soil,j + T c,j where S is the change in soil water storage for layer, j .q p is the percolation or infiltration (q i ) in case of the top soil layer, which is positive from the overlying store, negative towards the underlying store; C r is the capillary rise, which is negative towards the overlying store, positive when coming from the underlying store.E soil is the actual bare soil evaporation (limited to S 1 ), T c is the actual transpiration, Q is the lateral drainage (direct runoff, interflow or baseflow), and C is any sink due to human water consumption (all negative) [all in m d −1 ].
Local runoff, Q local [m d −1 ], is calculated as the sum of the direct runoff (Q direct ), the interflow (Q int ), and the baseflow (Q base ).
Local runoff is routed along the river network (see Sect.The parameterization of the vegetation and soil properties relies primarily on the GLCC 2.0, the FAO Digital Soil Map of the World (DSMW; FAO, 2003), and the WISE data set of global soil properties (ISRIC-WISE; Batjes, 2005).From the GLCC 2.0, the maximum rooting depth was used to obtain root content, the shape parameter b of the improved Arno scheme (Eq.A2), and the fractional vegetation cover and corresponding maximum interception storage capacity.From the DSMF and the ISRIC-WISE, soil properties including saturated hydraulic conductivity, saturated and residual (volumetric) water contents, porosity, air entry value, and coefficient β of the soil water retention curve were derived for each soil class for two different depths, i.e., from 0 to 30 cm (S 1 ) and from 30 to 150 cm (S 2 ).These values were first aggregated at the pedon level, where up to 8 soil classes and their fractional cover were specified per pedon at the spatial resolution of 0.5 • .The two soil layers represent the first and second store of the model except in those areas where soil formation is limited by bedrock or impeding layers, in which the two layers were reduced proportionally.For the third store of infinite capacity, the recession constant (J ) was estimated on the basis of the lithology and distance to the drainage network derived from the HYDRO1k, which was also used to determine the slope length (L s ) and slope tan(α s ).
In addition, the wilting point (θ E wp,j ) for each soil layer (j ) was calculated with matric suction [m] at wilting point (ψ wp ), matric suction [m] at air entry value (ψ ae,j ) according to Clapp and Hornberger (1978), and pore size distribution parameter (β j ) (varies on average between ∼ 4 and ∼ 11 over the range from sand to clay) according to Clapp and Hornberger (1978): )

Fig. 3 .
Fig. 3. Comparison of simulated total water withdrawals and water withdrawals per water source (surface water and groundwater) to reported values [km 3 yr −1 ] for the year 2005 over (a) the globe per country (N = 100), (b) Europe per country (N = 34), (c) the USA per state (N = 50), and (d) Mexico per state (N = 32) in log-log plots.Simulated water use at 0.5 • was spatially aggregated to country and state.Simulated value indicates the mean of the simulation with the ERA-Interim and MERRA climate.Error bars show standard deviation (σ )among the simulation with the ERA-Interim and MERRA climate.The dashed lines represent the 1:1 line.The reported water withdrawal per source was obtained from the FAO AQUASTAT database for the globe, from the Eurostat database (http://epp.eurostat.ec.europa.eu/portal/page/portal/environment/data/database) for Europe, from the US Geological Survey (Water Use in the United States; http://water.usgs.gov/watuse/) for the USA, and from the CONAGUA (Statistics on Water in Mexico; http://www.conagua.gob.mx/english07/publications/Statistics_Water_Mexico_2008.pdf) for Mexico.
from the global change in consumptive water use during the period 1979-2010.Siebert et al. (2010), Kummu et al. (2010), and Wada et al. (2012a) also report an increasing dependency of consumptive water use on groundwater resources in recent decades.

Fig. 4 .
Fig. 4. Comparison of simulated and reported trends of groundwater withdrawal per country over the period 1980-2005 (N = 19).The comparison is given in 5 yr interval according to the reported values including missing values for some years.Countries are identified with their ISO country codes: (a) the USA (USA) and China (CHN); (b) Mexico (MEX); (c) France (FRA), the UK (GBR), Poland (POL), Greece (GRC), and Spain (ESP); (d) Austria (AUT), Belgium (BEL), the Czech Republic (CZE), Finland (FIN), Israel (ISL), Luxemburg (LUX), Namibia (NAM), the Netherlands (NLD), Puerto Rico (PRI), Slovakia (SVK), and Sweden (SWE).Reported groundwater withdrawal was obtained from the FAO AQUASTAT database.Simulated value indicates the mean of the simulation with the ERA-Interim and MERRA climate.Error bars show standard deviation (σ ) among the simulation with the ERA-Interim and MERRA climate.The dashed line represents the 1 : 1 slope.

Fig. 5 .
Fig. 5. Regional trends of water withdrawal per source (desalination water, surface water, and groundwater) over the period 1979-2010.The results were obtained from the mean of the simulation with the ERA-Interim and MERRA climate.The global figure is shown in the left corner.

Fig. 6 .
Fig. 6.Regional trends of consumptive water use per source (desalination water, surface water, and groundwater) over the period 1979-2010.The results were obtained from the mean of the simulation with the ERA-Interim and MERRA climate.The global figure is shown in the left corner.

Fig. 7 .
Fig. 7. Comparison of simulated monthly TWS anomalies with those of the GRACE observations [m] for selected major basins over the period 2003-2010.The results were obtained from the mean of the simulation with the ERA-Interim and MERRA climate.Black solid line, blue dashed line, and red dashed line indicate the GRACE observation, pristine condition (natural climate variability only), and human-induced change (water use and reservoir operations), respectively.Monthly GRACE terrestrial water storage anomaly data were obtained from the DEOS Mass Transport release 1/1b (DMT-1) model(Liu et al., 2010).

Fig. 8 .
Fig. 8. Sensitivity of water allocation scheme in simulated global (a) groundwater and (b) surface water withdrawals [km 3 yr −1 ] over the period 1979-2010 based on three scenarios.Scenario 1 correspondsto the water allocation algorithm of this study (Eq.13).Under Scenario 2, water demand is first met from the groundwater storage (S 3 ) and the remaining water demand is withdrawn from surface water (including reservoirs).Under Scenario 3, water demand is first met from surface water and the remaining water demand is withdrawn from the groundwater storage.The results were obtained from the mean of the simulation with the ERA-Interim and MERRA climate.
2.4).Water can be lost from the drainage network by evaporation or human water consumption.

Wada et al.: Surface water and groundwater resourcesTable 1 .
Previous global studies to simulate irrigation water requirement (IWR); "corr" stands for corrected.

Table 2 .
Correlation of simulated IWR to reported statistics per country for the year 2000 (N = 212).IWR was simulated with the CRU TS2.1 (C), ERA-Interim (E), and MERRA (M) climate, respectively.Average indicates the mean of the two or three results.Reported statistics were obtained from the FAO AQUASTAT data base (globe: 2434 km 3 yr −1 ).R 2 and α denote the coefficient of determination and the slope of regression line respectively.R 2 was derived from the comparisons between normal values.The value with the CRU TS2.1 climate is provided for a reference and is not included in our overall analysis.The values of irrigation water consumption (IWC) is also provided under each climate.

Table 3 .
Comparison of simulated IWR to reported statistics [km 3 yr −1 ] for major irrigated countries of the world for the year 2000 (N = 212).IWR was simulated with the CRU TS2.1 (C), ERA-Interim (E) and MERRA (M) climate respectively.Reported statistics was obtained from the FAO AQUASTAT database (http://www.fao.org/nr/water/aquastat/main/index.stm).Simulated IWR with the CRU TS2.1 climate is provided for a reference and is not included in our overall analysis.

Table 4 .
Correlation between simulated and reported water withdrawals per source (TWW: total water withdrawal, SWW: surface water withdrawal, GWW: groundwater withdrawal) for the year 2005 over the globe per country (N = 100), Europe per country (N = 34), the USA per state (N = 50), and Mexico per state (N = 32) in log-log plots.R 2 and α denote the coefficient of determination and the slope of regression line respectively.R 2 was derived from the comparisons between normal values.
The reported water withdrawal per source was obtained from the FAO AQUASTAT database for the globe, from the Eurostat database (http://epp.eurostat.ec.europa.eu/portal/page/portal/environment/data/database) for Europe, from the US Geological Survey (Water Use in the United States; http://water.usgs.gov/watuse/) for the USA, and from the CONAGUA (Statistics on Water in Mexico; http://www.conagua.gob.mx/english07/publications/Statistics_Water_Mexico_2008.pdf) for Mexico.
. This is due to the fact that the low flow periods coincide with the growing season of irrigated crops (spring), which require large amounts of water.Irrigation water use thus decreases both surface water and groundwater storage during the low flow season.This improves R 2 from 0.85 to 0.90 (p value < 0.001) for the Ganges Basin.The impact of reservoir operations is less obvious over the Ganges.Over the Syr Darya and the Euphrates basins, similar to most of the basins, human water use decreases the seasonal amplitude of TWS change, but does not substantially improve the correlation between the simulated and observed TWS.
hibits very different interannual trends over the years, which are captured reasonably well by our model.Major reservoirs are mostly located in upstream regions of the basin and release water during the crop growing season in spring and summer.Over the Ganges Basin, contrary to the other basins, human water use increases the seasonal amplitude of TWS change.

Table 6 .
Global estimates of groundwater abstraction.
). Vertical water fluxes [m d −1 ] between S 1 and S 2 are driven by the effective degree of saturation of both layers, θ E,S1 = S S1 /SC S1 and θ E,S2 = S S2 /SC S2 or θ E,S1 = θ S1 /θ sat,S1 and θ E,S2 = θ S2 /θ sat,S2 .θ is the effective moisture content defined as the fraction of water storage to soil depth (θS1 = S S1 /Z S1 and θ S2 = S S2 /Z S2 ).Percolation, q p,S1→S2 [m d −1 ],from S 1 to S 2 is governed by the rate of unsaturated hydraulic conductivity of S 1 , k(θ E ) S1 [m d −1 ] if sufficient soil moisture is available.However, when θ E,S1 < θ E,S2 , upward capillary (capillary rise) can occur at the rate of unsaturated hydraulic conductivity of the second soil layer, k(θ E ) S2 [m d −1 ], but is limited to the portion of soil moisture deficit of the overlaying layer (1-θ E,S1