The understanding of water resource changes and a proper projection of their future availability are necessary elements of sustainable water planning. Monitoring GWS change and future water resource availability are crucial, especially under changing climatic conditions. Traditional methods for

The arid and semi-arid regions of the world have historically suffered from depleting freshwater resources, including TWS (Terrestrial Water Storage) and GWS (Ground Water Storage), apart from low rainfall and rising water demands. The depletion of these water resources primarily depends on climatic parameters (precipitation and temperature) and the rising water demand for municipal, agricultural and industrial purposes. The population of Saudi Arabia increased by 15.48% from 27.4 million in 2010 to 32.4 million in 2016, and the volume of extracted groundwater also rose by 27% from 15.8 × 109 cubic meters in 2010 to 21.6 × 109 cubic meters in 2016 [

The GRACE satellite missions provide an excellent opportunity to monitor TWS and GWS [

The present investigation also contributes to two major aspects: forecasting the future status of GWSC using LSTM modeling from July 2020 to December 2025 and analyzing the correlation between groundwater change and rainfall. There are four main objectives of the present investigation: (1) to analyze the changes in TWSC and GWSC from 2003–2020, (2) to analyze the relationships among TWSC, GWSC, satellite-based rainfall and rainfall from meteorological stations from 2009 to 2016; (3) to estimate the net budget modeling of groundwater based on GWSC and rainfall; and (4) to forecast TWSC and GWSC changes from July 2020 to June 2025 using developed and tuned LSTM models. All these objectives were accomplished for entire Saudi Arabia and five basins within Saudi Arabia. A total of 364 tiles of SRTMGL3 DEM were downloaded from LPDAAC (

LSTM networks are a type of RNN that use special units (cells) and standard units to overcome the limitation of traditional RNN [

Input node

Recent studies have used LSTM models [

Terrestrial water storage (TWS) was obtained from GRACE and GRACE-FO satellites. These satellite missions were jointly launched by NASA (

The monthly rainfall dataset of Saudi Arabia from January 2009 to December 2018 was procured from the General Authority for Statistics [

TWS was estimated using GRACE and GRACE-FO satellites from Jan 2003 to June 2020 at 0.5° spatial and monthly temporal resolutions. GRACE, on its own, is incapable of separating anomalies from several elements of TWS (e.g., surface water storage, canopy water, and soil moisture content). Therefore, it is essential to subtract non-groundwater components. When this process is carried out for TWS data, the GWS can be obtained. The current study area being the arid and semi-arid part of Saudi Arabia, the contribution of surface water storage and canopy water was likely to be insignificant in the overall calculation. However, the soil moisture content (SMC) had to be subtracted from the change in TWS (TWSC) to get the change in GWS (GWSC), see

In the equation, GWSC and SMCC showed changes in groundwater storage and changes in soil moisture content, respectively. Information related to soil moisture was obtained using the Global Land Data Assimilation System (GLDAS) [

The calculation of total TWSC uncertainties (∂TWSC) was performed using the methodology described by [

Mann-Kendall tests [

Here,

An attempt was made to study the correlation analysis among GWSC, TWSC, SatRain and METRain for five basins from Jan 2009 to Dec 2016 as per the availability of a common temporal dataset [^{o} × 0.5^{o}).

Data transformation is crucial before implementing the LSTM model. Three data transformations were applied in the current investigation. First, a lag variable of 1 was applied to remove the decreasing trend in the dataset. The second step was the transformation of time series data into input and output so that the output of a step becomes the input of the next step to forecast the value of the current time step. As described earlier, total data in the time series were 186 monthly values. The first 150 months’ dataset for all five basins were taken for the training (126 months) and testing (24 months) of the LSTM model; the remaining 36 months of data (July 2015-June 2020, data gaps of 2016 and 2017) were kept separate from the training process for the unbiased external validation of the LSTM prediction. The third transformation was the scaling of time series data from −1 to 1. All these three transformations were inverted after the prediction step to get the values at the original scale so that the uncertainty calculation could be adequately assessed.

Keras library version 2.3.0 with TensorFlow version 2.0 at the backend and Python version 3.8.0 was used to build the LSTM models in the current study. The libraries used in the current investigation were NumPy, Pandas, Matplotlib, and scikit-learn. A 4 step procedure was applied to implement LSTM.

The first step was to define the LSTM network, which is organized as a sequence of layers contained by a sequential model. LSTM layer requires that the number of inputs in a three-dimensional shape consists of samples, time steps and features; therefore, reshape function was used to convert our data into a three-dimensional shape based on 130 samples of training, monthly time-step, and features represented by five basins. Two LSTM layers were used in the current investigation. In the first LSTM layer, the return_sequences were set to true, and it indicated that the output of each neuron's hidden state was used as an input to the next LSTM layer. Several hyper-parameters such as optimizer, learning rate, number of units, momentum, and activation functions have to be chosen

The second step was compiling the network. It required several parameters, such as an optimization algorithm to train the network and the loss function to evaluate the network. Several optimizers were tested based on their performances. The third step was the fitting of the LSTM model. The objective of model fitting is to adapt the weights based on a training dataset. It requires the training data to be specified for inputs and corresponding outputs. The initial value of epochs was given as 100 with a validation split value of 0.2 (i.e., 30 values for internal validation of the model) and input and output values.

The fourth and essential step was the prediction using the LSTM model. We forecasted the output step by step for the test data, and the model fed the current forecasted value back into the input window by moving it one step forward to aid forecasting at the next time step using the moving-forward window technique [

Tuning the hyper-parameters of any neural network model is essential for evaluating and appraising the model. The major hyper-parameters of LSTM tuned in the present investigation were (1) number of nodes, (2) number of training epochs, (3) choice of optimizer, and (4) learning rate. A walk-forward model evaluation using the hyper-parameters mentioned above was attempted for different configuration values. The prediction on the test dataset was evaluated based on root mean squared error (RMSE).

The first configuration tuned was based on the number of nodes, which influences the learning capability of the LSTM model. Generally, more nodes can learn more complex mapping at the cost of training time and sometimes cause overfitting. Different nodes (1, 2, 4, 6, and 8) were tested for different configurations. By using four numbers of nodes, a lesser average RMSE value of 1.4 and the least variance based on 20 experimental runs were obtained. However, since it could be an indication of overfitting, dropout was applied to prevent overfitting. The box and whisker plot (

The optimization algorithm was the third tuned hyper-parameter of the present LSTM model. Various optimization algorithms such as Stochastic Gradient Descent (SDG), Adagrad, Adam, AdaDelta, and RMSProp, were tested. The Adam (Adaptive Moment Estimation) algorithm is an extension of SGD that calculates each parameter's learning rate, such as alpha, beta1, beta2 and epsilon [

Learning rate (LR) was the fourth tuned parameter in the current study. The effects of different learning rates for the ADAM optimizer were evaluated; the default value of parameter LR in Keras is 0.001. Similarly, the default value for β1 is 0.9. In the same way, β2 is specified to be 0.999 and ε is specified to be 1e-07. An attempt was made to tune the learning rate with values of 0.1, 0.01, 0.001, 0.0001, and 0.00001. The box plot indicated that LR of 0.1, 0.01, and 0.001 showed good RMSE values of 0.9 and less variance without any significant difference. On the contrary, LR of 0.0001 and 0.00001 showed slightly lower RMSE values with higher variances (

An AutoRegression (AR) based model was used in the present investigation to compare the forecasting performance of the LSTM model (^{2} values for the autoregression model were 0.65 and 0.76, respectively for all the 5 basins, which was significantly good.

The Moment Correlation Coefficient (MCC), Root Mean Square Error (RMSE), Mean Absolute Percentage Error (MAPE) and Nash-Sutcliffe coefficient (NSE) were utilized to evaluate the uncertainty of the LSTM model output. The Mean Absolute Deviation (MAD) was also considered to analyze the LSTM model's accuracy between measured and predicted values.

The MCC summarizes the direction and degree of linear relations between actual and modeled datasets. The correlation coefficient can take values between −1 (perfectly negative correlation) through 0 (no correlation) to +1 (perfectly positive correlation). The MCC formula to compute the correlation coefficient is given in

Here, N represents the number of pairs of data. The terms X and Y are parameters.

RMSE is a method to calculate the error or accuracy in predicting models based on standard deviation

Here, P_i is the i^{th} LSTM predicted value, and Oi is the i^{th} original value.

The MAPE method was used to calculate the prediction accuracy of the LSTM forecast. The calculation was based on the difference between the original values and values forecasted by the LSTM and dividing the original value difference. It was then multiplied by the number of observations and 100 to obtain the percentage error,

Here,

The NSE or efficiency coefficient test determines the magnitude between the variance of residual time series and variance of actual data, and its value ranges from -∞ to 1, see

_{f}

As described earlier, the total uncertainty based on the first residual, second residual and the standard deviation of the second residual was computed for TWSC (±1.2 mm/year). The uncertainty in SMC based on the mean monthly standard deviation from the three GLDAS models was computed as ±0.2 (mm/year). Therefore, the total uncertainty for GWSC was 1.5 (mm/year). Results with a confidence factor ≥ of 90% indicate decreasing annual TWSC and GWSC change for all the five sub-basins (

Basins | TWSC (±1.2 mm/year) | GWSC (±1.5 mm/year) | SNS (TWSC) | SNS (GWSC) | Rainfall change (mm/year) | SNS (Rainfall) | Annual rainfall |
---|---|---|---|---|---|---|---|

1 | −5.88 | −3.5 | −0.53 | −0.33 | +0.15 | +0.01 | 95.60 |

2 | −8.24 | −6.7 | −1.17 | −0.85 | NT | 0.001 | 86.81 |

3 | −8.06 | −5.9 | −1.07 | −0.83 | +0.4 | 0.03 | 79.40 |

4 | −14.12 | −10.7 | −1.45 | −1.22 | NT | 0.001 | 44.18 |

5 | −10.59 | −7.4 | −1.12 | −1.02 | +0.32 | 0.03 | 60.73 |

The analysis of rainfall data from 26 MET stations from 2009 to 2018 showed an interesting trend. Areas 1, 3 and 5 showed a slight increase in rainfall from 2009–2018; however, areas 2 and 4 showed no rainfall change for the same observation years. The rainfall period showed two peaks, one in winter [November-January] and another in summer [March-April]. A positive increasing trend of 0.2 mm/year was observed with a confidence factor ≥ of 80% when combined areas were taken into account. Based on METRain data, for Saudi Arabia, the average annual rainfall value increased from 52 mm for the year 2000 to 60 mm for the year 2015 [

The MCC was performed to understand the relationship between GWSC, TWSC, SatRain and METRain for the five basins from January 2009 to December 2016 as per the common temporal dataset availability (

We have processed GRACE TWSC and GWSC for the KSA; therefore, we have compared the current investigation results with other studies (

Studies | [ |
Current |
[ |
Current |
[ |
Current |
[ |
Current |
---|---|---|---|---|---|---|---|---|

Observation period | 2002–2016 | 2003–2020 | 2002–2016 | 2003–2020 | 2002–2016 | 2003–2020 | 2007–2016 | 2003–2020 |

TWSC | RAK = −3.0 ± 1.1 | RAK = −5.9 ± 1.2 | LMAS = −9.1 ± 1.3 | LMAS = −10.34 ± 1.2 | SAQ = −8.55 ± 0.22 | SAQ = −9.06 ± 1.2 | KSA = 7.85 ± 0.44 mm/year | KSA = −8.01 ± 1.2 |

MAS = −7.2 ± 1.0 | MAS = −8.05 ± 1.2 | WASB = −9.4 ± 1.4 | WASB = −8.05 ± 1.2 | |||||

GWSC | RAK = −1.76 ± 1.4 | RAK = −3.5 ± 1.5 | LMAS = −7.8 ± 1.3 | LMAS = −7.58 ± 1.5 | SAQ = −6.34 ± 0.22 | SAQ = −6.79 ± 1.5 | ||

MAS = −4.8 ± 1.4 | MAS = −5.9 ± 1.5 | WASB = −8.2 ± 1.4 | WASB = −7.44 ± 1.5 |

Studies | Current study | Current study | [ |
[ |
[ |
[ |
[ |
[ |
[ |
[ |
[ |
---|---|---|---|---|---|---|---|---|---|---|---|

Model | LSTM | AR | LSTM | LSTM | SVR | RF | XGBL | XGBT | CESA | Stat mod | ANN |

R2 | 0.81 | 0.76 | 0.36 | 0.71 | 0.8 | 0.74 | 0.8 | 0.74 | 0.31 | 0.30 | 0.81 |

RMSE | 0.5 | 0.65 | 0.7 | 0.26 | 0.28 | 0.33 | 0.27 | 0.32 | 0.56 | 0.23 | 1.5 |

There was a possibility while predicting the future values that the LSTM model output may be uncertain as the model's output was fed back into it as input. Therefore, we initially forecasted TWS and GWS from July 2015 to June 2020 and compared them with the actual values based on the values of coefficient of determination, RMSE, MSE, MAPE, and MAD (

Sub basins | R^{2} |
MSE | RMSE | MAPE (%) | MAD | NSE |
---|---|---|---|---|---|---|

1 | 0.8 | 0.17 | 0.42 | −6.80 | 0.6 | 0.82 |

2 | 0.9 | 0.18 | 0.43 | −2.70 | 1.17 | 0.91 |

3 | 0.84 | 0.37 | 0.61 | −4.25 | 1.21 | 0.87 |

4 | 0.79 | 0.65 | 0.8 | −3.75 | 1.29 | 0.83 |

5 | 0.71 | 0.27 | 0.52 | −3.57 | 0.59 | 0.74 |

The area of the basins with the respective volume of TWSC and GWSC from January 2003 to June 2020 (^{3}/year) and −3.206 (km^{3}/year), respectively, from 2003 to 2020. While comparing the rate of historical (2003–2020) extracted groundwater (GWS) with the forecasted (2020–2025) rate of extracted groundwater (GWS), it was observed that basin three and basin one had shown higher rates −29.68% and −29.58%, while the other three basins had shown rates ranging from −8.89% to 13.66%.

In the absence of the data on groundwater wells in Saudi Arabia, the annual extracted groundwater dataset [

Basins | TWSC July 2020– | GWSC July 2020– |
---|---|---|

Dec 2025 (±1.2 mm/year) | Dec 2025 (±1.5 mm/year) | |

1 | −7.78 | −4.97 |

2 | −9.28 | −7.37 |

3 | −11.6 | −8.39 |

4 | −15.6 | −12.21 |

5 | −12.2 | −8.57 |

Basins | Area (km^{2}) |
TWS vol (km^{3}/year) (03–20) |
GWS vol (km^{3}/year) (03–20) |
TWS vol (km^{3}/year) (20–25) |
GWS vol (km^{3}/year) (20–25) |
---|---|---|---|---|---|

1 | 916100 | −5.387 ± 1.099 | −3.206 ± 1.374 | −7.127 ± 1.099 | −4.553 ± 1.374 |

2 | 24430 | −0.201 ± 0.029 | −0.164 ± 0.037 | −0.227 ± 0.029 | −0.180 ± 0.037 |

3 | 288000 | −2.321 ± 0.346 | −1.699 ± 0.432 | −3.341 ± 0.346 | −2.416 ± 0.432 |

4 | 152700 | −2.156 ± 0.183 | −1.634 ± 0.229 | −2.382 ± 0.183 | −1.864 ± 0.229 |

5 | 193900 | −2.053 ± 0.233 | −1.435 ± 0.291 | −2.366 ± 0.233 | −1.662 ± 0.291 |

However, the present investigation found that the depletion rate did not show steady or increasing trends until June 2020. LSTM-based forecasting until 2025 also suggested that there is a requirement for sustainable water resources planning (

Year | Groundwater extracted (km^{3}) |
GWS(Grace) (km^{3}) |
---|---|---|

2010 | −15.791 | −9.745 |

2011 | −17.299 | −13.3707 |

2012 | −18.94 | −16.812 |

2013 | −20.395 | −18.6753 |

2014 | −21.352 | −20.2442 |

2015 | −22.648 | −21.8787 |

2016 | −21.595 | −25.7241 |

Basin | Extracted volume of groundwater (GWS) (km^{3}/year) |
Recharge estimate based on rainfall (km^{3}/year) |
Water-saving required to achieve sustainability (km^{3}/year) |
---|---|---|---|

1 | 8.081 ± 1.374 | 7.006 ± 1.374 | 1.074 ± 1.374 |

2 | 0.302 ± 0.037 | 0.170 ± 0.037 | 0.132 ± 0.037 |

3 | 3.482 ± 0.432 | 1.829 ± 0.432 | 1.652 ± 0.432 |

4 | 3.234 ± 0.229 | 0.540 ± 0.229 | 2.694 ± 0.229 |

5 | 3.080 ± 0.291 | 0.942 ± 0.291 | 2.137 ± 0.291 |

Total | 18.177 ± 2.363 | 10.487 ± 2.363 | 7.690 ± 2.363 |

During the model tuning, there were clear differences between the number of iterations and associated computational time required. It was often found that the best suited LSTM model in the present investigation used 3 sec/epochs, whereas the LSTM model developed by [

The present investigation provides an understanding of the historical and projected GWS changes at a large-scale using GRACE data for Saudi Arabia. Deep Learning LSTM models were developed based on rigorous hyper-parameters tuning to forecast the water storage changes based on GRACE-derived GWS. The major hyper-parameters of LSTM tuned in the present investigation were (1) number of nodes, (2) number of training epochs, (3) choice of optimizer, and (4) learning rate. A walk-forward model evaluation using the hyper-parameters mentioned above was attempted for different configuration values. The optimized model took 25% computational time compared to the model, which took 8000 epochs in 400 min. The correlation coefficient, MAPE, MAD, NSE and RMSE values were applied to test the LSTM models’ outcomes. GRACE-derived GWS indicated that all the five investigated basins were experiencing groundwater depletion rates of 18.17 km^{3}/year. Future forecasting based on the developed LSTM model indicates that, from July 2020 to Dec 2025, the investigated five basins are likely to experience depletion of water at rates ranging from −7.78 ± 1.2 to −15.6 ± 1.2 for TWSC and −4.97 ± 1.5 to −12.21 ± 1.5 for GWS. The future scope of the present investigation is to add

The significant limitations of the present study can be discussed with respect to the following points: (1) Data continuity and temporal availability of different datasets. GRACE data was available from 2002 to 2020 with data gaps in 2002 and 2017–2018 due to no data or no coverage. (2) Different resolutions of the dataset, e.g., GLDAS resolution was 1.0^{o}, and GRACE-derived GWS and TWS were gridded to 1.0^{o}. Similarly, SatRain resolutions were 0.25^{o} and 0.10^{o} for TRMM and GPM, respectively. (3) GWS and TWS were the change rate values, whereas SatRain and METRain were the time series data. (4) In the absence of data about groundwater wells for entire Saudi Arabia, the annual extracted groundwater dataset [

We thank the anonymous referees for their useful suggestions.