Investigating Spatio-temporal behaviour for Groundwater in North-West India: A Deep Learning Approach

doi:10.21203/rs.3.rs-3638857/v1

Sustainability and conservation of natural resources require amalgamation of the novel approaches with the current usages and weather conditions. Groundwater resources is one of the natural resources which is varying annually and requires regular attention for the prediction. Our proposed deep learning (DL) approach, namely Convolution Neural Network-Long short term memory (ConvLSTM) has been implemented for the groundwater level (GWL) prediction model. The model is designed based on U-Net framework with up-sampling and down-sampling modules and also induces non-linearity using the ReLU activation function. Each module in the LSTM unit is responsible for pattern recognition based on the temporal information of GWL. The assessment of the groundwater in North-West India (NWI) has been carried out using several fundamental factors such as precipitation, soil moisture, evapotranspiration and satellite-based groundwater storage. In addition, in-situ groundwater has been used to get groundwater fluctuation scenarios (i.e., categorised into four cycles PrePre, PrePost, PostPre, and PostPost) w.r.t monsoon season to understand the difference (Δh) in GWL. The proposed model has been tested with other DL frameworks such as; Artificial neural network (ANN) and Convolution neural network (CNN). The model has been trained using the stochastic gradient method to optimise the internal parameter and validated using several geo-locations information of NWI, where ConvLSTM outperformed compared to the benchmark method. The proposed model has shown consistent least error in terms of root mean square root (RMSE) and mean square error (MAE) for the year 2014-17 with an overall score of 0.0957 and 0.0520, respectively.

Groundwater level prediction

deep learning

convolution neural network

long short term memory

spatio-temporal modelling

The presence of groundwater on the Earth is due to the interaction of the hydrological, climatic, physiographical, geological, and ecological factors. Due to the effortless accessibility of freshwater resources around the globe, overexploitation of aquifers leads to a higher degree of groundwater (GW) level fluctuation. GW manages 50% of the global irrigation water demand and 97% of the freshwater demand (Mallick et al. 2021). The largest freshwater resource available in the form of groundwater plays a vital role in human sustenance and global food security via irrigated agriculture wherever the surface water resource is inadequate (Aeschbach-Hertig & Gleeson 2012). The Indian subcontinent (ISC) hosts around 24% of the world’s population, however, it comprises merely ∼4% of the land area with more than 30% of the global irrigated land (Bhanja et al. 2018). In total, 89% of the groundwater in India is consumed for irrigation, the rest 9% and 2% are used for domestic and industrial use respectively (Kumar et al. 2021). Due to prompt changes in economic/lifestyle and growth in urbanised areas, ISC has witnessed a sharp escalation in freshwater demand that leads to an acute shortage of groundwater resources available for drinking and other purposes (Bhanja et al. 2018). Moreover, the densely populated area around the globe is facing a rapid loss of GW resources (Rodell et al. 2009) which may soon lead to socio-economic stress to billions of people. The large population growth (or population density) in ISC renders the increasing demand of groundwater. As a result, people are indirectly extracting more and more groundwater. Thus, the GW resources are getting projected to a disquieting state all over the country (CGWB 2014).

Groundwater totally feeds on irrigated land in India, which has drastically expanded from 30% to ∼50% within the time span of 1960 to 1995 (Scanlon et al. 2010), and extending gradually in northwest India (NWI). The highest volume of groundwater is extracted in the agriculture sector in Punjab (lies in NW India), also known as agriculture-based economy (Sidhu et al. 2021). The world’s second-most extracted groundwater resources lie under the Indus drainage basin lies beneath Punjab (Richey et al. 2015). Similarly, National Aeronautics and Space Administration’s (NASA’s) based Gravity Recovery and Climate Experiment (GRACE) observations have also shown a temporal pattern of declining drainage basins of NW Indian aquifers, which trigger the high risk region for GW depletion. The observation primarily specifies that the NW Indian states of Delhi-national capital region (NCR), Punjab, Rajasthan, and Haryana, are depleting at a rate of 109 km³/year (Rodell et al. 2009). GRACE mission observes the spatio-temporal variation in the Earth’s gravity field that assesses the Earth’s total integral mass variations, which can be used in hydrological studies more effectively at regional-scale.

In modelling perspective, previous researches have already proved that the deep learning approach has a broader prospect in the field of groundwater level studies and its management (Huang et al. 2019). Various models formulated on machine learning (Alizadeh et al. 2018) and deep learning (Ye et al. 2019) approach are confirmed to be used in numerous fields, including climate change forecasting (Scher 2018), predicting groundwater potential (Kumar et al. 2021), projection of forest cover (Ye et al., 2019), typhoon and flood forecasting (Jiang et al., 2018). The standard regression methods have been adopted in GW studies and improve the performance in groundwater dynamics prediction (Huang et al. 2019).

In this study, a deep learning approach has been used for the first time to predict the GWL to learn the spatial as well as the temporal behaviour of the NWI. In the deep learning approach, CNN-LSTM model is constructed from the scratch, where the seasonal cycle has been used to render the temporal behaviour on the geospatial (or image) datasets. For the dataset creation, in-situ groundwater level data and hydrological parameters are included that actively influence the GW scenario in NWI. The proposed method is tested with standard models such as ANN and CNN, for qualitative and quantitative comparisons. The deep learning architecture and data set creation with detailed explanation of hydrological parameters, are discussed in the section 4.2 and 3 respectively.

NW part of ISC has been selected for this study, as it is primarily among the most populated and highly GW extracted regions around the globe (Rodell et al. 2009; Bhanja et al. 2018; Pranjal et al. 2021). The NWI has an area of 11,70,553 km² (shown in Fig. 1), covering the Delhi-National Capital Region (NCR) and other parts of the states (namely, Punjab, Haryana, Uttar Pradesh, Madhya Pradesh, Rajasthan, Gujarat, and Uttarakhand). These parts of the states lie within the five major drainage basins i.e.; Chambal, Luni-Ghaggar, Narmada, Upper Ganga and Yamuna. This region receives moderate to low rainfall during the monsoon months of June-August, while most of the region faces aridness for the rest of the year.

