The study area is located in the OER River Basin, in the north-central part of Morocco, between latitudes 31° 00' and 33° 00' North and longitudes 5° 00' and 9° 00' West (Figure 1A). The study area comprises two main regions with very different hydrogeological, land use and relief contexts: the phosphate plateau and the Tadla plain (Figure 1B). The climate in the study area is semi-arid to arid with a continental character. Precipitation is irregularly distributed in time and space within the basin. There is a dry season from May to October and a wet season from November to April. During the dry season the minimum inter-annual monthly rainfall averages are as low as 0.8 mm, however the average maximum is over 72 mm during the winter. From north to south, the increase in altitude is accompanied by an increase in precipitation amount and a decrease in temperature.

In the OER River Basin, groundwater aquifers are subject to constant monitoring, encompassing both qualitative and quantitative variables. This region benefits from a particularly extensive historical groundwater dataset, unlike other parts of Morocco where data is sparse. The comprehensive database available in the OER basin owes its central role as one of Morocco's agricultural hubs (Ouatiki et al., 2019; Khellouk et al., 2021; Lionboui et al., 2021). Five wells are investigated in this study, wells P1 (Figure 2A) and P2 (Figure 2B) are located in the Tadla Plain which is characterized by a generally regular topography at an average altitude of 400 m and a gentle slope toward the west (Figures 1B, C). The Tadla plain is divided by the OER wadi into two hydrogeological zones: the Beni Amir aquifer to the north and the Beni Moussa aquifer to the south, in which the P1 and P2 wells are located (Figure 1D). The water table of the Tadla plain is formed by Villafranchian and Quaternary lacustrine limestone arranged in discontinuous lenses separated by low permeability formations. The agricultural areas in this plain cover 320 000 ha that include 124 000 ha as irrigated areas (Lionboui et al., 2021).

Wells P3, P4 (Figure 2C) and P5 (Figure 2D) are located in the phosphate plateau (Figure 1B). These wells belong to the Turonian deep water table which is the most important one in the Tadla (Figure 1D), circulating in limestone and dolomitic limestone formations. However, these three wells are located in distinct hydrological contexts. Wells P3 and P4 are located in the floodplain of the Wadi Zem, which dissipates further south into the alluvium without merging with the OER River. Well P5, shows the greatest depth to the water table (average depth = 72 m), is located further north in the center of the phosphate plateau in a semi-arid environment near the phosphate mines (Figure 2D). Table 1 presents the descriptive statistics of groundwater level observed for each well.

Predicting groundwater levels using remote sensing data from several climate variables for several wells with very different dynamics necessitates an understanding of the data that goes into the model. As a result, a variable importance analysis can help us understand the relationship between the descriptive variables and the target variable, and determine which features are not relevant to the model in each well (Sahu et al., 2020; Sharafati et al., 2020; Khaire and Dhanalakshmi, 2022).

Traditionally, important variables are chosen using correlation coefficient analysis, with the selected variables being those with the highest correlations with groundwater levels (Derbela and Nouiri, 2020; Band et al., 2021; El Bilali et al., 2021; Vu et al., 2021). However, strongly correlated input variables can have a negative impact on the performance of a modeling algorithm. When using multiple correlated variables, the effect of collinearity occurs, resulting in unstable and unreliable model results. To avoid this, choosing the best single variable eliminates the need to use multiple variables and their potential interactions, preventing over-learning and allowing the model to generalize better to new, unknown data (Dormann et al., 2013; Anh et al., 2023).

Thus, before building the model, it is critical to examine the correlations between the input variables and the variable to be predicted in order to assess their relationships in a fairly straightforward and descriptive manner. However, it is recommended to use input feature selection techniques, which calculate a score for all input features of a given model that simply represent the importance of each variable (Sahu et al., 2020). A higher score means that the variable will have a greater effect on the model that is used to predict the groundwater level. In this study we used the Extreme Gradient Boosting (XGBoost) method, which constructs boosted trees to obtain the scores of the variables, indicating their importance for the modeling model (Hastie et al., 2009; Breiman et al., 2017; Zheng et al., 2017). The scores provided by XGBoost indicate the usefulness or value of each variable in building the boosted decision trees in the model. The more a variable is used to make key decisions with the decision trees, the higher its relative importance (Friedman, 2002; Zheng et al., 2017). XGBoost counts the importance by “gain”, which is the main reference factor of the variables importance in the tree branches (Hastie et al., 2009; Breiman et al., 2017; Zheng et al., 2017):

With T as a decision tree, ω_l importance score for each predictor variable, M is the additive number of trees et τt2 is the estimated maximum improvement of the variable in the risk of squared error over that of a constant fit.

Specifically designed for handling sequential data, the LSTM model is highly suitable for forecasting both qualitative and quantitative groundwater variables, which often exhibit temporal dependencies. One of the main advantages of the LSTM model is its ability to retain information on events that occurred at the start of the time series, even as it processes new data. This feature is crucial for capturing long-term dependencies and patterns present in time series data. This is especially relevant for groundwater variables, as they are influenced by climatic seasonality and other phenomena with return times spanning several years or even longer periods.

The mathematical representation of a RNN is defined by the equation (Hochreiter and Schmidhuber, 1997; Gers et al., 2000):

Where, h_t is the initial state at time t, h_t−1 is the hidden state of the cell at time t-1, X is the input vector with the sequence, σ and the activation function of TanH, b_h and b_o bias vector, U is the weight (vector) for the hidden layer, V is the weight (vector) for the output layer, W is the same weight vector for different time steps, and y the output vector with the sequence.

The RNN uses the time backpropagation approach, which is a type of gradient descent algorithm that is used to update the weight and minimize errors using a derived chain rule. However, the way the RNN is structured does not allow for effective long-term dependency treatment, as its learning process leads to the vanishing gradients (Hochreiter, 1998). This problem is due to the fact that during the backpropagation the weights are the same for all time steps, as the weight is updated continuously, so that the gradient becomes either too weak or too strong with the updates. One of the solutions to this problem is the use of long-term memory (LSTM) (Graves, 2012). The LSTM was built to avoid the problem of long-term dependencies, by having an additional feature compared to RNNs, which is called memory or the internal stat, which is specifically designed to store information over long periods of time.

This memory deals with the evanescent gradient problem by implementing a structure consisting of four elements, namely the CEC cell (the Constant Error Carrousel) and three types of gates: input gates, output gates and forget gates, which ensure the preservation of previous information with a stable gradient calculation. The information is processed by a sequential calculation using the following equations iteratively along the time series:

where w_i, w_f and w_o denote the matrix of weights of the input, forget and output gates at the input, respectively. Similarly, U_i, U_f, and U_o denote the weights matrix of the input, forget and output gates to the hidden, respectively. b_i, b_f, and b_o denote the bias vectors of the input, oblivion and output gates, respectively. σ is a non-linear activation function by elements: logistic sigmoid. i_t, f_t, o_t and C_t are the input, forget, output gates and state vectors of the cell at time t, respectively, which are all the same size as the output vector of the cell ht. The element-wise multiplication of two vectors is denoted with ⊗ (Zhang et al., 2018).

The LSTM is designed to handle large databases. However, with a large number of parameters, LSTM can easily be over-fitted, especially when data is limited. To solve this problem, the 'Dropouts' method is used where randomly selected neurons are ignored during training. They are randomly 'dropped'. This avoids over-fitting the model.

The three LSTM models developed in this study consist of two layers connected with the dense layer. It should be noted that the exact number of LSTM layers to be used as hidden layers in order to obtain the optimized results remained unknown. However, in general, two layers have proven to be sufficient to detect the complex interrelationships between the variable to be predicted and the input data (Zhang et al., 2018).

In deep neural networks models there are a set of parameters called hyperparameters that must be defined before the training process. Hyperparameters are frequently set manually or at random; however, this type of setting remains difficult to perform and becomes extremely tedious, when network architecture is complex. Several hyperparameter tuning algorithm are therefore used to identify the optimal combination of model hyperparameters in an efficient, fast and reliable manner (Bowes et al., 2019). The GridsearchCV method, a traditional method of hyperparameter optimization, was used in this study, which is simply an exhaustive search in a specified subset of a learning algorithm's hyperparameter space (Chang et al., 2022; Elzain et al., 2022; Anh et al., 2023). In other words, it iterates through all of the hyperparameters entered into the parameter grid to find the best parameter combination. Furthermore, GridsearchCV performance is often measured using cross validation on the training set (Ghawi and Pfeffer, 2019). In this study, the LSTM hyperparameters tuned are the number of neurones, batsh size, learning rate, dropout rate and number of epochs (Table 3). These parameters have a significant impact on prediction performance and are considered the most important to optimize (Bowes et al., 2019). The dropout rate helps avoid overfitting during the training process by skipping randomly chosen neurons. The batsh size is the number of samples introduced to the model in order to distinguish common features in the training data set and prior to performing a weight update. The number of epochs defines the number of complete iterations of the training data set to be executed.

Three prediction performance evaluation indices were used, namely, mean absolute error (MAE), root mean squared error (RMSE) and the R². The MAE measures the average magnitude of errors in a set of predictions, regardless of their direction. It is the average, over the training or test sample, of the absolute differences between the prediction and the actual observation, with all individual differences having the same weight. It is less sensitive to outliers, making it more suitable for cases where extreme values might be present in the dataset. The RMSE is a quadratic scoring rule that also measures the average magnitude of the error. It is the square root of the average of the squared differences between the prediction and the actual observation. RMSE ensures that outliers and large errors have a substantial impact on the evaluation, providing a more robust assessment of the model's performance. R² is an indicator used in data analysis to judge the quality of a model and to measure the degree of accuracy in reproducing observed values, R² values are between (0, 1) where an R² score close to 1 represents optimal model prediction. The selection of performance indices in this study is based on their widespread use and established importance in hydrogeological modeling, including groundwater level prediction (Rajaee et al., 2019).

where y_j is the measured value at time j, yj̄ is the average of y_j (j= 1,..., N) and i yj^ is the predicted value at time j.

The groundwater level shows very distinct correlation values with the input variables and within the five wells (Figure 4). In well P1, the groundwater level shows a weak correlation with NDVI and evapotranspiration; on the other hand it shows a slight positive correlation with precipitation (0.3) and negative correlation with land surface temperature (−0.3). The groundwater level in well P2 has low correlations with the four input variables, the highest correlation is 0.21 obtained with evapotranspiration. Well P3 shows relatively high positive correlation values of groundwater level with NDVI (0.43) and evapotranspiration (0.39), while precipitation and land surface temperature show no significant correlations. The level of groundwater in well P4 shows slight correlation values with NDVI (0.28) and evapotranspiration (0.3). In well P5, the level of groundwater has weak correlations with all four input variables, and the correlation value that is slightly high (0.28) is observed with evapotranspiration (Figure 4). Moreover, it can be observed that within the five wells. the input variables exhibit strong correlations (both positive and negative) among themselves.

The benefit of using XGBoost is that once the boosted trees are constructed, it is quite simple to retrieve the importance scores for each variable. The XGBoost result shows that land surface temperature is the most important input variable for groundwater level in the wells P1 and P5. Furthermore, precipitation, NDVI, and evapotranspiration are the important input variables, respectively, for P2, P3, and P4 (Figure 5). Thus, based on these results, the most important input variable for each well is used as a predictor variable in the LSTM-XGB-GS model.

The outcomes of hyperparameter tuning reveal a notable distinction between the LSTM-XGB-GS model and the LSTM-GS model. Conversely, for each individual model, the differences are relatively minor across the wells (Table 4). In the case of the LSTM-XGB-GS model, the optimal number of neurons is 70, except for well P2 where it is 50. The preferred batch size is 48 for wells P1, P2, and P3, 24 for well P4, and 12 for well P5. Across most wells, the favored learning rate is 0.001, with the exception of well P5 where it is 0.01. The highest dropout rate (0.3) is utilized for wells P1, P4, and P5, while for wells P2 and P3, the optimal dropout rate is 0.2. The optimal number of epochs for all wells stands at 200.

For the LSTM-GS model, the most effective number of neurons is 120 for wells P1 and P2, and 70 for wells P3, P4, and P5. The optimal batch size across all wells, except for well P2, is 48. Well P2 features a batch size of 24. Learning rates of 0.001 are consistent across all wells. The highest dropout rate (0.3) is employed for wells P3 and P4, while for wells P2 and P5, the ideal dropout rate is 0.2, and for P1 it's 0.1. The optimum number of epochs remains consistent at 200 for all wells.

These results highlight the sensitivity of the LSTM model to the quantity and selection of input data, a factor that can have an impact on the results of individual models and wells.

The performance of the models, in terms of three statistical indices during both the training and test periods, is presented in Table 4. The optimized model with the best input feature (LSTM-XGB-GS) demonstrated significantly superior performance compared to both the optimized model with all input features (LSTM-GS) and the standalone model (LSTM) across all five wells.

The R² values of the LSTM-XGB-GS model during the training period ranged from 0.95 to 0.99, while the corresponding values for LSTM-GS and LSTM were found to be in the range of 0.89 to 0.99 and 0.67 to 0.84, respectively. Upon analyzing the computed RMSE and MAE values, it was observed that the LSTM-XGB-GS model achieved the minimum values for all five studied wells. The RMSE ranged from 0.03 to 0.68 meters, while the MAE ranged from 0.03 to 0.51 meters. Specifically, the wells in the irrigated areas P1 and P2, as well as those in the floodplain P3 and P4, exhibited acceptable values of RMSE (0.33 m, 0.57 m, 0.68 m, and 0.25 m, respectively) and MAE (0.22 m, 0.45 m, 0.51 m, and 0.20 m, respectively). Notably, the P5 well located on the phosphate plateau demonstrated the lowest values of RMSE (0.03 m) and MAE (0.02 m). Furthermore, it is worth mentioning that the LSTM model showed the lowest modeling performance in terms of RMSE and MAE during the training period.

During the test period, the LSTM-XGB-GS model demonstrated good to excellent performance, achieving an R² between 0.72 and 0.99. On the other hand, the LSTM-GS model exhibited varying performance, ranging from poor to excellent depending on the well location, with an R² between 0.27 and 0.97. Meanwhile, the LSTM standalone model showed the poorest performance across all five wells, attaining an R² between 0.09 and 0.96.

As for the RMSE and MAE indices, the LSTM-XGB-GS model consistently displayed the lowest values among the three models for all five wells. In particular, when utilizing MODIS LST as the input variable for well P1, the LSTM-XGB-GS model yielded RMSE and MAE values of 0.20 m and 0.15 m, respectively (Table 5). For well P2, employing IMERG-GPM precipitation as the input, the LSTM-XGB-GS model resulted in RMSE and MAE values of 0.66 m and 0.58 m, respectively. For well P3, using MODIS NDVI as the input, the model provided acceptable RMSE and MAE values of 0.42 m and 0.39 m, respectively. Similarly, well P4, with MODIS evapotranspiration as the input, demonstrated acceptable RMSE and MAE values of 0.35 m and 0.29 m, respectively. Notably, the phosphate plateau well P5, utilizing MODIS LST, exhibited the best prediction performance during the test period, with both RMSE and MAE values at 0.02 m.

The scatter plots of groundwater levels observed and predicted using the three models during the training and test periods are shown in Figure 6. In general, the results show that the groundwater level predicted by the three models closely reflect, or even align with, the trend in the observed data (the black line) during the training period. This alignment is also supported by the R² values, which exceed 0.7 for all models over the five wells (Table 5). It should be noted, however, that the scatterplot associated with the LSTM-XGB-GS model adheres more closely to the trend in the observed data than the scatterplots for the LSTM-GS model and the stand-alone LSTM model. The latter two show slightly more dispersed scatterplots, with some values far from the trend of the observed data.

The results from the test period give a clearer overview of the prediction performance of the models compared with the results from the training period. Most notably, the scatter plot for the LSTM-XGB-GS model shows better alignment with the trend in the observed data, with points showing reduced dispersion and closer proximity to the trend line. Conversely, the scatterplots of the LSTM-GS and standalone LSTM models are less closely aligned with the trend of the observed data, with points showing increased scatter. These observations are also supported by the R² values, where the LSTM-XGB-GS model has the highest values for all five wells (Table 5).

Time series plots of monthly variations in groundwater levels observed and predicted over the test period by the three models is shown in Figure 7. The plots show that all three models effectively capture seasonal fluctuations and variations in groundwater levels in each well. However, the LSTM-XGB-GS model shows superior performance in reproducing the groundwater levels observed. It shows remarkable alignment with minimal overestimation or underestimation. In particular, for wells P1 and P4, the LSTM-XGB-GS model almost perfectly reproduces the observed data. On the other hand, the LSTM-GS and standalone LSTM models, although capable of capturing seasonal variations in groundwater levels, show significant overestimates and underestimates. This discrepancy is most pronounced, although slightly attenuated, in well P4, where the values of the latter two models are very close to the observed data.

Finally, the ability of the three models to reproduce the probability distribution of observed groundwater level data during the training and test periods for all five wells was visually evaluated by preparing violin plots. The violin plots displaying the observed and predicted groundwater level data by the models during the training period are presented in Figure 8. The violin plots demonstrate that all three models have a similar overall shape to the observed groundwater level data, with their median values closely aligned, suggesting similarity in the distribution of predicted and observed groundwater level data. However, it was also observed that the violin plots of the LSTM-XGB-GS model appear slightly narrower compared to the other models and the observed data for wells P1, P2, P3, and P5. This indicates a slightly lower variance in the predictions. Additionally, the LSTM-XGB-GS model shows a slight tendency to underestimate maximum values and overestimate minimum values for these wells. On the other hand, the LSTM-GS model exhibits a high degree of similarity with the observed data for wells P1, P2, and P4, without underestimating or overestimating extreme values. In contrast, the LSTM model shows some divergence from the observed data, with underestimation and overestimation of maximum and minimum values, except for well P4, where the LSTM model exhibits a distribution similar to the observed values.

The violin plots of the test period show that the LSTM-XGB-GS model exhibits the highest similarity to the observed data compared to the other two models (Figure 8). For wells P1, P4, and P5, the LSTM-XGB-GS model displays a similar shape and range to the observed data, with a median value at the same level and less underestimation and overestimation of the maximum and minimum values compared to the other models. For well P2, all three models overestimate the median and extreme values compared to the observed data, but the LSTM-XGB-GS model shows the lowest divergence and a shape similar to the observed data. For well P3, all three models display a wider shape than the observed data; however, the LSTM-XGB-GS model exhibits the least divergence by showing a median value at the same level as the observed data, with the lowest overestimation and underestimation of extreme values.

The performances of the different models were also assessed using the Taylor diagram. The Taylor diagram graphically compares two statistics, namely the standard deviation and correlation, to provide a reliable and more direct evaluation of the relative performance of the different models. The results for the five wells are presented in Figure 9. In the Taylor diagram, the dashed black circle represents the observed groundwater level, and the black star denotes the reference point for the observed data. A model is considered proficient if it is close to the reference point. During the training period (Figure 9), the LSTM-XGB-GS model exhibited the best performance across all five wells. On the other hand, the standalone LSTM model showed the weakest performance in predicting the observed data. However, for well P5, the performances of the LSTM-XGB-GS and LSTM-GS models were roughly equal. For the results of the test period (Figure 9), the LSTM-XGB-GS model also displayed the best performance in predicting the observed data during the test period. Conversely, the standalone LSTM model demonstrated better performance than the LSTM-GS model for wells P2 and P3. However, for the P5 well, the performance of the LSTM-XGB-GS and LSTM-GS models was fairly equal compared to the standalone LSTM model, whose performance was much lower.

The good performance achieved by the LSTM-XGB-GS model in this study clearly emphasizes the importance and impact of identifying the best input feature and fine-tuning the hyperparameters for groundwater level prediction. The hydrodynamic characteristics of aquifers and the groundwater levels fluctuation in monitoring wells exhibit significant variations influenced by several factors (Rinderer et al., 2016; Iqbal et al., 2020; Band et al., 2021), primarily encompassing hydro-geodynamic and climatic conditions, as well as land use and land cover. To achieve enhanced model performance, it becomes imperative to pinpoint the variable that reflects the dominant factor that explains most the variation of the groundwater level in each well. Moreover, the optimization of the LSTM model's hyperparameters using the GridSearchCV method enables a more effective capture of patterns and relationships between groundwater level and the selected remote sensing variable.

However, it is worth noting that the number of observed values in the five wells averages 134, which represents a relatively small dataset for ensuring the LSTM model's optimal performance. LSTM model initially designed for learning on large datasets, commonly known as “big data,” shows remarkable performance. However, the model's formidable learning capacity can present difficulties when applied to a limited dataset of groundwater level observations. In such cases, the LSTM model becomes susceptible to over-fitting, which can compromise its predictive accuracy (Baek and Kim, 2018). This overfitting issue is evident in the case of the LSTM-GS model. By employing all available remote sensing variables without specifically selecting the most informative ones for each well, the LSTM-GS model becomes overly specialized in capturing the training data, leading to decreased performance when confronted with new, unseen data during the testing phase (Table 5). Furthermore, the poor performance of the LSTM standalone model was expected, since it uses all available input variables without the hyperparameters optimization. As a result, this model is prone to overfitting and its ability to accurately capture fluctuations in groundwater levels is limited. This has been observed in several other studies, where standalone ANN models exhibited the weakest performance in predicting groundwater variables (Rahman et al., 2020; Cui et al., 2022; Zhang et al., 2023).

Despite the high fluctuation of groundwater levels in well P1 (Standard deviation > 2 m) (Table 2), the integration of the best input feature (MODIS LST) using XGBoost and the hyperparameter optimization of the LSTM-XGB-GS model led to a remarkably accurate prediction of its groundwater levels during the test period (Figure 7). This good prediction performance can be attributed to the relationship between irrigation and land surface temperature, as evidenced by several studies (Bright et al., 2017; Liu et al., 2019; Yang et al., 2020; Zhang Z. et al., 2022). These studies, demonstrated that in arid regions, irrigation has been shown to substantially alter the land surface temperature, resulting in lower daytime land surface temperature in irrigated areas compared to adjacent non-irrigated regions (Yang et al., 2020). Given that irrigation in the Beni Moussa aquifer heavily relies on pumping water from the aquifer (Lionboui et al., 2021), a direct relationship may exist between the dynamics of groundwater levels and the variations in land surface temperature. Consequently, this interplay between land surface temperature and groundwater levels could explain LSTM-XGB-GS model's good performance in simulating groundwater levels for well P1.

The LSTM-XGB-GS prediction results for wells P2 and P3 are generally acceptable and reliable. However, a noticeable decrease in model prediction performance is observed when compared to their training period or to the other wells (Table 5). This decrease in performance can be attributed to several factors. To begin with, the correlation coefficient for groundwater level in well P2 does not exhibit a distinct correlation with any of the four input variables (Figure 4). XGBoost indicates that IMERG-GPM precipitation is ranked first, but it is not isolated from the other variables since they all hold roughly the same level of importance (Figure 5). This observation might be explained by the significant unsaturated zone at well P2, where the depth to the water table is 36.78 m (Table 2). At certain depths, surface climatic conditions such as land surface temperature and evapotranspiration may have only a minor influence on short-term groundwater fluctuations. Conversely, changes in the precipitation regime can have a long-term impact on groundwater levels. Nevertheless, the effect of precipitation on groundwater is not always immediate; it can vary from month to month and is heavily dependent on the aquifer's hydrological conditions (Jan et al., 2013; Grinevskii, 2014; Cai and Ofterdinger, 2016; Qi et al., 2018; Mogaji and Lim, 2020). Moreover, the IMERG-GPM product's low spatial resolution (10 km) might hinder its ability to accurately capture the fine-scale spatial variability in precipitation, which is crucial in this semi-arid irrigated plain (Ouatiki et al., 2019). Additionally, numerous evaluations of the IMERG-GPM product have revealed certain biases and estimation errors, leading to significant over- and underestimations of the duration and intensity of precipitation (Ramsauer et al., 2018; Ouatiki et al., 2019; Kazamias et al., 2022; Pradhan et al., 2022).

In addition, NDVI exhibits a strong positive correlation with well P3 in the floodplain (Figure 4), supporting the findings of XGBoost, where NDVI is identified as the most important input feature. This good correlation/importance can be attributed to the temporal variation of NDVI in floodplain areas, where it has been shown to effectively capture seasonal and inter-annual variability in vegetation dynamics and the hydro-sedimentological regime along and across the floodplain (Shrestha et al., 2013; Powell et al., 2014; Marchetti et al., 2016). Moreover, NDVI values provide insights into the relationship between groundwater and plant growth. In cases where a high NDVI area is experiencing drought, it indicates that vegetation is consuming more groundwater, which potentially impact groundwater levels variation (Hassan et al., 2019; Hadri et al., 2021; Lebrini et al., 2021; Moumane et al., 2021; Sajjad et al., 2022; Zhang Z. et al., 2022). However, given the positioning of well P3 in the floodplain, it is plausible that the influence of the nearby wadi contributes significantly to the observed decline in model performance between the training and test periods. The dynamic interactions that occur between the wadi and the aquifer have the potential to elucidate a substantial part of the fluctuations observed in the water table, as highlighted by previous studies (May et al., 2011; Wang et al., 2014; Kazamias et al., 2022). To establish this hypothesis with greater certainty, a more conclusive assessment of the modeling performance would require the integration of flow data from the wadi. Unfortunately, such data were not available in the present case.

The strong relationship between groundwater levels in well P4 and MODIS evapotranspiration (Figures 4, 5) clearly enabled accurate prediction of groundwater level fluctuations in this well. This can be explained by the fact that the floodplains in this region serve as designated irrigation areas, used to extract water from surface and groundwater reservoirs. This is particularly important during prolonged periods of limited rainfall, when the extracted water supports plant growth, resulting in increased evapotranspiration (Zhang et al., 2018; Hssaisoune et al., 2020; Ahmed et al., 2021; Malakar et al., 2021). In this context, the measurement of evapotranspiration is of key importance and is considered an essential climatic parameter for vigilant monitoring and effective management of irrigation operations and floodplain ecosystems. Many research studies have seamlessly integrated evapotranspiration data acquired from remote sensing resources, including MODIS data, to assess irrigation needs and decipher hydrological models (Droogers et al., 2010; Van Houdt et al., 2020; Vogels et al., 2020).

It's important to acknowledge that groundwater resources in irrigated and floodplain areas bear a substantial impact from intensive pumping activities (Zhang et al., 2018; Malakar et al., 2021). This phenomenon accounts for the substantial fluctuations in groundwater levels within irrigated (wells P1 and P2) and floodplain (wells P3 and P4) regions, as evidenced by their high standard deviations (> 2m) (Table 2). The influence of farming practices, particularly irrigation and pumping, is even more pronounced in semi-arid and arid climates, where high pressure and reliance on groundwater are prevalent during water scarcity periods (Wang et al., 2014; Kirby et al., 2015; Mao et al., 2017; Pulido-Bosch et al., 2018; Salem et al., 2018; Tweed et al., 2018; Cavelan et al., 2022). In Morocco, the combination of meteorological drought driven by reduced rainfall and elevated temperatures has escalated irrigation water demand. Consequently, irrigation sources often shift from surface water to groundwater during these periods. Such challenges have been closely monitored and evaluated across the OER Basin and its irrigated areas, where human activities have exerted significant pressure on water resources, leading to notable declines in groundwater levels averaging from 71.9 to 148.8 cm/year (Ahmed et al., 2021) and ~30 m over the past three decades (Hssaisoune et al., 2020). As a consequence, the model prediction performances highlight the importance of precise knowledge regarding groundwater pumping rates and volumes. Such information proves critical for enhancing future groundwater modeling and forecasting in arid and semi-arid climates characterized by intense irrigation and pumping activities. These insights enable effective training of models to adeptly capture the abrupt shifts in hydraulic head that result from these activities. However, the challenge often lies in obtaining accurate data concerning the quantity of water pumped from individual wells (Trichakis et al., 2011).

Well P5, located in the phosphate plateau area, achieves the best modeling results during the training and testing periods in both LSTM-XGB-GS and LSTM-GS models (Table 5), and the simulated groundwater level maintains the same trend, distribution and quasi-stationary aspect of the observed groundwater level variation. This result was anticipated given that well P5 have the lowest groundwater level fluctuation (standard deviation of 0.11 m) and the deepest water table depth (mean depth of 72 m) (Table 2). These aspects suggest that groundwater levels in well P5 remain relatively unaffected by climatic conditions or surface activities such as irrigation or pumping, as demonstrated by the low correlation values of the groundwater level with the four input variables (Figure 4). In general, groundwater level prediction studies that have shown excellent performance (RMSEs and MAEs < 0.1 m), have often been applied to groundwater levels that did not show large groundwater fluctuations, or they showed a long-term trend, over the whole time series (Bowes et al., 2019; Band et al., 2021; Ghazi et al., 2021; Wunsch et al., 2021). This near-stationary or regular dynamics allow ANN models, and in particular recurrent neural networks such as the LSTM, to adequately predict the time series, as the output variable becomes increasingly easy to predict as the time series advances. in such conditions, it would be possible to predict the groundwater level of the well P5 by adopting only a univariate prediction approach, which relies solely on the groundwater level observed data as input to the LSTM optimized model (Raghavendra and Deka, 2016; Mohanasundaram et al., 2019; Roy et al., 2021; Sarma and Singh, 2022).