1. Introduction
Global warming is threatening surface water supply in many parts of the word, particularly in arid regions. In these regions, groundwater, a viable alternative and an important source of freshwater, is often limited [
1,
2]. Considering that groundwater level (GWL) is an indicator of groundwater availability at any given time, monitoring GWLs provides significant insights into the dynamics of recharge and withdrawals and how they influence the long-term availability of groundwater. In arid regions, this can be challenging due to the inadequate distribution of groundwater wells and the presence of spatial and temporal data gaps in monitoring records [
3]. Therefore, accurate and reliable predictive tools are essential for supporting the sustainable management of groundwater in these areas [
3].
The relationships between GWL fluctuations and explanatory variables are generally complex and nonlinear [
4,
5,
6,
7]. However, machine learning (ML) algorithms can effectively learn these relationships. One of such algorithms is the support vector machine (SVM) for regression purposes (SVR), especially when embedded with the radial basis function (RBF) kernel [
5,
8]. Another algorithm capable of learning these relationships is random forest (RF) [
9], which is the most employed ML technique for GWL predictions [
10]. Both SVMs and RF are known to give accurate results [
10,
11].
Several studies have shown the RBF–SVR model to outperform various other techniques, such as artificial neural network (ANN) [
12,
13,
14], radial basis function neural network (RBF–NN) [
15], the autoregressive integrated moving average (ARIMA) model [
16], RF [
8], and the gradient boosting mechanism (GBM) [
8]. These studies also attribute the success of SVMs to their strong prediction capability and the ability to generalize well to unseen data.
Likewise, various studies have reported success using RF for GWL prediction. For example, it outperformed K-nearest neighbor (KNN), ANN, and SVR based on root mean square error (RMSE) values during testing [
17]; ANN and SVR based on R-squared (R
2), mean absolute error (MAE), and RMSE values in training and R
2 and MAE values in testing [
11]; multilinear regression (MLR) based on R
2, MAE, and RMSE values in both training and testing [
18]; decision trees (DTs) and SVRs based on its R
2 and RMSE values in testing [
19]; and the XGBoost regressor based on its MAE and RMSE values in testing [
20].
In addition to traditional ML methods, geostatistical interpolation (GI) techniques are also commonly used for GWL prediction. Kriging is the most utilized GI technique [
21,
22], where the measured values of a variable (GWL, in this case) at specific locations are used to make predictions at unmeasured areas. It relies on the correlation between the measured values as a function of distance, depicted by a semivariogram, to assign weights that describe the contribution of each measured point to the prediction at unmeasured locations [
23]. Kriging presents an important advantage over other interpolation methods in being able to quantify and minimize prediction uncertainties [
23].
Classical kriging (CK), the traditional form of kriging, relies on a single semivariogram assumed to be the true representation of the measured data. In contrast, empirical Bayesian kriging (EBK), an advanced kriging approach, incorporates multiple semivariograms to account for the uncertainties associated with estimating a single semivariogram [
23]. Thus, EBK is a more robust kriging algorithm [
23], and studies such as Bouhout et al. (2022) [
24] and Hussain et al. (2016) [
25] have demonstrated its superior performance in GWL prediction applications.
Our review of the existing literature suggests that RBF–SVR, RF, and EBK models are some of the most effective GWL prediction tools. To enhance the prediction accuracy, we propose an approach that integrates all three techniques to predict the monthly GWLAs across the state of Arizona (representative of arid/semi-arid systems) between January 2010 and December 2019, using remotely sensed predictor variables.
2. Materials and Methods
Natural groundwater recharge in arid regions is often limited, raising concerns about the sustained availability of freshwater in these regions. Climate change exacerbates these issues by intensifying the hydrologic cycle, resulting in increased evapotranspiration rates and a reduction in the soil moisture available to recharge groundwater systems [
26,
27]. In a recent study, Zowam et al. (2023) [
27] quantified terrestrial water cycle intensity (WCI) changes across the contiguous United States (CONUS) attributable to climate change and showed that the state of Arizona might be experiencing much higher relative WCI rates on average than other arid regions in the CONUS. These factors underscore the need for accurate prediction of the GWL in such regions to continue to effectively manage the potentially limited freshwater resources therein.
2.1. The Study Area and Target Variable
The study area, covering about 114,000 mi
2, is located in the southwest U.S. (
Figure 1). Its surficial geology is characterized primarily by unconsolidated deposits in the south/southwest–northwest corner of the state, with various rocks dominating the other regions [
28]. The consolidated rocks are mainly sedimentary and extrusive igneous (volcanic) rocks and constitute the mountain ranges that border the basins filled with unconsolidated materials [
29]. These rocks are the main sources of sedimentary materials that fill the basin and have very low permeability and groundwater flow rates [
29]. Many communities depend solely on groundwater to meet their water needs, which has led to a long history of over-extraction in many parts of the state [
30]. In the late 1980s, the Colorado river arrived in Arizona and eased some of the pressure on groundwater to meet these needs, but the prolonged drought in the Colorado river basin, coupled with projected warming temperatures, is expected to not only reduce the surface water availability in the state but also further stress the aquifers in the region [
30].
Daily GWL data from 59 monitoring wells were downloaded from the National Groundwater Monitoring Network (NGWMN) portal (
https://cida.usgs.gov/ngwmn/index.jsp, accessed on 29 January 2023). Among these wells, 38 were drilled into unconsolidated sand and gravel aquifers and are managed by the Arizona Department of Water Resources (ADWR), while the remaining 21 were drilled into consolidated rock aquifers (
Figure 1). Three of the rock aquifer wells used in this study are maintained by the U.S. Geological Survey (USGS), and the rest are managed by the ADWR. The depths of these wells varied and ranged from 89 to 1600 feet below ground surface (bgs) for the sand and gravel aquifers, and about 25 to 851 feet bgs for the rock aquifers, and the aquifers were predominantly unconfined. Daily GWL measurements were aggregated into monthly averages from January 2010 to December 2019.
Missing data is the most common challenge in real-world ML applications [
31], and various methods exist to address this issue. The simplest of these methods is mean imputation (MI) [
31]. In MI, the mean values of available observations are used to fill in missing observations, which has proven to work well with small variance distributions [
32,
33], i.e., distributions with a coefficient of variation (CV) less than 10% [
34]. Therefore, with CV values ranging from 0.01 to 0.91% for the unconsolidated aquifers and 0.001 to 0.21% for the rock aquifers, the MI method was ideal for the target variable (monthly GWL). Missing monthly values were replaced with the annual average GWL for the given year.
To compute monthly GWLA, we calculated the mean GWL at each well and subtracted the monthly measurements from this mean value.
2.2. Predictor Variables
The selection of the input variables was informed by the Seyoum et al. (2019) [
35] study as well as established hydrogeological principles. The initial variables included precipitation, soil moisture, evapotranspiration, land surface temperature, vegetation index, curve number, saturated hydraulic conductivity, and groundwater storage anomalies (
Figure 2).
2.2.1. Precipitation (P)
Precipitation is the most important hydrological variable for predicting groundwater recharge [
36], thus playing a crucial role in determining GWL. In the last two decades, satellite-based precipitation measurement techniques have seen significant advancements, and the Global Precipitation Measurement Mission (GPM) using the Integrated Multi-Satellite Retrievals for the GPM (IMERG) algorithm stands out as one of the best alternatives to ground-based measurements [
27,
37]. In particular, the GPM mission demonstrates significant potential to mitigate the challenges associated with estimating precipitation in arid regions [
38]. In addition, under light rainfall conditions (typical of these regions), IMERG tends to produce lower detection errors and generally more accurate estimates [
39]. The final run of the IMERG system provides the most accurate precipitation measurements, making it ideal for research purposes [
40]. A monthly IMERG (final run) dataset with a 0.1° × 0.1° grid resolution was downloaded from the National Aeronautics and Space Administration (NASA) data portal for January 2010 to December 2019 (
https://gpm.nasa.gov/data/directory, accessed on 11 May 2021).
2.2.2. Soil Moisture (SM)
SM and GWL tend to demonstrate a negative relationship [
41,
42,
43], which can be much stronger for shallow groundwater [
41]. For this study, we utilized a research product that integrated measurement efforts from both the European Space Agency (ESA) and NASA (
https://doi.org/10.1594/PANGAEA.940409, accessed on 12 February 2023). The dataset was generated by downscaling ESA’s Climate Change Initiative (CCI) data using NASA’s Soil Moisture Active and Passive (SMAP) data [
44]. The global dataset had a grid size of 9 km (~0.08° × 0.08°) and a daily temporal resolution covering a 43-year period from 1978 to 2020 [
44]. We converted the daily rasters into monthly averages spanning our study period and extracted the monthly data for the study area.
2.2.3. Evapotranspiration (ET)
The combined effects of liquid water losses from soil surfaces (evaporation) and its uptake by plants (transpiration) limit the amount that infiltrates the ground, thereby affecting GWLs. We utilized a global dataset with a fine grid size of 1 km (~0.01° × 0.01°) and a monthly temporal resolution (
https://doi.org/10.7910/DVN/ZGOUED, accessed on 3 July 2021). The dataset was obtained by synthesizing the best-performing satellite ET products following validation against flux eddy covariance ET and performed better than local products across the United States, China, and the continent of Africa [
45].
2.2.4. Land Surface Temperature (LST)
LST tends to exhibit a positive relationship with GWL, which is more pronounced for shallow groundwater [
41]. This study utilized a gap-filled, continuous LST dataset generated by filling in missing pixels in the Moderate-Resolution Imaging Spectroradiometer (MODIS) 1 km resolution LST daily product [
46]. Daytime (1:30 PM) and nighttime (1:30 AM) global datasets were downloaded from the Iowa State University research repository (
https://doi.org/10.25380/iastate.c.5078492.v3, accessed on 12 February 2023). The downloaded rasters were converted into monthly averages, and the final LST dataset was obtained by averaging the daytime and nighttime monthly estimates.
2.2.5. Vegetation Index (VI)
The presence of vegetation may affect groundwater recharge (and GWLs) in various ways, including slowing down runoff and enhancing ET through transpiration. VI values are unitless and help to visualize the locations and relative abundance of green vegetation. A monthly dataset with a 0.1° × 0.1° grid resolution was accessed and downloaded from the NASA Earth Observations (NEO) data portal (
https://neo.gsfc.nasa.gov/view.php?datasetId=MOD_NDVI_M&date=2010-12-01, accessed on 12 February 2023).
2.2.6. Curve Number (CN) and Runoff Depth (R)
CN is a dimensionless parameter that characterizes the runoff potential of a surface. It is influenced by land use, soil characteristics, and the antecedent moisture conditions of soils. Lower numbers typically correspond to permeable soils with high infiltration rates, while higher numbers are associated with impervious surfaces and limited infiltration capacities. The CNs utilized in this study were generated using a 250 m hydrological soil group dataset (HYSOG250m) and the 2015 ESA–CCI 300 m land cover dataset, and available at a 250 m grid resolution (
https://doi.org/10.6084/m9.figshare.7756202.v1, accessed on 12 February 2023).
With the CN data, we generated a monthly dataset of R for the period of study using Grove et al. (1998)’s [
47] equations:
Given that
where R = runoff depth, P = precipitation, S = potential maximum retention, and CN is the curve number [
47].
2.2.7. Soil Saturated Hydraulic Conductivity (Ks and PKs)
Soil saturated hydraulic conductivity describes the ability of soils to transmit water under saturated conditions [
48]. Using remotely sensed environmental variables and the RF ML technique, Gupta et al. (2021) [
49] generated four global K
s maps representing four different soil depths at a 1 km (~0.01° × 0.01°) resolution. All four maps were utilized in this study (
https://doi.org/10.5281/zenodo.3935359, accessed on 11 October 2023).
To obtain a continuous and dynamic K
s dataset, we multiplied each K
s value by the monthly precipitation estimates. We selected precipitation as it is a primary driver of water input into the soil and significantly influences K
s. By integrating the precipitation data with K
s, we added temporal variability to the K
s variable, making it suitable for our ML models.
Here, P is the precipitation dataset, Ks is the soil saturated hydraulic conductivity dataset at a depth of i, and i ranges from 1 to 4. The four resulting outputs from Equation (3) were averaged into a single comprehensive PKs dataset for the study area.
2.2.8. Groundwater Storage Percentile (GWSP)
Using terrestrial water storage (TWS) observations from the Gravity Recovery and Climate Experiment (GRACE) satellite mission and a numerical model representing the interactions between water and energy across the Earth’s surface, scientists at NASA are able to determine weekly groundwater conditions, expressed as percentiles, based on comparison with historical data [
50]. These percentiles indicate the probability of occurrence within the 1948 to 2014 period of record and are generated at a spatial resolution of 0.125° × 0.125° over North America from April 2002 to the present [
50]. We downloaded the monthly averages of the GWSP from the Giovanni data portal (
https://giovanni.gsfc.nasa.gov/giovanni/, accessed on 10 October 2023) for our period of study. Considering all the input datasets, the GWSP had the largest spatial resolution (0.125° × 0.125°), requiring all the other datasets to be resampled (upscaled) to the spatial resolution of the GWSP rasters.
Groundwater recharge is typically computed using a simple water balance approach assuming negligible changes in soil water storage in the unsaturated zone [
51]. Based on this idea, we introduced a secondary input variable, called the recharge index (RI), to represent the balance between the water inflows and outflows within each 0.125° × 0.125° grid and the amount potentially available to recharge groundwater.
where P is precipitation; ET is evapotranspiration; R is runoff depth; and RI is the recharge index for a given grid.
GWLs from monitoring wells located within a grid were assumed to be representative of the entire grid, as were the predictor variables. Ultimately, this study utilized six predictor variables, namely soil moisture (SM), land surface temperature (LST), the vegetation index (VI), saturated hydraulic conductivity (K
s), the groundwater storage percentile (GWSP), and the recharge index (RI) (
Table 1).
2.3. Model Algorithms
Developing an acceptable ML model with a monthly temporal resolution requires approximately 10–12 years of data [
52]. Here, we utilize 10 years of monthly data (2010–2019), meeting the acceptable sample size threshold. Selecting algorithms or models that are most appropriate for the data is another important aspect of ML. Getting this step wrong could result in unreliable predictions, leading to a disappointing predictive performance and misleading conclusions [
53]. In this study, we selected two ML algorithms (and a GI technique) from a pool of candidates. The selected techniques are discussed in detail below:
2.3.1. The SVM and SVR
SVM was introduced by Vladimir Vapnik based on the idea of nonlinear mapping of input vectors to a high-dimensional feature space and constructing an optimal hyperplane to effectively separate the different groups or classes within that space [
54]. The generalization capabilities of SVMs led to the development of the less popular SVR for real-value (regression) problems [
55].
First proposed in 1996 by Harris Drucker and his colleagues [
56], the SVR has become an effective tool for prediction problems, demonstrating excellent generalization abilities and a high prediction accuracy [
55]. It works by incorporating a loss function (known as the epsilon insensitive margin of error, ϵ) in the form of a flexible tube formed symmetrically above and below the estimated function, where the prediction errors (ζ*) within the tube are accepted and those that fall outside are penalized (
Figure 3). The objective of the SVR is to find the narrowest possible tube around the estimated function in a way that minimizes the prediction errors [
55,
57].
In general, the performance of the SVR model depends on the tube size (epsilon, ϵ), the regularization constant (C), and the choice of kernel function [
58,
59,
60]. The C hyperparameter controls the complexity of the model, where large values may lead to overfitting [
58,
61]. By overfitting, the model learns the training data well but generalizes poorly, i.e., predicts poorly on new unseen data. Kernel functions are used to transform the data into the higher-dimensional feature space, enabling linear machine learning to improve the representation of the nonlinear relationships that exist in the original input space [
62]. While there is no guide to the appropriate kernel functions for specific datasets, the most commonly used are the RBF and polynomial functions [
63]. RBFs are versatile kernels used when there is a lack of prior knowledge about the data [
64]. Such models (RBF-SVR) require an additional hyperparameter, gamma (γ), in addition to C and ϵ, which controls the width of the RBF [
60,
65,
66]. The C hyperparameter must be a positive number, while ϵ and γ can be positive or zero [
66].
2.3.2. RF
The RF ML algorithm is an ensemble algorithm of multiple trees that improves the prediction accuracies of the single DT algorithm [
19,
67]. The different DTs are trained with subsets of the input variables and bootstrapped samples of the original training data such that each DT is unique, resulting in reduced variance [
11]. By bootstrapping, samples are randomly drawn (with replacement) from the original training data, maintaining the sample size of the training data. Because the sampling is carried out with replacement, a particular observation may appear multiple times in a bootstrapped sample.
Decision points in the DT structure are called nodes. At the nodes, tree branches are created based on the splitting criteria (
Figure 4). The first node (without prior branching) is the root [
68]. From the root, each node is split using the best variable among the subset of input variables chosen at that node [
69]. The leaf is the final node (with no further branching) associated with an output value [
68].
RF model hyperparameters are the number of trees (DTs) and the number of input variables in the random subset at each node [
11,
69]. The final predictions are either determined by majority votes from individual DTs (classification) or by averaging the predictions from all the trees (regression) [
69,
70].
2.3.3. EBK
EBK differs from CK in the way it optimizes the parameter uncertainty associated with creating a single semivariogram and the way it automates the optimization process. In a single semivariogram, the semivariance (the y-axis) measures the spatial dependency between pairs of observations or samples, and the lag (the x-axis) is their separation distance (
Figure 5a). Depending on the characteristics of the data, a semivariogram may display three important components: a sill, a range, and a nugget (
Figure 5a). The range is the distance (or lag) beyond which samples are not spatially autocorrelated, and the sill is the semivariance at that distance. The nugget is the y-intercept of the semivariogram and represents variability at distances much smaller than the minimum spacing between pairs of sample points.
With EBK, the input data are divided into subsets, specifying parameters such as the size of the subsets (subset size) and the degree of overlap between them (overlap factor). Within each subset, a semivariogram distribution is produced (
Figure 5b), and predictions are made for each location using the distribution from one or more subsets [
23].
2.4. Model Design
The initial phase of the analysis aimed to assess the feasibility of ML to capture spatiotemporal patterns of GWLAs across the study area. For spatial patterns, we trained 12 different SVR models (each corresponding to a specific month) on sixty percent of the dataset using manually tuned values of γ, C, and ϵ and tested each of the trained models with the remaining forty percent (
Figure 6a). Each model was trained using three predictor variables—LST, RI, and the previous month’s GWLA (except for the first month). Incorporating GWL from a prior time step is common practice. It is the most employed type of input data in GWL prediction [
52]. To evaluate ML feasibility for temporal patterns, we utilized the same dataset as the test set but used all six predictor variables and the predicted output—replacing the “previous month GWLA” variable to ensure consistency in the number of predictors for each observation (
Figure 6b). Following another random train/test split procedure, an RF model was trained on eighty-five percent of the data and tested on the remaining fifteen percent. The models were evaluated for their performance and predictive accuracy using the Nash–Sutcliffe efficiency (NSE) and the coefficient of determination (R
2), respectively. The NSE ranges from –∞ to 1, while R
2 ranges from 0 to 1. For both metrics, a perfect prediction would yield a value of 1.
Based on the results of the preliminary assessments, it was evident that incorporating some approximation of the GWL as a predictor variable would significantly enhance the model performance and predictive accuracy. Therefore, we sought to use EBK for those approximations. To ensure that the EBK predictions covered the study area, we obtained additional GWL data from 18 nearby monitoring wells outside the study area. A total of 120 EBK models were constructed using the Geostatistical Wizard interpolation tool in ArcGIS PRO, and monthly predictions from those models were converted into anomalies and incorporated into the dataset. A single RF ML model (
Figure 7) was then developed to learn the patterns in substantial portions of the augmented dataset.
This approach was first validated at three monitoring well locations. The measurements corresponding to those locations were removed from the dataset before performing the EBK. Because the accuracy of kriging is significantly influenced by the number and density of kriging points, we could not afford to eliminate additional wells (and their corresponding data) from the dataset. Subsequently, the EBK process was repeated using the entire dataset, and a final RF model was trained on a significant portion of this dataset, tested on a smaller subset, and eventually deployed at locations without monitoring wells.
4. Discussion
Although ML is able to understand complex relationships between GWLs and contributing factors, this study revealed much better predictive performance for the unconsolidated material aquifers, where the relationships are generally more straightforward. Validation wells 1 and 3 were in close proximity to other wells (
Figure 10) and could have benefited from the spatial autocorrelation of GWLs between them, but well 3 did not. In fact, validation well 2 performed better than well 3, despite not having that advantage (
Figure 10,
Table 4). This discrepancy may have been due to several factors, including the intricate heterogeneities and geologic structures in rock aquifers. The PDP for the PK
s variable (
Figure 11) also showed a distinction between the two aquifer types. This suggests that even with increased precipitation, changes in aquifer properties (such as reduced permeability in rock aquifers) can restrict groundwater flow and contribute to a declining GWL, considering all other factors at play. Both aquifers showed a negative relationship with PK
s, but GWLs were relatively higher for the unconsolidated material aquifers (
Figure 11). This shows that under similar hydrological conditions and assuming all other factors are kept constant, the properties of unconsolidated aquifer materials may allow them to maintain higher GWLs compared to rock aquifers.
The discrepancy in model performance could also have been due to both model and data limitations. RF, the most employed ML algorithm for GWL prediction [
10], failed to effectively capture GWL trends in a dolomite rock aquifer in a semi-arid region [
73] and was also outperformed in an aquifer with fractured hydrogeology in a similar climate [
74], as shown by two separate studies. In both cases, deep learning (DL) models demonstrated the best performance. DL is a branch of ML that is based on the concept of deep neural networks and is especially known to outperform traditional (shallow) ML techniques in applications involving large amounts of data [
75,
76] and high dimensionality [
75]. But drilling a large number of wells into hard rock formations in dry regions might not be practical for various reasons, including cost and limited water availability [
77], which can pose significant challenges for ML-based GWL predictions in rock aquifers in arid regions. Specifically, fewer wells were drilled into the consolidated rocks in the study area compared to the basin fill unconsolidated materials [
29], which potentially limited the adequate representation of the geologic complexities of these rock aquifers. However, the relatively weak performance of the RF models for the rock aquifers presents opportunities for future research to investigate and enhance the model performance in such complex geologic and climate settings. These efforts should begin with acquiring the maximum amount of good-quality data for a comprehensive analysis.
The prediction errors from the EBK models (
Figure 12b) underscore the importance of the spatial density and distribution of the kriging data in ensuring the reliability of predictions. The largest uncertainties were seen around the boundaries of the study area, and the predictions in the vicinity of the monitoring wells were relatively more accurate. Two of the wells excluded in the validation phase (
Figure 10) were included in the final model as test wells (
Figure 13) and showed improved predictions (
Table 4,
Table 5), further underscoring the importance of data density and quality in kriging.
Based on the percent increases in the mean squared error (MSE) when important predictor variables are left out, the RI and EBK variables were the most important in the validation RF model, where both variables showed similar levels of importance. However, in the final deployment model, EBK was the most important variable (by a significant margin). This suggests that incorporating spatial interpolation techniques such as EBK can substantially enhance the performance of ML models. Although the kriging process can be tedious and challenging, the model improvements they offer make these efforts worth it.
As shown, the average GWLA for the period of study was predominantly negative. In fact, only about twenty-eight percent of the study area showed a positive average anomaly during this period. This trend reflects the challenges in a dry/arid region with high groundwater demand and withdrawal rates possibly exceeding natural recharge. Historically, groundwater in Arizona has been pumped out faster than it has been replenished by natural means [
78,
79], resulting in overdraft in many agricultural and urban areas [
79]. But as the quest to exploit deeper aquifers continues, the costs of drilling to these depths are much higher, as are the energy costs of pumping water from them [
79]. This study therefore could be useful for optimizing the drilling process by identifying locations for new wells and increasing the likelihood of accessing groundwater at optimal depths. Efforts to manage the groundwater overdraft issue began with identifying regions with a high reliance on groundwater (known as Active Management Areas (AMAs)) and subsequently empowering the ADWR to monitor compliance within the AMAs with the regulatory frameworks in place [
79,
80]. Within these AMAs and beyond, this study can also aid in the monitoring and allocation of groundwater resources by identifying groundwater-deficient areas based on the average predicted GWLA values, and offer data-driven support and recommendations towards the effective management of groundwater resources in vulnerable areas.
5. Conclusions
Groundwater is the largest reservoir of available freshwater in the world and a critically important resource. Its global relevance is amplified by the direct impacts of climate change on surface water sources, particularly in arid regions. In this study, we demonstrated the effectiveness of ML in predicting monthly GWLAs when combined with reliable spatial interpolation models and developed the first statewide GWLA prediction model for the state of Arizona. Following satisfactory performance based on average NSE/R2 values of 0.62/0.63 and 0.72/0.76 during the validation and testing phases, respectively, monthly GWLA rasters were produced for January 2010 to December 2019. Moving forward, future studies may focus on addressing some of the challenges of applying traditional ML techniques to rock aquifers in dry regions discussed in this study, in terms of leveraging the available data and reducing prediction uncertainties in such complex settings.
With well depths ranging from 25 to 1600 feet, this study demonstrated effectiveness for both shallow and deep aquifers. The model design utilized remotely sensed datasets from satellites with global coverage, enabling replicability for similar climates across the globe. Our remote sensing approach ensures that data-sparse regions of the world, where field-based hydrological variables are limited or largely inaccessible, are not left out. It is our hope that this work contributes substantially to the science of monitoring groundwater resources in the face of global warming and climate change threats, ensuring the availability of groundwater to meet domestic, agricultural, and industrial water needs.