The study area has diverse geomorphological processes that lead to the advancement of different landforms, such as structural landforms, fluvial landforms, volcanic landforms, and aeolian landforms (Chatterjee et al. 2020). Also due to physiographic conditions and diverse climate, it has varied vegetation types such as, cultivation of crops primarily wheat, rice, sugarcane, besides swamp vegetation occur tropical dry deciduous forest, tropical semi-evergreen forest, subtropical pine forest and littoral tropical moist deciduous forest, and tropical thorn forest in the region.

3.1 Groundwater storage estimates using in-situ GW level (GWL) data

Central Ground Water Board (CGWB; a Govt. of India body) maintains a strong network for regularly monitoring GW wells across the country, which was started in 1969 ((CGWB) 2014). There are around 22,730 monitoring wells operating across India, in which 16,375 (~ 72%) are dugwells (as of March, 2017; (CGWB) 2017). However, there are around 8,219 wells available in NWI, out of which 3,176 (~ 39%) are tubewells. CGWB collects in-situ GWL data for four seasons (namely, late post-monsoon [January], pre-monsoon [May], monsoon [August] and early post-monsoon [November]), along with all the well locations across NW India (shown in Fig. 2). Groundwater wells offer numerous benefits, as it’s sustainable use requires careful management to prevent over-extraction, drought resilience, GW quality control, contamination, and other environmental concerns. Local regulations, monitoring programs, are some of the essential components of maintaining a strong network of groundwater wells.

For this research, GWL data between January 2005 and November 2017 have been considered, mostly from unconfined aquifers (~ 87%) of the study area. The GWL fluctuation has been calculated based on their respective seasonal change for every consecutive season. The annual GWL changes have been calculated using pre- and post-monsoon data collected from CGWB. The difference in GWL height (Δh) between pre- to post- monsoon, pre- to pre- monsoon, post- to pre- monsoon and post- to post- monsoon represented as pre-post (P_rP_o), pre-pre (P_rP_r), post-pre (P_oP_r), and post-post (P_oP_o) monsoon head change respectively for the years 2005–2017 (illustrated in Fig. 3). The seasonal variability (or head change) of GWL fluctuation data have been spatially interpolated (integrating the statistical tools, namely, Kriging and IDW). The interpolated GWL Δh value provided the spatial picture of the seasonal GWL change. Spatial layers generated using Kriging and IDW tools were combined to a form a single interpolated layer by applying weighted average. Finally, the layers were downscaled to homogeneous 44 × 64 gridded datasets for the study area. The temporal continuity of the data was taken into account to finalise the data for further processing with at least three seasons for the respective year.

The trend observed out-of-the data collected from observation wells (i.e., water-table and piezometric well) reveals that the short-term or local rainfall can only improve the soil moisture of the area or serve farmers. The rainfall plays a major role in influencing the groundwater fluctuation in this region. The GWL change for various season in relation with the annual rainfall suggest groundwater extraction/recharge all through the year (Kadiyan et al. 2021). However, only ample rainfall for a longer period of time can help in groundwater recharge. As a result, all the modelled recharge takes place in the monsoon (Kharif) season (i.e, P_rP_o), whereas most of the discharge occurs during the non-monsoon (Rabi) seasons.

3.2 Satellite-based groundwater storage estimates

In a collaborative mission instigated by NASA and DRL in 2002, the launch of GRACE has been continuously observing the variations in Earth’s gravity field with exceptional accuracy (Tapley et al. 2004). GRACE uses the variation in the mass of the Earth’s fluidic layers (particularly air and water), to recognize the temporal changes in the Earth’s gravity field. Various algorithms are being adopted to estimate the terrestrial water storage (TWS) from GRACE gravity anomaly. These algorithms also remove the inferences present due to solid earth systems, oceanic waves and atmospheric inundations (Flechtner 2007). Many studies have been proven globally that GRACE satellite mission is extensively used in computing groundwater storage change (GWSC), even at places where there is inaccessibility or limitations in collecting data (Bhanja & Mukherjee 2019). Besides, it has also been reported that satellite-based GRACE mission deliberates the GWSC (either at basin/aquifer scale or continental scale) (Rodell et al. 2009; Reager & Famiglietti 2013; Richey et al. 2015). The observed GRACE gravity data deliberates the TWS change over the glacier as well as non-glacier area (Chen et al., 2016). The study assessed the quantitative measures of signal attenuation and leakage error resulting from post-processing filters. This evaluation was conducted in relation to simulated monthly TWS anomalies derived from data generated by the GLDAS-Noah land surface model (Chatterjee et al. 2019). Additionally, variations in GRACE gravity are attributed to alterations in mass distribution within both the solid Earth and the Earth's surface. Key contributors to gravity variation from solid Earth processes encompass phenomena like glacial isostatic adjustment (GIA) and earthquakes, particularly at a local scale. Moreover, as the study area mostly covers the non-glacier area i.e., land area. Consequently, by eliminating the hydrological parameters (like; surface water storage, soil moisture storage, evapotranspiration) effecting the mass redistribution over the area quantifies the GWSC.

The original product of GRACE (or spherical harmonic products) was basically intended for the regional study with spatial resolution of 3°×3° (Tapley et al. 2004). Several algorithms have made it possible to achieve 1°×1° of spatial resolutions (Watkins et al. 2015), which help users to work on local-scale to understand the GWS change impacted from anthropogenic activities (Bhanja & Mukherjee 2019). The improved the spatial resolution of GRACE products up to 0.25°, is obtained through the statistical regression downscaling method (Ning et al. 2014). Subsequently the mascon solution of GRACE provides data at a spatial resolution of 0.5°×0.5° and even 0.25°×0.25°. At each grid point GRACE mass rate is estimated using unweighted least square to fit of a linear trend, plus annual and semi-annual sinusoids to GRACE-derived TWS time series. Mascon data offers benefits over spherical harmonics, such as (1) minimizing signal loss from land to ocean, resulting in stronger signals, and (2) incorporating geophysical data constraints in processing with minimal need for additional empirical adjustments. This makes it more accessible for non-geodetic users. Comparing the outcomes of this product assessment should assist hydrologists in making more informed choices regarding the most appropriate GRACE solutions for hydrological applications (Scanlon et al. 2016).

