Integrating groundwater irrigation into hydrological simulation of India: Case of improving model representation of anthropogenic water use impact using GRACE

Study region: India Study focus: India boasts the largest irrigated agricultural system in the world relying on groundwater. To address the strong linkages between the natural groundwater and the anthro- pogenic irrigated system requires innovative hydrological modeling geared at informing national policies on groundwater management and future development of irrigated agriculture. For this, we developed a predictive, integrated hydrological and groundwater use model and evaluated the model using total water storage (TWS) data from the Gravity Recovery and Climate Experiment (GRACE). The utility of the model was demonstrated in a case study in which the model was applied to project the groundwater balance in northwest India under four RCP (Representative Concentration Pathways) scenarios. New hydrological insights for the region: The model shows good identi ﬁ ability to GRACE data in northwest India and with incorporated groundwater irrigation simulation module the model can adequately replicate the declining trend in TWS over this region. It is concluded that by assuming a unchanged pattern of agricultural water use climate change is likely to help reduce the magnitude of the groundwater de ﬁ cit, but the bene ﬁ cial e ﬀ ect is insu ﬃ cient to halt the trend of groundwater depletion. This result provides new evidence for the importance of groundwater conservation through changes in cropping patterns and improved groundwater governance.


Study region: India
Study focus: India boasts the largest irrigated agricultural system in the world relying on groundwater. To address the strong linkages between the natural groundwater and the anthropogenic irrigated system requires innovative hydrological modeling geared at informing national policies on groundwater management and future development of irrigated agriculture. For this, we developed a predictive, integrated hydrological and groundwater use model and evaluated the model using total water storage (TWS) data from the Gravity Recovery and Climate Experiment (GRACE). The utility of the model was demonstrated in a case study in which the model was applied to project the groundwater balance in northwest India under four RCP (Representative Concentration Pathways) scenarios. New hydrological insights for the region: The model shows good identifiability to GRACE data in northwest India and with incorporated groundwater irrigation simulation module the model can adequately replicate the declining trend in TWS over this region. It is concluded that by assuming a unchanged pattern of agricultural water use climate change is likely to help reduce the magnitude of the groundwater deficit, but the beneficial effect is insufficient to halt the trend of groundwater depletion. This result provides new evidence for the importance of groundwater conservation through changes in cropping patterns and improved groundwater governance.

Introduction
Hydrological research traditionally focuses on physical processes, but the importance of investigating linkages between the anthroposphere and the hydrosphere has been well recognized (Wagener et al., 2010). In this paper, we developed a model that integrates groundwater use for irrigation into a national-scale hydrological simulation system of India, where the anthropogenic influence of groundwater irrigation has dramatically changed the country's water resource base. India boasts one of the largest irrigated agricultural systems in the world, with more than 60 million hectares equipped for irrigation (FAO AQUASTAT). Furthermore, more than 60 % of the irrigated area is fed by groundwater (FAO AQUASTAT). The heavy use of groundwater for irrigation has created extensive concerns about overdraft and depletion and has made India a global hotspot for water resources studies (Kumar et al., 2007;Foster et al., 2008). Any hydrological model of India that is developed or used for decision-support will therefore need to be able to quantify the impact of groundwater irrigation. We developed the model using the Soil and Water Assessment Tool (SWAT) (Arnold et al., 1998) and publicly available data, in particular total water storage (TWS) variation data from the Gravity Recovery and Climate Experiment (GRACE) (Tapley et al., 2004) against which model performance was evaluated. GRACE TWS data provide observations of terrestrial total water storage anomalies that integrate variations of all vertical components of water storages and have been applied for groundwater resources assessments in India at a national scale by Rodell et al. (2009) ;Tiwari et al. (2009) ;Chen et al. (2014); Long et al. (2016) and Asoka et al. (2017) etc. They were also used by Barik et al. (2017) to examine the food-water-energy nexus in India over recent time periods. What distinguishes our study from these studies is that past applications of GRACE TWS data in India focused on using trends derived from the GRACE TWS data to estimate historical groundwater depletion rates while our study developed a predictive model with enhanced representation of linkages between groundwater irrigation and groundwater storage. Modules representing anthropogenic water use have been developed and incorporated into a number of large-scale hydrologic simulation models (Cai and Rosegrant, 2002;Alcamo et al., 2003;Hanasaki et al., 2008;Wada et al., 2012;Döll et al., 2014;Pokhrel et al., 2015;Wada et al., 2016), and GRACE data have been used in the validation of some of these models (f. ex. Döll 2012 and. However, these models have global coverage and did not include detailed results, discussion and information on model uncertainty for individual countries. In this study, we also tested the feasibility of applying process-based crop simulation techniques in the development of a largescale, integrated hydrological-water use model in a country with varying climate and intensive agricultural water use. Simulation of irrigation water demand and irrigation water use in global models is primarily based on the FAO-56 crop coefficient approach (Allen et al., 1998;Döll and Siebert, 2002;Wada et al., 2014), which links consumptive water use of crops to climate variables through a series of coefficients. SWAT, originally designed as a meso-scale river basin model, incorporates algorithms from the EPIC (Erosion-Productivity Impact Calculator) crop simulation model (Williams et al., 1989) and therefore is capable of simulating plant water uptake and irrigation water application processes in much greater detail. The size of study river basins in SWAT applications typically range from tens to several thousand square kilometers, but in recent years the model has also been applied at national and continental levels (Schuol et al., 2008;Xie et al., 2012).
Moreover, as described in section 3, with the incorporation of the groundwater irrigation simulation module, the model can adequately replicate the declining trend in TWS in northwest India. In section 4, we applied the model developed in this study to provide a "partial" analysis of the impact of climate change on groundwater balances in northwest India as a demonstration of the utility of the model. The prospects of extending the model to provide fuller analysis to inform policy making for groundwater management and groundwater-fed irrigated agriculture under changing climate and socioeconomic conditions are also discussed.

SWAT and watershed delineation
SWAT is a semi-distributed river basin model. Using digital elevation data from the HydroSHEDS (Hydrological data and maps based on SHuttle Elevation Derivatives at multiple Scales) (Lehner et al., 2008) database, the study region was partitioned into 677 subbasins (Fig. 1). The area modeled is 3.9 million km 2 and, in addition to India, covers Nepal, Bangladesh, Bhutan as well as parts of China, Myanmar and Pakistan.

HRU definition
Within a subbasin, hydrologic response units (HRUs) can be defined. An HRU in a subbasin is a land segment consisting of areas with the same soil, land cover/use and land management practices. The definition of HRUs helps reflect the heterogeneity of soil, land cover/use and land management practices within a subbasin (Winchell et al., 2007).
For this study, we defined HRUs for the main crops under groundwater irrigation. To this end, we harmonized land use/cover data from the GLC (Global Land Cover) 2000 database (Bartholome and Belward, 2005) and downscaled crop harvested area data by crop and by production system from the SPAM (Spatial Production Allocation Model) database (You et al., 2006). GLC is a comprehensive global land cover and land use database containing land cover/land use data classified under the FAO Land Cover Classification System that covers the complete range of land cover and land use from natural vegetation to anthropogenic land uses. In GLC2000, pixels with dominant land use for agriculture are classified as irrigated agriculture, intensive irrigated agriculture and rainfed agriculture. More precise and detailed estimates of irrigated agricultural area by crop are obtained from the SPAM database. In the production of the SPAM database, national or sub-national statistics are downscaled to a 5 arc-minute (approximately 10 km × 10 km) grid. The SPAM 2000, version used in this study, contains estimated harvested areas for 20 major crops around 2000 on a 5 arcminute grid in four production systems: irrigated, rainfed-high input, rainfed-low input and rainfed-subsistence. The irrigated harvested area by crop in the study region reported by SPAM is summarized in Table 1. We defined groundwater irrigated HRUs for the five key groundwater irrigated crops in the country: rice, wheat, cotton, sugarcane and maize (Fig. 2). They account for 84 % of total irrigated harvested area in the study region. The harvested area of each of the five crops under groundwater irrigation in each subbasin was further estimated using data on area equipped with groundwater irrigation from a global groundwater irrigation inventory developed by Siebert et al. (2010) (Fig. 3). We calculated: where A GW C , is groundwater irrigated harvested area for crop c in a subbasin, A Total C , is total irrigated harvested area of crop c in the  is the share of irrigated area equipped with groundwater irrigation in the subbasin (0-1) derived from the groundwater irrigation inventory data developed by Siebert et al. (2010). Moreover, among the five crops, sugarcane is a perennial crop, maize and cotton are mainly Kharif (summer) season crops while the cultivation of wheat mainly occurs in the Rabi (winter) season. For these crops, a single growing season is assumed in each year in the model with planting dates and length of the growing  season estimated by Sacks et al. (2010) and the physical irrigated areas in each growing season are equal to the annual harvested irrigated areas. As for rice, India features up to three rice-growing seasons: Kharif, Rabi and pre-Kharif. The physical areas of irrigated rice in each season were estimated using data on rice cropping patterns compiled by Frolking et al. (2006). It is worth noting that in this study we limited our effort to simulating groundwater irrigation and viewed this as a step toward developing a more full-fledged national hydrological model of India. Simulation of irrigation from surface water is linked to modeling storage variations and water conveyance and allocation processes in man-made surface water infrastructure. However, there is a lack of data to support such simulations. Although global hydrological modeling studies simulate assumed water infrastructure operation rules for modeling human water uses including surface water irrigation (f. ex. Hanasaki et al., 2006;Döll et al., 2012), such approaches are less justifiable at the national level. The unaccounted surface water irrigation has implications for modeling groundwater balances (see the discussion for Eq. 6), and this constitutes a major limitation of our study. Omitting the water storage variation from man-made water infrastructure, which is a component of the TWS variation detected in GRACE, also implies greater uncertainty in GRACE data-based model evaluation. In fact, there is uncertainty originating from unvalidated model performance in simulating whole surface water system. However, national river discharge data are difficult to obtain for India. In this study, we left the model uncalibrated/unvalidated in terms of its performance in modeling surface water hydrology.
The soil data used in this study were obtained from the WaterBase (George and Leon, 2008). This data set is developed based on the FAO/UNESCO Soil Map of the World and contains derived soil properties required for SWAT modeling. A dominant soil type was identified in each subbasin and used for HRU definition.

Irrigation and groundwater storage simulation
SWAT operates at a daily time step. The simulation algorithms in the model are described in detail in (Neitsch et al., 2011). Only those components that are relevant to simulating irrigation and groundwater balance are discussed here. Modifications were made to address limitations in the SWAT model in simulating irrigation for paddy rice and groundwater storage variation.

Irrigation simulation in SWAT
The crop simulation algorithm in SWAT originates from the EPIC (Erosion-Productivity Impact Calculator) model (Williams et al., 1989). On each day during a growing season, the model keeps track of root, leaf and biomass development of crops and calculates crop water demand as well as the amount of water the plant can actually uptake from the soil. SWAT provides functions to simulate irrigation events. For this study, the auto-irrigation scheme was used to simulate irrigation activities related to the production of cotton, maize, sugarcane and wheat. Under this auto-irrigation scheme, irrigation can be triggered automatically based on water stress experienced by crops or irrigation water is applied whenever the water stress of the crop exceeds a specified threshold. Crop water stress is calculated as: where wstr is water stress on a given day with a value ranging from 0 to 1. E t is the potential transpiration or water demand of a crop on a given day (mmH 2 O), and E a is the amount of water actually taken up or transpired by the plant (mmH 2 O). Once auto-irrigation is invoked, soil moisture will be restored to soil field capacity.
While the auto-irrigation algorithm in SWAT is suitable to simulate irrigation for upland crops, the model has limitations in simulating irrigated rice production. Pothole module in SWAT is often used as an option in modeling rice paddy, which is designed to simulate hydrology in closed depressions with ponded water. However, the module cannot adequately simulate infiltration in saturated soils underlying paddy fields and tends to underestimate percolation losses (Sakaguchi et al., 2014). Furthermore, the autoirrigation algorithm with crop water stress as a trigger is no longer applicable as the key decision variable for irrigated rice is ponded water depth. In this study, we adopted the modifications proposed by Kang et al. (2006); Xie and Cui (2011) and Sakaguchi et al. (2014) to improve the capacity of the SWAT pothole module in simulating water balances and irrigation for irrigated rice production. We assumed a constant daily percolation rate for the water balance calculation when the field is flooded (Kang et al., 2006, andSakaguchi et al., 2014). Second, we implemented the irrigation simulation algorithm for paddy proposed by Xie and Cui (2011). In this algorithm, the water depth in paddy varies between a predefined minimum fitting depth and a maximum ponding depth. Irrigation occurs when the water depth drops below the minimum fitting depth, and the objective of the irrigation operation is to fill the paddy field with water to the maximum fitting depth (minimum fitting depth < maximum fitting depth < maximum ponding depth). The values of the three critical water depths vary by the growth stage of rice (Table 1 in Xie and Cui, 2011).
Using the simulation algorithms described above, the pumping rate of groundwater for irrigation was first estimated as soil moisture deficit with respect to soil field capacity (for crops other than paddy rice) or ponded water depth deficit with respect to maximum fitting depth of paddy fields and then adjusted for water losses in irrigation application, which are unaccounted for in the simulation (e.g., water conveyance losses), by an irrigation efficiency coefficient (Santhi et al., 2005).
where Q irr is the estimated groundwater pumping rate for irrigation (mmH 2 O), ′ Q irr is the irrigation water demand estimated without consideration of unaccounted-for-water losses in irrigation (mmH 2 O), and γ is the irrigation efficiency coefficient. The volume of pumped groundwater for irrigation was also constrained by the water storage in aquifers. Irrigation stops when aquifer water is depleted.

Groundwater balance simulation in SWAT
SWAT uses the Soil Conservation Service curve number method (Soil Conservation Service (SCS, 1972) to calculate surface runoff, and the model simulates infiltration, evaporation, plant water uptake and interflow generation processes in the soil profile. Groundwater aquifers are recharged by water infiltrating through the soil.
In SWAT, for each HRU groundwater storage is partitioned into two components: one represents shallow aquifer and the other deep aquifer. Either component can be assigned as the water source for irrigation from which groundwater is abstracted. Moreover, both components can contribute to baseflow. Groundwater baseflow from each component is calculated as where Q gw i , is the groundwater flow on a given day i (mm H 2 O), − Q gw i , 1 is the groundwater flow on the preceding day i-1(mmH 2 O), α gw is the baseflow recession coefficient (day −1 ), t Δ is the time step (1 day), and Q rchrg i , is the amount of recharge entering aquifer on day i (mmH 2 O).
Water in shallow aquifers can also be removed by the "revap" process, which is defined in SWAT to account for water removed from shallow aquifers through the capillary fringe that separates saturated from unseparated zones. "Revap" occurs when water storage in shallow aquifers exceeds a certain threshold, and the amount of water removed by the "revap" process on a given day where β revap is the revap coefficient, and E 0 is the potential evapotransipration for the day (mm H 2 O).
In this study, we substituted a new groundwater module for the module in SWAT. The water balance for groundwater storage in the new module is: where S Δ is groundwater storage change on a given day (mmH 2 O), Q rchrg is the recharge (mmH 2 O), Q gw is groundwater flow out of the storage on a given day (mmH 2 O), and Q irr is the amount of water abstracted for irrigation, or pumping rate (mmH 2 O). We chose to omit the "revap" process since the effect of the "revap" process tends to be confounded with the effect of irrigation water withdrawal: a large "revap" coefficient can lead to similar declines of groundwater storage as irrigation water withdrawals, and data to describe the "revap" process is lacking.
Aquifer formations in India generally fall into two groups (Central Ground Water Board (CGWB, 2012): unconsolidated and consolidated/semi-consolidated. Unconsolidated aquifers underlie the Indo-Gangetic and Brahmaputra plains in northern India (MacDonald et al., 2016) ( Fig. 1). They serve as large underground reservoirs to support the intensive irrigated agriculture in the region. By contrast, aquifers in central and southern India are predominantly consolidated/semi-consolidated with limited water storage capacity. In the new module, we distinguish between these two types of aquifers. For subbasins with unconsolidated aquifers as the dominant formation, the groundwater storage is still partitioned into two components although the two components are referred to as "active" and "inactive" instead. In the simulation, only the active storage contributes to groundwater baseflow and when irrigation occurs water is first abstracted from the active storage. As water in the active storage component is depleted, no further baseflow will be discharged from the aquifer to streams and water will start to be abstracted from the inactive component. The initial storage in the active component is established through a spin-up (see section 2.5). Unlimited groundwater storage is assumed in unconsolidated aquifers. A large positive value was thus assigned as initial storage of the inactive component of groundwater storage.
The baseflow from the active storage of the aquifer is calculated as: where k is the groundwater flow recession constant, and Q active is the volume of water in the active component of groundwater storage (mmH 2 O). This equation has been extensively applied in hydrologic modeling to estimate groundwater flow (f.ex. Chapman, 1999;Fenicia et al., 2006;Beck et al., 2013). The motivation of replacing Eq. (4) with Eq. (7) is to better reflect the impact of irrigation water withdrawal on groundwater flow. There is a lack of response in groundwater baseflow calculated using Eq. (4) to changes in groundwater storage caused by irrigation (Zeng and Cai, 2014). For river basins with consolidated/semi-consolidated aquifers, the entire storage was assumed to be active with baseflow being calculated using Eq. (7).
The recharge term Q rchrg in Eq. (6) includes contributions from precipitation and "return flows" from groundwater irrigation in our model. In reality, groundwater can also be replenished by water from surface water irrigation which includes not only percolation from fields irrigated with surface water but also seepage from irrigation canals during water conveyance (MacDonald et al., 2016;Singh, 2011). Irrigation canals are typically built as components of large-scale surface water irrigation projects. The seepage rate of water from irrigation canals is a function of canal length, lining materials, water residence time and wetted perimeter, and effective hydraulic conductivity of the channel alluvium (Santhi et al., 2005). As noted above, considering limited data availability, recharge from surface water irrigation was omitted, and this is a primary source of uncertainty in our study.
The pumping rate, Q irr , was estimated using the irrigation simulation algorithm described in section 2.3.1.

GRACE data and comparison of SWAT-and GRACE-based TWS variation data
The GRACE TWS data used for this study were obtained from CNES-GRGS (Centre National d'Etudes Spatiales-Groupe de Recherches de Géodésie Spatiale), RL3-v3 product, which provides an improved temporal resolution (10-day means) with high spatial resolution. The Stokes coefficients were truncated at degree 50 (i.e. ∼400 km spatial resolution) to address the issue of decreasing signal to noise ratio with increasing spatial resolution. Spherical harmonic data were recombined and projected on a 0.5°l atitude/longitude grid following Wahr et al. (1998).
For model identification, we compared the GRACE-based and SWAT-simulated TWS variations during 306 10-day intervals between July 29, 2002 and June 2, 2011 (GRACE TWS variation data are missing for 18 10-day intervals during the periods from November 26, 2002to February 23, 2003, May 25, 2003to July 3, 2003, January 20, 2004to January 29, 2004, and December 24, 2010to February 1, 2011. The strategy applied to compare GRACE-and SWAT-based TWS variations is discussed in detail in Xie et al. (2012) and is summarized below: Firstly, we computed basin-wide TWS variations for each 10-day interval and disaggregated them on a 0.5 degree grid based on daily SWAT simulation results. Secondly, regional model outputs were "nested" in a global model. Simulated TWS variations from the Noah land surface model (Ek et al., 2003) in NASA's Global Land Data Assimilation System (GLDAS) (Rodell et al., 2004) were used as prior information to fill the regions outside the SWAT simulated area. The global TWS data were then projected on spherical harmonics using SHTOOLS (https://shtools.oca.eu/shtools/). The projection converted the model-based TWS variation fields to GRACE resolution and allowed a comparison of the two data sets (Güntner, 2008).

Climate data and spin-up
SWAT provides three options for estimating potential evapotranspiration (Neitsch et al., 2011): the Priestley-Taylor method (Priestley and Taylor, 1972), the Penman-Monteith method (Monteith, 1965), and the Hargreaves method (Hargreaves and Samani, 1985). The Penman-Monteith method was selected. For this option, climate data for precipitation, temperature, solar radiation, relative humidity and wind speed data are required in the SWAT simulation. The gridded One-Degree Daily (1DD) data from 1997 through 2011 of these variables covering the GRACE TWS data period were obtained through the agroclimatology portal of the NASA Prediction of Worldwide Energy Resource (POWER) Project; original data sources include the Global Precipitation Climate Project (GPCP), the Goddard Earth Observing System model version 4 and version 5 (GEOS-4 & GEOS-5), the Release 3 of the NASA/GEWEX Surface Radiation Budget (GEWEX SRB 3.0) project, and the NASA's Fast Longwave And SHortwave Radiative Fluxes (FLASHFlux) project. During model development, the subbasin-wide estimates of these climatic variables were calculated by averaging griddedbased climate data using fractions of area covered by different climate data grid cells as weights (Schuol and Abbaspour, 2007).
In SWAT modeling, initial conditions of hydrologic state variables, including active groundwater storages, are unknown and are established through spin-up simulations. Given the need to initialize active groundwater storage, a long spin-up period, 1910-2001, was used. Climate data obtained from ISI-MIP, which were downscaled from the NorESM1-M model, were used for simulation during the spin-up period. Groundwater based irrigated agriculture in India started to take off during the 1960s as part of the Green Revolution. Groundwater irrigation activities during 1950-2001 were simulated in the spin-up simulation, in which sizes of groundwater irrigated crop HRUs in each subbasin were adjusted dynamically by a ratio derived from historical data on state-wise groundwater irrigated area reflecting the actual pathways of groundwater-driven irrigation expansion (Bhaduri et al., 2006;Gandhi and Namboodiri, 2009).

Results of model identification
Figs. 4-6 present simulation results using a set of selected/default values of model parameters (water percolation rate in paddy rice = 8 mm H 2 O/day, threshold of upland crop water stress triggering irrigation = 0.7, Irrigation efficiency = 0.7, soil evaporation compensation factor = 0.95, groundwater recession coefficient = 0.05). The map in Fig. 4(a) shows the model quality on a grid basis obtained from simulation that does not consider groundwater irrigation. The agreement between the observed and simulated TWS data series was measured by the Nash-Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970). NSE is a coefficient extensively used in hydrologic modeling for model evaluation. Its value ranges from − ∞ to 1 and equals unity when model output and observations are identical. Fig. 4 partitions India into three regions: the northwest (region I), the northeast and east (region II) and the central and southern region (region III). The northwest region defined in this study covers an area of 0.9 million km 2 and stretches across seven Indian states: Punjab, Haryana, Delhi, Himachal Pradesh, Uttar Pradesh, Rajasthan and Jammu and Kashmir. This region accounts for 40 % of irrigated area of the country (the northeast and east region for 28 % and the central and southern region for 32 %). Moreover, approximately 70 % of the irrigated area in the northwest region is irrigated by groundwater (the shares are 61 % and 63 % in the northeast and east region and the central and southern region, respectively). As evident in the map, the model fit is very good in the northeast and east region of India and good in the central and southern region of the country. However, poor values of NSE can be found in northwest India indicating that the simulated TWS values deviate substantially from GRACE observations. The reasons  behind the displayed model misfit over the region are revealed in Fig. 5(a) and (b) and the aggregated TWS time series plots in Fig. 6(a). Fig. 5(a, b) shows the trends in gridded TWS data from GRACE and the model simulation (in filtered space without simulated groundwater use for irrigation), respectively. Large-scale depletion trends (up to 70 mm H 2 O/yr) in northwest India can be seen in the trend map based on GRACE data (Fig. 5a), whereas such patterns cannot be seen when groundwater irrigation simulation is turned off (Fig. 5b). Fig. 6(a) shows the aggregated GRACE-based and simulated TWS data series over northwest India. The map also shows a downward trend in the GRACE-based TWS data series while such a trend is lacking in the simulated TWS data series.
When simulations include groundwater irrigation activities, model fit improves significantly in the northwest region of India and the model successfully replicates TWS trends in GRACE data, as indicated by the gridded NSE map in Fig. 4(b) and the corresponding gridded TWS trend maps in Fig. 5(c) (at GRACE resolution) and Fig. 5(d) (at full resolution) and the aggregated model-based TWS data series in Fig. 6(b).
Encouraged by the promising results obtained from the initial model evaluation, further analysis was performed with quantification of parametric uncertainty in the simulation. The analysis was conducted using the conceptual framework of generalized likelihood uncertainty estimation (GLUE) (Beven and Binley, 1992). Five model parameters were included in the analysis. These parameters and their value ranges used for the analysis are shown in Table 2. In addition to those parameters defined to characterize groundwater use processes for irrigation and groundwater flow (PRRP-water percolation rate in rice paddy, WSTR-threshold of upland crop water stress triggering irrigation, γ -irrigation efficiency, and k -groundwater recession coefficient), we also took the soil evaporation compensation factor (ESCO) into account. Soil moisture variability is a major source of the simulated TWS variation, and ESCO is a key parameter defined in SWAT for modeling soil moisture variation. The implementation of GLUE involved Monte Carlo sampling of uncertain parameters. Two hundred random samples were drawn uniformly from the specified parameter value space. Plots in Fig. 7 provide a comparison of pair-wise NSE values calculated using aggregated GRACE-based and simulated TWS data series over three regions within India in cases with and without simulated groundwater irrigation. Consistent with observations from Figs. 4-6, groundwater use for irrigation is the main process influencing TWS trends in northwest India in the SWAT simulation; by contrast, NSE values for the other two regions only vary slightly after groundwater irrigation is incorporated into the simulation.
GLUE is designed to accommodate parametric uncertainty in the model simulation, referred to as equifinality, arising from the existence of multiple "behavioral" parameter sets. All "behavioral" parameter sets enable the model to fit the observation data acceptably well. Fig. 8 presents dotted plots showing correspondence between sampled parameter values and the resultant NSE values which indicates model fit over the northwest region of India, obtained from Monte Carlo sampling. Based on these results and following the spirit of the GLUE method, an NSE cutoff value of 0.8 was selected, and 83 out of 200 sampled parameter sets associated with NSE values greater than 0.8 were accepted as behavioral. The strong equifinality revealed in the GLUE analysis leads to large uncertainty in the simulated water balance and TWS. The ensemble of TWS variations over northwest India simulated using the identified behavioral parameter sets are shown in Fig. 6(b). The distributions of estimated groundwater recharge, groundwater baseflow, irrigation water withdrawal and annual average groundwater depletion are displayed in Fig. 9. Modeled groundwater depletion over the region ranges from 14 km 3 H 2 O/yr to 29 km 3 H 2 O/yr, with a mean value of 23.9 km 3 H 2 O/yr and the mode at 24.3 km 3 H 2 O/yr.
Note that the uncertainty CNES-GRGS GRACE TWS trend is -27 ± 2 km 3 from weighted least-square regression that lays within the uncertainty interval derived from the GLUE analysis. We also considered capturing GRACE uncertainty using a wider range of GRACE products (Farinotti et al., 2015). When using the full range of GRACE products, namely official products (CSR, GFZ, JPL, release 5) and alternative products processed by TU GRAZ release 2016 (Mayer-Gürr et al., 2016), HuaZhong University of Science and Technology release 16 (Zhou et al., 2017), Tongji University release 2  and the Astronomical Institute from Univ. Bern release 2, the estimated trend is -29 +/-2.5 km 3 , which still partially overlaps with the GLUE interval. In the GLUE analysis, model fit to GRACE TWS data was measured by the Nash-Sutcliffe efficiency, which incorporates the effect of seasonal variation of TWS. It might be possible to use the information on uncertainty of the GRACE TWS trend to further narrow down the selection of the "behavioral" parameters in future work when more "certain" estimates for uncertainty bounds of GRACE TWS trends become available. As the variability in the different estimates for GRACE TWS trend uncertainty is not fully understood, we did not use such information in the GLUE analysis in this study.

Assessing impacts of climate change on groundwater storage in northwest India
The validity of the developed hydrological-agricultural groundwater use model in simulating groundwater storage variation in northwest India was established through the results of validation analysis presented in the preceding section. As a demonstration of the predictive skill of the model, we extended the simulation to 2050 to simulate future groundwater storage/groundwater irrigation Table 2 Feasible value ranges in GLUE parameters.

Range of feasible values
Water percolation rate in rice paddy (PRRP, mm H 2 O/day) 2 ∼ 12 Threshold of upland crop water stress triggering irrigation (WSTR) 0.5 ∼ 0.9 Irrigation efficiency (γ ) 0.5 ∼ 0.9 Soil evaporation compensation factor (ESCO) 0 ∼ 1 Groundwater recession coefficient (k) 0 ∼ 1 H. Xie, et al. Journal of Hydrology: Regional Studies 29 (2020) 100681 in northwest India under changed climate conditions. The model developed in this study takes data on crop mix in groundwater-fed irrigated agriculture as an input. Groundwater irrigation is expected to increase in the future in line with income and population growth, and because groundwater is considered to be the most productive component of the agricultural production system. At the same time, concerns about groundwater depletion are counteracting factors and may lead to policies and regulations to curb the development of groundwater-based irrigated agriculture. India is also a major player in global agriculture markets. As a result, policy options and interventions, which may result in major changes in food production, might affect global food security and need to be carefully evaluated in the context of global agriculture trade. For these reasons, any future crop pattern changes would be part of a complex decision-making process. A tradeoff analysis that considers food security and groundwater sustainability as objectives and includes future cropping pattern as a decision variable would therefore be needed to inform a thorough discussion about policies for groundwater management and groundwater-fed irrigated agriculture under changing climate and socioeconomic conditions, and more integrated modeling system that couples the hydrological-groundwater use model with food production-agricultural economic/ trade modeling tool needs to be applied to support the tradeoff analysis. The paper focuses on presenting the development of a hydrological model with improved representation of agricultural groundwater use processes. Instead of providing a full tradeoff analysis, in this case study we considered a special case in which the cropping pattern of irrigated agriculture, and irrigation efficiencies, were kept constant. The case study provides insights that could be used to inform policy discussions qualitatively on the direction and urgency of shifting agricultural groundwater use patterns.
The projection analysis used climate change data downscaled by the Inter-Sectoral Impact Model Intercomparison Project (ISI-MIP) (Hempel et al., 2013). ISI-MIP offers a framework for consistently assessing the impacts of climate change across sectors. A large collection of food and agricultural trade models used these data (Nelson et al., 2014). The ISI-MIP climate data are downscaled from outputs of five general circulation models (GCMs) from the Coupled Model Intercomparison Project Phase 5 (CMIP5): GFDL-ESM2M, HadGEM2-ES, IPSL-CM5A-LR, MIROC-ESM−CHEM and NorESM1-M generated under four Representative Concentration Pathways (RCP) scenarios: RCP2.6, RCP4.5, RCP6.0, and RCP8.5 (Moss et al., 2010). The five GCMs were selected by ISIMIP because projections from these GCMs span a wide range of the space of projected global mean temperature change and relative precipitation changes (Warszawski et al., 2014). Studies evaluating the performance of GCMs over India, for example, Saha et al. (2014), revealed H. Xie, et al. Journal of Hydrology: Regional Studies 29 (2020) 100681 that the majority of CMIP5 GCMs failed to simulate the decreasing trend of the Indian monsoon in historical experiments. In a more recent study (Gusain et al., 2020), improvement was found to be made in CMIP6 models, which can capture Indian monsoon characteristics reasonably well. At the time of this study, the full set results from CMIP6 were not available. We therefore choose to only report results from NorESM1-M scenarios. NorESM1-M is a top-rated model among a collection of CMPI5 GCMs based on its performance in simulating present-climate characteristics of Indian summer monsoon (Sharmila et al., 2015). Changes in annual precipitation and average daily temperature between base period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) and the last decade of the projection period (2041-2050) indicative of trends of climate change over northwest India projected by the NorESM1-M are shown in Table 3. The extent to which climate change will influence groundwater balances also depends on the changes of other climate variables (solar radiation, humidity and wind speed) and temporal and spatial variability of these variables within the time horizon of the projection analysis and within the spatial extent of the study region. The model developed in this study helps translate these trends and variabilities associated with climate change into quantitative estimates for the magnitude of impact on groundwater storage under the presence of intensive of agricultural groundwater water use.
Box plots indicating distributions of projected annual average groundwater depletion rates during the projection period (2012-2050) are shown in Fig. 10 (a). The model behavioral parameter sets identified through the GLUE analysis were used to generate these projections. The projected groundwater depletion rates vary over a wide range. Compared to the estimated groundwater storage depletion rates in the base period, slightly accelerated groundwater depletion rates were projected under the RCP 4.5 scenario while decreased depletion rates were found under the other three scenarios. The histogram in Fig. 10 (b) further summarizes the differences between groundwater depletion rates in the projection and base periods. Model realizations associated with all identified behavioral parameter sets and under all RCP scenarios were combined to construct the histogram, and all realizations were assigned equal probability. According to the histogram, there is a 84 % chance that the groundwater depletion rate would be lower, with a reduction of up to 10 km 3 H 2 O/yr. The projected increases in precipitation under climate change (Table 3) may lead to an anticipation that there could be amelioration in groundwater depletion. The simulation results in this projection analysis confirmed the likely existence of the ameliorating effect. However, in all cases the depleting trend will continue. Even with the most optimistic forecast, there will still be a deficit of 9 km 3 H 2 O in the annual groundwater storage budget (Fig.10a). The model results are therefore certainly not a call to sit back and relax but instead a call for more systematic efforts to measure, monitor and manage the important groundwater resources of India with a focus on arresting depletion trends. The variability of results across RCP scenarios furthermore highlights the challenges of managing groundwater resources in the face of uncertainty brought about by climate change.

Conclusions and discussions
India is a country with intensive groundwater-fed irrigated agriculture; and the large extraction of groundwater for food production has created concerns about the sustainability of groundwater resources for both current and future citizens. In this paper, we incorporated groundwater use processes for irrigation into a national-scale SWAT model of India. The development of the model was based on publicly available data, and the model was evaluated and calibrated using GRACE TWS data. The results of the model identification confirm that the integrated hydrologic-groundwater use model can capture groundwater storage depletion in northwest India and therefore can serve as a promising predictive modeling tool to quantify the anthropogenic impact of irrigation on groundwater resources in the region.
The utility of the developed model is demonstrated in a case study presented in section 4, in which the model was applied to assess the impact of climate change on groundwater storage in northwest India under four RCP scenarios. The assessment results suggest that, without consideration of the effect of changes in cropping patterns and irrigation water use efficiency, it is well possible that climate change will ameliorate the groundwater deficit in northwest India. However, the beneficial effect is not strong enough to reverse depletion trends.
Mathematical models provide simplified representations of real-world systems and processes and the modeling exercises inevitably involve uncertainties. Discussions of some key uncertainties and limitations arising from model development and applications have been provided in the previous sections. They are summarized here with a few more added.
Firstly, in this study we used the GLUE approach to treat parametric uncertainty in the modeling. As an additional remark from the GLUE analysis, while GRACE TWS data prove to provide highly-informative information on spatial patterns of vertically-integrated water storage variations and are of great value for understanding agricultural groundwater use processes in northwest India, the large remaining uncertainty in the model parameters, which were indicated by the wide ranges of estimated groundwater depletion rates, suggests that it is desirable, when possible, to introduce observations on other hydrologic variables to provide further constraints in model identification.
In addition to the parametric uncertainty, there are also uncertainties arising from other sources. As already noted, a major limitation in our study is the omission of surface water irrigation, which may bias estimates for those variables in the groundwater balance. Additional uncertainty related to the CO 2 effect of plant evaporation exists in projections of the impact of climate change on groundwater balances. It is well known that plant growth can benefit from CO 2 enrichment that helps improve the water use efficiency of plants. SWAT incorporates an algorithm to simulate the effects of atmospheric CO 2 enrichment. However, the algorithm was formulated based on data obtained in a controlled experimental environment (Stockle et al., 1992) on which experiments from open-air field conditions (Long et al., 2006) cast uncertainty. Given this uncertainty, the CO 2 effect was not simulated in this study. Lastly, there is also uncertainty associated with input data (f. ex. on irrigated crop area and climate) and GRACE TWS data processing, which are exemplified by the existence of multiple data products in each category (Sakumura et al., 2014;Anderson et al., 2015;Gehne et al., 2016). Choosing different data products may have implications for model identification and simulation results.
Finally, as indicated in section 4, the projection analysis presented in this paper is partial because potential cropping pattern changes in groundwater-fed irrigated agriculture were not considered. There is a need to build more integrated modeling systems to support more comprehensive analysis for better national polices on groundwater management and development of groundwater-fed irrigated agriculture. We will endeavor to investigate these uncertainties and limitations in future work.

Declaration of Competing Interest
The authors declare that they have no conflict of interest.