Satellite gravity observation and hydrological modelling-based integrated groundwater storage change in Northwestern India

Satellite-based Gravity Recovery and Climate Exchange (GRACE) studies by Rodell et al. (2009) and Tiwari et al. (2009) identified severe groundwater depletion in the northwestern part of India (NWI) and Indo-Gangetic basin in northern India, respectively. Later, several studies were conducted to evaluate groundwater storage (GWS) change, groundwater depletion, and recharge scenarios in NWI and Indo-Gangetic basin using GRACE data, land surface model (LSM), and in situ groundwater level (GWL) data (Chen et al. 2014; Dasgupta et al. 2014; Bhanja et al. 2016; Long et al. 2016; Shamsudduha & Panda 2019; Chatterjee et al. 2020) showing an improved estimate of GWS change with more detailed spatiotemporal variations. GRACE gravity observation on GWS change has also been successfully studied in other parts of the globe, such as in the Mississippi River Basin, USA (Rodell et al. 2007), California Central Valley, USA (Famiglietti et al. 2011), Middle East and Arabian countries (Voss et al. 2013; Richey et al. 2015), Canning Basin, Australia (Richey et al. 2015), Central Asia (Forootan et al. 2014; Li et al. 2019), and in parts of China (Feng et al. 2018; Gong et al. 2018; Chen et al. 2019). It may be noted that monsoon rainfall is the primary source of groundwater recharge in India (Sinha & Sharma 1988). Spaceborne GRACE gravity-based monthly terrestrial water storage (TWS) change was used to derive monthly average GWS change. Besides, GWS change of a drainage basin can be directly estimated from GWL change during the observation period and indirectly by hydrological modelling of incoming and outgoing water fluxes from each of the spatial grids of the drainage basin. Using pre-monsoon and post-monsoon GWL data, GWL fluctuation-based monsoonal GWS change was estimated in the present study. However, due to the lack of adequate in situ observations, inherent complexity in the hydrogeological properties of subsurface sediments and rocks, and complex groundwater recharging processes (Döll et al. 2012), the accurate quantification of GWS change and its spatiotemporal variability are challenging problems (Chen et al. 2014). Presently, GWS change is estimated from spaceborne GRACE observations on a few hundred kilometres (∼300–400 km) spatial grids (Rodell et al. 2007; Landerer & Swenson 2012) and extrapolated into the finer scale of 1° × 1° or 0.5° × 0.5° spatial grids based on the simulation of TWS variations by land-hydrology models. However, quantifying GWS change based on the land-hydrology model is not straightforward and it is difficult to attain a high level of accuracy (Rodell et al. 2004; Döll et al. 2012; Chen et al. 2014) due to adverse terrain and geological factors governing groundwater infiltration and base flow. Thus, it is unable to reflect finer spatial details in GWS change. Alternatively, hydrological modelling estimates the balance surface water volume undergoing infiltration during the monsoon by a simple water balance (SWB) model (Schaake et al. 1996). It represents groundwater recharge vis-à-vis GWS change in the monsoon. Based on hydrological model-based groundwater recharge estimates, groundwater withdrawal may be adjusted (Zareian & Eslamian 2021), and conjunctive use of surface water and groundwater may be promoted for groundwater sustainability. A hydrological model duly considers terrain conditions and hydrogeological properties of surface materials of an area. However, it does not consider groundwater withdrawal from the aquifers and therefore represents underestimated GWS change.

Spaceborne GRACE observation and hydrological model-based GWS change may be integrated to reflect finer spatial details with due consideration to terrain, hydrogeological properties of surface materials, hydrometeorological factors, and groundwater withdrawal conditions. The integrated GWS change product may be obtained by combining two TWS change data layers derived from GRACE observation and hydrological modelling, followed by removing significant hydrological components (e.g., soil moisture change) other than GWS change. In multiple linear (multi-linear) regression analysis of individual TWS change products with significant hydrological components, R² values (coefficient of determination) of regression represent how well the governing hydrological components (independent variables) fit the regression models (the goodness of fit) with individual TWS change products. Finally, using the R² values obtained from two sets of multi-linear regressions, an integrated TWS change product can be generated. After removing the significant hydrological components other than GWS change (e.g., soil moisture change) from the integrated TWS change product, an integrated GWS change product may be obtained. By spatial correlation analysis of derived GWS change products with in situ GWS change, the evaluation of the individual products and the improvement in the integrated product can be assessed. This paper presents an approach for operational but robust groundwater recharge estimates, which is immensely important to formulate a strategy to address severe groundwater depletion in the northwestern part of India, covering the states of NCR Delhi, Haryana, Punjab, Rajasthan, Gujarat, and Uttar Pradesh. This paper is organised as follows: Section 2 describes the study area and data preparation. Section 3 deals with the detailed procedure for estimating satellite observation, hydrological models, and in situ GWL fluctuation-based GWS change products. Finally, the procedure for integrating satellite observation and hydrological model-based GWS changes was elaborated. Section 4 provides an in-depth discussion of the results. Finally, Section 5 concludes with the salient findings of the study.

