Abstract
In the present study, a new hybridization strategy for predicting the groundwater table (GWT) and drought analysis is presented. Therefore, a hybrid of the bi-directional long short-term model (BLSTM) and the Harris hawk optimization (HHO) algorithm, namely the BLSTM–HHO algorithm, is applied. In this algorithm, the lagged data of the GWT are used as the input, whereas the current GWT data are used as the output. Additionally, the standalone BLSTM, the long short-term model (LSTM), artificial neural networks (ANN), Seasonal Autoregressive Integrated Moving Average (SARIMA), and the Autoregressive Integrated Moving Average (ARIMA) are employed as benchmark simulating algorithms. The results show that the BLSTM–HHO algorithm has more accuracy than the other investigated simulating algorithms based on the different evaluation criteria such as relative root mean squared error (RRMSE), Nash–Sutcliffe coefficient (NSE), and refined Willmott index (dr). The prediction results (from 2018 to 2022) in all three investigated aquifers show the decline of the GWT (−5.40 m for Brojen aquifer, −7.23 m for Javanmardi aquifer, and −5.81 m for Shahrekord aquifer). Accordingly, the drought analysis by the ground resource index (GRI) in the investigated areas shows that drought is expected to be continued for the next 5 years with an increasing magnitude of severity.
HIGHLIGHTS
A new hybridization strategy (BLSTM–HHO) has been introduced to predict the GWT and drought analysis.
The BLSTM–HHO algorithm has a high performance compared with other well-known algorithms.
The HHO algorithm has enhanced the accuracy of the BLSTM in simulating the GWT.
The BLSTM–HHO algorithm has the potential to predict various hydrological and geohydrological parameters.
Graphical Abstract
INTRODUCTION
Groundwater is one of the essential freshwater resources for domestic, irrigation, and industrial uses, especially in arid and semiarid aquifers. It supplies more than the required water for 1.5 billion people and almost 40% of the irrigation water demand (Melaku & Wang 2019). According to previous studies, hydrological parameters such as precipitation, temperature, and transpiration evaporation in the coming years will be associated with changes that can reduce the inflow to aquifers. For example, the rainfall trend will be declining in some parts of Serbia (Amiri & Gocić 2021; Arab Amiri & Gocić 2021). It was reported that evapotranspiration in Iran's Urmia and Sefidrood basins would increase in the future (Kadkhodazadeh et al. 2022). On the other hand, land-use changes, growing population, and overexploitation of resources lead to depleted groundwater. This issue can lead to a failure in sustainable development. In these conditions, modeling and prediction of groundwater can be beneficial. There are different modeling methods in which time series and artificial intelligence are competitive compared to physical and mathematical ones; due to no requirement of initial and boundary conditions, compatible accuracy, and lower computation.
Time series analysis methods and their hybrids with adaptive neuro-fuzzy inference systems were used for groundwater table (GWT) prediction (Shirmohammadi et al. 2013). Mirzavand & Ghazavi (2014) employed different time series methods for GWT prediction in arid aquifers, including Autoregressive (AR), Moving Average (MA), Autoregressive Moving Average (ARMA), and Autoregressive Integrated Moving Average (ARIMA). In this study, it was shown that the AR method had better accuracy than other the investigated methods. Applying ARIMA and artificial neural networks (ANN) for groundwater depletion prediction in semiarid aquifers (Choubin & Malekian 2017) showed that ARIMA had better applications than the ANN. GWT time series prediction by the ARIMA method indicated a slightly increasing trend in most wells (Gibrilla et al. 2018). The mathematical models, the Bayesian network (BN), and the ANN employed for GWT prediction have shown that the BN was superior to mathematical models and ANN (Moghaddam et al. 2019). The application of the ANN, genetic programming (GP), support vector machine (SVM), and extreme learning machine (ELM) showed that ELM had better outcomes than the other investigated algorithms (Natarajan & Sudheer 2020).
However, GWT's fluctuation and the nonlinear nature of GWT time series cause difficulties in modeling traditional time series and artificial intelligence methods. Consequently, deep learning (DL) methods can be an alternative for enhanced modeling of the GWT. In recent years, the DL method has been applied successfully in hydrology and hydrogeology. The long short-term memory recurrent neural network (LSTM-RNN) was used for low-flow hydrological time series forecasting (Sahoo et al. 2019). The DL neural network was used for predicting flash flood susceptibility (Tien Bui et al. 2020). This method had better performance than the ANN and SVM. The DL method, ELM, and the Gaussian process regression (GPR) were compared in modeling the GWT, and the results revealed DL to be superior (Kumar et al. 2020). The LSTM-based methods and the hydrological simulated program-fortran (HSPF) were compared for sub-surface flow estimation and it was found that LSTM-based methods had better accuracy than the HSPF (Abbas et al. 2020). Thus, DL methods have been employed in other fields such as water quality (Barzegar et al. 2020; Liang et al. 2020; Ma et al. 2020), remote sensing (Li et al. 2019; Shen et al. 2020), and landslide prediction (Bui et al. 2020; Dao et al. 2020).
Note that as DL methods have numerous parameters, changes in any part could have a significant impact on the ultimate accuracy. However, there is no particular method for estimating these parameters except trial and error (Gocić & Arab Amiri 2021), which requires a lot of computation time, so using metaheuristic optimization algorithms is an alternative to determine the parameters of DL methods with an accurate precision. The metaheuristic optimization algorithms have received the attention of many researchers in many fields such as water engineering by genetic algorithm (GA), particle swarm optimization (PSO), bat algorithm (BA), shark smell optimization algorithm (SMO), hybrid BA, and PSO (Valikhan-Anaraki et al. 2019), hydraulic engineering by PSO, BA, and hybrid of BA and PSO (Farzin & Valikhan Anaraki 2021a), hydrology by GA, PSO, and grasshopper optimization algorithm (GOA) (Farrokhi et al. 2020), and optimization of machine learning parameters by ant colony optimization for continuous domains (ACOR), GA, PSO, differential evolution (DE), firefly algorithm (FFA), and whale optimization -algorithm (WOA) (Azad et al. 2019a, 2019b; Farzin et al. 2020; Anaraki et al. 2021). One of the state-of-the-art metaheuristic optimization algorithms is the Harris hawk optimization (HHO) algorithm introduced by Heidari et al. (2019). The HHO algorithm has been practiced in various fields such as optimization of the ANN for predicting the slope stability (Moayedi et al. 2019), optimization of least square support vector machine (LSSVM), ANN, and multiple linear regression (MLR) for rainfall–runoff modeling (Tikhamarine et al. 2020), designing optimal water distribution network for Homashahr city of Iran (Khalifeh et al. 2020), and finding optimal parameters of the ANN for predicting scour depth (Sammen et al. 2020).
One of the main aims of GWT simulation and prediction is the drought assessment. Drought can impact groundwater resources and different aspects in the development of a country. Among various groundwater drought indices, the groundwater resource index (GRI) (Mendicino et al. 2008) has more popularity owing to its greater simplicity and accuracy (Azimi et al. 2019; Azimi et al. 2020).
To the best of the authors’ knowledge, there is no study on the GWT time series modeling and prediction by a combination of state-of-the-art DL and metaheuristics optimization algorithms. Hence, in the present study, the hybrid of the BLSTM and HHO algorithm will be employed to model and predict the GWT in the selected aquifers of Chaharmahal and Bakhtiari province in Iran. It is worth mentioning that BLSTM and HHO are the newest and most potent DL and optimization methods, respectively. The hybrid of bi-direction DL methods and optimization algorithms has not been used to model and predict groundwater levels or other fields. Therefore, the algorithm developed in the present study (BLSTM–HHO) can be used as a new tool in groundwater resource management. Continuously, the results will be compared with ARIMA, SARIMA, ANN, LSTM, and standalone BLSTM as benchmark methods. Also, to achieve less uncertain outcomes, the GWT simulation and prediction will be done based on the multi-run of the BLSTM–HHO algorithm. Furthermore, the drought based on the GRI will be expected for the next 5 years.
MATERIALS AND METHODS
Presented strategy for simulation and prediction of the GWT and drought analysis
In the present study, a new hybrid algorithm is introduced based on the hybrid of BLSTM and HHO for GWT prediction and analysis of following years’ drought. This algorithm consists of preparing data, simulating the GWT by different time series methods, selecting the best algorithm, predicting the GWT for the next 5 years, and calculating the GRI as drought criteria. The mentioned sections are stated as follows:
The GWT data from different piezometer wells in various aquifers were collected from IRAN Water Resources Management Company (http://www.wrm.ir/). Then, these data were prepared by other lag times (1, 3, 6, and 12 months).
The time series algorithms including ARIMA, SARIMA, ANN, LSTM, BLSTM, and BLSTM–HHO were employed for the GWT time series simulation. This section determines the parameters of ARIMA, SARIMA, ANN, LSTM, and BLSTM by sensitivity analysis such as in Farzin & Valikhan Anaraki (2021a). Besides, the BLSTM has a complex structure and essential parameters that can significantly impact the final accuracy of results. Hence, the HHO algorithm was used to find the best parameters. The hybrid of BLSTM and HHO is called the BLSTM–HHO in the present study.
After simulation, the assessment criteria were used for choosing the best algorithm. The best algorithm ran 200 times, and confidence intervals were calculated.
The most selected fitted algorithm was practiced for the prediction of the GWT for 200 times. Next, a variation of the GWT was assessed, and confidence intervals for the prediction were regulated.
The GRI is estimated for different aquifers in the observed and future periods.
Figure 1 shows the presented strategy for simulation and prediction of the GWT and drought analysis.
Case study and data source
The case study is Chaharmahal and Bakhtiari province in Iran as it consists of 10% of the country's water resources. It covers 16,421 km2 with the population equal to 947,763. Groundwater resources supply 70% of water usage of Chaharmahal and Bakhtiari province. Therefore, decreasing the GWT could cause harmful impacts on the environment as well as economic crises. In the present study, the GWT of several piezometer wells in three aquifers, namely Javanmardi, Brojen, and Shahrekord, were used for GWT simulation and prediction. Figure 2 shows the location of the study area and its piezometer wells.
The data used in the present study include the GWT for Brojen from October 1985 to September 2017, for Javanmardi from October 1986 to September 2017, and for Shahrekord from October 1985 October to September 2017. The statistical criteria of GWT data for these three aquifers are presented in Table 1.
Statistical criteria of data used
Aquifers . | Average (m) . | Standard deviation (m) . | Coefficient of variation (%) . | Range (m) . |
---|---|---|---|---|
Brojen | 2,193 | 0.1414 | 0.0064 | 13 |
Javanmardi | 1,783 | 1.4142 | 0.0793 | 20 |
Shahrekord | 2,067 | 0.0707 | 0.0034 | 21 |
Aquifers . | Average (m) . | Standard deviation (m) . | Coefficient of variation (%) . | Range (m) . |
---|---|---|---|---|
Brojen | 2,193 | 0.1414 | 0.0064 | 13 |
Javanmardi | 1,783 | 1.4142 | 0.0793 | 20 |
Shahrekord | 2,067 | 0.0707 | 0.0034 | 21 |
Table 2 shows the input combination for simulating the GWT. It is observed that four different input combinations include the GWT with 1-month (Lag = 1); 1- and 3-month (Lag = 3); 1-, 3-, and 6-month (Lag = 6); and 1-, 3-, 6-, and 12-month (Lag = 12) lag times are considered as input data.
Statistical criteria of data used
Models . | Input combination . |
---|---|
Lag = 1 | GWTt−1 |
Lag = 3 | GWTt−1, GWTt−3 |
Lag = 6 | GWTt−1, GWTt−3, GWTt−6 |
Lag = 12 | GWTt−1, GWTt−3, GWTt−6, GWTt−12 |
Models . | Input combination . |
---|---|
Lag = 1 | GWTt−1 |
Lag = 3 | GWTt−1, GWTt−3 |
Lag = 6 | GWTt−1, GWTt−3, GWTt−6 |
Lag = 12 | GWTt−1, GWTt−3, GWTt−6, GWTt−12 |
Autoregressive Integrated Moving Average
Here, GWTt is the GWT for tth month, is a random number with zero mean and a variance equal to
, c is a constant value, and
is the parameter of AR.
More information about ARIMA is presented in Appendix A.
Seasonal Autoregressive Integrated Moving Average
Other explanations of SARIMA are presented in Appendix B.
Artificial neural network
The ANN is one of the well-known ML methods inspired by human brain neural networks. It consists of one input layer, multiple layers, and one output layer. Each of the layers in the ANN is comprised of one or more neurons. In the input layer, input data are imposed on the network, and each input is equivalent to one neuron. Then, input data are processed in hidden layers by passing from hidden layers neurons. Each hidden layer neuron in the ANN consists of one transform function, where the input after multiplying in corresponding weights and summing with bias coefficient are imposed on this transform function. Then, the output of this neuron is used as the input for other neurons in the next layers. During this process, the corresponding weights and biases for each neuron are determined by training algorithms (Valikhan Anaraki et al. 2022). There are different training algorithms, and one of the best ones is adam as previously reported by Farzin et al. (2020). The diagram of the ANN structure is shown in Figure C.1 (Appendix C).
Bi-directional long short-term memory
The BLSTM is one of the state-of-the-art simulation algorithms proposed by Graves & Schmidhuber (2005). To understand the BLSTM method, first, the LSTM method is stated. The LSTM is presented by Hochreiter & Schmidhuber (1997). Similar to the RNN, the LSTM has several memory cells.
One of the new extensions of the LSTM is the BLSTM, which learns the relationship between the time series of inputs and outputs based on the complete time series. The BLSTM has two blocks in each layer that work in opposite directions (feedforward and backward). Therefore, the BLSTM uses the past and future information of the time series at any point in the time series. Figures D.1 (Appendix D) and Figure 3 show the diagram of the LSTM and BLSTM, respectively. It is worth mentioning that the structure of the LSTM, such as the ANN, is composed of one input layer, one or more hidden layers, and one output layer. Thus, each layer consists of multi neurons. The formulation of the BLSTM is given in Appendix D.
HHO algorithm















Hybrid of LSTM and HHO (BLSTM–HHO)


Groundwater drought index




Assessment criteria
Finally, different criteria are employed for assessing the performance of the introduced model in simulating the GWT including mean absolute error (MAE), relative root mean squared error (RRMSE), Nash–Sutcliffe coefficient (NSE), and refined Willmott index (dr), Pearson correlation coefficient (R). The MAE, RRMSE, and R are formulated in G.1–G.5 relations of Appendix G. The mentioned criteria were previously studied (Valikhan-Anaraki et al. 2019; Farzin & Valikhan Anaraki 2021b; Kadkhodazadeh & Farzin 2021).
RESULTS
Analyzing the GWT time series
Changes in the yearly GWT and the 5-year MA in the Brojen, Javanmardi and, Shahrekord aquifers are shown in Figure 4. Note that positive and negative values indicate the decreasing and increasing GWT, respectively. As shown in Figure 4, the GWT was oscillating and nonlinear, so there is a need for robust methods to predict future changes.
Changes of GWT in (a) Brojen, (b) Javanmardi, and (c) Shahrekord aquifers.
ARIMA parameter setting
The ARIMA parameters included p, d, and q, which are determined by automatic sensitivity analysis in the python software. Based on the mentioned method, p, d, and q parameters were set to 6, 1, and 0 in the Brojen aquifer; 12, 1, and 1 in the Javanmardi aquifer; and 6, 1, and 0, respectively, for the Shahrekord aquifer.
SARIMA parameter setting
In the SARIMA method, the important parameters p, d, q, P, D, Q, and m are determined by automatic sensitivity analysis such as the ARIMA method. The p, d, q, P, D, Q, and m parameters in the Brojen aquifer were determined by the values of 0, 1, 0, 2, 0, 1, and 12, respectively. In the Javanmardi aquifer, these parameters were 0, 1, 0, 2, 0, 1, and 12, respectively. In the Shahrekord aquifer, mentioned parameters were 1, 1, 0, 2, 0, 1, and 12, respectively.
ANN parameter setting
Figure 5 shows the results of the sensitivity analysis of ANN parameters in different lag times and aquifers. Variations of MSE (MSE of the observed and simulated GWT) by changing ANN parameters, including the number of neurons in the first and second layer, are shown in Figure 5. Thus, the most and most minor sensitivities of the ANN were for the Shahrekord and the Brojen aquifers, respectively, and thus the dependency on lag time can be concluded.
The best lag times (lag time with minimum MSE) for the ANN algorithm in different aquifers are shown in Figure 6. It is evident that in the Brojen aquifer, the best lag time was 6 months, whereas in the Javanmardi and Shahrekord aquifers, the best lag time was 12 months.
The best number of neurons in the first and second layers of ANN are, respectively, shown in Figure 7. As seen, the best number of neurons of first and second layers for best lag times in the Brojen aquifer (lag = 6 months) were 3 and 4, in the Javanmardi aquifer (lag = 12 months) were 1 and 2, and in the Shahrekord aquifer (lag = 12 months) was 5.
LSTM parameter setting
The LSTM algorithm sets the number of hidden layers and neurons in both layers as 2 and 100, respectively. The results of MSE for different lag times in different aquifers are shown in Figure 8. The best lag time was 1 month for all areas.
BLSTM parameter setting
In the BLSTM, the number of hidden layers and the number of neurons in them are 2 and 100, respectively. Also, the MSE results for different lag times in different aquifers are compared. As shown in Figure 9, the minimum MSE in all aquifers was related to 1 month of the lag time.
BLSTM–HHO parameter setting
Figure 10 demonstrates the results of the optimization of the BLSTM by the HHO. It is evident that for the Brojen aquifer, the best results are for 1-month lag time. The minimum value of MSE for this lag time was 0.1898 m2. In the Javanmardi aquifer, the minimum value for MSE was for a12-month lag time by a value of 1.0964 m2. Based on the other results of this figure, the minimum MSE in the Shahrekord aquifer was 0.3016 m2 and was for a 12-month lag time.
Results of optimization BLSTM by the HHO for the different lag time.
In Figure 11, the best parameters of BLSTM–HHO for Brojen, Javanmardi, and Shahrekord aquifers are shown. As shown in this figure, the best number of hidden layers and number of neurons in the first and second layers at the Brojen aquifer were 2, 255, and 65; at Javanmardi aquifer they were 2, 120, and 74; and at Shahrekord aquifer they were 1 and 51, respectively.
Comparing time series algorithms
The results of the evaluation criteria of investigated algorithms in the Brojen aquifer including ARIMA, SARIMA, ANN, LSTM, BLSTM, and BLSTM–HHO are compared and presented in Table 3. As seen, all algorithms had good results in the training period. However, during the testing period, ARIMA and SARIMA failed, and the ANN had moderate accuracy. In contrast, the accuracy of LSTM, BLSTM, and BLSTM–HHO was significantly better than other algorithms. Thus, the results of LSTM, BLSTM, and BLSTM–HHO were very close.
Performance of investigated algorithms in the Brojen aquifer
Training period . | MAE (m) . | RMSE (m) . | RRMSE . | dr . | NSE . |
---|---|---|---|---|---|
ARIMA | 0.3976 | 0.4614 | 0.1655 | 0.9725 | 0.9110 |
SARIMA | 0.3155 | 0.3609 | 0.1295 | 0.9832 | 0.9294 |
ANN | 0.6805 | 0.8727 | 0.3131 | 0.9016 | 0.8476 |
LSTM | 0.5371 | 0.6691 | 0.2412 | 0.9416 | 0.8792 |
BLSTM | 0.5372 | 0.6690 | 0.2412 | 0.9417 | 0.8792 |
BLSTM–HHO | 0.5400 | 0.6695 | 0.2413 | 0.9416 | 0.8786 |
ARIMA | 0.6040 | 0.7908 | 0.9476 | 0.0902 | 0.5326 |
SARIMA | 2.8399 | 3.1807 | 3.8113 | −13.7171 | −0.5449 |
ANN | 0.4576 | 0.6156 | 0.7377 | 0.4486 | 0.6460 |
LSTM | 0.3109 | 0.4224 | 0.5111 | 0.7354 | 0.7571 |
BLSTM | 0.3109 | 0.4224 | 0.5111 | 0.7355 | 0.7571 |
BLSTM–HHO | 0.3127 | 0.4237 | 0.5126 | 0.7339 | 0.7557 |
Training period . | MAE (m) . | RMSE (m) . | RRMSE . | dr . | NSE . |
---|---|---|---|---|---|
ARIMA | 0.3976 | 0.4614 | 0.1655 | 0.9725 | 0.9110 |
SARIMA | 0.3155 | 0.3609 | 0.1295 | 0.9832 | 0.9294 |
ANN | 0.6805 | 0.8727 | 0.3131 | 0.9016 | 0.8476 |
LSTM | 0.5371 | 0.6691 | 0.2412 | 0.9416 | 0.8792 |
BLSTM | 0.5372 | 0.6690 | 0.2412 | 0.9417 | 0.8792 |
BLSTM–HHO | 0.5400 | 0.6695 | 0.2413 | 0.9416 | 0.8786 |
ARIMA | 0.6040 | 0.7908 | 0.9476 | 0.0902 | 0.5326 |
SARIMA | 2.8399 | 3.1807 | 3.8113 | −13.7171 | −0.5449 |
ANN | 0.4576 | 0.6156 | 0.7377 | 0.4486 | 0.6460 |
LSTM | 0.3109 | 0.4224 | 0.5111 | 0.7354 | 0.7571 |
BLSTM | 0.3109 | 0.4224 | 0.5111 | 0.7355 | 0.7571 |
BLSTM–HHO | 0.3127 | 0.4237 | 0.5126 | 0.7339 | 0.7557 |
The values of assessment criteria of investigated algorithms for simulating the GWT in the Javanmardi aquifer are tabulated in Table 4. All algorithms were trained with reasonable accuracy, but, there was a significant difference between the results of BLSTM–HHO and other algorithms in the testing period. The BLSTM–HHO decreased the MAE, RMSE, and RRMSE by 73, 72, and 73%, respectively, and increased the dr and NSE by 260 and 110%, respectively, in comparison to other algorithms. Also, BLSTM and LSTM were placed in the second and third ranks after the BLSTM–HHO.
Performance of investigated algorithms in the Javanmardi aquifer
Training period . | MAE (m) . | RMSE (m) . | RRMSE . | dr . | NSE . |
---|---|---|---|---|---|
ARIMA | 0.8064 | 0.9440 | 0.1705 | 0.9708 | 0.9068 |
SARIMA | 0.4940 | 0.6735 | 0.1208 | 0.9854 | 0.9434 |
ANN | 0.8497 | 1.0967 | 0.1937 | 0.9623 | 0.9043 |
LSTM | 0.9484 | 1.1653 | 0.2116 | 0.9551 | 0.8896 |
BLSTM | 0.9502 | 1.1654 | 0.2116 | 0.9551 | 0.8893 |
BLSTM–HHO | 0.7842 | 1.0457 | 0.1531 | 0.9765 | 0.9273 |
ARIMA | 3.1385 | 3.9364 | 1.2524 | −0.5895 | 0.4046 |
SARIMA | 3.3024 | 3.9599 | 1.2384 | −0.5537 | 0.3836 |
ANN | 0.8314 | 1.0955 | 0.3452 | 0.8792 | 0.8427 |
LSTM | 1.1850 | 1.3440 | 0.4203 | 0.8210 | 0.7788 |
BLSTM | 1.1849 | 1.3430 | 0.4200 | 0.8213 | 0.7789 |
BLSTM–HHO | 0.8557 | 1.1163 | 0.3338 | 0.8865 | 0.8480 |
Training period . | MAE (m) . | RMSE (m) . | RRMSE . | dr . | NSE . |
---|---|---|---|---|---|
ARIMA | 0.8064 | 0.9440 | 0.1705 | 0.9708 | 0.9068 |
SARIMA | 0.4940 | 0.6735 | 0.1208 | 0.9854 | 0.9434 |
ANN | 0.8497 | 1.0967 | 0.1937 | 0.9623 | 0.9043 |
LSTM | 0.9484 | 1.1653 | 0.2116 | 0.9551 | 0.8896 |
BLSTM | 0.9502 | 1.1654 | 0.2116 | 0.9551 | 0.8893 |
BLSTM–HHO | 0.7842 | 1.0457 | 0.1531 | 0.9765 | 0.9273 |
ARIMA | 3.1385 | 3.9364 | 1.2524 | −0.5895 | 0.4046 |
SARIMA | 3.3024 | 3.9599 | 1.2384 | −0.5537 | 0.3836 |
ANN | 0.8314 | 1.0955 | 0.3452 | 0.8792 | 0.8427 |
LSTM | 1.1850 | 1.3440 | 0.4203 | 0.8210 | 0.7788 |
BLSTM | 1.1849 | 1.3430 | 0.4200 | 0.8213 | 0.7789 |
BLSTM–HHO | 0.8557 | 1.1163 | 0.3338 | 0.8865 | 0.8480 |
Table 5 presents the results of six investigated algorithms in the Shahrekord aquifer based on the assessment criteria. All algorithms had good accuracy in the training period in this aquifer, similar to other investigated aquifers. However, ARIMA and SARIMA were the weakest algorithms in the testing period, and BLSTM–HHO was the best algorithm. The BLSTM and LSTM had better accuracy than ANN, SARIMA, and ARIMA. In comparison, the BLSTM performed slightly better than the LSTM.
Performance of investigated algorithms in the Shahrekord aquifer
Training period . | MAE (m) . | RMSE (m) . | RRMSE . | dr . | NSE . |
---|---|---|---|---|---|
ARIMA | 0.3173 | 0.0915 | 0.9916 | 0.9588 | 0.3173 |
SARIMA | 0.2953 | 0.0844 | 0.9929 | 0.9635 | 0.2953 |
ANN | 0.5248 | 0.1497 | 0.9775 | 0.9327 | 0.5248 |
LSTM | 0.5261 | 0.1525 | 0.9767 | 0.9274 | 0.5261 |
BLSTM | 0.5262 | 0.1526 | 0.9766 | 0.9274 | 0.5262 |
BLSTM–HHO | 0.4336 | 0.1194 | 0.9857 | 0.9461 | 0.4336 |
ARIMA | 3.9162 | 4.4846 | 1.9115 | −2.7000 | −0.0212 |
SARIMA | 4.1953 | 4.5813 | 1.9527 | −2.8613 | −0.0863 |
ANN | 0.3703 | 0.4702 | 0.2185 | 0.9516 | 0.8961 |
LSTM | 0.5167 | 0.5843 | 0.2546 | 0.9343 | 0.8624 |
BLSTM | 0.5167 | 0.5836 | 0.2543 | 0.9345 | 0.8624 |
BLSTM–HHO | 0.3055 | 0.3881 | 0.1931 | 0.9622 | 0.9075 |
Training period . | MAE (m) . | RMSE (m) . | RRMSE . | dr . | NSE . |
---|---|---|---|---|---|
ARIMA | 0.3173 | 0.0915 | 0.9916 | 0.9588 | 0.3173 |
SARIMA | 0.2953 | 0.0844 | 0.9929 | 0.9635 | 0.2953 |
ANN | 0.5248 | 0.1497 | 0.9775 | 0.9327 | 0.5248 |
LSTM | 0.5261 | 0.1525 | 0.9767 | 0.9274 | 0.5261 |
BLSTM | 0.5262 | 0.1526 | 0.9766 | 0.9274 | 0.5262 |
BLSTM–HHO | 0.4336 | 0.1194 | 0.9857 | 0.9461 | 0.4336 |
ARIMA | 3.9162 | 4.4846 | 1.9115 | −2.7000 | −0.0212 |
SARIMA | 4.1953 | 4.5813 | 1.9527 | −2.8613 | −0.0863 |
ANN | 0.3703 | 0.4702 | 0.2185 | 0.9516 | 0.8961 |
LSTM | 0.5167 | 0.5843 | 0.2546 | 0.9343 | 0.8624 |
BLSTM | 0.5167 | 0.5836 | 0.2543 | 0.9345 | 0.8624 |
BLSTM–HHO | 0.3055 | 0.3881 | 0.1931 | 0.9622 | 0.9075 |
Figures A.1–Figure A.6 in Appendix H show the distribution patterns of the observed and simulated GWT and R2 values for investigated algorithms with the best lag times. All algorithms showed fitting proficiency in the training period (Figure A.1–Figure A.3). However, the DL methods can estimate the GWT more accurately in the testing period than the ANN, SARIMA, and ARIMA. The value of R2 for BLSTM–HHO, BLSTM, LSTM, ANN, SARIMA, ARIMA approximately ranged in (0.75–0.98), (0.75–0.94), (0.75–0.94), (0.57–0.96), (0.04–0.68), and (0.06–0.47), respectively. Also, it can be seen that the estimated GWT by DL methods was less scattered and closer to straight-line than those of the other investigated algorithms. It is clear from Figures A.1–A.3 that the BLSTM–HHO algorithm was superior to the BLSTM, LSTM, ANN, SARIMA, and ARIMA. The performances of BLSTM and LSTM are similar, with small differences. Also, Figure A.1–Figure A.6 were in accordance with Tables 3–5.
Simulated GWT time series by the best algorithm
The results of 200 simulated GWT time series by the best algorithm, BLSTM–HHO, in different aquifers are shown in Figure 12. The simulated GWT is fitted to observe the GWT in training and testing periods in all areas. Also, the 95% confidence interval for the three aquifers shows that the Brojen aquifer has a narrower 95% confidence interval band and the Javanmardi aquifer has a wider 95% confidence interval band. The average 95% confidence interval widths for Brojen, Javanmardi, and Shahrekord were 0.68, 6.41, and 3.04 m.
Predicting the GWT by the best algorithm
Figure 13 demonstrates the predicted GWT for the future 60 months from 23 September 2017 to 23 September 2022. The average of predictions in all aquifers showed a decreasing trend. In the Brojen aquifer, the average predicted GWT decreased from 2,189.10 to 2,183.70 m (−5.40 m). For the Javanmardi aquifer, the average predicted GWT reduced from 1,819.20 to 1,811.97 m (−7.23 m). Also, in the Shahrekord aquifer, the average predicted GWT decreased from 2,055.30 m to 2,049.49 m (−5.81 m). Moreover, the maximum width of band 95% confidence interval was for the Javanmardi aquifer by 24.14 m, and the minimum width of band 95% confidence interval was for the Brojen aquifer by the value of 3.03 m. The width of the band 95% confidence interval for the Shahrekord aquifer was also 10.66 m. This issue shows more uncertainty in simulating and predicting the GWT in the Javanmardi aquifer.
Groundwater resource index
Finally, the GRI is calculated for the observed and predicted GWT in Brojen, Javanmardi, and Shahrekord aquifers. Figure 14 shows the GRI in all areas. As seen, the GRI was lower than the base-line (GRI = 0) in the Brojen aquifer from October 2008 to the end of the observed period (October 2017), in the Javanmardi aquifer from July 2007 to the end of the observed period, and in the Shahrekord aquifer after October 2008 to the end of the observed period. Therefore, there was a drought indication in the mentioned periods in Brojen, Javanmardi, and Shahrekord aquifers. The drought was observed in all areas in the prediction period. It is worth noting that the magnitude of the severity of this drought is increased based on the average results in all investigating aquifers.
GRI variation in (a) Brojen aquifer, (b) Javanmardi aquifer, and (c) Shahrekord aquifer.
GRI variation in (a) Brojen aquifer, (b) Javanmardi aquifer, and (c) Shahrekord aquifer.
DISCUSSION
According to the results of GWT simulation by the investigated algorithms, BLSTM–HHO, BLSTM, LSTM, SARIMA, and ARIMA were placed in first to fifth order according to the accuracy. Improved accuracy of BLSTM–HHO compared to BLSTM and LSTM can be explained by the optimal determining structure of this algorithm by HHO. Indeed, HHO, by finding the optimal configuration of the BLSTM, leads to simulating the GWT with a maximum accuracy. The BLSTM has also shown improved closed results in comparison to the LSTM. The closed results for these two algorithms happened in a 1-month lag time, whereas the BLSTM provides more accurate results than the LSTM in other lag times. This issue can be discussed through the simplicity of the structure of two algorithms in a 1-month lag time, while, by increasing lag times, the complexity of the structure of two algorithms increased, and the results changed. In addition, the slightly better results of the BLSTM than the LSTM were due to the use of past and future information of time series in the BLSTM, whereas the LSTM only used the past knowledge of time series. The greater accuracy of the BLSTM and LSTM than the ANN was because of utilizing the past and future information of time series in the BLSTM and past data of time series of the LSTM, whereas the ANN only used data on the current state of time series. According to other results, it can be seen that SARIMA and ARIMA had less accuracy than other algorithms. The reason could be that other algorithms process data in multilayers, whereas SARIMA and ARIMA do not have such an ability. The better accuracy of DL-based methods (BLSTM–HHO, BLSTM, and LSTM) than the ANN, SARIMA, and ARIMA is shown in other studies (Sufi Karimi et al. 2019; Apaydin et al. 2020; Dhakal et al. 2020; Gao et al. 2020). However, more accuracy of ANN to SARIMA and ARIMA is similar, as mentioned by Dhakal et al. (2020).
According to prediction results, the GWT for 5 future years decreased in all aquifers based on the 200 prediction results of the best algorithm (BLSTM–HHO ). Thus, the variations in the different areas while predicting the GWT were due to the different characteristics of the observed GWT (Table 1) and uncertainty while simulating and predicting the GWT. In addition, the decrease in the GWT of Chaharmahal and Bakhtiari aquifer, Karaj plain up to 19.089 m (Yousefi et al. 2019), and Tasuj plain up to 5 m (Ghazi et al. 2021) was also reported, which is in parallel with the findings of the present study. The decreasing trend of the GWT can be due to precipitation reduction, increase in temperature and evapotranspiration, and increasing GWT exploitation. Note that a decrease in the GWT can increase drought and shortage of household, agricultural, and industrial water needs and damage to these sectors. The decrease in groundwater level causes shortage of household, agricultural, and industrial water demands.
According to the 95% confidence interval results in simulating and predicting the GWT, it is evident that Javanmardi and Shahrekord aquifers have more uncertainty than the Brojen aquifer. A wider range of the GWT can explain more uncertainty in the Javanmardi and Shahrekord aquifers considering the Brojen aquifer (Table 1). Indeed, a wider range and outliers can lead to more uncertainty and difficulty in simulating and predicting one time series.
The GRI results also show a drought for 5 future years in all aquifers. This drought is for decreasing the GWT in all investigated areas. Indeed, population growth and land-use changes lead to more water needs and more depletion of groundwater resources. Besides, the recharge of groundwater resources is decreased because of climate change, consequently reducing precipitation and increasing temperature. These two issues mainly resulted in decreasing GWT causing drought in different aquifers, as shown in the present study.
CONCLUSION AND REMARKS
In the present study, a new algorithm based on the hybrid of DL and the metaheuristic algorithm was introduced for GWT simulation, prediction, and drought assessment in Brojen, Javanmardi, and Shahrekord aquifers in Chahrmahal and Bakhtiari provinces in Iran. The BLSTM was used as a powerful DL method in this algorithm. The HHO algorithm, one of the novel optimization algorithms, was used as a metaheuristic algorithm. The hybrid of these two algorithms is called the BLSTM–HHO algorithm. The standalone BLSTM, LSTM, ANN, SARIMA, and ARIMA were used as the benchmark algorithms. The results showed that BLSTM–HHO in simulated time series of the GWT had superiority over other investigated algorithms. Also, the results of 200 prediction time series of the GWT based on the BLSTM–HHO showed that the GWT decreases up to −5.40, −7.23, and −5.81 m in Brojen, Javanmardi, and Shahrekord aquifers, respectively, for the years 2018–2022. The analysis of groundwater drought based on the GRI showed that drought would exist in all three investigated areas in the mentioned period. Thus, the severity of drought would increase. This issue hinders the sustainable development of the aquifer and can lead to water shortages for stakeholder consumption. Therefore, there is a need for the optimal operation of groundwater resources in this province. By managing groundwater aquifers such as flood spreading or artificial recharging in rainy seasons, the conditions of groundwater resources in this province can be improved.
The BLSTM–HHO has a high potential for modeling and predicting various parameters of water resources such as rainfall, temperature, humidity, discharge. The effect of data decomposition methods on the accuracy of this method can also be tested.
DECLARATION OF INTEREST STATEMENT
The authors declare that there are no known conflicts of interest associated with this article.
DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.