In this research, 133 newly released RL05 mascon monthly gridded (0.5°×0.5°) Jet Propulsion Laboratory (JPL) solution has been used for analysis from 2005–2017. These GRACE mascon products have been derived through liquid water equivalent thickness (LWET, which is also recognized as terrestrial water storage (TWS)), archived from the NASA JPL. Total terrestrial water storage anomaly (TWSA) is relative to a long-term mean TWS (e.g., the average of 2004.00–2009.99 used in this study), which is expressed as an equivalent water thickness (EWT; unit: mm) (Wahr et al. 1998; Yang et al. 2020). Several reports in India suggest that GRACE mascon products have resulted to be a better solution than GRACE spherical harmonic products for estimating groundwater storage using in-situ groundwater storage measurement. The difference in this mascon product is that the entire globe is divided into 30 spherical mass concentration units with equal area (Watkins et al. 2015) and several techniques applied on the product to estimate TWS. To make the land product of TWS further representative, the provided scale-factor with the data has been multiplied with the TWS solution. Finally, to achieve the common spatial resolution, post-processed data have been downscaled to 44 × 64 gridded data for the study.

3.3 Precipitation

Tropical Rainfall Measuring Mission (TRMM) mission is a collaborative venture amongst NASA and JAXA, with several other space agencies to provide Earth’s frequent precipitation measurements. The TRMM mission is a part of NASA’s Earth systematic missions program that provides full global coverage using a series of satellite constellations. It also provides TRMM multi-satellite precipitation analysis (TMPA) data used in this study. The TMPA comes up to provide optimal rainfall estimation around the globe (Huffman et al. 2010). Moreover, the Global Precipitation Measurement (GPM) mission is also a continuation of the TRMM mission. The spatial coverage of both the products are provided with a gridded (0.25°×0.25°) solution and having a cycle-based product (i.e., 3-hourly, daily, monthly and yearly) across the globe. Accurate monitoring of precipitation, both in terms of its quantity and distribution over time and space, remains crucial for assessing its adequacy or scarcity. These datasets enable the tracking of precipitation patterns and its broader impact on hydro-ecological systems and the management of water resources. Globally, reliable precipitation information is vital for fostering sustainable socio-economic development, particularly in regions with lower rainfall levels and in some developing nations where rainfall is pivotal for economic activities like, groundwater recharge and agriculture. Given the highly variable nature of precipitation in terms of both location and time, establishing a robust network of gauge installations is imperative for effectively monitoring the occurrence and intensity of rainfall. However, in numerous countries, precipitation data has traditionally been gathered from a sparsely distributed set of gauge stations (Nwachukwu et al. 2020). These products include the integration of several satellite products with the rain gauge data products collected from Global Precipitation Climatology Centre (GPCC) archive (Huffman et al. 2007; Huffman et al. 2010).

Precipitation is the primary element that influence recharging the groundwater reservoirs. Therefore, in that respect precipitation has been taken into account for the prediction of groundwater level. Its primary objective is to enhance the comprehension of precipitation distribution and variability within the tropics, a crucial component of the current climate system's water cycle. By encompassing the tropical and sub-tropical regions, TRMM furnished vital data on rainfall patterns and the accompanying release of heat, which drives the global atmospheric circulation governing both weather and climate. In conjunction with other satellites within NASA's Earth observing missions, TRMM utilized various space-based instruments to provide essential insights into the interactions among water vapour, clouds, and precipitation. These interactions play a central role in the regulation of Earth's climate. Precipitation holds a pivotal position in the GWL prediction model, exerting influence on both the groundwater recharge process and subsequent variations in groundwater levels. Its roles encompass pivotal aspects including groundwater recharge, temporal fluctuations in GWL, seasonal variations, hydrological connectivity, groundwater-surface water interactions, and groundwater management. Alterations in rainfall patterns can modify the interconnectivity between surface water bodies like rivers and lakes and the underlying groundwater systems. Heightened rainfall may trigger amplified surface runoff, potentially recharging neighbouring groundwater systems or fostering interactions between groundwater and surface water. As CGWB collects the in-situ GW fluctuation data four times a year, therefore, the rainfall products are arranged based on the four time-periods (Bhanja & Mukherjee 2019), and the data is further post-processed to achieve 44 × 64 gridded data for analysis.

3.4 Soil moisture change

Remote sensing-based Global Land Data Assimilation System (GLDAS) satellite mission has been used to prepare soil moisture using the Noah land surface model (Rodell & KH 2007). By assimilating observed soil moisture data and other relevant information, GLDAS refines its estimates of soil moisture conditions at various depths within the soil profile. These estimates are expressed as soil moisture equivalents. The monthly product has a spatial extent of 0.25°×0.25° and provides the soil moisture data in 4 different layers with depth: 0–10 cm, 10–40 cm, 40–100 cm, and 100–200 cm. Henceforth, the monthly soil moisture data has been prepared by combining the four layers of each month and further rescaled to 44 × 64 gridded data. As soil moisture is a state function, therefore, monthly soil moisture change was assessed based on the long-term average (Chatterjee et al. 2020).

3.5 Evapotranspiration

GLDAS follows the assistance for land surface modelling activities (ALMA) definitions and conventions for the land surface model output variables outlined in the form of various products. The Land Data Assimilation System (LDAS) mission aims to provide the land surface states by combining ground and satellite-based data products, using the advanced data assimilation technique and land surface modelling. In recent progress on satellite GRACE gravity measurements which has facilitated the quantification of a sea level rise up to 0.71 ± 0.20 mm/year from 2002 to 2014 (Zan et al. 2019). This deceleration is attributed to the retention of water in land due to climate variability (Wada et al., 2017). This assessment is pivotal in reducing uncertainties in estimates of terrestrial water budgets. Moreover, accounting for changes in water flux induced by ET significantly enhance the accuracy of uncertainty estimates in regional or basin-scale hydrological modeling (Melesse et al., 2009). Among various other GLDAS products, evapotranspiration (ET) data is an average monthly 0.25°×0.25° solution (provided in kg¹ m^-2 s^-1), which is further processed to common 44 × 64 gridded data products.

4.1 Data processing