IRS Resourcesat-2 AWiFS multispectral data with a spatial resolution of 56 m and a combined swath of 740 km acquired during November–December 2013 was used. Resourcesat AWiFS data have four 10-bit spectral bands such as green (0.52–0.59 μm), red (0.62–0.68 μm), near infra-red (0.77–0.86 μm), and short-wave infra-red (1.55–1.70 μm). After pre-processing the data for radiometric and geometric corrections, the false colour composite (FCC) images for the study area were generated to visualise and interpret land use land cover classes, geomorphology, and broad litho-contacts. The paper presented the standard FCC image (Bands432 in a red-green-blue colour scheme) of the study area from AWiFS data.

The paper used Shuttle Radar Topography Mission (SRTM) 1 arc-second DEM available in 1-degree × 1-degree tiles. The DEM is available in a geographic (lat/long) projection system with the WGS84 horizontal and EGM96 vertical datum.

Daily and monthly aggregated rainfall data for the study area were generated by combining satellite-based and ground-based observations. Satellite-based daily rainfall measurements with spatial resolution 0.25° × 0.25° were obtained from the Tropical Rainfall Measuring Mission (TRMM), a joint space mission of NASA, USA and JAXA, Japan. Ground-based rainfall data at 0.25° × 0.25° grid was obtained from the Indian Meteorological Department (IMD 2014). IMD generated the gridded daily rainfall data for the entire country by interpolating in situ rainfall measurements available from 6,995 rainfall stations. The present study used a combined rainfall product of TRMM's gridded rainfall data and IMD's gridded and point rainfall data. A 3 × 3 low-pass filter filtered the combined product to minimise the edge effect between the grids of combined rainfall products. The SCS-Curve Number method used gridded daily rainfall data for surface runoff estimation. Also, the monthly aggregated gridded rainfall data was obtained by aggregating daily rainfall data during the observation period.

Surface runoff flows over the surface of a drainage basin and through the river channels to go out of the system as outgoing water flux. The present study adopted the SCS runoff curve number method (Kent 1973) to calculate runoff from daily rainfall data at 0.25° × 0.25° spatial grids. The curve number depends on the hydrologic soil group, land use/treatment, and hydrologic condition. The hydrologic soil groups, land use, and slope maps were prepared from the existing data such as the NBSS-LUP soil map (NBSS-LUP 2002), land use land cover map (modified after Agrawal et al. (2003)), and SRTM 1 arc-second DEM.

It is important to note that surface runoff appears in the hydrograph after the initial demands of infiltration and surface storage are satisfied. The initial demand for surface storage is satisfied at the onset of the monsoon. During the monsoon months, a large number of rainy days, substantial daily rainfall, and high antecedent soil moisture content help reduce runoff estimation uncertainties.

Terrestrial evapotranspiration (ET) includes water vapour fluxes from soil evaporation, wet canopy evaporation, and plant transpiration at dry canopy surface. Evapotranspiration plays a crucial role in hydrological water balance. The MOD 16 dataset provides global evapotranspiration (ET) and potential evapotranspiration (PET) at 1 km × 1 km regular grids for the 109.03 million km² global vegetated land areas at 8 days, monthly and annual intervals (Mu et al. 2013). MOD 16 monthly ET data were used after upscaling into 0.25° × 0.25° spatial grids by pixel aggregated averaging. The MOD16 ET data products are generated by the algorithm of Mu et al. (2007, 2011) using MODIS land cover and fraction of photosynthetically active radiation/leaf area index (FPAR/LAI) data, and global surface meteorology data available from the Global Modeling and Assimilation Office (GMAO v. 4.0.0).

