An integrated model for the assessment of global water resources – Part 1: Model description and input meteorological forcing

To assess global water availability and use at a subannual timescale, an integrated global water resources model was developed consisting of six modules: land surface hydrology, river routing, crop growth, reservoir operation, environmental flow requirement estimation, and anthropogenic water withdrawal. The model simulates both natural and anthropogenic water flow globally (excluding Antarctica) on a daily basis at a spatial resolution of 1 × ◦ (longitude and latitude). This first part of the two-feature report describes the six modules and the input meteorological forcing. The input meteorological forcing was provided by the second Global Soil Wetness Project (GSWP2), an international land surface modeling project. Several reported shortcomings of the forcing component were improved. The land surface hydrology module was developed based on a bucket type model that simulates energy and water balance on land surfaces. The crop growth module is a relatively simple model based on concepts of heat unit theory, potential biomass, and a harvest index. In the reservoir operation module, 452 major reservoirs with >1 km3 each of storage capacity store and release water according to their own rules of operation. Operating rules were determined for each reservoir by an algorithm that used currently available global data such as reservoir storage capacity, intended purposes, simulated inflow, and water demand in the lower reaches. The environmental flow requirement module was newly developed based on case studies from around the world. Simulated Correspondence to: N. Hanasaki (hanasaki@nies.go.jp) runoff was compared and validated with observation-based global runoff data sets and observed streamflow records at 32 major river gauging stations around the world. Mean annual runoff agreed well with earlier studies at global and continental scales, and in individual basins, the mean bias was less than±20% in 14 of the 32 river basins and less than ±50% in 24 basins. The error in the peak was less than ±1 mo in 19 of the 27 basins and less than ±2 mo in 25 basins. The performance was similar to the best available precedent studies with closure of energy and water. The input meteorological forcing component and the integrated model provide a framework with which to assess global water resources, with the potential application to investigate the subannual variability in water resources.


Abstract.
To assess global water availability and use at a subannual timescale, an integrated global water resources model was developed consisting of six modules: land surface hydrology, river routing, crop growth, reservoir operation, environmental flow requirement estimation, and anthropogenic water withdrawal. The model simulates both natural and anthropogenic water flow globally (excluding Antarctica) on a daily basis at a spatial resolution of 1 • ×1 • (longitude and latitude). This first part of the two-feature report describes the six modules and the input meteorological forcing. The input meteorological forcing was provided by the second Global Soil Wetness Project (GSWP2), an international land surface modeling project. Several reported shortcomings of the forcing component were improved. The land surface hydrology module was developed based on a bucket type model that simulates energy and water balance on land surfaces. The crop growth module is a relatively simple model based on concepts of heat unit theory, potential biomass, and a harvest index. In the reservoir operation module, 452 major reservoirs with >1 km 3 each of storage capacity store and release water according to their own rules of operation. Operating rules were determined for each reservoir by an algorithm that used currently available global data such as reservoir storage capacity, intended purposes, simulated inflow, and water demand in the lower reaches. The environmental flow requirement module was newly developed based on case studies from around the world. Simulated Correspondence to: N. Hanasaki (hanasaki@nies.go.jp) runoff was compared and validated with observation-based global runoff data sets and observed streamflow records at 32 major river gauging stations around the world. Mean annual runoff agreed well with earlier studies at global and continental scales, and in individual basins, the mean bias was less than ±20% in 14 of the 32 river basins and less than ±50% in 24 basins. The error in the peak was less than ±1 mo in 19 of the 27 basins and less than ±2 mo in 25 basins. The performance was similar to the best available precedent studies with closure of energy and water. The input meteorological forcing component and the integrated model provide a framework with which to assess global water resources, with the potential application to investigate the subannual variability in water resources.