There are basically three sets of approach that can be considered for GWL fluctuation prediction using spatio-temporal/time-series models, such as; taking GWL as an input component, taking rainfall as an input component for GWL prediction, and direct and recursive prediction strategies using hydrological components (Yoon et al., 2016). In the GWL perspective, several hydrological parameters have been considered for the study that actively influences the GW scenario. For instance, GLDAS simulation products have been used to estimate the soil moisture (SM) equivalent and surface water (SW) equivalent from NASA’s archive (Rodell et al. 2004). Surface water equivalents represent the amount of water stored in surface water bodies (e.g., rivers, lakes) within the model's grid cells. This information is also influenced by precipitation, and other hydrological processes simulated by the LSMs. Integrating GLDAS soil moisture data into hydrological models can help in understanding how soil moisture influences processes like runoff, evapotranspiration, and groundwater recharge. In India, particularly the northwestern subcontinents do not fall under the snowfall zone (except the minor sections of the Himalayas) and are negligible among other hydrological factors. Therefore, snow water equivalent has not been considered for this study (Bhanja & Mukherjee 2019). Surface runoff has been derived from GLDAS, which acts as an alternative for SW equivalent (as measurements are not available for NW India) (Asoka et al. 2017). Noah-based land surface model (LSM) has been implemented for the calculation of SM and SW (Koren et al. 1999). Noah LSM models the movement of water within the land surface, including precipitation, infiltration into the soil, runoff generation, groundwater recharge, and evapotranspiration. The model calculates soil moisture at multiple soil layers, representing the vertical profile of soil properties. It considers processes such as infiltration, percolation, and capillary rise to estimate soil moisture content. In addition to soil moisture, Noah LSM accounts for surface water processes. This includes the calculation of surface runoff, which is the water that flows over the land surface in response to excess precipitation and may contribute to streamflow.

Time series anomalies (Δ) of hydrological components have been used after filtering the all-time mean from the individual datasets. Space-borne GRACE gravity anomalies are being used worldwide for calculating TWS change, especially for the regions with limited or very few observations well. Seasonality in hydrological modeling has been refers to the variations in hydrological processes (i.e., precipitation, et, and rainfall) that occur in response to changing seasons throughout the year. It is a crucial aspect to consider in hydrological modeling because it directly affects the availability of water resources and can have significant implications for various applications like water resource management etc.

All these aforementioned parameters have been collected from various sources (see section 3), containing the cyclic data for different years (i.e., 2005–2017). To maintain the temporal continuity of the datasets, all the monthly data were collated for the mentioned duration and rescaled to similar spatial resolution. In general, the spatio-temporal model requires the homogeneous size and rectangular shape of the input feature for training and testing. Therefore, to begin with the processing steps, all the images were initially masked from the raw data using the shape-file of the NWI region using ArcGIS 10.5 software. Thereafter, all the images were resized to a quadrangular shape of 44×64 pixels. The rectangular input has several benefit such as it preserves the spatial arrangement of feature, several operations like pooling, down sampling will relatively easy to operate and memory efficient. In addition, the setup demonstrates the design to predict the difference in the GWL based on the different factor and cycle (w.r.t the monsoon season). Input features of the images were stacked w.r.t cycle (see Fig. 4), and further restacked w.r.t factor. Consequently, the resultant matrix (or features) becomes 44 × 64 × k × z, where k and z are the cycle and factor affecting the end result, respectively.

After arranging the data, the internal (or brightness) values need to be calibrated. The resultant matrices of each selected hydrological layer were used for further processing to best fit the deep learning model. The choice and amount of dataset are significant aspects of any deep learning method. The size of the training dataset is limited, therefore, the diversity of the datasets has been increased by several transformations using data augmentation, including random flipping, colour enhancement rotation, zooming, etc. The processing steps are summarised in Fig. 4.

In the deep learning approach, the internal parameters (such as weights and biases) ranges between [0,1], and input features can vary differently (see Table 1). These mismatched ranges between input feature and internal parameter cause abnormal gradient behaviours during the training and lead toward poor model performance (Soni et al. 2020). In order to calibrate the values, min-max normalisation has been used with the following function:

$$\underset{\_}{x}=\frac{x-{x}_{min}}{{x}_{max}-{x}_{min}}$$

1

Where, ${x}_{min}$ and ${x}_{max}$ are the maximum and minimum value within the same feature vector, x is the pixel value in the feature map and $\underset{\_}{x}$ is the normalised pixel value.

4.2 Model parameter

In the previous work, spatial and temporal information was handled differently by the conventional approach. Even in the deep learning model, it has been addressed separately by convolution neural networks and LSTM. Designing a model with separate modalities could not incorporate (or take advantage) the spatial and temporal properties while learning the parameter in the model. The CNN architecture is neural network especially deal with image information and extract feature from the input image by applying the several operations (or layer) and develop the correlation between spatial feature in the low-level feature (i.e edges). The layers include convolution layer, pooling layer, etc. to enhance the quality of extracted feature, however, it doesn’t provide the adequate operation to extract the information from an image in case of temporal changes. The temporal information has been handled by different DL architecture such as LSTM. The LSTM is a type of recurrent neural network (RNN) architecture designed to overcome the limitations of traditional RNNs in capturing long-range dependencies in sequential (or temporal) data. At its core, LSTM uses hidden cores (or layer) to remember the sequence and it pattern. Therefore, we have addressed these issues for the prediction of groundwater level using the composite module CNN and LSTM, known as ConvLSTM (Xingjian et al. 2015). There are some challenges and difficulties associated with their integration such as vanishing gradient, overfitting, etc, which were carefully addressed during the parameter setting and pre-processing (see section 4.1).

ConvLSTM extracts the feature information using the convolution operation in input to state and state to state transition (shown in Fig. 5). In the preliminary operation, the input feature should comprise of 5 dimensions: sample size, row, column, channel, timestep (i.e. S×44 × 64 × k × z, where S is the sample size). ConvLSTM layer produces the output of a similar size to the input image. The internal structure of ConvLSTM is shown in Fig. 5, with interoperation is expressed below:

${i}_{t}=\sigma ({W}_{xi}{X}_{t}+{W}_{hi}{H}_{t-1}+{W}_{ci}o{C}_{t-1}+{b}_{i})$	(2)
${f}_{t}=\sigma ({W}_{xf}{X}_{t}+{W}_{hf}{H}_{t-1}+{W}_{cf}o{C}_{t-1}+{b}_{f})$	(3)
${C}_{t}=tanhtanh ( {W}_{xc}{W}_{t}+{W}_{hc}{H}_{t-1}+{b}_{c})$	(4)
${C}_{t}={f}_{t}o{C}_{t-1}+{i}_{t}o{C}_{t}$	(5)
${o}_{t}=\sigma ({W}_{x0}{X}_{t}+{W}_{ho}{H}_{t-1}+{W}_{co}o{C}_{t}+{b}_{o})$	(6)
${H}_{t}={o}_{t}\left({C}_{t}\right)$	(7)