In this study, soil moisture data derived by the GLDAS Noah LSM on spatial resolution 0.25° × 0.25° and temporal resolution 1-month (Rodell & Beaudoing 2007) were used. The data were produced by advanced land surface modelling using satellite-based meteorological forcing data and limited ground-based measurements. Soil moisture is available in four layers: the first layer for 0–10 cm, the second layer for 10–40 cm, the third layer for 40–100 cm, and the fourth layer for 100–200 cm. Summing up the values in four layers, the volumetric soil moisture (volumetric soil water content) was obtained. Soil moisture represents a land surface state. For obtaining monthly average soil moisture change corresponding to GRACE TWS change, a long-term monthly average similar to the GRACE TWS time-mean baseline was derived, from which the deviation for each month was calculated. Finally, GLDAS 2.0 monthly average monsoonal soil moisture change was estimated for the observation period.

It is observed that GRACE GWS change is low in the Luni-Ghaggar drainage basin and the northwestern part of the Yamuna drainage basin covering the desert/arid region of Rajasthan, the whole of Haryana and NCR Delhi, the southern part of Punjab, the western part of Uttar Pradesh (U.P.), the hilly northern part of Uttarakhand and Kutch region of Gujarat states. On the other hand, it is moderate to high in the Narmada drainage basin, the southern part of the Chambal and Yamuna drainage basins, and the eastern part of the Upper Ganga drainage basin covering major parts of Madhya Pradesh (M.P.), Gujarat, the northern part of Maharashtra, and the eastern part of U.P. states, respectively.

The GWS change is conventionally monitored by in situ GWL change data. It gives precise point data, which is useful for local and regional-scale GWL monitoring, provided the monitoring wells are spatially dense and evenly distributed, representing a statistically random spatial distribution. The dug wells, which provide regular GWL (e.g., water table) data from shallow aquifers, are dense and evenly distributed (Figure 6(a)). In contrast, the piezometric wells, which provide GWL (e.g., piezometric head) data from deeper aquifers, are not evenly distributed, and their spatial density is not adequate at places (e.g., in parts of U.P., M.P., and Rajasthan states) (Figure 6(b)). Hence, lateral and vertical variations of groundwater resources are not adequately accounted for to achieve reliable spatial continuity. Moreover, it is cumbersome and labour-intensive to replicate the exercise on a regular basis. Furthermore, the availability of in situ GWL data on a regular temporal interval is a matter of concern. It was found that the observations were not temporally consistent for all the monitoring wells. Year-wise point maps on GWL change (Δh) were prepared where both pre-monsoon and post-monsoon observations are available. It overcomes the problem of temporal inconsistency of GWL data. Finally, from these point maps, year-wise interpolated Δh surface maps were prepared, and a stack of year-wise GWL change surface maps was prepared for further analysis. Recently, artificial intelligence adopting artificial neural networks, deep learning algorithms, advanced geospatial analysis techniques (e.g., wavelet and Gaussian process regression), and modelling approaches, several studies were made for predicting groundwater levels (Taormina et al. 2012; Gholami et al. 2015; Rajaee et al. 2019; Pandey et al. 2020; Afan et al. 2021; Band et al. 2021). A combination of in situ and predicted GWL data would hopefully address the spatial and temporal consistency of GWL data.

Presently, GWS change can be estimated from spaceborne GRACE observations on temporal gravity field variations of the Earth's surface on monthly time scales and a few hundred kilometres (∼300–400 km) spatial grids (Rodell et al. 2007; Landerer & Swenson 2012). However, based on the simulation of TWS variations by land-hydrology models, independent of the actual GRACE data, GRACE TWS estimates could be extrapolated into a finer scale of 1° × 1° or 0.5° × 0.5° spatial grids. But, it is difficult to model and quantify GWS changes accurately by land-hydrology models (Rodell et al. 2004; Döll et al. 2012; Chen et al. 2014). It is observed that 1° × 1° or 0.5° × 0.5° GRACE TWS change products do not show finer details on the spatial variation of GWS change, but they provide meaningful information on regional-scale groundwater withdrawal. Sometimes, adverse terrain and geological factors governing groundwater infiltration and base flow are not duly reflected in in situ GWS change observed in GRACE-derived regional estimates. For example, at locations G1–G8 (Figure 8(a)), the differences occur due to adverse terrain, geological, and/or base flow conditions (in G1–G5 due to terrain and geological conditions and G6–G8 due to high base flow conditions). On the other hand, in and around G9 and G10 (western U.P. and northern Rajasthan), GRACE GWS change is largely different from in situ GWS change over extended regions. In these areas, excessive groundwater withdrawal has occurred from the deeper aquifers, but very few piezometric wells are available to monitor the GWL of the rapidly depleting deeper aquifers (Figure 6(b)).

