Satellite-based estimates of groundwater depletion in India

Groundwater is a primary source of fresh water in many parts of the world. Some regions are becoming overly dependent on it, consuming groundwater faster than it is naturally replenished and causing water tables to decline unremittingly. Indirect evidence suggests that this is the case in northwest India, but there has been no regional assessment of the rate of groundwater depletion. Here we use terrestrial water storage-change observations from the NASA Gravity Recovery and Climate Experiment satellites and simulated soil-water variations from a data-integrating hydrological modelling system to show that groundwater is being depleted at a mean rate of 4.0 ± 1.0 cm yr-1 equivalent height of water (17.7 ± 4.5 km3 yr-1) over the Indian states of Rajasthan, Punjab and Haryana (including Delhi). During our study period of August 2002 to October 2008, groundwater depletion was equivalent to a net loss of 109 km3 of water, which is double the capacity of India’s largest surface-water reservoir. Annual rainfall was close to normal throughout the period and we demonstrate that the other terrestrial water storage components (soil moisture, surface waters, snow, glaciers and biomass) did not contribute significantly to the observed decline in total water levels. Although our observational record is brief, the available evidence suggests that unsustainable consumption of groundwater for irrigation and other anthropogenic uses is likely to be the cause. If measures are not taken soon to ensure sustainable groundwater usage, the consequences for the 114,000,000 residents of the region may include a reduction of agricultural output and shortages of potable water, leading to extensive socioeconomic stresses.