Where X_t denotes the input feature of the current cell, C_t-1 and H_t-1 are the state of output feature of the previous cell. i is input gate, f is forget gate, o is output gate, c is cell state, h is hidden state, W is weight, ⌠ is activation function (i.e. Relu), o means Hadamard product and * is the convolution operation.

The network architecture is designed in encoder-decoder fashion (Sari et al. 2020). It consists of two major parts: 1. Encoder part; 2. Decoder part. The encoder part consists of three sets of ConvLSTM, batch normalisation, and maxpooling operation. This part extracts the spatial feature and adjusts the parameter based on the temporal feature. In each set of operations, the size of the feature is downsampled to enhance the field of view. The decoder part is similar to the encoder part, where upsampling operation is used instead of maxpooling operation. This part re-generate the feature using upsampling operation and spatial feature extracted from the other part. This repetitive upsampling recovers the spatial properties (such as, height and width) similar to the input feature.

The architecture was trained from the scratch using the stochastic gradient descent (SGD) optimizer with learning rate of 0.001 and momentum of 0.9. The iteration contains only 10 samples (and 5x samples using augmentation) with different time-series. The loss function is based on the root mean square error (RMSE) and mean square error (MAE), which evaluate the training and testing performance of the network (see section 4.3). To avoid the overfitting during the training, validation (equivalent to testing) set has also been used to consistently monitor the model performance. In addition, early stopping function with patience of 5 was used, which eventually check five continuous iterations and its improvement.

4.3 Evaluation metrics

This section explains the metric that has been used to perform quantitative analysis between the models. All the models follow the regression approach, therefore, root mean square error and mean absolute error has been used.

Root means square error is the standard deviation of the prediction errors. The error is defined as the distance between the actual and prediction values. The mathematical expression is shown in Eq. 8:

$$RMSE=\sqrt{\frac{{\sum }_{i=0}^{n}{({x}_{i}-{\widehat{x}}_{i})}^{2}}{N}}$$

8

Means Absolute Error is an arithmetic means between absolute error produced by the model. The absolute error is the difference between actual and prediction values. The mathematical expression is shown in Eq. 9.

$$MAE=\sqrt{\frac{{\sum }_{i=0}^{n}\left|{x}_{i}-{\widehat{x}}_{i}\right|}{N}}$$

9

Where, ${x}_{i}$ and $\widehat{{x}_{i}}$ are the actual and prediction observation of the data point i, respectively. N denotes total number of data points.

Seasonality plays an important role in deriving values from the hydrological parameters. Table 1 manifests the seasonal variation of different hydrological parameters analysed for the NWI. It can be perceive that the range of P_oP_r cycle is inclined toward the negative region, as a result scarcity of water level has been observed during this season. In addition, the other parameters also indicate the reduced values for all the parameters, which appears to be the decline in terrestrial water storage in the region is highest for this season. The monsoon season abruptly affects the GWL and mostly increases the water level. Such a scenario has been observed during P_rP_o season (in Table 1), where the rainfall directly affects the GWL. In the inverse scenario i.e., PoPr, the rainfall observed during this season is relatively low and directly reflects the changes in the GWL. The annual cycles i.e., P_oP_o and P_rP_r indicate the variations in the monthly average values of the parameters. The values are having high fluctuations due to high spatial distributions. These fluctuations in data (in Table 1) are managed and the predictions have been done by applying deep learning models. The values within the parameter were rescaled using the min-max normalisation, which also preserves the data distribution properties. It is noted that the influence of the observed GWL values in neural networks modelling were essential, including the range of the GRACE, soil moisture, evapotranspiration, rainfall, and surface water, therefore the values were normalised.

Table 1

Seasonal variation in the values of various hydrological parameters during the study period (2005–2017) in NWI (values in cm)
Parameter/ Cycle	GWL		GRACE		SM		ET		Rainfall		SW
Parameter/ Cycle	Min	Max	Min	Max	Min	Max	Min	Max	Min	Max	Min	Max
P_rP_o	-16.4	275.4	-37.3	30.1	34.9	73.9	0.3	15.9	0.01	71.5	0	59.2
P_rP_r	-260.2	293.0	-37.9	30.1	34.8	73.9	0.2	15.9	0.03	71.5	0	59.2
P_oP_r	-393.2	-7.1	-35.7	21.5	35.4	58.2	0.3	6.3	0.03	12.8	0	8.7
P_oP_o	-212.1	193.3	-37.9	30.1	34.8	73.9	0.3	15.9	0.03	71.5	0	59.2

A complex mapping function has been defined such as; ConvLSTM, CNN and ANN for the GWL modelling. The model uses datasets between the period 2005–2013 for the training and 2014–2017 for the testing. The overall performance of the prediction models are found to be acceptable. However, it has been observed that the predicted GWL data between 2014-17, the ConvLSTM based-prediction shows less RMSE and MAE as compared to the CNN and ANN models (Table 3). In the similar manner, Fig. 6 clearly represents the visual comparison of ConvLSTM, CNN and ANN for different cycles, where ConvLSTM shows homogeneity in the predicted results to the ground reference compared to CNN and ANN based-predictions. The spatial variation in the values (change in GWL fluctuation i.e., Δh) among the models and their accuracies have been represented using various box plots and graphs.

As Fig. 6 indicates the visual comparison of the prediction model with the ground references for the three consecutive years. In comparison, it is observed in the various cases that ConvLSTM based-prediction gives almost similar results as compared to the ground reference data. Whereas, in Fig. 7, it demonstrates the statistical representation of the distribution of datapoints from the modelled prediction using box plot. In this plot, the horizontal red lines show the central lines of the reference data for the comparison. In several scenarios, most of datapoints of ConvLSTM is concentrated around the reference lines. Thus, ConvLSTM has displayed similar distribution with relatively less error as compared to CNN and ANN (see Table 3).