Alternatively, the hydrological model duly considers terrain conditions and hydrogeological properties of surface materials for estimating the GWS change. However, groundwater withdrawal for irrigating vast agricultural fields in NWI and domestic and industrial water supply are unaccounted for in the hydrological model-based GWS change. In and around the locations M1–M6 (Figure 8(b)) e.g., Patiala-Ambala urban area of Punjab-Haryana states, Lucknow-Kanpur urban area of U.P. state, Gwalior city and surrounding region of M.P. state, Ajmer-Vijainagar and Jalore-Jodhpur-Jaisalmer urban areas of Rajasthan state, and Banaskantha-Palanpur agricultural belt of Gujarat state, respectively, excessive groundwater pumping has been taking place to cater to domestic and industrial water supply in the urban and peri-urban areas and agricultural irrigation in the intensive agricultural belts. Excessive groundwater depletion in such areas is unaccounted for in hydrological model-based GWS change estimates. Besides, in and around M7 (Kutch wetland in Gujarat state), model-based GWS change is largely different from in situ GWS change estimate because the surface materials (saltpan, mud flat area – often inundated) possess extremely low infiltration capacity, and therefore, model-based GWS change is utterly low whereas in situ GWS change was essentially interpolated from the southern meadow area due to non-availability of GWL monitoring wells in wetland areas (saltpan, mud flat area – often inundated) (Figure 2(d)) which led to fairly anomalous in situ GWS change estimate.

The integrated GWS change estimate is significantly improved compared with the individual estimates. The integrated product considers terrain, geological, and hydrometeorological factors, groundwater leakage by base flow, and groundwater withdrawal from shallow and deep aquifers. It also compensates for inadequate spatial continuity of lateral and vertical variations in groundwater resources estimated from in situ GWL data due to the limited spatial density of piezometric wells in some places (in parts of U.P., M.P., and Rajasthan states).

This study elaborated the procedure for estimating GWS change individually from satellite-based GRACE gravity observation, hydrological modelling, and in situ GWL fluctuation data. Finally, a novel approach was presented to produce an improved regional estimate of GWS change by integrating satellite-based GRACE gravity observation, hydrological modelling of hydrometeorological data, and in situ GWL fluctuation data. By combining GRACE gravity observation and hydrological model-based TWS change products using multi-linear regression coefficients as weight coefficients, the approach integrates the advantages of individual datasets and overcomes their limitations. The integrated product on 0.25° × 0.25° spatial grid overcomes the limitations of GRACE observation (spatial resolution, local geology, terrain, and hydrometeorological factors), hydrological model (groundwater withdrawal conditions), and conventional method (lack of spatial continuity at places and cumbersome, labour-intensive exercise). The integrated GWS change product is better than any of the three individual products in terms of the operability and quality of the product. In Pearson product-moment correlation analysis, the hydrological model and GRACE observation-based individual GWS change products and the integrated GWS change product show positive correlation in ∼59, ∼69, and ∼73% of the study area with in situ GWS change. Overall, the integrated GWS change product is substantially improved. It may be noted that due to inadequate spatial density of piezometric data in parts of U.P., M.P., and Rajasthan states, the lateral and vertical variations are not adequately accounted for achieving a reliable spatial continuity in GWS change which led to low/negative correlation coefficients at such places. The integrated product overcomes such problems of inadequate spatial continuity. Altogether, the approach presents an operational but robust method that is immensely useful to formulate a strategy for ensuring groundwater sustainability and minimising adverse impacts of groundwater depletion.

The authors gratefully acknowledge the Indian Institute of Remote Sensing (Indian Space Research Organization) for providing institutional support for conducting the study. The authors are also grateful to the Central Ground Water Board, India for providing groundwater level data of the study area.

Most of the data are publicly available. For any specific data or clarification, readers should contact the corresponding author.

The authors declare there is no conflict.