Satellite-based estimates of groundwater depletion in India
Matthew Rodell 1 , Isabella Velicogna 2,3,4 & James S. Famiglietti 2 Groundwater is a primary source of fresh water in many parts of the world. Some regions are becoming overly dependent on it, consuming groundwater faster than it is naturally replenished and causing water tables to decline unremittingly 1 . Indirect evidence suggests that this is the case in northwest India 2 , but there has been no regional assessment of the rate of groundwater depletion. Here we use terrestrial water storage-change observations from the NASA Gravity Recovery and Climate Experiment satellites 3 and simulated soil-water variations from a dataintegrating hydrological modelling system 4 to show that groundwater is being depleted at a mean rate of 4.0 6 1.0 cm yr 21 equivalent height of water (17.7 6 4.5 km 3 yr 21 ) over the Indian states of Rajasthan, Punjab and Haryana (including Delhi). During our study period of August 2002 to October 2008, groundwater depletion was equivalent to a net loss of 109 km 3 of water, which is double the capacity of India's largest surface-water reservoir. Annual rainfall was close to normal throughout the period and we demonstrate that the other terrestrial water storage components (soil moisture, surface waters, snow, glaciers and biomass) did not contribute significantly to the observed decline in total water levels. Although our observational record is brief, the available evidence suggests that unsustainable consumption of groundwater for irrigation and other anthropogenic uses is likely to be the cause. If measures are not taken soon to ensure sustainable groundwater usage, the consequences for the 114,000,000 residents of the region may include a reduction of agricultural output and shortages of potable water, leading to extensive socioeconomic stresses.
Groundwater responds more slowly to meteorological conditions than the near-surface components of the terrestrial water cycle 5 . Its residence time (the ratio of quantity in storage to average rate of recharge or discharge) ranges from months in shallow aquifers to a million or more years in deep desert aquifers 6 . Hence, groundwater can be slow to recover from perturbations to its state of dynamic equilibrium. In particular, withdrawals can easily surpass net recharge in arid and semi-arid regions where people depend on fresh groundwater for domestic needs and irrigation 1 . Despite the increasing pressure placed on water resources by population growth and economic development, the laws governing groundwater rights have not changed accordingly, even in developed nations 7 . Nor is groundwater depletion limited to dry climates: pollution and mismanagement of surface waters can cause over-reliance on groundwater in regions where annual rainfall is abundant.
India now suffers severe water shortages in many of its states. It averages about 120 cm yr 21 of precipitation, which is more than any other country of comparable size 8 , but the rain is unevenly distributed. In New Delhi, India's richest city, most middle-class residents do not have a dependable source of clean water (Sengupta, S., 'In India, water crisis means foul sludge', New York Times, 29 September 2006). The World Bank has warned that India is on the brink of a severe water crisis 9 . Nationally, groundwater accounts for about 50-80% of domestic water use and 45-50% of irrigation 8,10 . Total irrigated area in India nearly tripled to 33,100,000 ha between 1970 and 1999 11 . In neighbouring Pakistan, which is largely arid, groundwater is essential for much of the country's agriculture. Competition for precious water in transboundary aquifers is likely to exacerbate already strained relations between the two nations.
India's government is aware that groundwater is being withdrawn at unsustainable rates in some areas, and in 1986 it established a Central Ground Water Authority with the power to regulate groundwater development 12 . However, as in other nations composed of smaller sovereignties and encompassing competing interests that have become dependent on a certain level of water availability, it is difficult to implement a coordinated and appropriately stringent response. Political and aquifer boundaries bear no resemblance to each other, and aquifers themselves are interconnected, so that one state's (or country's) groundwater management practices are likely to affect its neighbour. Holistic regional groundwater assessments would be valuable in promoting appropriate policies and for hydrologic research, but such assessments are difficult to generate on the basis of well surveys, which are typically unsystematic.
The Gravity Recovery and Climate Experiment (GRACE) satellite mission, launched by NASA and the German Aerospace Centre (DLR) in 2002, measures temporal variations in the gravity field, which can be used to estimate changes in terrestrial water storage (TWS) 3 . Although its spatial resolution (no better than ,160,000 km 2 ) and temporal resolution (ten day to monthly) are low in comparison with those of other Earth-observing satellites, GRACE has the major advantage that it senses water stored at all levels, including groundwater. Unlike radars and radiometers, it is not limited to measurement of atmospheric and near-surface phenomena.
Groundwater storage variations can be isolated from GRACE data given auxiliary information on the other components of TWS, from either in situ observations 13 or land-surface models 14 . We used the second approach to produce a time series of groundwater storage anomalies (deviations from the mean state) averaged over the area encompassed by Rajasthan (342,239 km 2 ), Punjab (50,362 km 2 ) and Haryana (45,695 km 2 including the National Capital Territory of Delhi) between August 2002 and October 2008. This region was chosen because the Indian Ministry of Water Resources estimates that groundwater withdrawals in each of the three states exceed recharge 2 (Fig. 1). Figure 2 maps the averaging function used to retrieve regional TWS time series from the GRACE data.
Rajasthan, Punjab and Haryana are semi-arid to arid, averaging about 50 cm of annual rainfall overall [15][16][17] , and encompass the eastern part of the Thar Desert. The 114,000,000 residents of the region have benefitted from India's 'green revolution', a massive agricultural expansion fuelled largely by increased production of groundwater for irrigation, which began in the 1960s. Wheat, rice and barley are the major crops. The region is underlain by the Indus River plain aquifer, a 560,000 km 2 unconfined-to-semiconfined porous alluvial formation that straddles the border between India and Pakistan 11 . TWS variations observed by GRACE include the combined contributions of groundwater, soil water, surface water (lakes, rivers, canals and rice paddies), snow, ice and biomass. Interannual variations in biomass have been shown to be well below the detection limits of GRACE 18 . The climate of the study region is warm and the terrain shallow; hence, snow is uncommon. However, glaciers abound in the Himalayas, as close as 100 km northeast of Punjab, so we assessed the potential for leakage into our region of the gravity signal associated with non-seasonal glacier melt. We determined that even if 13.4 km 3 yr 21 of glacier loss (equal to the estimated rate of melt for all of the Himalayas during 1962-2004 19 ) occurred along the 150-km stretch of mountains (6% of the Himalayan range) closest to Punjab and Haryana, only 2.8 km 3 yr 21 (15%) of the estimated TWS trend could be attributed to glacier retreat. We assessed potential leakage from other adjacent trends and included the effects in our error budget.
Surface-water storage also merits consideration. The 203 main reservoirs 20 within or on the border of the three states, plus a large salt lake in Rajasthan, occupy a total area of 4,320 km 2 and have a gross storage capacity of 39.5 km 3 . On the basis of reservoir storage reports from the Indian Central Water Commission, seven of the largest reservoirs, totalling 34.4 km 3 of gross storage capacity, experienced a net increasing trend of 0.5 km 3 yr 21 during the study period. The total area of rice paddies in the three states is 38,061 km 2 (8.7% of the region), but all are seasonally flooded, shallow-water (,1-m) paddies 21 , which dry annually. Hence, if anything, surface-water storage changes probably reduced the rate of TWS depletion observed by GRACE.
Root-zone soil water is often the dominant contributor to TWS variability in temperate regions 22 . To remove its effect and thus isolate groundwater storage changes, we estimated soil-water storage variations by averaging results from five simulations of the Global Land Data Assimilation System (GLDAS) 4 . These simulations used combinations of three land-surface models and three meteorological-forcing (input) data sets. Total soil-column depths in the models ranged from 200 to 340 cm; thus, they did not account for water storage variations in deep unsaturated soil. However, sub-root-zone soil dries only by gravity drainage or by diffusion to drier layers above. The lack of a drying trend in the root zone (as described below) indicates that deep soil-water storage was likewise stable. Figure 3 shows the resulting 6-yr time series of monthly groundwater storage anomalies as equivalent heights of water averaged over the three-state region. It compares favourably with the TWS and soilwater time series. Its seasonal cycle (Fig. 3    is 5.9 cm, compared with 10.3 cm for soil water and 13.8 cm for TWS. These relationships are consistent with those reported in previous studies of soil-water/groundwater covariability 5,14,22 . The soil-water time series reflects rainfall anomalies during the period (discussed below) and exhibits no significant trend. On the other hand, TWS and groundwater decline steadily from 2003 onwards. We calculate the rate of depletion of groundwater to be 4.0 6 1.0 cm yr 21 . Assuming 2 a specific yield of 0.12, the regional mean rate of water Groundwater levels also appear to be declining quickly in western Uttar Pradesh, to the east of Haryana. If there is groundwater depletion in Pakistan, to the northwest, it seems to be much less severe. Although six years is a short period from which to assess a longterm trend with confidence, two pieces of evidence support our conclusion that severe groundwater depletion is occurring as a result of human consumption rather than natural variability. First, the Indian Ministry of Water Resources reports that groundwater withdrawals exceed recharge in the three states we studied 2 . Irrigation accounts for about 95% of the consumption 2 ; about 28% of the area is irrigated 23  GLDAS modelled soil-water fields integrate the effects of precipitation, solar radiation, air temperature and other meteorological factors that directly or indirectly influence groundwater storage 4 . The trend in simulated soil-water storage during the period of study was 0.4 cm yr 21 . This supports the notion that groundwater declines were not caused by natural climate variability. It also confirms that the computed groundwater trend is not a mathematical artefact caused by the subtraction of a large positive soil-water trend from the GRACE-derived TWS trend.
We conclude that withdrawals for irrigation and other uses are depleting the groundwater reserves of Rajasthan, Punjab and Haryana at a rate of 4.0 6 1.0 cm yr 21 equivalent height of water, or 17.7 6 4.5 km 3 yr 21 . The Indian Ministry of Water Resources reports that the difference between annual available recharge and annual withdrawals in the region is a 13.2 km 3 yr 21 deficit 2 . Our result implies that the portion of irrigated water that replenishes the aquifers is less and/or the rate of withdrawal is more than the Indian government has estimated. Apparently, most of the groundwater withdrawn subsequently is lost from the region as a result of increases in run-off and/or evapotranspiration. Between August 2002 and October 2008, the region lost 109 km 3 of groundwater, which is double the capacity of India's largest reservoir, the Upper Wainganga, and almost triple the capacity of the largest man-made reservoir in the United States, Lake Mead. Depletion is likely to continue until effective measures are taken to curb groundwater demand or until the supply or quality of the resource is diminished to the point at which farmers and residents of the region are forced to react. Severe shortages of potable water, reduced agricultural productivity, conflict and suffering surely would accompany the supplylimited solution.

METHODS SUMMARY
We used 73 monthly GRACE gravity solutions generated by the Center for Space Research at the University of Texas at Austin, consisting of sets of spherical harmonic (Stokes) coefficients, C lm and S lm , to degree l and order m, both #60. After removing the temporal mean to obtain gravitational anomalies, we filtered 24 each field to reduce noise and converted it to mass in units of equivalent water thickness. We then computed the uncertainty associated with the GRACE measurements 25,26 .
To compute the groundwater storage time series, we removed GLDAS 4 estimates of soil-water storage variations from the GRACE TWS 14 . For consistency, we applied degree-60 truncation and the same filtering process to the soil-water time series. Uncertainty in the GLDAS modelled soil-water time series was computed as the standard deviation of results from the five contributing simulations 27 .
We then generated time series of monthly TWS, soil-water and groundwater storage in the study region by applying an averaging kernel (Fig. 2) to the Stokes coefficients. The kernel is spatially smoothed primarily as a consequence of three factors: the truncation of the Stokes coefficients to degree 60, the filtering process and convolution with a Gaussian function. When mass variations inside and outside the region differ, such smoothing causes the amplitude of the retrieved signal to be damped 26,28 . To offset this effect, we computed and applied a scaling factor of 1.95. We estimated leakage from the surrounding regions using various scenarios to determine all possible impacts. The total uncertainty estimate for the groundwater trend, 1.0 cm yr 21 , combines errors associated with the GRACE measurements (0.54 cm yr 21 ), GLDAS soil water (0.38 cm yr 21 ), scaling (0.50 cm yr 21 ) and leakage (0.60 cm yr 21 ). An expanded description of the methods used is provided in the online Methods.
Full Methods and any associated references are available in the online version of the paper at www.nature.com/nature.

METHODS
We used 73 monthly GRACE gravity field solutions generated by the Center for Space Research at the University of Texas at Austin 3 . Each monthly gravity field consisted of a set of spherical harmonic (Stokes) coefficients, C lm and S lm , to degree l and order m, both #60. The C 20 coefficients showed unreasonable variability, so we replaced them with values derived from satellite laser ranging 29 . After removing the temporal mean to obtain gravitational anomalies, we filtered 24 each monthly field to reduce noise ('striping') and converted it to mass in units of equivalent water thickness.
We corrected the GRACE mass anomalies for the solid-Earth contributions generated by the high-latitude Pleistocene deglaciation using an independent model 30 . However, the contribution was negligible in this analysis.
We used GLDAS 4 time series of soil-water storage to isolate groundwater storage variations from the GRACE TWS anomalies 14 . For consistency, before we calculated Stokes coefficients for the soil water storage we truncated them to degree 60 and applied the same filtering process as used for the GRACE fields.
We used the filtered, monthly coefficients to generate time series for TWS, soil water and groundwater storage in the study region. To do this we constructed an averaging kernel by convolving a Gaussian function with a gridded map of the combined area of Punjab, Haryana (including Delhi) and Rajasthan (438,296 km 2 ). After testing various radii for the Gaussian function, from 300 km to 0 km, we selected an exact (0-km) averaging function (that is, no Gaussian averaging) because it minimized uncertainty and was generally best suited for our region (see below).
We applied the kernel to the Stokes coefficients to estimate mass changes averaged over the region. Ideally, the kernel should take the value 1 within a study region and 0 outside. In practice, the kernel is smoothed as a consequence of three factors: the truncation to degree 60 (the maximum degree for the GRACE Stokes coefficients), the filtering process and the convolution with the Gaussian function. Such smoothing causes an amplitude damping of the recovered mass. To obtain a realistic estimate of the mass changes, the damping of the signal amplitude caused by smoothing must be mitigated 26,28 . Also, we must account for potential leakage of mass changes from adjacent regions. In particular, the leakage effects may or may not be significant. In cases where mass changes are uniform across the region's boundaries, bias is minimal 25,26 and the signal recovered using the averaging kernel does not need to be rescaled because the smoothing effects are compensated by leakage of identical variations from the outside regions. However, if mass variability immediately outside the region is uncorrelated with that inside or they are negatively correlated or have different amplitudes, the recovered signal will be biased 28 . This must be evaluated on a case-by-case basis. We determined, on the basis of the GLDAS soil-water fields and the GRACE fields, that seasonal variations in TWS inside and outside the region are well correlated. In fact, the exact average of the GLDAS soil water over the study region and the average calculated when applying the truncation, filtering and Gaussian averaging differed only by few per cent.
On the other hand, the secular change (trend) varies between the inside and the outside of the study region, with the negative anomaly mainly concentrated inside the region. Therefore, to recover the total mass trend we rescaled the time series to account for amplitude damping 26,28 . The scaling factor was estimated by distributing a synthetic mass-change signal uniformly over the study region, processing it in the same manner as the GRACE data (that is, converting it to the spectral domain, truncating it to degree 60 and filtering and spatially averaging it) and comparing the retrieved signal with the original synthetic signal. We then multiplied the GRACE time series trend by this factor, which was determined to be 1.95. The scaling depends on how mass change is distributed within the region; therefore, to account for the uncertainty associated with our representation of that distribution, we calculated the sensitivity to different mass-change distributions within the region. We found that mass distribution effects caused an uncertainty of 0.5 cm yr 21 , which we included in our final error budget.
We also evaluated the leakage from mass-change signals in the surrounding regions. To estimate the leakage from mass loss associated with non-seasonal Himalayan glacier melt, we distributed a trend of 13.4 km 3 yr 21 , equal to the estimated rate of melt for the all of the Himalayas during 1962-2004 19 , along the 150-km stretch of glaciers closest to Punjab and Haryana, and applied our averaging function and scaling. This is a highly conservative approach, given that 150 km is only 6% of the length of the mountain range. The portion of that signal that leaked into our region was 2.8 km 3 yr 21 , or about 15% of the estimated groundwater trend. Because that is almost surely a significant overestimate, we did not use it to adjust the estimated groundwater trend, but continued with a more careful analysis of leakage effects. In that analysis, we included leakage from all TWS trends in the surrounding region that were detected by GRACE. We investigated different leakage scenarios to determine all possible effects on our final estimate. In the most likely scenario, the leakage effects nearly cancel. Nevertheless, uncertainty due to leakage was estimated to be 0.6 cm yr 21 , which is included in the final error budget.
The decision to use an exact (radius, R 5 0 km) Gaussian function was based on the following reasoning. We first evaluated the effects of filtering in the studied region by comparing trend maps obtained by applying 300-, 250-, 200-, 150-, 100-, 50-and 0-km Gaussian averaging. Even if R 5 0 km and there is no Gaussian smoothing, the resulting map will still be smoothed as a consequence of the necessary truncation of the GRACE solutions to degree 60. We found that the seven smoothing radii produced insignificant differences in our region; thus, none of them was eliminated from consideration.
We constructed seven different averaging kernels by convolving Gaussian functions for each of the seven radii with a gridded map of the region of interest. We calculated TWS time series using the seven kernels, and we estimated the corresponding scaling factors and measurement errors. Both the scaling factor and the error are a function of the specific kernel (note that in the case of R 5 0 km, use of the scaling factor is necessitated only by the truncation to degree 60 and filtering of the GRACE solutions). For example, the averaging kernel produced by convolving a Gaussian function with R 5 250 km with a gridded map will deliver results that are more smoothed than in the case in which R 5 200 km (that is, the associated scaling factor will be larger for R 5 250 km than for R 5 200 km) but will have smaller measurements errors.
A priori, we did not know which would be the best averaging function because of this trade-off between smoothing and errors. Also, it is possible for the net effect of striping to nearly cancel itself out, depending on the specific location. We compared the TWS time series, scaling factor and errors for each of the seven averaging functions and found that the GRACE measurement errors (defined by the scatter of the monthly values about their seasonal cycle) were not significantly larger in the case of R 5 0 km than in the case of any of the greater radii (again, truncation to degree 60 produces some smoothing, even in the case of R 5 0 km).
Uncertainties in the Stokes coefficients were determined by assuming that the scatter of the monthly values about their seasonal cycle is due entirely to errors 25 . This represents the upper bound on the random component of the error. The estimate is conservative, because intra-annual variations in the signal will be interpreted as error. The 1s error estimates in the spatially averaged GRACE time series were then calculated from the uncertainty in the individual Stokes coefficients 26 . The scaled root-mean-squared error in the trend is 0.54 cm yr 21 .
Uncertainty in the GLDAS modelled soil-water time series was computed as the standard deviation of results from the five contributing simulations 27 . This calculation was done separately for the seasonal and secular components. The uncertainty in each monthly soil-water estimate for the study region is 1.71 cm and uncertainty in the trend is 0.38 cm yr 21 . The total error estimate for the ground water trend, 1.0 cm yr 21 , combines the GRACE measurement errors, the GLDAS error, the scaling error and the leakage error, under the assumption that the individual errors are uncorrelated.