A box plot, also known as a whisker plot, is a graphical representation used here to display the distribution, central tendency, and variability of the datapoints (shown in Table 2). It is particularly useful for visualizing summary statistics and identifying potential outliers. A boxplot graph has been used to indicate the insight values of the distributed datasets spatially spread throughout the region. It is the convenient way to visually display the data distribution throughout the quartiles. The distribution of datapoints have been shown in Fig. 7, also a random set of points were selected from the dataset to perform the comparison between the ConvLSTM and ground reference (shown in Fig. 8). The red and blue line denotes the fluctuation of ground reference and ConvLSTM based predictions values, respectively. The graph interprets that the values of red and blue lines mostly overlap with each other with relatively fewer disagreement points.

Table 2

Pixel wise error statistics w.r.t ground reference for the year 2014-17.
		2014-15				2015-16				2016-17
		P_rP_o	P_rP_r	P_oP_r	P_oP_o	P_rP_o	P_rP_r	P_oP_r	P_oP_o	P_rP_o	P_rP_r	P_oP_r	P_oP_o
Ground Reference	Median	-0.139	-0.126	-0.200	-	0.152	-0.111	-0.105	-0.051	-0.041	0.263	-0.676	-0.051
Ground Reference	SD	0.222	0.309	0.319	-	0.216	0.270	0.319	0.177	0.226	0.204	0.356	0.177
ANN	Median	-0.47	-0.157	0.417	-	0.123	0.336	-0.163	-0.054	-0.033	0.135	-0.174	0.392
ANN	SD	0.103	0.165	0.159	-	0.135	0.194	0.466	0.132	0.135	0.197	0.199	0.109
CNN	Median	-0.45	-0.104	0.296	-	0.218	0.247	-1.351	-0.115	-0.217	0.505	-0.256	0.164
CNN	SD	0.194	0.206	0.262	-	0.204	0.143	0.631	0.189	0.170	0.158	0.356	0.174
LSTM	Median	-0.95	-0.135	-0.150	-	0.234	-0.129	-0.154	-0.133	-0.047	0.211	-`0.541	0.270
LSTM	SD	0.190	0.210	0.315	-	0.215	0.199	0.471	0.164	0.194	0.224	0.110	0.089

Consequently, we interpreted from the results that the comparison among the observed and predicted models tested for the data sets of groundwater using ANN, CNN, and ConvLSTM model found to be well-fitted in case of ConvLSTM. Further, it was attempted to compare the simulation results of ANN, CNN and ConvLSTM models using various error evaluation metrics, to estimate the accuracy of the models. It has also been observed that by increasing the number of image layers, accuracy improves. Likewise, to evaluate the accuracy via increasing the input parameters of the dataset’s size, the accuracy gets ameliorated (the results have been shown in Fig. 6). The RMSE and MAE statistics of the three identified deep learning models for the predicted period of 3 years (2014–2017) are summarised below (in Table 3). It has also been observed that all the models have performed well, however, the error in ConvLSTM predictions are relatively less than CNN and ANN. In general, the ANN and CNN have, on average, comparable prediction performance despite their basically dissimilar learning scheme. Based on the deep learning models, hydrological datasets processed for the NWI and their accuracies were evaluated (see Table 3).

Table 3

Relationship among the predicted models and regression methods with variation in the values of various seasons during 2014-17 in NW India
Regression methods / Model(Cycle)		2014-15		2015-16		2016-17		Overall
Regression methods / Model(Cycle)		MAE	RMSE	MAE	RMSE	MAE	RMSE	MAE	RMSE
ANN	PrPr	0.1047	0.1626	0.1052	0.1638	0.1063	0.1663	0.1063	0.1662
	PrPo	0.0781	0.1268	0.0783	0.1269	0.0782	0.1271	0.0781	0.1268
	PoPr	0.1052	0.1604	0.1050	0.1604	0.1054	0.1304	0.1053	0.1640
	PoPo	0.0979	0.1352	0.0980	0.1648	0.0991	0.1672	0.0990	0.1652
CNN	PrPr	0.0805	0.1341	0.0803	0.1335	0.0800	0.1329	0.0804	0.1336
	PrPo	0.0672	0.1159	0.0673	0.1161	0.0674	0.1163	0.0677	0.1189
	PoPr	0.0850	0.1223	0.0892	0.1235	0.0894	0.1233	0.0890	0.1231
	PoPo	0.0962	0.1089	0.0662	0.1171	0.0685	0.1250	0.0679	0.0119
RNN-CNN	PrPr	0.0600	0.1093	0.0602	0.1091	0.0604	0.1103	0.0604	0.1099
	PrPo	0.0596	0.1083	0.0573	0.1081	0.0570	0.1082	0.0571	0.1082
	PoPr	0.0631	0.1001	0.0628	0.1015	0.0625	0.1059	0.0626	0.1005
	PoPo	0.0756	0.0961	0.0522	0.0952	0.1063	0.0966	0.0520	0.0957

The groundwater resource is declining significantly at an exponential rate around the globe, especially some of the major cities in India are already out of groundwater and several are in the queue. The groundwater is greatly influenced by several factors such as rainfall, soil moisture etc. In this study, information of the driving factor of groundwater has been extracted by different sources of remote sensing images to analyse the spatio-temporal behaviours of GWL in NW India. The deep learning model has been devised for the first time for GWL prediction. The deep learning model such as ANN, CNN, and ConvLSTM have been used, where ConvLSTM has outperformed among other models. The predicted values based on ConvLSTM for GW fluctuation (for the period 2014–2017), found to be almost equal to the observed/ in-situ groundwater fluctuation change. Due to drastic improvement of high-resolution space-borne and ground observation along with observation wells around the NWI region, it has become vital to design a model to forecast the GWL and unearth its pattern. Therefore, deep learning methods have been proposed with different functionality to train and test the pattern observed from the actual result and its robustness.

The result indicates the capability of neural networks in modelling GWL using GRACE, rainfall, ET, and soil moisture data with the trained value of the groundwater level. The deep learning techniques manifest the models to have performed well while incorporating the different cycle information. The proposed ConvLSTM-based model in this study demonstrated very promising results for improving water supply planning in semiarid areas. The results were analysed by referring to the correlation coefficient analysis, coefficient of determination, i.e., MAE and RMSE. It was also noted that the influence of the temporal information for GWL was very significant in neural networks modelling along with the value of the GRACE, rainfall, ET, and soil moisture data. A substantial advantage of the ConvLSTM model is that it has provided good predictions with the limitations in the GW level. As our approach towards prediction of GW change using deep learning-based algorithms provided encouraging results across NW India. These are the fundamental reasons behind this study to apply both the spatial and time series modelling methods using deep learning techniques.

Furthermore, the performance of the model can likely be improved by taking into consideration several additional information, such as discharge rates, water withdrawal, demographic factors etc. More interestingly, the model is less costly and can be used for regional studies for GWL fluctuation prediction. However, the prediction results obtained for NW India are relatively good depending upon the choice of model. We understand that the result of this study would help the scientific society and policymakers for decision making for the sustainable management of groundwater resources.

Acknowledgement

The authors would like to acknowledge IIT (ISM) Dhanbad and IIRS-ISRO, Dehradun for providing a working facility to execute this research. We would also like to appreciate CGWB for providing groundwater data, and NASA’s GRACE along with other missions, for processing and maintaining satellite dataset. We also want to give our deepest sense of gratitude to the editors and all the anonymous reviewers for their valuable suggestions and meticulous comments, which will help in improving the manuscript.

Data and code availability statement

The data generated and analysed in the current study are not publically available for the study area but can be available from the corresponding author on the reasonable request. The method describe in this article are openly available in: https://github.com/Ashy3874/Ground_Water_Prediction .

Author Contribution

Conceptualisation and methodology has been done by Dr. Pranshu Pranjal and Dr. Ashish Soni, whereas investigation, data preparation, formal analysis, writing and preparation of the original draft has been done by Dr. Pranshu Pranjal. Dr. Ashish Soni has contributed in data curation and processing, along with its validation and visualisation. Supervision, administration and procurement of the resources has been done by Prof. Dheeraj Kumar. While, writing, reviewing and editing has been done by Prof. Dheeraj Kumar and Dr. R.S. Chatterjee. All authors have read and agreed to the published version of the manuscript.

Funding Details

This research received no specific grant from any funding agency

Conflict of Interest

The authors have no conflicts of interest.

(CGWB) CGWB. 2014. Dynamic groundwater resources of India (as on March, 2011).
(CGWB) CGWB. 2017. Dynamic ground water resources of India (as on March, 20017).
Aeschbach-Hertig W, Gleeson T. 2012. Regional strategies for the accelerating global problem of groundwater depletion. Nat Geosci. 5(12):853–861.
Alizadeh MJ, Kavianpour MR, Danesh M, Adolf J, Shamshirband S, Chau K-W. 2018. Effect of river flow on the quality of estuarine and coastal waters using machine learning models. Eng Appl Comput Fluid Mech. 12(1):810–823.
Asoka A, Gleeson T, Wada Y, Mishra V. 2017. Relative contribution of monsoon precipitation and pumping to changes in groundwater storage in India. Nat Geosci. 10(2):109–117.
Bhanja SN, Mukherjee A. 2019. In situ and satellite-based estimates of usable groundwater storage across India: Implications for drinking water supply and food security. Adv Water Resour. 126:15–23.
Bhanja SN, Mukherjee A, Rodell M. 2018. Groundwater storage variations in India. In: Groundw South Asia. [place unknown]: Springer; p. 49–59.
Chatterjee RS, Pranjal P, Jally S, Kumar B, Dadhwal VK, Srivastav SK, Kumar D. 2019. Potential groundwater recharge in northwestern India vs spaceborne GRACE gravity anomaly based monsoonal groundwater storage change for evaluation of groundwater potential and sustainability. Groundw Sustain Dev. 10:100307.
Chen, J., Famigliett, J.S., Scanlon, B.R., Rodell, M. 2016. Groundwater Storage Changes: Present Status from GRACE Observations. Remote Sensing and Water Resources. 55.
Flechtner F. 2007. GFZ Level-2 Processing Standards Document For Level-2 Product Release 0004 (Rev. 1.0, February 19, 2007).
Huang F, Zhang Y, Zhang D, Chen X. 2019. Environmental groundwater depth for groundwater-dependent terrestrial ecosystems in arid/semiarid regions: A review. Int J Environ Res Public Health. 16(5):763.
Huffman GJ, Adler RF, Bolvin DT, Nelkin EJ. 2010. The TRMM multi-satellite precipitation analysis (TMPA). In: Satell rainfall Appl Surf Hydrol. [place unknown]: Springer; p. 3–22.
Huffman GJ, Bolvin DT, Nelkin EJ, Wolff DB, Adler RF, Gu G, Hong Y, Bowman KP, Stocker EF. 2007. The TRMM Multisatellite Precipitation Analysis (TMPA): Quasi-global, multiyear, combined-sensor precipitation estimates at fine scales. J Hydrometeorol. 8(1):38–55.
Kadiyan, N., Chatterjee, R.S., Pranjal, P., Agrawal, P., Jain, S.K., Angurala, M.L., Biyani, A.K., Sati, M.S., Kumar, D., Bhardwaj, A. and Champati Ray, P.K., 2021. Assessment of groundwater depletion–induced land subsidence and characterisation of damaging cracks on houses: a case study in Mohali-Chandigarh area, India. Bulletin of Engineering Geology & the Environment, 80(4).
Koren V, Schaake J, Mitchell K, Duan Q-Y, Chen F, Baker JM. 1999. A parameterization of snowpack and frozen ground intended for NCEP weather and climate models. J Geophys Res Atmos. 104(D16):19569–19585.
Kumar, R., Dwivedi, S.B. and Gaur, S., 2021. A comparative study of machine learning and Fuzzy-AHP technique to groundwater potential mapping in the data-scarce region. Computers & Geosciences, 155, p.104855.
Melesse, A. M. , Abtew, W. , & Dessalegne, T. 2009. Evaporation estimation of rift valley lakes: comparison of models. Sensors (Basel, Switzerland). 9(12):9603–9615.
Ning, S.; Ishidaira, H.; Wang, J. 2014. Statistical Downscaling of GRACE-derived Terrestrial Water Storage Using Satellite and GLDAS Products. J. Jpn. Soc. Civ. Eng. Ser. B1 Hydraul. Eng. 70:133-138.
Nwachukwu, PN, Satge, F.; Yacoubi, S.E.; Pinel, S.; Bonnet, M.P. 2020. From TRMM to GPM: How Reliable Are Satellite-Based Precipitation Data across Nigeria? Remote Sens. 12, 3964.
Pranjal P, Kadiyan N, Chatterjee RS, Kumar D, Sati MS. 2021. Interpreting land subsidence impacts due to groundwater depletion using remote sensing-based GRACE gravity anomaly and DInSAR technique: a study on northwestern parts of India. Environ Earth Sci. 80(17):1–12.
Reager JT, Famiglietti JS. 2013. Characteristic mega-basin water storage behavior using GRACE. Water Resour Res. 49(6):3314–3329.
Richey AS, Thomas BF, Lo M-H, Reager JT, Famiglietti JS, Voss K, Swenson S, Rodell M. 2015. Quantifying renewable groundwater stress with GRACE. Water Resour Res. 51(7):5217–5238.
Rodell M, Houser PR, Jambor, UEA, Gottschalck J, Mitchell K, Meng C-J, Arsenault K, Cosgrove B, Radakovich J, Bosilovich M, others. 2004. The global land data assimilation system. Bull Am Meteorol Soc. 85(3):381–394.
Rodell M, KH B. 2007. GLDAS CLM land surface model L4 monthly 1.0$\times$ 1.0 degree..
Rodell M, Velicogna I, Famiglietti JS. 2009. Satellite-based estimates of groundwater depletion in India. Nature. 460(7258):999–1002.
Sari AP, Suzuki H, Kitajima T, Yasuno T, Prasetya DA, Rabi A. 2020. Prediction of wind speed and direction using encoding-forecasting network with convolutional long short-term memory. In: 2020 59th Annu Conf Soc Instrum Control Eng Japan. [place unknown]; p. 958–963.
Scanlon BR, Mukherjee A, Gates J, Reedy RC, Sinha AK. 2010. Groundwater recharge in natural dune systems and agricultural ecosystems in the Thar Desert region, Rajasthan, India. Hydrogeol J. 18(4):959–972.
Scanlon, B. R., Zhang, Z., Save, H., Wiese, D. N., Landerer, F. W., Long, D., Longuevergne, L., and Chen, J. 2016. Global evaluation of new GRACE mascon products for hydrologic applications, Water Resour. Res., 52:9412–9429.
Scher S. 2018. Toward data-driven weather and climate forecasting: Approximating a simple general circulation model with deep learning. Geophys Res Lett. 45(22):12–616.
Sidhu BS, Sharda R, Singh S. 2021. Spatio-temporal assessment of groundwater depletion in Punjab, India. Groundw Sustain Dev. 12:100498.
Soni A, Koner R, Villuri VGK. 2020. M-unet: Modified u-net segmentation framework with satellite imagery. In: Proc Glob AI Congr 2019. 47–59.
Tapley BD, Bettadpur S, Ries JC, Thompson PF, Watkins MM. 2004. GRACE measurements of mass variability in the Earth system. Science (80- ). 305(5683):503–505.
Wada, Y. , Reager, J. T. , Chao, B. F. , Wang, J. , Lo, M.‐H. , Song, C. , Li, Y. , & Gardner, AS. 2017. Recent changes in land water storage and its contribution to sea level variations. Surveys in Geophysics, 38(1):131–152.
Wahr, J.; Molenaar, M.; Bryan, F. 1998. Time variability of the Earth’s gravity field: Hydrological and oceanic effects and their possible detection using GRACE. J. Geophys. Res. Solid Earth. 103:30205–30229.
Watkins MM, Wiese DN, Yuan D-N, Boening C, Landerer FW. 2015. Improved methods for observing Earth’s time variable mass distribution with GRACE using spherical cap mascons. J Geophys Res Solid Earth. 120(4):2648–2671.
Xingjian SHI, Chen Z, Wang H, Yeung D-Y, Wong W-K, Woo W. 2015. Convolutional LSTM network: A machine learning approach for precipitation nowcasting. In: Adv Neural Inf Process Syst. [place unknown]; 802–810.
Yang, X.; Tian, S.; Feng, W.; Ran, J.; You, W.; Jiang, Z.; Gong, X. 2020. Spatio-Temporal Evaluation of Water Storage Trends from Hydrological Models over Australia Using GRACE Mascon Solutions. Remote Sens. 12:3578.
Ye L, Gao L, Marcos-Martinez R, Mallants D, Bryan BA. 2019. Projecting Australia’s forest cover dynamics and exploring influential factors using deep learning. Environ Model & Softw. 119:407–417.
Yoon, H., Hyun, Y., Ha, K., Lee, K.K. and Kim, G.B., 2016. A method to improve the stability and accuracy of ANN-and SVM-based time series models for long-term groundwater level predictions. Computers & geosciences. 90:144-155.

No competing interests reported.

Investigating Spatio-temporal behaviour for Groundwater in North-West India: A Deep Learning Approach

Status:

Journal Publication

Version 1

Abstract

Figures

1. Introduction

2. Study Area

3. Dataset Used

3.1 Groundwater storage estimates using in-situ GW level (GWL) data

3.2 Satellite-based groundwater storage estimates

3.3 Precipitation

3.4 Soil moisture change

3.5 Evapotranspiration

4. Experimental Setup

4.1 Data processing

4.2 Model parameter

4.3 Evaluation metrics

5. Results and discussion

6. Conclusion

Declarations

References

Additional Declarations

Status:

Journal Publication

Version 1

\({i}_{t}=\sigma ({W}_{xi}{X}_{t}+{W}_{hi}{H}_{t-1}+{W}_{ci}o{C}_{t-1}+{b}_{i})\)	(2)
\({f}_{t}=\sigma ({W}_{xf}{X}_{t}+{W}_{hf}{H}_{t-1}+{W}_{cf}o{C}_{t-1}+{b}_{f})\)	(3)
\({C}_{t}=tanhtanh ( {W}_{xc}{W}_{t}+{W}_{hc}{H}_{t-1}+{b}_{c})\)	(4)
\({C}_{t}={f}_{t}o{C}_{t-1}+{i}_{t}o{C}_{t}\)	(5)
\({o}_{t}=\sigma ({W}_{x0}{X}_{t}+{W}_{ho}{H}_{t-1}+{W}_{co}o{C}_{t}+{b}_{o})\)	(6)
\({H}_{t}={o}_{t}\left({C}_{t}\right)\)	(7)