Introduction
Water is one of the most fundamental resources for humans and society. Rapid growth of the world population and economy has brought major increases in water demand during the 20th century, and this trend is projected to continue in the 21st century (Shiklomanov, 2000). Several global water resource assessments were released to project the current and future distributions of water-deficient areas worldwide (Arnell, 1999(Arnell, , 2004Vörösmarty et al., 2000;Oki et al., 2001Oki et al., , 2003Alcamo et al., 2003aAlcamo et al., , 2003bAlcamo et al., , 2007Oki and Kanae, 2006). These assessments used global hydrological models to estimate the distribution of runoff (renewable freshwater) and various world statistics to estimate water use. A typical approach is to display the global distribution of per capita Published by Copernicus Publications on behalf of the European Geosciences Union. annual water availability or the ratio of withdrawal to availability on an annual basis. However, extreme seasonality in both water availability and water use occurs in some parts of the world. For example, in the Asian monsoon region, conditions change dramatically between the rainy and dry seasons. Moreover, global warming is projected to alter future temperature and precipitation patterns and consequently affect both the amount and timing of water availability and use (Kundzewicz et al., 2007). Therefore, subannual variability must be taken into account.
A model suitable for such assessments requires the following three functions. First, it must simulate both renewable freshwater resources and water use at a subannual timescale. Second, it must deal with major interactions between the natural hydrological cycle and anthropogenic activities. For example, withdrawal from the upper stream affects availability in the lower stream, and reservoir operation may contribute to increased water availability in the lower stream. Third, it must explain key mechanisms regarding the effects of global warming on water availability and water use for future projections.
Several integrated global water resources models that can simulate not only the natural water cycle, but also anthropogenic water flow, have been published. Alcamo et al. (2003aAlcamo et al. ( , 2003b developed a global water assessment model called "WaterGAP 2" which consists of a global water use model and a global hydrology model, and assessed the current and future water resources globally. Haddeland et al. (2006) developed and implemented a reservoir model and an irrigation model in the Variable Infiltration Capacity (VIC) land surface model (Liang et al., 1994) and studied the effects of reservoirs and irrigation water withdrawal on continental surface water fluxes for part of North America and for Asia. Jachner et al. (2007) enhanced the LPJmL (Lund-Potsdam-Jena managed land) dynamic global vegetation model  with a river routing model, including lakes and reservoirs, and withdrawals for households and industry, and assessed how much water is consumed by global irrigated and rain-fed agriculture and by natural ecosystems. Several other global hydrological models have been developed, but most have focused on the natural hydrological cycle, with less emphasis on anthropogenic aspects.
In contrast to earlier works, we set three basic policies. First, our primary purpose was to assess global water availability and use at a subannual timescale. No previous studies set this as their primary purpose. Second, both water and energy balances on the land surface are closed in our model. This is not only the most fundamental consideration of hydrology, but is also one of the key requirements for the interdisciplinary coupling of submodules (e.g., hydrological models and crop growth models). Recently, several advanced earth system modeling efforts have been reported such as coupling a land surface model (LSM) with a crop model (Gervois et al., 2004;Mo et al., 2005) and coupling a global climate model (GCM) with a crop model (Osborne et al., 2007). In these systems, soil moisture, evaporation, and other variables are shared by more than one submodel; therefore, to maintain consistency among submodels, energy and water balances should be conserved. In particular, the closure of the energy balance is a fundamental requirement of GCM and LSM approaches. Third, as much as possible, we tried to avoid model calibration involving the fit of simulated results to available observation records. Only two hydrological parameters were tuned by climatic zones, not individual basins (Sect. 3.1). It is well established that hydrological models do not reproduce observed hydrographs very well without model calibration (or model parameter tuning). However, in global-scale hydrological modeling, model calibration is a difficult issue. There are a few reasons for this. First, it is virtually impossible to calibrate the model worldwide because of the limited availability of observations, especially in developing countries. Second, both models and input meteorological forcing and validation data contain considerable uncertainty (Oki et al., 1999), and it is not always easy to attribute errors in simulations to improper settings of model parameters. Moreover, we intended to apply the model to future projection under climate change. Thus, the transparency and physical validity of the model are quite important because the simulated results are highly model dependent. Therefore, we extensively examined the simulated results of the model using model inherent parameters; even this sometimes produces large errors.
We developed an integrated global water resources model consisting of six modules: land surface hydrology, river routing, crop growth, reservoir operation, environmental flow requirements, and anthropogenic water withdrawal. The model simulates both natural and anthropogenic water flow globally (excluding Antarctica) at a spatial resolution of 1 • ×1 • (longitude and latitude) at a daily time interval.
Six modules were selected for the following reasons. First, to estimate renewable freshwater resources, models for land surface hydrology and river routing are indispensable. The estimation of agricultural water demand is particularly important in the assessment of global water resources because 85% of consumptive water use is used for agriculture (Shiklomanov, 2000), with considerable seasonality because the water is needed only during cropping periods. The crop growth module estimates the timing and amount of the irrigation water requirement. Reservoirs buffer the seasonal variation in streamflow and fill the gap between the timing of streamflow and water demand. Therefore, the reservoir operation module plays an important role in a subannual assessment. There is no doubt about the need for the anthropogenic withdrawal module. One of the most basic dynamics in water resources problems is that withdrawal from the upper stream limits that from the lower stream. Global modeling of environmental flow requirements is a relatively new topic in global hydrology. However, one cannot withdraw water from a channel until it completely runs out. Although there are other issues to consider, we started with these six modules that we judged to be most essential for our goal.
The modeling of anthropogenic activities is in its very initial stage in terms of global hydrological modeling. We do not expect that our model can reproduce individual events in the real world. What we introduced into our model was the minimum basic anthropogenic activities in current global hydrological models; such basic anthropogenic activities are indispensable in analyzing the seasonality of water availability and water use. Here, we assume that humans act rationally. We do not expect that our model can reproduce the daily operation of individual reservoirs or daily irrigation practices at individual farms in the real world. However, reservoir operators seldom release water in floods, and farmers seldom sow in periods that are unsuitable for cropping. Historical reservoir operations were fairly reproduced using a simplistic model that assumes rational actions by humans .
There are two potential beneficiaries of this model: the climate change impact assessment community and the earth system modeling community. A number of global water resources models contributed climate change impact projections to the fourth assessment report (AR4) of the Intergovernmental Panel on Climate Change (IPCC; Kundzewicz et al., 2007). The global water resources assessments of AR4 were on an annual basis and projected the effects of annual to decadal changes in precipitation and temperature on water resources. However, the effects of subannual change (e.g., decrease in snowfall, earlier snowmelt, increase in the intensity and frequency of heavy precipitation, and change in the timing of monsoon onset) were not examined explicitly. Our model calculates both water availability and use at a daily interval and has few tuning parameters (i.e., parameter tuning is impossible for future projection). This model will thus provide impact assessment for a currently missing time range. Also, the model has the potential to benefit the earth system modeling community (atmosphere-ocean-land-carbon coupled model). Earth system researchers have recently started to make anthropogenic activities such as irrigation and reservoir operation an important component of their earth system models (Boucher et al., 2004;Lobell et al., 2006;Osborne et al., 2007). Our model closes water and energy balances on the land surface, which is a fundamental basis of the earth system models. Thus, our methodology and results will be directly applicable and comparable to the results of earth system models.
In this two-part report, we introduce the integrated global water resources model and use the model to assess global water resources. Here, we describe the input meteorological forcing and the six modules (i.e., land surface hydrology, river routing, crop growth, reservoir operation, environmental flow requirement estimation, and anthropogenic water withdrawal). In modeling and simulations, the preparation of reliable meteorological forcing inputs is essential. First, we revisited the original meteorological forcing inputs of the second Global Soil Wetness Project (GSWP2). We traced some of its shortcomings and developed improved meteorological forcings, as described in Sect. 2. In Sect. 3, we present the six modules. Finally, in Sect. 4, we discuss the validation of the simulated runoff and streamflow, confirming that the global hydrological cycle was properly reproduced (Sect. 4). In a forthcoming paper , we present the results of the model application and global water resources assessments, which focused on subannual variation in water availability and water use.
Here, "runoff" indicates the water that drains from surfaces and subsurfaces of a certain area of land [mm yr −1 or mm mo −1 ]. "Streamflow" indicates the flow of water in river channels [m 3 s −1 ].

Meteorological forcing input
The simulation was conducted using the framework of the second Global Soil Wetness Project (GSWP2; Dirmeyer et al., 2006), which is an international project that estimated the global energy and water balance over land, with emphasis on variation in soil moisture. This framework has two significant benefits. First, it provides quality-checked meteorological forcing input (e.g., air temperature and precipitation) and consistent surface boundary conditions (e.g., land-sea mask and albedo) with which to simulate energy and water balances globally. Second, it allows for the comparison of our model with the 15 state-of-the-art land surface models involved in the GSWP2.
Some preceding studies have validated the output of the GSWP2 using field observations and have evaluated the meteorological forcing inputs and the performance of the participating LSMs. Decharme and Douville (2006) validated the precipitation forcing input and runoff output of the GSWP2 at the Rhône River basin, where reliable observations are available, and at 80 gauging stations distributed over the globe. They showed that the runoff of the GSWP2 was overestimated at middle and high latitudes, which was commonly observed in all of the LSMs participating in GSWP2, and attributed this overestimation to the precipitation forcing input. These findings indicate that the meteorological forcing input of GSWP2 needs to be revisited. First, we revisit the original meteorological forcing input of the GSWP2 (hereafter F-GSWP2-B0; F stands for forcing, B0 for baseline experiment version zero). In the GSWP2, seven input meteorological components were provided to drive the LSMs: downward longwave radiation, downward shortwave radiation, wind speed, surface air pressure, specific humidity, air temperature, and precipitation. All of these components were provided for 10 years (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995) at 3-h intervals, covering all land excluding Antarctica at a spatial resolution of 1 • ×1 • . F-GSWP2-B0 is a hybrid product of the NCEP-DOE reanalysis (Kanamitsu et al., 2002) and various observation-based monthly meteorological data on a global grid (see Table 1). GPCC was used for grids where rain gauges were densely located, whereas GPCP was used for grids where they were sparsely located. b GPCP was not used. c The algorithm of Motoya et al. (2002) and NCEP-DOE wind speed at the height of 10 m were used. d The algorithm of Motoya et al. (2002) and ERA40 wind speed at the height of 2 m (originally 10 m) were used. e Daily maximum and minimum temperature changes were scaled linearly by CRU data. f Adjusted to corrected air temperature so that the relative humidity of the original reanalysis was conserved. g Adjusted to ISLSCP2 elevation. The NCEP-DOE reanalysis was corrected linearly to match the monthly mean values to the observation-based data. For the precipitation data, an algorithm for the gauge correction of wind-caused undercatch was applied to the rainfall and snowfall input data (Motoya et al., 2002). The methodology for producing these components has been described in detail by Zhao and Dirmeyer (2003), and a short description is provided in Appendix A.
To revisit the findings of Decharme and Douville (2006), the global zonal mean distributions of wind speed and precipitation of F-GSWP2-B0 are provided (Fig. 1). As a yardstick, the mean 1961-1990 global observation-based data of the Climate Research Unit at the University of East Anglia (CRU; New et al., 1999) are provided for wind speed and precipitation. The precipitation of F-GSWP2-B0 is clearly larger than that of the CRU at middle to high latitudes, but smaller at low latitudes. One possible cause of the large precipitation of F-GSWP2-B0 is its wind speed, which is much higher than that of the CRU except at low latitudes. Motoya et al.'s (2002) undercatch correction is correlated with wind speed (see Appendix A), especially in regions at high latitudes in which precipitation is dominated by snow. Moreover, we obtained the original source program code that was used to produce F-GSWP2-B0 precipitation data and found that wind speed at a height of 10 m was used, whereas Motoya et al.'s (2002) algorithm expects a height of 2 m. This is another cause of overcorrection in F-GSWP2-B0 precipitation data because wind speed is stronger at higher altitudes.
To correct and improve the accuracy of F-GSWP2-B0, a new input meteorological data set (F-GSWP2-B1) was developed. Table 1 lists the differences in each variable between the two data sets. The major differences included a change in reanalysis data from NCEP-DOE (Kanamitsu et al., 2002) to ERA40 (Betts and Beljaars, 2003) and the proper application of Motoya et al.'s (2002) wind correction to precipitation data. Tanaka et al. (2005) compared air temperature, vapor pressure, wind speed, and precipitation of NCEP-DOE and ERA40 with daily ground meteorological observations collected by the National Climatic Data Center (http://www.ncdc.noaa.gov) for 1994-1995 at 2349 stations around the world and found better agreement with ERA40 than with NCEP-DOE. Because daily and diurnal variation in the GSWP2 meteorological forcing inputs is dependent on the reanalysis data, we substituted the ERA40 data for the NCEP-DOE data.
The wind speed of F-GSWP2-B1 is more similar to that of the CRU than the F-GSWP2-B0, but it is smaller than that of the CRU at southern low latitudes (Fig. 1). For precipitation, F-GSWP2-B1 has greater precipitation at latitudes higher than 50 • N in the Northern Hemisphere and higher than 35 • S in the Southern Hemisphere because of the undercatch correction, but the difference from the CRU is much smaller than that from the F-GSWP2-B0 (Fig. 1).

Land surface hydrology module
A land surface hydrology module calculates the energy and water budget, including snow, on the land surface from the forcing data. This module is based on a bucket model (Manabe, 1969;Robock et al., 1995), but differs from the original formulation in the following three aspects. First, soil temperature is calculated using the force restore method (Bhumralkar, 1975;Deardorff, 1978) to simulate the diurnal cycle of surface temperature reasonably using three-hourly meteorological forcing inputs. Second, a simple subsurface runoff parameterization is added to the model. Third, two independent land surface conditions can be simulated within a single grid that is intended to separate irrigated cropland from other land types. The bucket model is simple, but is still widely used in current global water hydrological studies. Soil moisture is expressed as a single-layer bucket 15 cm deep for all soil and vegetation types. When the bucket is empty, soil moisture is at the wilting point; when the bucket is full, soil moisture is at field capacity. Evapotranspiration is expressed as a function of potential evapotranspiration and soil moisture. In the original bucket model, runoff is generated only when the bucket is overfilled, but we used a "leaky bucket" formulation in which soil moisture drains continuously. Potential evapotranspiration and snowmelt are calculated from the surface energy balance. A detailed description of this module can be found in Appendix B.
At first, the parameters of the land surface hydrology module were set as globally uniform (Appendix B). However, there was substantial regional bias in runoff and streamflow simulations. Therefore, we set the parameters of the land surface hydrology modules for four climatic zones from an analysis of energy and water constraint (Appendix C). Hereafter, only the results obtained using the modified parameters are shown for clarity of discussion.

River routing module
The river module is identical to the Total Runoff Integrating Pathways model (TRIP; Oki and Sud, 1998;Oki et al., 1999). The module has a digital river map with a spatial resolution of 1 • ×1 • (the land-sea mask is identical to the GSWP2 meteorological forcing input) and a flow velocity fixed at 0.5 m s −1 . The module accumulates runoff calculated by the land surface hydrology module and outputs streamflow. This module does not deal with lakes or swamps, human-made reservoir operation, diversion or withdrawal, or evaporation loss from water surfaces. Human-made reservoir operation and withdrawal are modeled in the respective modules.

Crop growth module
We coded a crop growth module with reference to the Soil Water Integrated Model (SWIM; Krysanova et al., 1998Krysanova et al., , 2000. The SWIM model is an eco-hydrological model for regional impact assessments in mesoscale watersheds (100-10 000 km 2 ). The model integrates hydrology, vegetation, erosion, and nitrogen dynamics at the watershed scale. We used only the formulation and parameters of crop vegetation. The SWIM model can deal with more than 50 types of crops. For 18 of these crop types, Leff et al. (2004) provided the global distribution of the areal proportion at a 0.5 • ×0.5 • spatial resolution. The remaining crops were simulated using a generic parameter set. Because the SWIM model is a descendant of the Erosion Productivity Impact Calculator (EPIC) model (Williams, 1995) and the Soil Water Assessment Tool (SWAT) model (Neitsch et al., 2002), the formulation and parameters of the crop module are quite similar to those of the earlier models (Appendix D).
The agricultural water demand was assumed to be equal to the volume of water needed to maintain soil moisture at 75% of field capacity in irrigated fields. Above this threshold, the land surface hydrology module's evaporation is identical to the potential evaporation, and consequently, the water stress factor of the crop growth module is zero. Below this threshold, the water stress prevents the optimal growth of plant biomass and crop yield (see Appendix D). In the case of rice, soil moisture was maintained at saturation water content to express paddy inundation. Irrigation began 30 days prior to the planting date, increasing the soil moisture content linearly from 0% to 75%. Otherwise, heavy irrigation was required at the planting date to increase the soil moisture to at least 75%, especially in arid areas. If the soil moisture was above the target content, irrigation water was not applied. To identify the proportion of irrigated area in each grid cell, the global map of Döll and Siebert (2000) was used. This map has a spatial resolution of 0.5 • ×0.5 • with its original landsea mask. We re-gridded it to 1 • ×1 • spatial resolution with the common GSWP2 land-sea mask (Dirmeyer et al., 2006) used in this model. Table 2. Environmental flow requirements of the environmental flow requirement module of the integrated model; q indicates monthly runoff.

Regional classification
Monthly environmental flow requirement (q env ) Description Minimum Maximum Condition Monthly monthly streamflow monthly streamflow requirement Wet (wet throughout a year) 10≤q min 100≤q max q env =0.3q+q flood Stable (stable throughout a year) 1≤q min q max <100 q env =0.1q Variable (dramatic difference Other than above 0≤q<1 q env =0 between rainy and dry seasons) 1≤q<10 q env =0.1q 10≤q q env =0.3q+q flood A simplified cropping pattern was assumed because of limited detailed information on global cropping practices. We assumed that the same crop species was planted in both irrigated and nonirrigated croplands. The global distribution of crop species was obtained from Leff et al. (2004). To further simplify the simulation, only information on primary and secondary crop types in terms of the cultivated area reported by Leff et al. (2004) was used. We then assumed that the primary crop was cultivated as the first crop, and the secondary crop was cultivated as the second crop. The crop intensity of irrigated cropland (the areal proportion of cultivated area to the total irrigated area) was obtained from Döll and Siebert (2002). According to their estimation, crop intensity varied from 0.8 (i.e., on average, 80% of the total irrigated cropland is used) in parts of the former USSR, Baltic states, and Belarus to 1.5 in eastern Asia, Oceania, and Japan. For the former group of countries, we assumed that only 80% of the irrigated land was cultivated for the first crop and that no second crop was planted. For the latter group, we assumed 100% cultivation for the first crop and 50% for the second crop.

Reservoir operation module
In the global river map of the river routing module, the 452 largest reservoirs with storage capacity >10 9 m 3 each worldwide were geo-referenced, and available reservoir information (e.g., name, purposes in order of priority, and storage capacity) was compiled . The total storage capacity of these 452 reservoirs was 4140 km 3 , accounting for approximately 60% of the total reservoir storage capacity in the world (ICOLD, 1998). The reservoir operation module set operating rules for individual reservoirs. For reservoirs for which the primary purpose was not irrigation water supply, the reservoir operating rule was set to minimize interannual and subannual streamflow variation (i.e., storage capacity and inflow). For reservoirs for which irrigation water supply was the primary purpose, daily release from the reservoir was proportional to the irrigation water requirement in the lower reaches (Appendix E).

Environmental flow requirement module
We estimated a monthly environmental flow requirement using the algorithm of Shirakawa (2004Shirakawa ( , 2005. Because both reports by Shirakawa were published in Japanese, we first describe Shirakawa's methodology. First, using the 10-year mean monthly gridded streamflow data simulated by the land surface hydrology module and the river routing module, all grids were classified into four regions following specific criteria ( Table 2). The environmental flow requirement consisted of two factors: the base requirement, which was the minimum streamflow in the channel; and the perturbation requirement, which allowed flush streamflow in the rainy season. The perturbation requirement was 10% of the mean monthly streamflow and should occur for 2-3 days. However, considering the spatiotemporal resolution of our study, the perturbation requirement was not implemented; instead, an allocated amount was simply added to the base requirement (Appendix F).
We did not require that the environmental flow requirement estimated by this algorithm be sufficient for aquatic ecosystems. Unless withdrawal from and pollution of rivers is prohibited, an ecosystem will be more or less affected and changed by human activities. Whether the change or damage is tolerable is highly dependent on societal value judgments. The algorithm was based on case studies and fieldwork in semi-arid to arid regions or in heavily populated regions, which provided actual examples of value judgments in relatively water-scarce regions. There is no universally accepted theory regarding environmental flow requirements. Indeed, value judgments are largely influenced by regional welfare and culture. However, because of limited information regarding global applicability, the algorithm only accounted for natural hydrological conditions; cultural and economic perspectives were not considered.

Anthropogenic water withdrawal module
The anthropogenic water withdrawal module withdraws the amount of consumptive water use for domestic, industrial, and agricultural purposes from river channels in that order at each simulation grid cell. This module plays an important role in coupling water fluxes among the land surface hydrology, river routing, crop growth, and environmental flow requirement modules.
Domestic and industrial water use was not estimated by the integrated model. Instead, this information was obtained from the AQUASTAT database (Food and Agriculture Organization, 2007). The AQUASTAT database provides statistics-based national water withdrawal data for domestic, industrial, and agricultural sectors. These data were converted to gridded data by weighting the population distribution and national boundary information provided by the Center for International Earth Science Information Network (CIESIN) of Columbia University and Centro Internacional de Agricultura Tropical (CIAT;2005). The values were then converted to the consumptive amount, which is the evaporated portion of the total withdrawal. We used 0.10 for domestic water and 0.15 for industrial water, from the study of Shiklomanov (2000). Seasonal variation was not taken into account for these water uses. For future simulations, projections of other studies such as that by Shen et al. (2008) will be used.
When streamflow was less than the total water demand, streamflow except for the share of environmental flow, was withdrawn. Withdrawn irrigation water was added to the soil moisture in irrigated areas, and domestic and industrial waters were removed from the system. This latter process was an exception to the closure of the energy and water balances; however, the sum of consumptive domestic and industrial water was 132.4 km 3 yr −1 in 1995−1 in (Shiklomanov, 2000. This amount was two orders of magnitude less than global runoff and evaporation (40 000 km 3 yr −1 and 71 000 km 3 yr −1 , respectively; Baumgartner and Reichel, 1975) and was thus judged to be negligible.

Validation methodology
To confirm that the global hydrological cycle was properly simulated using our meteorological forcing data sets and the land surface hydrology module, simulated runoff and streamflow were validated by comparison with observations and earlier studies. First, the land surface hydrology module and the river routing module were coupled, and a global natural hydrological simulation was conducted using two meteorological forcing data sets. Hereafter the runoff/streamflow product obtained using F-GSWP2-B1 is called R-GSWP2-B1 (R stands for runoff; the forcing data were F-GSWP2-B1) and that using F-GSWP2-B0 is called R-GSWP2-B0.
First, observed streamflow data were obtained from the Global Runoff Data Centre (GRDC). From the 3045 gauging stations available, the 37 stations used here have catchment areas totaling >200 000 km 2 , monthly streamflow records for more than 60 months between January 1986 and December 1995, and are the farthest downstream gauging stations in their respective basins (Fig. 2). Of the 37 river gauging stations, five river basins in arid zones (Niger River, Darling River, Orange River, Colorado River, and Cooper Creek) were excluded because most previous studies significantly overestimated observation records from these basins. Thus, we used 32 stations and their catchment areas. To compare simulated runoff with that in previous studies, four published runoff data sets, namely those reported by Baumgartner and Reichel (1975 (Table 3). Of these four data sets, R-BR75, R-D03, and R-F02 are regarded as observation-based runoff products. R-BR75 is a compilation of the statistical record. R-F02 and R-D03 are simulation results in which the simulated runoff was corrected at every interstation basin using an adjustment multiplier defined as the ratio of observed interstation runoff to simulated runoff. As far as we know, the global runoff data of R-BR75, R-F02, and R-D03 are   (Baumgartner and Reichel, 1975;Fekete, 2002;Döll et al., 2003). If the bar height lies within the shaded region, then runoff can be considered plausible generally regarded as the best available data and have been cited extensively by earlier studies. Our focus here was to examine how closely our results fit these previous data. Three gridded global runoff data sets of earlier studies were also routed using the river routing module. R-D03, R-F02, and R-N01 were re-gridded so that they matched both the land/sea distribution and the spatial resolution of F-GSWP2-B1. For R-N01 data, which have 2.0 • ×2.0 • resolution, identical runoff was allocated to four 1.0 • ×1.0 • grids. For R-F02 and R-D03 data, which have 0.5 • ×0.5 • resolution, the runoff of four grids was aggregated into one grid. The Antarctic, Greenland, and lake grids (e.g., Great Lakes in USA and Canada) were excluded from analysis. Finally, simulated streamflow was obtained at 32 river gauging stations. Figure 3 shows simulated runoff at the continental scale for three runoff products (R-GSWP2-B0, R-GSWP2-B1, R-N01). R-GSWP2-B1 was within the plausible range (i.e., the runoff range of the observation-based data sets R-F02, R-D03, and R-BR75) of runoff in Asia, North America, South America, Oceania, and the globe (Fig. 3). In Europe and Africa, it exceeded plausible values by 27% and 10%, respectively. In contrast, R-GSWP2-B0 was the largest among the data sets. In this case, simulated runoff greatly exceeded the range of the three earlier projections in Europe and North America.

Continental runoff
Earlier studies showed a linear relationship between precipitation and runoff (Fig. 4). The observation-based data sets, namely R-F02, R-D03, and R-BR75, gave similar results, indicating that they are consistent in precipitation and runoff. This linear relationship emphasizes precipitation as the dominant factor in the production of continental runoff. Except for Europe and Africa, R-GSWP2-B1 precipitation is similar to that of R-F02, R-D03, and R-BR75; consequently, the projected runoff is also similar. Except for a few cases, R-GSWP2-B0 projects into the upper right, which means larger input precipitation was used and larger runoff was produced compared to earlier studies. The hydrological simulation of the land surface hydrology module shows better performance with F-GSWP2-B1 than with F-GSWP2-B0, mainly because its precipitation input is plausible.
The large precipitation in Europe from F-GSWP2-B1 produced large runoff (Fig. 4). This large precipitation was caused by the rain gauge undercatch correction applied to F-GSWP2-B1 precipitation forcing inputs. Even larger precipitation (similar to F-GSWP2-B1) was given; R-N01 is similar to that of R-F02, R-D03, and R-BR75. This is primarily a result of the basin-level hydrological parameter tuning applied to the model used for R-N01.
In contrast with other continents, for Africa, R-GSWP2-B1 has a large precipitation input when a linear relationship is assumed between precipitation and runoff. This result is attributable to factors other than precipitation, mainly low wind speed. The zonal mean runoff distribution in Africa indicates that the runoff at lower latitudes exceeds that of previous studies, although there is no significant difference in zonal mean precipitation among the earlier studies and F-GSWP2-B1 (not shown). The F-GSWP2-B1 wind speed is low at low latitudes, and evaporation is considered to be restricted (Fig. 1). In the land surface hydrology module, evaporation is basically correlated with wind speed (Appendix B).

Runoff in individual basins
The simulated streamflow data sets were validated at 32 major river gauging stations. Using the simulated and observed data, the normalized bias of mean annual streamflow (NBIAS), the difference in the arrival of peak streamflow (PEAK), and the correlation coefficient (CC) of variation in annual streamflow were calculated as follows: where o is the mean annual observed streamflow (calculated from available records between 1986 and 1995), s is the mean annual simulated streamflow (calculated for months in which o was available), m y,sim is the month in which the simulated monthly hydrograph recorded the maximum streamflow, m y,obs is the month of observation, s y is the monthly simulated streamflow, and o y is the monthly observed streamflow. The subscript y indicates the year. NBIAS is calculated to evaluate the simulated water balance in basins, PEAK to evaluate the timing of streamflow peaks in basins, and CC to evaluate the interannual variation in streamflow. Because lakes and reservoirs affect PEAK and CC considerably, five river basins, namely, the Don, Parana, Sao Fransisco, St. Lawrence, and Nelson, were excluded from calculations of PEAK and CC. R-BR75 was excluded because it reports only zonal mean runoff data. First, we examine NBIAS (Fig. 5a). For R-GSWP2-B1, NBIAS is less than ±50%, except for some river basins in Africa and northeastern South America. These basins are located in semi-arid to arid climatic zones, and runoff was significantly overestimated in these basins (>50% and sometimes >100% of the mean annual difference). The runoff of these basins was commonly overestimated in most of the earlier studies. R-N01 reproduced runoff fairly well, especially for river basins in Siberia. R-F02 and R-D03 showed even better agreement, although these two results are not surprising because the data sets were scaled so that simulated streamflow matched longterm mean annual streamflow. There were some basins with errors >20% because the period selected for scaling in these studies may have differed from ours. In the case of R-D03, there might be another reason for this disagreement, namely that a maximum change of 100% was allowed for the correction factor. In contrast, for R-GSWP2-B0, the runoff in a large number of basins in North America, Europe, and western Siberia was significantly overestimated (>50% of observed); NBIAS in eastern Siberia (e.g., the Amur River and Lena River) was an exception because the simulated runoff was well reproduced. Second, we examine PEAK (Fig. 5b). R-GSWP2-B1 reproduced the timing of long-term mean monthly streamflow well. In most of the basins, the error was within ±2 mo yr −1 . The results of the earlier studies for R-N01, R-F02, and R-D03 are not shown because their runoff data are at monthly intervals and thus are too coarse for a discussion of monthly peak flow.
Third, we examine CC (Fig. 5c). Because the simulation period of R-F02 involved 1 year of climatological information, the CCs for R-F02 are not shown. As in Fig. 4, simulated runoff (or streamflow) was strongly correlated with input precipitation. Because precipitation in the earlier studies was based on ground observations, it seems that the annual variation in simulated runoff agrees well with the observations.
We set arbitrary thresholds for NBIAS, PEAK, and CC. The first set of thresholds was ±20% for NBIAS, ±1 mo for PEAK, and 0.8 for CC. The threshold for NBIAS was derived from Fig. 5a; the R-F02 and R-D03 data supported these criteria for most of the basins. The number of basins that met these criteria was counted for each study (Table 4). The number of basins below the threshold of NBIAS clearly differed among the data sets. R-GSWP2-B1 showed fair performance   -B0  16  25  21  7  19  10  R-GSWP2-B1  24  25  22  14  16  13  R-N01  22  -22  14  -15  R-F02  30  --25  --R-D03  27  -22  19  -13  Total validation basins  32  27  27  32  27  27 and was similar to R-N01; however, only 14 of the 32 river basins met the criteria. R-D03 was generated so that longterm mean annual discharges matched, but not all river basins agreed with observations within ±20% error. The number of basins that met the criteria of PEAK and CC were 16 and 13 of the 32 river basins, respectively. The performance of CC was quite similar among R-GSWP2-B1, R-N01, and R-D03. Our results indicate that it is still a challenge for global hydrology models to simulate annual river discharge year by year. We changed the set of thresholds to ±50% for NBIAS, ±2 mo for PEAK, and 0.6 for CC (Table 4). In this case, NBIAS, PEAK, and CC of R-GSWP2-B1 agreed with the criteria for 24, 25, and 22 of 32 river basins, respectively. Because 70-80% of the validation basins fell within the criteria, we can state that these criteria indicate the simulation performance of the R-GSWP2-B1. In other words, water resource assessments should take the limited ability of the land surface hydrological module and the river routing model into account. The error is not negligible; however, considering that site-specific parameter tuning was not used in the series of simulations, the results indicate the validity of the model and of the input meteorological forcing. The goal of the integrated model was to assess the subannual distribution of water availability and water use. Figure 6 shows the normalized monthly streamflow of R-GSWP2-B1 and observations at 32 validation basins from 1986 to 1995. The monthly streamflow was normalized so that the mean annual streamflow from 1986 to 1995 equaled one. It is clear that significant seasonality occurs in many basins; these exhibited more than three times the mean annual streamflow for only a few months per year, and in the remaining months, streamflow was far less than one. The results indicate that we can move forward to assess the seasonal variability in global water resources.

Summary
To assess global water resources taking into account subannual variability, we developed an integrated model that comprised six modules. We revisited the original GSWP2 meteorological forcing input (F-GSWP2-B0) and developed an improved meteorological forcing data set (F-GSWP2-B1). We then presented the six modules: land surface hydrology, river routing, crop growth, reservoir operation, environmental flow requirement estimation, and anthropogenic water withdrawal. Finally, simulated runoff and streamflow were validated by comparison with observations and earlier works. The performance was similar to the best available precedent studies with closure of energy and water. This result indicates the validity of the model and input meteorological forcing. In our companion paper, we present the model application. Using the daily simulation outputs, we conduct global water resource assessments, focusing on subannual variation in water availability and water use.

Appendix A F-GSWP2-B0 meteorological forcing input
The F-GSWP2-B0 meteorological forcing input is a hybrid product of NCEP-DOE reanalysis (Kanamitsu et al., 2002) and various observation-based monthly gridded meteorological data. The air temperature input is a hybrid product of NCEP-DOE reanalysis (Kanamitsu et al., 2002) and observation-based monthly temperature data of the Climate Research Unit at University of East Anglia (CRU; New et al., 2000). First, the air temperature from the original NCEP-DOE reanalysis was re-gridded from the native 1.9 • ×1.9 • resolution to the ISLSCP2 required 1 • ×1 • resolution and processed from hourly data to three-hourly data (Zhao and Dirmeyer, 2003). The air temperature of the CRU was scaled to adjust for the altitude difference between the CRU grid and the ISLSCP2 mean altitude. The NCEP-DOE reanalyses were linearly scaled so that the monthly mean values were identical to the CRU values. The air temperature data were linearly scaled again so that the diurnal temperature range for each month was identical to that of the CRU. In this way, the air temperature of NCEP-DOE was linearly scaled so that the monthly maximum, minimum, and mean air temperatures were identical to those of the CRU. The daily and three-hourly variations were not corrected; they were determined by the NCEP-DOE reanalysis.
For specific humidity and air pressure, the original NCEP-DOE data were corrected so that they were consistent with the altitude of ISLSCP2 and air temperature. For wind speed, NCEP-DOE data were used without any correction. For longwave and shortwave downward radiation, threehourly Surface Radiation Budget (SRB) data produced at the NASA/Langley Research Center were used (Stackhouse et al., 2000).
The precipitation forcing input is also a hybrid product of the NCEP-DOE reanalysis, the global observational data set GPCC (Rudolf et al., 1994), and the GPCP (Huffman et al., 1997). GPCC data were used for grids where rain gauges were densely located, and GPCP data were used for grids where they were sparsely located. For precipitation data, an algorithm for rain gauge correction for wind-caused undercatch was applied to the rainfall and snowfall input data set (Motoya et al., 2002). In Motoya et al.'s (2002) algorithm, the catchment ratio of snowfall (CR snow ) and that of rainfall (CR rain ) are expressed as follows: where U is wind speed at a height of 2 m, and a and b are parameters from Sevruk (1982) and Sevruk and Hamon (1984) that depend on the rain gauge type. The corrected rainfall (Rainf) and snowfall (Snowf) are expressed as Snowf=Snowf org /CR snow ×100 Rainf=Rainf org /CR rain ×100.
where Rainf org and Snowf org are the original rainfall and snowfall, respectively. These equations indicate that stronger wind requires a stronger correction and that snowfall requires a stronger correction than rainfall. For example, if the wind speed is 5 m s −1 and the rain gauge type is unknown, CR snow is 48.7% and CR rain is 87.2%. The revised meteorological forcing F-GSWP2-B1 was prepared using the same methodology, but different data sets (see Table 1).

Appendix B
The land surface hydrology module

B1 Albedo
The albedo scheme was identical to that of Robock et al. (1995). The snow-free albedo (α soil ) was taken from the GSWP2 standard monthly land use data set and included plant phenological aspects. The snow surface albedo (α snow ) varied according to the surface temperature (T s ) as follows (Robock et al., 1995): where α max is the maximum snow albedo, fixed at 0.60; α min is the minimum snow albedo, fixed at 0.45; T crit is the critical temperature (263.15 K); and T melt is the melting point of ice (273.15 K). The surface albedo was expressed as where SWE is the snow water equivalent [kg m −2 ].

B2 Sensible heat and latent heat
Potential evaporation E P [kg m −2 s −1 ] was expressed as where ρ is the density of air [kg m −3 ], C D is the bulk transfer coefficient (0.003), U is the wind speed [m s −1 ], q SAT (T S ) is the saturated specific humidity at surface temperature [kg kg −1 ], and q a is the specific humidity [kg kg −1 ]. Evaporation from a surface (E) was expressed as where W is the soil water content [kg m −2 ] and W f is the soil water content at field capacity, which was fixed at 150 [kg m −2 ]. The sensible heat flux (H ) is expressed as where C * p is the specific heat of air [1005 J kg −1 K −1 ] and T a is the air temperature [K].

B3 Energy balance
The energy balance was expressed as where SW ↓ is the downward shortwave radiation, LW ↓ is the downward longwave radiation, σ is the Stefan-Boltzmann constant, ι is the latent heat [2.45×10 6 J kg −1 ], and G is the soil heat flux. The original bucket model (Manabe, 1969) does not have soil heat capacity or soil heat flux because it was not designed to simulate diurnal cycles; however, the meteorological forcing input of GSWP2 is threehourly. We added the force restore method (Bhumralkar, 1975;Deardorff, 1978) to simulate the surface temperature: where C s is the surface heat capacity, C d is the deep soil heat capacity and ω is the angular velocity (ω=2π 24·60·60[s −1 ]).

B4 Snow and soil water balance
The snow balance was expressed as where Q s is the surface runoff and Q sb is the subsurface runoff [kg m −2 s −1 ].

B5 Runoff
Surface runoff (Q s ) was generated if the soil water content exceeded the capacity of soil water (i.e., field capacity): Subsurface runoff (Q sb ), which was not in the original bucket model (Manabe, 1969;Robock et al., 1995), was incorporated to the model as where Q sb is the subsurface runoff [kg m −2 s −1 ] and τ is a time constant [s]. This equation is similar to the percolation rate of the LPJ model (Gerten et al., 2004). The γ was set at 2, and τ at 100 days=86400×100 s; both are globally constant.

B6 Mosaic
The module has two separate soil moisture regimes within a grid: one for irrigated areas and one for non-irrigated areas. Identical meteorological forcing input is given to both (a) Simulated results of the land surface hydrology module forced using F-GSWP2-B1. R, mean annual runoff of the basin; P , precipitation; R net , net radiation; solid line, energy and water constraint; curve, Budyko's semi-empirical curve of the relationship between (P − R) P and R net ιP . The colors indicate the Köppen climate classifications: red: tropical rainforest (Af); yellow: tropical monsoon (Am) and savanna (Aw); green: temperate (C), hot summer continental (Da), and warm summer continental (Db); cyan: continental subarctic (Dc), continental subarctic with extreme severe winters (Dd), and polar (E). (b) Same as in (a), but runoff (R) was substituted by GRDC observations. (c) Same as in (a), but the parameters of the land surface hydrology module were modified for four climatic zones.
of the land use types, but surface fluxes and state variables are calculated independently. This scheme is not used in natural hydrological cycle simulations, but is used in naturalanthropogenic coupled simulations when irrigation is taken into account. The soil moisture in irrigated croplands is distinguished from that in other areas.

Parameter setting of the land surface hydrology module
First, the land surface hydrology module (with globally uniform parameters) was driven using F-GSWP2-B1, and a global gridded runoff product and net radiation were obtained. The simulation period was 10 years (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995). The runoff product was routed using the river routing module. Second, for the 37 basins mentioned in Sect. 4.1, Budyko's aridity index (Budyko, 1974) and the evaporation to precipitation ratio were calculated and plotted in a Budyko diagram (Fig. C1a). Budyko's aridity index is expressed as where R net is net radiation, ι is the latent heat, and P is precipitation. The evaporation to precipitation ratio is expressed as where E is evaporation, R is runoff, and P =E+R is assumed. To investigate the relationship between Budyko's diagram and the climatic zone, all land grid cells were initially classified using Köppen's climate classification (McKnight and Hess, 2000), using the monthly temperature and precipitation data of F-GSWP2-B1. Köppen's climate classification was then integrated into four climatic groups: tropical rainforest (Af); tropical monsoon, savanna, and dry climates (Am, Aw, B); temperate and warmer continental climates (C, Da, Db); and cooler continental and polar climates (Dc, Dd, E). This grouping of Köppen's climate classification is similar to that proposed by Nijssen et al. (2001). Finally, a climatic group was assigned to each validation basin according to the majority of grid cells within the basin. The plots corresponded well to Budyko's semi-empirical curve (solid line), and there was no clear relationship regarding to which climate classification a basin belonged (Fig. C1a). This result confirms that the bucket model inherently reproduces the empirical energy and water balance relationship of the Budyko diagram. In some basins with a low aridity index, the plots exceeded the energy-constraint line. This resulted from the negative sensible heat flux in northern high latitudes and mountainous areas (Milly and Dunne, 2002). We assumed that the precipitation and simulated net radiation of GSWP2-B1 were correct and that only simulated runoff was biased. We then re-calculated the evaporation to precipitation ratio ((P −R) /P ) using the observed runoff of the GRDC and plotted the result (Fig. C1b). The distribution of the plot was clearly different from Budyko's curve, and there were relationships to the Köppen climatic classification. In tropical monsoon, savanna, and dry climates, the plots reached higher than Budyko's curve; in contrast, in polar and cooler continental climates, the plots were generally lower. The plots of temperate and warmer continental climates were distributed above Budyko's curve. The plots of basins in tropical rainforest climates were near the curve. These results indicate that the land surface processes in dryer climates and polar and cooler continental climates are not properly reproduced in the land surface hydrology module.
To correct the bias in runoff, the parameters of subsurface flow τ and γ in Eq. (B12) were modified for the four climatic groups. The parameter τ is a time constant that sets the daily maximum subsurface runoff. The parameter γ is a shape parameter that sets the relationship between subsurface flow and soil moisture. The global default values used to draw Fig. C1a were 100 days for τ and 2.0 for γ . A series of simulations was conducted, shifting parameter τ from 50 to 300 at 50-day intervals and parameter γ from 0.5 to 3 at an interval of 0.5. Of 36 combinations, the following combinations of θ and γ that minimized the sum of the root mean square errors of monthly streamflow in each climatic group were determined: (100, 2.0) for tropical forest; (300, 2.0) for tropical monsoon, savanna, and dry climates; (200, 2.0) for temperate and continental [warmer] climates; and (50.0, 1.0) for continental [cooler] and polar climates. Figure C1c shows the Budyko diagram with the modified parameter sets. The plots for basins in tropical monsoon, savanna and dry climates, and temperate and warmer continental climates shift upward, indicating a decrease in runoff. In contrast, the plots in polar and cooler continental climates shift downward, indicating an increase in runoff. Although the distribution in plots of observed runoff (Fig. C1b) still shows a large difference from that of simulated runoff (Fig. C1c), these modified parameter sets substantially improved the runoff simulation (shown in Sect. 4).
This parameter modification was based on findings that the bias in runoff is related to the climatic classifications proposed by Budyko and Köppen. These findings agree with those of Nijssen et al. (2001). Other than climate, soil and vegetation type might be tested to classify the hydrological parameters. Milly and Shmakin (2002) developed the LaD land surface model, which is based on the original bucket model (Manabe, 1969). In the LaD model, the parameters were set by soil and vegetation type, not by climatic classification. Their categorization and ours are similar in that the natural vegetation type is largely dependent on climate. plant transpiration in the second half of the growing season [kg ha −1 ], and SWP is the accumulated potential evapotranspiration in the second half of the growing season [kg ha −1 ].
If we use the formulation and parameters of the SWIM model globally, the cropping period of countries at low latitudes become unrealistically short (e.g., <100 days for the cropping period of grains). The SWIM model provides only one parameter set for each crop type, e.g., base temperature, optimal temperature, and potential heat units required for maturity. Heat units accumulate rapidly when the difference between the daily mean temperature and the base temperature is large, and the threshold of potential heat units required for maturity is attained quickly. Therefore, we set an upper limit for daily heat units. Doorenbos and Pruitt (1977) showed cropping periods for various crops planted in various places in the world; except for some vegetables, crops need at least 120 days to mature. The potential heat unit threshold for maturity is in many cases 1500 K in the SWIM model, so we set the daily maximum heat unit threshold to 12.5 K and excluded any excess heat units. In this case, at least 120 days are needed to reach maturity in all crops. This corresponds to altering the T B in Eq. (C1) to fit the local temperature or to planting different species that have a larger T B than that of SWIM. There is no upper limit of the cropping period, but it must be less than 365 days.

Reservoir operation module
This module sets operating rules for individual reservoirs. There are two types of operating rules: irrigation and nonirrigation. If the reservoir's primary purpose is not irrigation water supply, the "nonirrigation operating rule" is set as follows. This operation tries to reduce both the interannual and seasonal variation in streamflow; if conditions allow, the reservoir releases the mean annual inflow throughout the year. First, for every reservoir-georeferenced grid, each month is categorized as either a recharge month in which mean monthly inflow exceeds mean annual inflow or a release month. Second, we define the "operational year," which begins in the first month of a release period (longest continuous release months in a year). We assumed that the annual total release for an operational year is fixed at the beginning of the operational year. Thus, the annual total release for the operational year (R [kg yr −1 ]) is provisionally set as follows: where I mean is the mean annual inflow [kg yr −1 ], S ini is the storage at the beginning of the operational year [kg], and C is the storage capacity of the reservoir [kg]. The coefficient of 0.85 was set semi-empirically . If storage was smaller than 0.85C, the release for the next 12 months was smaller than the mean annual inflow to recover storage. In this way, the interannual variation in inflow is buffered by the reservoir. Once the annual release is fixed, the daily release from reservoirs (r [kg s −1 ]) is expressed as where i mean [kg s −1 ] is the mean annual inflow, i [kg s −1 ] is the daily inflow, and c is the ratio of the storage capacity (C [kg]) to the volume of mean annual inflow (I mean [kg]). The first equation is for reservoirs that have large storage capacity compared to annual inflow: release is independent of monthly inflow. The second equation is for reservoirs that have small storage capacity compared to annual inflow. To avoid overflow and storage depletion during the year, release is influenced by daily inflow. The squared exponent and criterion of 0.5 are set empirically. When c equals zero, reservoir operation is identical to run-of-the-river flow. If storage exceeds storage capacity even when allocated water has been released, the excess volume of water is also released. If the reservoir's primary purpose is irrigation water supply, the "irrigation operation rule" is applied. This operation tries to reduce the interannual variation in streamflow. Daily release is not constant, but is controlled to be proportional to the daily water demand in the lower reaches. The annual total release for an operational year is identical to that for the case of nonirrigation reservoir operation. The water demand of irrigated areas is calculated for 10 grids downstream of a reservoir. The delay from delivery is not taken into account.

Thresholds in the environmental flow requirement module
The 10% and 30% thresholds of base requirements are based on the work of Tennant (1976). He evaluated environmental flow as a proportion of mean annual streamflow and argued that 10% was the minimum requirement, 30% was good, and 60% was excellent. A number of studies have further supported the threshold of 30%. King et al. (2000) introduced some case studies of rivers in South Africa and argued that the annual total environmental flow was 37.3-45.7% of the annual streamflow for the Marite River and 30.2% for the Mhlatuze River. The water plan of the state of California, USA, projected the future share of water use by sector and suggested that the environmental flow requirements will be 32-46% of the total water demand in 2020 (Department of Water Resources State of California, 1998). In Japan, the environmental flow requirement was determined from several thresholds: