Groundwater rejuvenation in parts of India influenced by water-policy change implementation

The dwindling groundwater resource of India, supporting almost one fifth of the global population and also the largest groundwater user, has been of great concern in recent years. However, in contrary to the well documented Indian groundwater depletion due to rapid and unmanaged groundwater withdrawal, here for the first time, we report regional-scale groundwater storage (GWS) replenishment through long-term (1996–2014, using more than 19000 observation locations) in situ and decadal (2003–2014) satellite-based groundwater storage measurements in western and southern parts of India. In parts of western and southern India, in situ GWS (GWSobs) has been decreasing at the rate of −5.81 ± 0.38 km3/year (in 1996–2001) and −0.92 ± 0.12 km3/year (in 1996–2002), and reversed to replenish at the rate of 2.04 ± 0.20 km3/year (in 2002–2014) and 0.76 ± 0.08 km3/year (in 2003–2014), respectively. Here, using statistical analyses and simulation results of groundwater management policy change effect on groundwater storage in western and southern India, we show that paradigm shift in Indian groundwater withdrawal and management policies for sustainable water utilization appear to have started replenishing the aquifers in western and southern parts of India.


S2 GRACE mascon solution
S3 Satellite-based estimates using combination of 3 land surface models and using individual ones S4 PCR-GLOBWB simulation for detailed study area S5 Statistical analyses used

S6 Specific yield information
We used long term monthly mean precipitation data sets [Global Historical Climatological Network (GHCN)] from the year 1960 to 2010 for more than 120 locations spread over the study area. Continuous data were not available for several locations and hence the data were screened and restricted the selection of a location on the basis of availability of at least 70% of continuous time series data. Accordingly, we selected 37 locations that are evenly distributed over the entire study area. We used European Center for Medium Range Weather Forecasting (ECMWF) reanalysis 30 (ERA Interim) simulation output of specific humidity (SH) over the selected locations between 1979 and 2012 to constrain the boundary of the hydro-meteorological seasons. Precipitation as well as SH data exhibited higher values during monsoon and distinguishing the boundary between monsoon and non-monsoonal seasons. Pre-monsoon and post-monsoon seasons were defined by determining the rate of change in SH data, as precipitation is negligible during these seasons for most of the locations. The SH values showed decrease in trend immediately after monsoon, and then increasing trend from the advent of the monsoon. This sudden shift (decreasing to increasing) in SH values marked the boundary between postmonsoon and pre-monsoon in the present study. After defining hydrological seasons, we imposed inter-quartile range cut off for selecting the months that composed the seasons and accordingly, defined five hydro-meteorologic zones (HMZs). For example, location-wise longterm monthly mean of precipitation and SH values are shown in Supplementary Figure 1. Each two locations represent data from a HMZ in the following order, northern, western, central, eastern and southern. The seasons for the hydrologic year were defined for each of the HMZs as: Northern region [monsoon: June-September; post-monsoon: October-January; pre-monsoon:   (Watkins et al., 2015). Monthly gravity field variations were parameterized using 4,551 equal-area 3 0 surface spherical cap mascons. Inter-satellite range rate measurements between the two GRACE satellites are used for estimating variations in gravity field through application of partial derivatives. Degree 2 and order 0 coefficients are replaced in JPL processed level 1 GRACE observations by the coefficients derived from Satellite Laser Ranging reported by Cheng et al. (2011). Degree 1 coefficients are derived from a process described elsewhere (Swenson et al., 2008). Glacial isostatic adjustment (GIA) in the data has been removed following the process developed by Geruo and Wahr (2013). In order to minimize the measurement errors, a-priori constraints have been applied in gravity fields of equal-area spherical mascons.

S3 Satellite-based estimates using combination of 3 land surface models and using individual ones
We have performed the estimation of satellite-based GWS using a combination of 3 land surface models (LSMs) and also using estimates from individual models between 2005 and 2013, details can be found in Bhanja et al. (2016). The RMSE and Pearson's correlation coefficient are shown for the GWS obtained using 3 LSMs' combination and using each of the LSMs with the GWS obtained from in situ observation wells on 12 largest river basins in India (Bhanja et al., 2016). On the basis of selection of an individual model's estimate, the satellite-based GWS showing good match with in situ GWS in some of the basins, however, none of the models was showing consistent result while comparing the estimates in all of the river basins. In general, the combination of the models provides best satellite-based estimates (details can be found in Bhanja et al., 2016).

S4 PCR-GLOBWB simulation for detailed study area
In PCR-GLOBWB, total groundwater recharge (R sim ) is a combination of diffuse recharge through precipitation, and irrigational return flow; recharge was estimated by subtracting capillary rise from deep percolation [Wada et al., 2014]. The capillary rise is calculated by considering moisture flux in an upward direction under an existing upward gradient condition with soil moisture concentration lying below field capacity, and the amount could not exceed groundwater storage of underlying groundwater layer [Wada et al., 2014]. In order to produce a more realistic soil moisture, evaporation and transpiration, the irrigation scheme separately parameterizes paddy and non-paddy crops. The scheme also dynamically balance daily amount of surface and soil water by considering feedback from irrigation water. In the paddy fields, a 50 mm surface water depth [Wisser et al., 2010] was considered up to a late crop-growing stage and no irrigation phase was considered (~20 days) before the harvesting period [Wada et al., 2014].
Irrigated area was estimated from country-specific information provided by FAOSTAT (http://faostat3.fao.org/home/E). Domestic and industrial water requirements were also calculated on a daily basis. Domestic and livestock water requirement data in 2000 were collected from FAOSTAT and the daily water requirement was calculated as a function of temperature [Wada et al., 2011]. Change in water use intensity as a factor of economic and technological development also considered here [Wada et al., 2014]. Simulated Topological Networks [Vörösmarty et al., 2000] data were employed to simulate direct runoff, interflow and base flow. Irrigation water allocation was calculated from available surface water resources and readily available groundwater.
Tall and short vegetation, paddy rice, non-paddy crops, rain-fed crop, open water, different soil types are considered following FAO Digital Soil Map of the World (FAO, 2003) in each of the grid cell. The improved ARNO scheme are used to estimate saturated soil area fractions (Todini, 1996;Hagemann and Gates, 2003) and the groundwater depth frequency distribution are obtained from the HYDRO1k Elevation Derivative Database (HYDRO1k; U.S. Geological Survey Center for Earth Resources Observation and Science; https://lta.cr.usgs.gov/HYDRO1K/). The Hypothetical Test Scenario (HTS) considered in this study includes model simulation using all of the above mentioned data along with policy change related parameters in Gujarat. Global water use simulations were done using dynamic meteorological forcing parameters, i.e. temperature and precipitation data from ERA Interim and MERRA reanalysis products and reference evapotranspiration was calculated following FAO guidelines (with Penman-Monteith equation) in a grass surface [Wada et al., 2014].

Hodrick-Prescott (HP) filter analysis:
We have tried to estimate trend using linear regression, although non-linearity in the data restricts us to rely totally in linear trend analysis. We used Hodrick-Prescott (HP) filter [Hodrick and Prescott, 1997], a non-parametric, non-linear trend analysis widely used over the years. The HP filter separates trend (T t ) and cyclical (c t ) components from the data time series (y t ).
In order to separate the trend component, the HP filter solves the following equation: The smoothing parameter () is a positive number, which reduce the variability in the cyclical component [Hodrick and Prescott, 1997]. We have chosen the value of  for the quarterly data, 1600 [Hodrick andPrescott, 1997;Ravn and Uhlig, 2002]. Detailed descriptions on HP filter can be found in Hodrick and Prescott [1997] and Ravn and Uhlig [2002].

Bayesian VAR model analyses:
We have estimated Bayesian VAR and tested Granger causality for Gujarat and Andhra Pradesh. Lag values of surface water irrigation (SWI) has a significant (p value <0.05) positive impact on GWS anomaly in Andhra Pradesh. Lag values of groundwater irrigation (GWI) has a significant (p value <0.05) negative impact on GWS anomaly in Gujarat-this implies GWI Granger causes GWS anomaly in Gujarat, as a result, GWS increases after 2002-2003 due to reduction in electricity usage for groundwater pumping.

S6 Specific yield information
High-resolution specific yield map has been created using hydrogeology map of India (Fig. 1c). The map is constructed following CGWB (2012b) and Mukherjee et al. (2015). We segregate the S y information from the CGWB database (CGWB, 2012) and assigned the data according to aquifer characteristics shown in Fig. 1c. The range in S y values were determined from long-term pumping test data conducted over the years (MWR, 2009). On the basis of available information, we prepare a gridded (0.1 0 × 0.1 0 ) S y map for the Indian region. A unique specific yield value is assigned to each of the selected wells based on their geographic location. Observed GWS anomaly has been calculated by multiplying S y and ∆h values at all the in-situ well locations.