Spatial Mapping and Prediction of Groundwater Quality Using Ensemble Learning Models and SHapley Additive exPlanations with Spatial Uncertainty Analysis

Abstract

The spatial mapping and prediction of groundwater quality (GWQ) is important for sustainable groundwater management, but several research gaps remain unexplored, including the inaccuracy of spatial interpolation, limited consideration of the geological environment and human activity effects, limitation to specific pollutants, and unsystematic indicator selection. This study utilized the entropy-weighted water quality index (EWQI), the LightGBM model, the pressure-state-response (PSR) framework and SHapley Additive exPlanations (SHAP) analysis to address the above research gaps. The normalized importance (NI) shows that NO₃⁻ (0.208), Mg²⁺ (0.143), SO₄²⁻ (0.110), Cr⁶⁺ (0.109) and Na⁺ (0.095) should be prioritized as parameters for remediation, and the skewness EWQI distribution indicates that although most sampled locations have acceptable GWQ, a few areas suffer from severely poor GWQ. The PSR framework identifies 13 indicators from geological environments and human activities for the SMP of GWQ. Despite high AUROCs (0.9074, 0.8981, 0.8885, 0.9043) across four random training and testing sets, it was surprising that significant spatial uncertainty was observed, with Pearson correlation coefficients (PCCs) from 0.5365 to 0.8066. We addressed this issue by using the spatial-grid average probabilities of four maps. Additionally, population and nighttime light are key indicators, while net recharge, land use and cover (LULC), and the degree of urbanization have the lowest importance. SHAP analysis highlights both positive and negative impacts of human activities on GWQ, identifying point-source pollution as the main cause of the poor GWQ in the study area. Due to the limited research on this field, future studies should focus on six key aspects: multi-method GWQ assessment, quantitative relationships between indicators and GWQ, comparisons of various spatial mapping and prediction models, the application of the PSR framework for indicator selection, the development of methods to reduce spatial uncertainty, and the use of explainable machine learning techniques in groundwater management.

1. Introduction

Groundwater is an indispensable freshwater resource across many regions, supporting essential services such as water supply, agricultural irrigation, and industrial development [1]. However, in recent decades, groundwater quality (GWQ) has been significantly compromised due to factors such as urbanization, climate change, overexploitation and inadequate resource management [2,3,4,5,6]. In many areas, deteriorated groundwater is still used for drinking purposes due to the absence of alternative sources, which leads to significant health risks to local populations [7,8,9]. Moreover, as a crucial component of the global hydrological cycle, the degradation of GWQ directly impacts ecological stability and regional sustainability [10,11]. Under these circumstances, the precise spatial mapping and prediction of GWQ is essential for identifying pollution sources and informing comprehensive water management strategies for sustainable groundwater management.

Currently, numerous index-based methods are available for assessing GWQ, including the water quality index (WQI) and its adaptations, such as the entropy-weighted WQI (EWQI) [12], principal component analysis (PCA)—WQI [13], the integrated WQI (IWQI) [14], the drinking WQI (DWQI) [15], and the CCME-WQI [16]. Additionally, other methodologies, like the Nemerow index [17] and the comprehensive pollution index (CPI) [18], are also widely employed in GWQ assessment. However, while these methods enable us to indicate the overall GWQ based on samples, they fail to map the spatial distribution of GWQ across larger areas [19,20,21]. This research gap restricts decision-makers in managing groundwater spatially, thereby obstructing the implementation of effective groundwater management strategies across different regions.

Currently, for the spatial prediction of GWQ, spatial interpolation is primarily used (e.g., kriging and inverse distance weighting), which estimates values at unsampled locations by assuming that closer points are more similar than distant ones based on spatial autocorrelation [22,23,24,25,26,27]. However, while spatial interpolation offers a convenient means to estimate the spatial distribution of GWQ, its limitations are significant. Firstly, the accuracy of predictions heavily depends on the density and distribution of sampling points, with sparse data leading to unreliable results [28,29]. Secondly, these methods often assume spatial homogeneity, failing to account for complex environmental variability and external factors like human activities, both of which can significantly influence GWQ [30,31,32]. Lastly, the accuracy of these methods diminishes near the edges of the study area due to fewer data points, known as edge effects [33,34]. Therefore, there is a critical need for new approaches in GWQ mapping and prediction that address the limitations of these traditional interpolation methods, accommodate spatial heterogeneity, and ensure accuracy in data-sparse areas.

Recent studies have increasingly utilized machine learning (ML) models due to their ability to accurately assess and predict GWQ. For example, Singha et al. [35] conducted a comparative analysis of random forest (RF), extreme gradient boosting (XGBoost), artificial neural networks (ANNs), and deep learning (DL) methods in GWQ assessment, concluding that phosphates have a high relative importance. El Bilali et al. [36] found that while adaptive boosting (AdaBoost) and RF models performed better in forecasting GWQ for irrigation purposes, ANNs and support vector regression (SVR) models exhibited greater generalizability. Jeihouni et al. [37] used five decision-tree-based data mining algorithms to identify high-quality groundwater zones, finding RF to be the most accurate for creating reliable GIS-based GWQ maps. However, while these studies primarily focus on the physical and chemical characteristics of GWQ and their interrelations, they often overlook the critical impacts of geological environment settings (e.g., geology, topography and climate) and human activities (e.g., urbanization and pollution source) on GWQ. This oversight leads to the challenges in accurately assessing the factors contributing to groundwater pollution. Meanwhile, research on ML models for single pollutants like nitrate, arsenic, salinity, and fluoride is more detailed, focusing on their predictive accuracy in GWQ prediction. For example, the ML models, such as RF, ANN, XGBoost, CART, BRT, SVR, and KNN, and Bayesian-based methods are employed in mapping and the spatial prediction of nitrate distribution, incorporating geological environmental parameters [19,38,39,40,41,42]. Podgorski et al. [43] employed RF and multivariate logistic regression (MLR) to screen 25 indicators and predict the distribution of fluoride across India. Xia et al. [44] utilized four models (XGBoost, RF, AdaBoost, and SVM) to perform spatial predictions for fluoride, arsenic, and iodine in the Hetao Basin, China by considering different environmental factors. Tran et al. [45] identified ten environmental indicators to compare the performance of various ML models in predicting groundwater salinity in coastal areas and found that the CatBoost regression model provides the highest accuracy. Podgorski and Berg [46] used RF and eleven geological and climatic indicators to carry out spatial predictions of arsenic contamination globally. However, although these studies account for various geological environmental factors for GWQ spatial prediction, they primarily focus on specific pollutants. The analysis of all GWQ parameters and environmental stress indicators in these studies is insufficient, and the application of ML models for the spatial prediction of overall GWQ indexes, such as the EWQI, the IWQI, and the CPI, remains underexplored.

The pressure-state-response (PSR) framework, developed by the Organization for Economic Cooperation and Development (OECD), provides a structured framework to analyze the dynamic interactions between human activities, environmental conditions, and management responses [47,48,49]. This framework has been ingrained in regulatory approvals and development management plans across Australia, especially for mining applications and environmental assessments. Also, many studies, such as those focused on ecological security [50], ecological vulnerability [51], forest management [52], mine area contamination assessment [53], and urban resilience [54], also demonstrated its wide application. However, its application in GWQ prediction for indicator selection is still limited. In fact, in GWQ prediction, there is a lack of systematic approach in the selection of indicators, especially concerning those related to human activities. The PSR framework aptly provides a structured methodology to identify and analyze various indicators crucial for assessing GWQ. This framework is especially suitable in rapidly urbanizing areas, where the dual impact of human activities on GWQ (both detrimental and beneficial) presents a complex challenge that demands further exploration.

Building upon identified research gaps, this study aims to conduct spatial GWQ mapping and prediction using the Guanzhong Basin as the study area. This research integrates the PSR framework, the EWQI, the LightGBM model, and explainable machine learning techniques (EMLTs). In Guanzhong Basin, studies on GWQ primarily include contamination risk [55], human health risk [56], hydrogeochemical processes [57], and water quality assessment [58]. However, the spatial prediction of GWQ has not yet been explored. Therefore, this study include three innovations: (1) It pioneers the use of the advanced LightGBM model and the EWQI to perform spatial mapping and carry out predictions of GWQ, as this has not been explored; (2) It utilizes the PSR framework to systematically to select indicators for GWQ mapping and prediction, considering geological environment indicators, spatial uncertainty and dual impacts of human activities on GWQ; (3) It incorporates SHapley Additive exPlanations (SHAP), a widely used advanced explainable machine learning technique (EMLT), to visualize the influence of these indicators on GWQ distribution, thereby supporting decision-making in sustainable groundwater management.

2. Study Area

The Guanzhong Basin, located in the central part of Shaanxi Province, China, serves as an essential agricultural and industrial area; it is bordered by the Qinling Mountains to the south and the Bei Mountains to the north [59]. This basin covers an area of approximately 18,955.25 km² and is characterized by its longitudinal range of 107°–110°30′ E and latitudinal span of 34°00′–35°40′ N (Figure 1). Geologically, the basin is distinguished by a thick layer of Mesozoic sedimentary rocks and is underlain by complex hydrogeological structures formed from tertiary river-lake facies, heavily influenced by historical tectonic activities [60]. The region’s climate is classified as temperate, with four distinct seasons and an average annual temperature of 13.3 °C, and annual rainfall varies from 544 to 863 mm, predominantly occurring during the summer months [55]. However, the area is prone to droughts due to its high annual evaporation rate of 800–1200 mm [61]. Hydrologically, the Guanzhong Basin is dominated by the Weihe River, the largest tributary of the Yellow River, which plays a crucial role in the regional water system by linking surface water interactions with the groundwater [62]. Despite its natural water riches, the basin faces challenges related to water scarcity and the uneven seasonal distribution of rainfall, which can impact both agricultural productivity and urban water supply. Groundwater in the basin is found mainly in unconfined aquifers with thicknesses varying from 5 to 80 m, predominantly recharged by precipitation and lateral flows from adjacent mountainous regions [55]. The infiltration coefficients of the floodplains and terraces further highlight the complex interaction between surface and groundwater systems. The significant human activities, including dense urbanization and industrial operations in cities like Xi’an and Xianyang, intensify the demand for water and place additional pressures on the groundwater systems. Given these factors, this study aims to leverage advanced machine learning models to enhance the prediction and management of GWQ within the Guanzhong Basin, focusing on integrating environmental, climatic, and anthropogenic indicators to provide a comprehensive analysis of the region’s groundwater sustainability.

3. Methodology

Figure 2 illustrates the methodological framework of this study, which includes four parts: indicator determination based on the PSR framework; GWQ assessment based on the EWQI; GWQ mapping and spatial prediction based on LightGBM and the TPE; and indicator analysis based on SHAP values. It is crucial to note that GWQ assessment relies on solely groundwater sample analysis using the EWQI to obtain overall GWQ. Conversely, spatial GWQ mapping and prediction extend these EWQI calculations spatially through indicators identified by the PSR framework. Finally, future directions are proposed for sustainable groundwater management.

3.1. Groundwater Samples Descriptions

The groundwater samples used in this study were derived from the research conducted by Chengzhu et al. [63], which originally comprised 200 groundwater samples. After excluding 10 blank samples and 10 duplicates, a total of 180 samples remained for analysis. The samples were evaluated for various physicochemical parameters and selected for this study if they exceeded the Type III groundwater standards of China (GB/T 14848-2017) [64]. For parameters not covered by the GB/T 14848-2017 standards, the World Health Organization (WHO) drinking water guidelines and a related study [35] were used. A total of 16 indicators were selected, and a description of these data is provided in

Table 1. GWQ samples descriptions (mg/L).

Parameters	Min	Max	Mean	SD	Standard
pH	6.96	9.89	7.84	0.32	6.5–8.5
Total Hardness (TH)	9	1885	478.87	291.5	450
Total Dissolved Solids (TDS)	196	10570	1077.61	1098.60	1000
Calcium (Ca2+)	0.56	301	92.58	55.91	75
Magnesium (Mg2+)	1.8	352	60.24	53.36	30
Potassium (K+)	0.13	49.5	3.17	5.82	12
Sodium (Na+)	6.36	1160	140.62	155.88	200
Chloride (Cl−)	3.7	2135	106.54	195.99	250
Sulfate (SO42−)	1.33	4255	230.9	434.04	250
Bicarbonate (HCO3−)	117	1349	509.41	189.46	300
Nitrate (NO3−)	0	373	47.32	53.59	20
Fluoride (F−)	0.12	4.26	0.98	0.77	1
Zinc (Zn2+)	0.001	0.066	0.008	0.008	1
Hexavalent chromium (Cr6+)	0.001	0.45	0.033	0.055	0.05
Aluminum (Al3+)	0.003	0.1	0.008	0.01	0.2
Iron (Fe3+)	0	0.35	0.105	0.074	0.3

.

3.2. GWQ Assessment

3.2.1. EWQI Calculation

The EWQI method employs an entropy-based objective weighting system for each parameter, making it extensively applicable in the assessment of GWQ [12,65,66,67,68]. Compared to the traditional WQI, the EWQI provides a more objective evaluation by reducing the subjectivity in parameter weighting, thereby offering a more reliable assessment of GWQ across diverse environments. Given $m$ groundwater samples and $n$ parameters, an $m\times n$ matrix ($X$) can be constructed, as shown as Equation (1).

$$X=\left(\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1n}\\ x_{21}& x_{22}& \cdots & x_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ x_{m1}& x_{m2}& \cdots & x_{mn}\end{array}\right)$$

(1)

Due to variations in units and magnitudes across different parameters, the normalization step is essential, resulting in a standardized matrix $R$.

$$R_{ij}={\displaystyle \frac{x_{ij}-min(x_j)}{max(x_j)-\mathrm{m}\mathrm{i}\mathrm{n}(x_j)}}$$

(2)

where $R_{ij}$ denotes the normalized value of the $i^{th}$ groundwater sample for $j^{th}$ parameter, and $min(x_j)$ and $max(x_j)$ are the minimum and maximum values of the $j^{th}$ parameter, respectively.

The entropy value ($H_j$) is then calculated to determine its relative importance by Equation (3), as lower entropy indicates greater parameter significance.

$$H_j=-{\displaystyle \frac{1}{\mathrm{l}\mathrm{n}\left(m\right)}}{\sum }_{i=1}^m{P}_{ij}ln(P_{ij})$$

(3)

where $P_{ij}={\displaystyle \frac{R_{ij}}{{\sum }_{i=1}^m{R}_{ij}}}$ and $P_{ij}ln(P_{ij})=0$ if $P_{ij}=0$.

The weight $W_j$ for each parameter is subsequently derived from the entropy values using Equation (4).

$$W_j={\displaystyle \frac{1-H_j}{n-{\sum }_{j=1}^m{H}_j}}$$

(4)

To evaluate the quality index ($Q_{ij}$) of groundwater sample $i$ for parameter $j$, Equation (5) is employed.

$$Q_{ij}=\left\{\begin{array}{c}100\times {\displaystyle \frac{x_{ij}}{S_j}}\\ {\displaystyle \frac{C_{pH}-7}{S_{pH}-7}}\end{array}\right.$$

(5)

where $S_j$ refers to the reference value for parameter $j$ (See Table 1); $C_{pH}$ is the value of pH; and $S_{pH}$ is the permissible limit of pH, which can be either 6.5 or 8.5. When the pH is less than 7, the limit is set at 6.5; when it is greater than 7, the limit is set at 8.5.

Finally, the EWQI is computed the following equation:

$$EWQI={\sum }_{j=1}^n{W}_j\times Q_{ij}$$

(6)

3.2.2. Parameter Analysis

Similar to other WQIs, the EWQI calculates the overall GWQ of a sample point but cannot assess the conditions of individual parameters within the area. Moreover, existing methods for single-parameter water quality analysis do not fully consider the importance of the parameters and their exceedance rates. Therefore, this study defines a new term, normalized importance (NI), considering both importance and exceedance rates to reflect the priority level for management of a single GWQ parameter in the study area. Parameters with high NI values should be prioritized for management to reduce the impact of groundwater pollution in the study area. The formula for NI is

$$NI={\sum }_{j=1}^n{W}_j\times E_j$$

(7)

where $W_j$ is the weight value for parameter $j$, and $E_j$ is the exceedance rate of groundwater samples for parameter $j$ compared to the reference value.

3.3. GWQ Mapping and Prediction

3.3.1. Data Split

First, 180 groundwater samples were categorized based on their calculated EWQI values into two groups: 0 (90 samples), representing good GWQ, and 1 (90 samples), representing poor GWQ. To determine the optimal split between training and validation sets, we conducted a preliminary analysis by comparing the average model performance and standard deviation (SD) of six split ratios (65/35, 70/30, 75/25, 80/20, 85/15 and 90/10) based on four random selections. This preliminary analysis has been validated as an effective method for selecting training and validation sets [69,70,71,72].

3.3.2. Indicator Selection by PSR Framework

Based on the PSR framework, 13 indicators were selected, including potential pressure indicators, state indicators, and potential response indicators. The reason these are considered “potential” is that it is not yet clear whether the human responses they represent are beneficial or detrimental to GWQ, necessitating further investigation.

Table 2. Thirteen indicators determined by PSR framework.

Group	Indicators	Sources	Scale	Format
Pressure	Population	SEDAC	250 m	Raster
	PPSD	SPDEE	1:300,000	Point
	LULC	Yang and Huang [73]	30 m	Raster
State	Depth to groundwater	MWRPIC	1:300,000	Line
	Net recharge	Peng et al. [74]	1 km	Raster
	Aquifer water yield capacity	Hydrogeological map	1:300,000	Polygon
	Slope	NASADEM data	30 m	Raster
	Impact of the vadose zone	Zhang et al. [55]	1:300,000	Polygon
	Conductivity	Zhang et al. [55]	1:300,000	Polygon
Potential response	GDP2015	GRDC	1 km	Raster
	Ten years change of NDVI	Yang et al. [75]	30 m	Raster
	Degree of urbanization	SEDAC	250 m	Raster
	Nighttime lights	Elvidge et al. [76]	1 km	Raster

Note(s): SEDAC (Socioeconomic Data and Application Center); SPDEE (Shaanxi Provincial Department of ecology and environment); MWRPIC (Ministry of Water Resource of the People’s Republic of China); GRDC (Geographic Data Sharing Infrastructure, global resources data cloud).

displays information on the data sources, scales, and relevant details of these indicators. To ensure consistency across all indicators during modeling, we standardized the spatial resolution. For vector data such as the point, line and polygon, we utilized the vector-to-raster conversion tool in GIS to achieve a resolution of 30 m (Note: point density was used for PPSD). For raster data of varying resolutions, we applied resampling techniques to uniformize the resolution to 30 m. This standardization ensures uniformity and reliability in our data analysis.

Potential Pressure Indicators

Potential pressure indicators are the factors that many influence GWQ through external forces, such as agricultural activities, industrial emissions, and urban expansion. In this study, population, land use and cover (LULC) and potential pollution-source density (PPSD) are selected as potential pressure indicators for GWQ mapping and prediction (Figure 3). The population size may drive groundwater demand and contribute to wastewater and solid waste production, which can contaminate groundwater [77,78]. Higher population densities typically increase the risk of over-extraction and pollution. In terms of LULC, urbanization increases impervious surfaces, decreasing groundwater recharge and increasing runoff that may carry pollutants [79]. Agricultural practices often use fertilizers and pesticides, risking groundwater contamination through leaching and infiltration [80,81]. PPSD measures the concentration of potential contaminant sources like industrial areas, waste disposal sites, and chemically intensive agriculture [82]. Higher PPSD elevates the contamination risk to groundwater systems.

State Indicators

State indicators, referencing the DRASTIC model in groundwater vulnerability, are utilized to measure the fundamental hydrogeological conditions that determine GWQ. Groundwater vulnerability is used to measure capacity to resist pollutants entering the groundwater system, which reflects the current state and attributes of the area [83]. In this study, six indicators were selected, including depth to groundwater, net recharge, aquifer water yield capacity, slope, impact of the vadose zone, and conductivity (Figure 4). We removed the indicator of soil media because the net recharge is calculated by multiplying rainfall by an infiltration coefficient, which is determined based on soil type. The influence of these indicators on GWQ and the scores for each indicator can be found in many groundwater vulnerability studies [83,84,85,86,87,88,89]. A higher score indicates that the aquifer system in the area is more vulnerable to contamination, while a lower score suggests greater resistance to pollution.

Potential Response Indicators

Response indicators represent the measures and policies implemented to address or mitigate impacts on GWQ. Due to the dual impact of human activities on GWQ, we selected several potential indicators for the response category. These include GDP2015, ten-year changes in the NDVI, degree of urbanization, and nighttime lights (Figure 5). GDP reflects economic activity levels, where higher values are not only linked to greater environmental impacts from industrial and agricultural runoff affecting GWQ but also indicate increased potential for funding and implementing policies aimed at mitigating these impacts [51,90]. Ten-year changes in the NDVI, the degree of urbanization settlement, and nighttime lights collectively represent the dual impacts of human activities on GWQ. The changes in the NDVI indicate variations in vegetation cover that can either enhance groundwater recharge and pollutant filtration with increased greenery or reduce these capabilities through land degradation [91,92,93]. Urbanization increases impervious surfaces and pollution runoff, degrading natural water infiltration and quality, yet it also prompts opportunities for implementing advanced urban planning and sustainable infrastructure to protect groundwater [3,94,95]. Similarly, increased nighttime lights correlate with intensified urban and industrial activities that elevate contamination risks, but they also mark areas where targeted environmental regulations and remediation efforts can effectively mitigate these impacts [96,97]. Each indicator not only reflects the challenges posed by human activities but also underscores the potential for proactive groundwater management responses.

3.3.3. Correlation Analysis

To ensure the accuracy and effectiveness of our model, a correlation analysis was conducted on the dataset variables before modeling. We used the Pearson correlation coefficient (PCC) to identify linear relationships between variables. Any pair of variables with a correlation coefficient exceeding 0.7 was considered strongly correlated, and one variable from each pair was removed to avoid multicollinearity, which can impair model stability and interpretability. The equation of PCC is shown as follows:

$$r={\displaystyle \frac{\sum (x_i-\overline{x})(y_i-\overline{y})}{\sqrt{\sum {(x_i-\overline{x})}^2\sum {(y_i-\overline{y})}^2}}}$$

(8)

where $x_i$ and $y_i$ are the values of the two variables, $\overline{x}$ and $\overline{y}$ are the means of variables, respectively.

3.3.4. LightGBM Model

In our analysis, the LightGBM model, developed by Ke et al. [98], was employed due to its effectiveness in processing large-scale and high-dimensional datasets. This model integrates two key innovations: gradient-based one-side sampling (GOSS) and exclusive feature bundling (EFB) [99]. EFB reduces the dimensionality of the feature space by grouping features that are rarely active at the same time, treating the combination as a graph coloring problem where features are vertices connected by edges when not mutually exclusive. GOSS enhances the training process by focusing on instances with larger gradients by maintaining all instances in the top $a$ percentile (denoted as $G_{high}$) and sampling $a$ fraction $b$ from the lower gradients ($G_{low}$) [100]. The information gain from a feature $j$ at split $d$ is calculated using the formula:

$${\tilde{V}}_j\left(d\right)={\displaystyle \frac{1}{n}}({\displaystyle \frac{{\left({\sum }_{x_i\in A_l}{g}_i+\frac{1-a}{b}{\sum }_{x_i\in B_l}{g}_i\right)}^2}{n_l^j\left(d\right)}}+{\displaystyle \frac{{\left({\sum }_{x_i\in A_r}{g}_i+\frac{1-a}{b}{\sum }_{x_i\in B_r}{g}_i\right)}^2}{n_r^j\left(d\right)}})$$

(9)

where $A_l$ and $A_r$ are subsets of $G_{high}$, and $B_l$ and $B_r$ are subsets of the sampled $G_{low}$, with $n_l^j\left(d\right)$ and $n_r^j\left(d\right)$ representing the number of instances on the left and right sides of the split, respectively.

3.3.5. Hyperparameter Selection and Optimization

Hyperparameter selection and optimization are crucial in ML, as they significantly enhance model performance by fine-tuning settings to align precisely with specific data characteristics and learning objectives [101]. In this study, six hyperparameters and their corresponding parameter spaces are presented in

Table 3. Hyperparameter spaces and optimal hyperparameters.

Hyperparameters	Hyperparameter Spaces
bagging_fraction	hp.uniform(‘bagging_fraction’, 0.5, 0.9)
bagging_freq	hp.choice(“bagging_freq”, range(4, 7))
boosting_type	hp.choice(“boosting_type”, [‘gbdt’, ‘dart’, ‘rf’])
feature_fraction	hp.uniform(‘feature_fraction’, 0.5, 0.9)
learning_rate	hp.uniform(‘learning_rate’, 0.01, 0.5)
num_leaves	hp.choice(“num_leaves”, range(15, 128))

. The meanings of these hyperparameters are shown in Appendix A.

Also, we used the tree-structured parzen estimator (TPE) for hyperparameter optimization—a method within Bayesian optimization used for its efficiency in refining model parameters by leveraging historical performance data. Specifically, we utilized the “1-area under the receiver operating characteristic curve (AUROC)” as the objective function for minimization, and the iterations were set as 1000. The detailed information about the TPE approach can be found in Xiong et al. [102], Nguyen et al. [103], Rong et al. [104] and Tao et al. [105].

3.3.6. Model Performance Evaluation

Model performance evaluation is crucial in assessing the efficacy and reliability of ML models, providing insights into their predictive accuracy and guiding improvements to ensure robust real-world applications [106,107]. In this study, building on previous research [108,109], we utilized precision (Equation (10)), recall (Equation (11)), F₁ score (Equation (12)), overall accuracy (OA) (Equation (13)), and AUROC (Equation (14)) as metrics to comprehensively evaluate the performance of the model.

$$Precision={\displaystyle \frac{TP}{TP+FP}}$$

(10)

$$Recall={\displaystyle \frac{TP}{TP+FN}}$$

(11)

$$F_1 score=2\times {\displaystyle \frac{ 2TP}{2TP+FP+FN}}$$

(12)

$$OA={\displaystyle \frac{TP+TN}{TP+TN+FP+FN}}$$

(13)

$$AUROC={\sum }_{i=1}^{n-1}({FPR}_{i+1}-{FPR}_i)\times ({\displaystyle \frac{{TPR}_{i+1}+{TPR}_i}{2}})$$

(14)

In this analysis, $TP$ denotes instances where the model correctly identified poor GWQ samples, while $FP$ indicates cases where good GWQ samples were mistakenly classified as poor. $TN$ represents instances where good GWQ samples were correctly recognized, and $FN$ refers to cases where the model failed to identify poor GWQ samples. Additionally, ${FPR}_i$ measures the proportion of good GWQ samples incorrectly classified as poor at the $i^{th}$ threshold. ${TPR}_i$ quantifies the percentage of actual poor GWQ instances that were correctly identified at the same threshold.

3.4. Spatial Uncertainty Analysis

Some studies have indicated that different sampling methods can affect the uncertainty of the results [110], and similar model performance does not necessarily indicate similar spatial distributions [111]. Therefore, we conducted a spatial uncertainty analysis for four random selections. The spatial uncertainty of four GWQ maps was assessed by two methods. The first approach involves calculating the PCC for all grid cells between each pair of the four GWQ maps to obtain an overall correlation. Stronger correlations (closer to 1) in PCC analysis indicate lower spatial uncertainty

The second method is to spatially visualize the uncertainty between each pair of the four GWQ maps across the study area. For a specific grid cell at position $(i,j)$ in the study area, the spatial uncertainty can be calculated using Equation (15).

$$Uncertainty={\displaystyle \frac{2}{n(n-1)}}{\sum }_{k=1}^{n-1}{\sum }_{l=k+1}^n\left|x\left(k,j\right)-x(l,j)\right|$$

(15)

where $n$ is the number of GWQ maps, and $x\left(k,j\right)$ represents the probability of GWQ from the $k^{th}$ map at grid cell $\left(i,j\right).$ In this study, $n$ equals 4.

3.5. Indicator Importance Analysis and SHAP Analysis

SHAP analysis offers a systematic approach within an EMLT to quantitatively detail the contribution of each feature to a model’s predictions; it has been widely used in the groundwater field [112,113,114]. This approach is crucial for understanding the role of input features (indicators) in determining model outcomes [115]. SHAP analysis is frequently utilized with ensemble models like XGBoost [116], LightGBM [117], CatBoost [118], and RF [119] due to its code compatibility. In SHAP analysis, two key values are included: the SHAP value and the feature value. The former quantifies the impact of each feature (indicator) on the model’s prediction. Its positive or negative sign indicates contributions to the binary outcomes of 1 and 0, respectively. The latter refers to the actual value of the indicator itself, which is used as an input in the model. The SHAP method assigns a value to each feature based on its influence, calculated using Equation (16):

$${\varnothing }_j={\sum }_{S\subseteq N/j}{\displaystyle \frac{\left|S\right|!\left(n-\left|S\right|-1\right)!}{n!}}\left[f\left(S\cup j\right)-f\left(S\right)\right]$$

(16)

where ${\varnothing }_j$ represents the SHAP value for feature $j$, derived by summing contributions over all possible subsets of features excluding $j$; $\left|S\right|$ is the number of features in subset $S$; $n$ is the total number of features; $f\left(S\right)$ denotes the model’s output using subset $S$ without feature $j$; and $f\left(S\cup j\right)$ is the output when feature $j$ is included.

4. Results

4.1. GWQ Assessment Results

Table 4. Weights, exceedance percentage and NI of parameters.

Parameters	Weights	Exceedance Rate	NI
pH	0.008052	38.89%	0.015521
Total Hardness (TH)	0.021665	34.44%	0.036983
TDS	0.062166	1.67%	0.005146
Calcium (Ca2+)	0.024363	60.00%	0.072454
Magnesium (Mg2+)	0.041843	68.89%	0.142876
Potassium (K+)	0.104982	4.44%	0.023104
Sodium (Na+)	0.068902	27.78%	0.094873
Chloride (Cl−)	0.105335	7.22%	0.037696
Sulfate (SO42−)	0.107567	20.56%	0.109618
Bicarbonate (HCO3−)	0.014792	90.00%	0.065986
Nitrate (NO3−)	0.067453	62.22%	0.208023
Fluoride (F−)	0.042884	35.56%	0.075585
Zinc (Zn2+)	0.064638	0.00%	0
Hexavalent chromium (Cr6+)	0.124156	17.78%	0.109416
Aluminum (Al3+)	0.108349	0.00%	0
Iron (Fe3+)	0.032851	1.67%	0.002719

presents the weights, exceedance rate and NI of different parameters. The results revealed that the NI of NO₃⁻ is the highest (0.208), followed by Mg²⁺ (0.143), SO₄²⁻ (0.110), Cr⁶⁺ (0.109) and Na⁺ (0.095). Additionally, while the NI of HCO³⁻, Ca²⁺, and pH is not high, their exceedance percentages are significantly elevated, at 90%, 60%, and 38.89%, respectively. These findings should prompt significant attention from managers in the study area due to the potential impact on GWQ.

Figure 6 shows the boxplot and distribution of the EWQI based on 180 groundwater samples. As shown in the boxplot, the median EWQI value is at approximately 64.68, with the interquartile range extending from 48.72 to 100.26. Some outliers can be observed beyond the whiskers, indicating that the distribution is skewed. The histogram further illustrates the distribution of EWQI values, emphasizing a right-skewed trend (skewness = 3.26 and maximum EWQI = 534.72). A majority of the samples cluster towards the lower range of EWQI values (EWQI < 100), which indicates that most of the 180 groundwater samples from the study area are categorized as “Excellent”, “Good” and “Moderate,” using the standards outlined in many studies [67]. As a result, combined with Figure 6, the overall GWQ in the study area is generally good. However, a few areas suffer from severely poor quality, which warrants urgent attention.

4.2. Indicator Selection by Correlation Analysis

Figure 7 presents the PCCs between indicators. It is evident that all the PCCs are below 0.7, indicating that multicollinearity among these indicators is not significant. Notably, the coefficient between nighttime light and GDP is 0.63, which can be attributed to the fact that areas with higher economic output often have increased nighttime lighting due to urbanization and industrial activity [120,121]. Similarly, the PCC between aquifer media and net recharge is 0.66, and this is because the characteristics of aquifer media directly influence the rate and volume of groundwater recharge through their porosity and permeability properties [122]. Despite the correlations among the indicators, their correlations are not strong enough to undermine the models’ independence for modeling.

4.3. Optimal Hyperparameters, Model Performance, and Spatial GWQ Mapping

Based on the preliminary analysis, it was found that the 80/20 split ratio had high model performance (average AUROC = 0.8989) and the highest robustness (SD = 0.0083) (

Table 5. Data split determination for training and validation dataset by average AUROC and SD.

Selection No.	65/35	70/30	75/25	80/20	85/15	90/10
1	0.8291	0.823	0.8856	0.9074	0.8797	1
2	0.8291	0.8697	0.8573	0.8981	0.8571	1
3	0.7954	0.8011	0.8855	0.8858	0.8546	0.9383
4	0.7808	0.823	0.8601	0.9043	0.8731	0.9877
Average	0.8086	0.8292	0.8721	0.8989	0.8661	0.9815
SD	0.0211	0.0250	0.0135	0.0083	0.0106	0.0254

). Comparing other splits, 80/20 not only maximizes the predictive accuracy but also ensures the model’s stability across different test sets.

Table 6. Optimal hyperparameters for different random sample selections.

Hyperparameters	Selection 1	Selection 2	Selection 3	Selection 4
bagging_fraction	0.755129	0.578798	0.628146	0.817471
bagging_freq	5	4	5	5
boosting_type	gbdt	gbdt	gbdt	gbdt
feature_fraction	0.861003	0.561902	0.899878	0.677141
learning_rate	0.310296	0.069457	0.353949	0.314827
num_leaves	102	45	47	21

illustrates the optimal hyperparameters for four different random sample selections, which highlights that the best optimization results under the TPE approach vary with each sampling method.

Figure 8 displays the testing results of the model for four random splits with an 80/20 ratio. It is evident that the model exhibits excellent performance in all four cases, achieving AUROC values of 0.9074, 0.8981, 0.8858, and 0.9043, respectively. Other metrics, including precision, recall, F₁ score, and OA, further corroborate the results. The consistently high performance not only validates the appropriateness of the selected indicators for GWQ mapping and prediction but also reflects the effectiveness of the LightGBM model combined with the TPE, as stated in Li et al. [123], Guo et al. [124] and Li et al. [125]. However, as emphasized by Xiong et al. [102], similar model performances do not necessarily imply similar spatial distributions. We therefore mapped the spatial distribution of GWQ based on the optimal hyperparameters (Figure 9a–d). The natural break method was used to categorize spatial GWQ into five classes (very high, high, moderate, low and very low). The areas of the five classes in the four GWQ maps are displayed in

Table 7. Areas of the five classes in the four GWQ maps.

Areas	Selection 1	Selection 2	Selection 3	Selection 4
Very high (km2)	2625.53	1278.48	1984.58	2261.95
High (km2)	3736.37	4915.00	3261.57	3981.42
Moderate (km2)	3559.01	6026.96	4728.50	5181.86
Very low (km2)	5242.62	4353.83	4695.23	4305.54
Low (km2)	3791.72	2380.98	4285.36	3224.49

. It is found that despite similar model performances, there are significant differences in spatial distribution for the four maps. Therefore, we averaged the possibility of each grid cell from the spatial GWQ maps of the four random selections (Figure 9e). This may be a possible strategy to address the spatial uncertainty caused by selecting different datasets.

4.4. Spatial Uncertainty in GWQ Mapping

Figure 10 shows the spatial uncertainty analysis results with PCC values. It is found that PCCs among these four maps range from 0.5365 to 0.8066, confirming the observed differences in Figure 9. Particularly notable is the variance between selection 2 and 3, with a PCC of 0.5365. Even for the most similar pair of selections (1 and 4), the PCC is only 0.8066. Figure 10g displays the final results of spatial uncertainty with an average PCC of 0.6707. It indicates that the spatial uncertainty caused by four random selections of training and validation groundwater samples is pronounced. We have highlighted three typical areas with particularly high uncertainty on the map, providing a basis for supplementing groundwater samples in future groundwater management.

4.5. Indicator Analysis with Importance and SHAP Value

Figure 11 shows the results of feature importance for different randomly selected training and testing sets, based on the Python code “feature importance”.

Table 8. Accumulated importance of 13 indicators for four random selections.

Indicators	Accumulated Importance	Rank	Proportion
Population	74.22	1	18.55%
Nighttime light	70.59	2	17.65%
Aquifer media	41	3	10.25%
GDP2015	39.49	4	9.87%
Groundwater yield	34.37	5	8.59%
Conductivity	31.23	6	7.81%
Change of NDVI	28.85	7	7.21%
Depth to groundwater	25.21	8	6.30%
PPSD	21.32	9	5.33%
Topography	15.75	10	3.94%
LULC	9.29	11	2.32%
Degree of urbanization	5.99	12	1.50%
Net recharge	2.7	13	0.67%

shows the accumulated importance of 13 indicators for four random selections. It is evident that population (18.55%) and nighttime light (17.65%) are the most critical indicators, leading over aquifer media, GDP2015, and groundwater yield. After removing population and nighttime light and rerunning the model, a significant decrease in model performance was observed (Appendix B). In contrast, LULC (2.32%), degree of urbanization (1.50%) and net recharge (0.67%) have relatively low importance. Interestingly, both the most and least important indicators include those related to human activities.

Figure 12 shows the SHAP analysis results from four random selections, and high and low feature values (indicator values) are represented by red and blue, respectively. The broader the distribution of an indicator’s feature values, covering a wider range of SHAP values, the more important that indicator is considered. It is found that the indicator importance observed here is consistent with that shown in Figure 11 by “feature importance” code. By analyzing the relationship between feature values and SHAP values, we found that population density, nighttime lights, and GDP, which are theoretically positively correlated, exhibit diverse distributions (Figure 12). Specifically, the contributions of nighttime lights and population to GWQ prediction are opposing, with the areas of poorest GWQ (high SHAP values) being those with higher population density but not necessarily high GDP. Also, the areas with significant NDVI changes and high PPSD exhibit poor GWQ (positive SHAP values), which aligns with the expected outcomes. Additionally, the SHAP analysis results are generally consistent with the scores from the DRASTIC model for groundwater vulnerability, except for the conductivity indicator.

5. Discussion

5.1. Discussion on GWQ Assessment

The EWQI offers a refined approach compared to traditional composite GWQ assessment methods like the CPI (GB/T 14848-2017) and the Nemerow index [126]. While the CPI often penalizes overall water quality for a single poor parameter due to its averaging approach, the EWQI mitigates this by using entropy to weight parameters based on their variability and significance. This results in a more balanced and realistic assessment of GWQ [127]. In contrast to the Nemerow index method, which tends to emphasize the worst-case scenario, EWQI provides a broader perspective, integrating various indicators without letting a single outlier skew the overall results [128]. This makes the EWQI particularly useful for creating targeted and effective water management strategies. In fact, if multiple methods are used to evaluate groundwater samples, more comprehensive analysis results can be obtained, and such comparative studies can serve as a future research direction [129,130].

Unlike other studies based on the EWQI [68,126], this paper further considers the quality of individual parameters on the basis of the NI approach. For the single GWQ parameter, the NI value reflects the priority level of management for the petameter. NO₃⁻ has the highest priority, which is likely highly related to the extensive agricultural land and the use of fertilizers in the study area [56,131,132]. The relatively high NI of Mg²⁺, SO₄²⁻, and Na⁺ is primarily linked to specific hydrogeochemical processes such as rock weathering, cation exchange, and evaporation [57,133,134]. Additionally, the significant presence of SO₄²⁻ may be attributed to the oxidation of pyrite [135]. It is noted that the NI of Cr⁶⁺ reached 0.109, highlighting the serious health risks it poses when concentrations exceed safe levels in groundwater. The elevated Cr⁶⁺ concentrations in groundwater primarily stem from low groundwater velocity in the loess aquifer, cation exchange in alkaline environments, and industrial activities [136]. In terms of other common GWQ parameters such as HCO³⁻, Zn²⁺ and F⁻, Fe³⁺, Al³⁺, comparing NI values, weights and exceedance rates can provide managers with prioritized information for management. The managers can develop varying levels of macro-strategic groundwater management plans based on their circumstances, such as economic and policy factors.

5.2. Model Performance and Spatial Uncertainty

In this study, all four selections demonstrated high model performance, with AUROC values of around 0.9. However, spatial differences and uncertainties were evident (average PCC = 0.6707). In fact, some studies have demonstrated that similar model performances do not necessarily imply similar spatial distributions in groundwater potential mapping [102] and landslide susceptibility mapping [111]. This is primarily because of the training and testing sets being sourced from specific geographic locations. When these models are generalized spatially, the diversity of indicators and the variation in optimal hyperparameters can lead to inconsistencies in spatial distribution. However, in ML-based spatial predictions, the dataset used for training is inherently limited. It is challenging to guarantee that spatial samples (e.g., groundwater samples) are both sufficient and evenly distributed, which implies that spatial uncertainty cannot be fully eliminated. Therefore, considering that sampling is both time-consuming and expensive, balancing the number of sampling points with spatial uncertainty is an important direction for future research.

Many studies have employed k-fold cross validation to reduce uncertainty in ML models [137,138,139,140]. This technique enhances model reliability by ensuring robustness and consistency across different data subsets [141]. However, the effectiveness of k-fold cross-validation in addressing spatial uncertainty is still very limited. This study proposed a possible method that averaged the grid cells from the spatial GWQ maps of four selections. This averaging approach reduces variability and enhances the stability of spatial predictions by mitigating the effects of outliers and random sampling errors. Also, we have highlighted three typical areas with particularly high uncertainty in Figure 10, identifying potential locations for additional groundwater sampling. With more sampling data, the spatial uncertainty of the groundwater quality map is expected to decrease accordingly. However, these methods serve as a starting point, and the current research on spatial uncertainty is notably insufficient. Under this condition, we strongly encourage the development and discussion of more solutions to address spatial uncertainty. Such approaches should be extended beyond GWQ to a broader range of ML-based spatial prediction applications like landslide susceptibility [142], groundwater salinity [143], groundwater potential [72] and nitrate concentrations [19].

5.3. SHAP Observation and Discussion

SHAP analysis is one of the most important EMLTs, and this study further verifies its applicability in spatial GWQ mapping and prediction. In very poor GWQ areas (high SHAP value), the feature value trends of population density and GDP2015 show the variations rather than the expected similarity. The population is depicted in red in poor GWQ areas, while GDP2015 shows the different pattern. The red points for the GDP2015 indicator mainly appear at slightly positive SHAP values, indicating that regions with severe groundwater pollution are not necessarily high-GDP areas. This inconsistency further highlights the importance of considering both the positive and negative impacts of these human related indicators in the analysis. This is why we use the term “potential” in Section 3.3.2 for indicator selection. However, it is found that when it comes to population growth and economic development, many researchers assume their impact on groundwater is primarily groundwater deterioration or some other negative effect [144,145,146,147]. Here, we want to emphasize that the positive impacts of these indicators on GWQ should not be ignored. Given the current lack of detailed studies exploring the impacts of economic or population indicators on GWQ, we propose two viewpoints. First, we hypothesize that the impact of these indicators on GWQ may exhibit an inverse U-shape, similar to findings in studies on greenspace and economic growth [148] or population aging and economic growth [149] and carbon emissions and population size [150]. This may be because as economic or population development reaches a certain level, corresponding groundwater protection and remediation measures are likely to improve, and these indicators may shift from being “pressure” indicators to “response” indicators. Second, we also strongly recommend conducting more research to find evidence that supports the hypothesis or to further investigate the relationship between economic and population indicators and GWQ.

We observed that in areas with poor GWQ (positive SHAP value), the feature value of nighttime light is low while PPSD is high, whereas in areas with good GWQ, the situation is reversed. Nighttime light is often a complex factor that often correlates with GDP, LULC, population, PPSD, urbanization, and other socio-economic indicators (Figure 7). Considering the possible inverse U-shape relationship observed in the analysis of GDP and population, along with the low importance of LULC and degree of urbanization, we have made a reasonable inference regarding the causes of poor GWQ in the study area. Point-source pollution (e.g., industry, farms, mine exploration, hazardous waste disposal sites, and landfills) is the primary cause of GWQ deterioration. The pollution in these areas is characterized by concentrated, localized contamination, often confined to specific sites. Additionally, these areas have low nighttime light intensity and are less influenced by land use and urbanization, which further supports the positive correlation with the PPSD indicator. The impact of point-source pollution on GWQ explains why there is a substantial difference between the spatial distribution map of GWQ created in this study and the groundwater vulnerability map created in current study [55]. This also implies that when using specific pollutant parameter (NO₃⁻) to validate groundwater vulnerability models, it is crucial to consider and mitigate the influence of point-source pollution on the results. Therefore, for these potential pollution sources, conducting regular GWQ testing, enhancing wastewater treatment facilities, and providing education on best practices for pollution prevention are top priorities for protecting groundwater in the study area. Additionally, implementing stricter regulations on industrial discharges and monitoring land use changes can further mitigate the risks posed by point-source pollution.

The poorer GWQ in areas with significant NDVI changes indicates that human expansion activities over the past decade have generally degraded GWQ. The low importance of the degree of urbanization further suggests that it is the process of human activity or expansion, rather than the presence of constructed areas themselves, that has led to the deterioration. Human expansion activities often include deforestation, land conversion for agriculture, industrial development, and infrastructure construction [151,152]. These processes contribute to GWQ degradation by disturbing soil and vegetation during land clearing, increasing surface runoff, introducing pollutants before adequate infrastructure is in place, and disrupting natural water recharge areas. Although groundwater protection measures may improve following urban expansion, it is crucial to regulate human activities, particularly by controlling land use changes, implementing sustainable development practices, and enforcing stricter environmental regulations during the human expansion process to address this issue.

In terms of the “state” indicators group, the SHAP analysis results are generally consistent with the scores from the DRASTIC model for groundwater vulnerability [153,154,155], except for the conductivity indicator. This can also be explained by point-source pollution. In cases of point-source pollution, low conductivity indicates that the polluted groundwater does not easily disperse, thereby affecting regional GWQ. In areas with high conductivity, pollutants from the source diffuse with the water flow, diluting the contaminants and resulting in improved GWQ. Based on the SHAP analysis, the preliminary causes of GWQ condition in the study area have been identified. The next step is for managers to conduct detailed investigations according to these preliminary hypotheses and develop corresponding management strategies. However, it is important to emphasize that SHAP analysis only explains the association between indicators and outcomes from a statistical perspective and does not necessarily imply a definitive causal relationship, particularly for complex indicators like nighttime light. The inferences we make based on SHAP results require further validation through concrete evidence. Causal analysis, as explored by Jia et al. [156], is a direction worth pursuing in future research following SHAP analysis.

5.4. Limitations and Future Research

Although the findings of this study are interesting, there are four limitations to consider. Firstly, for a study area of nearly 20,000 km², the use of 180 sample points is somewhat insufficient for both training and validation sets, which may be a significant reason for spatial uncertainty. However, this issue is common in groundwater studies due to the limited economic and time resources available for extensive sampling. This also highlights the value of our research on spatial uncertainty and provides a foundation for locating where additional groundwater bores need to be drilled. Secondly, although SHAP analysis is an ELMT, the explanation provided is statistical. For some complex socio-economic indicators, the causal relationships between the indicators themselves and between the indicators and the outcomes have yet to be confirmed. Thirdly, this study only used the LightGBM model with an 80/20 split ratio of training and validation datasets for spatial uncertainty analysis. Introducing more models and a broader range of data-split comparisons could increase the stability of the results. Fourth, when determining GWQ standards, we used Chinese, international, and other literature standards, which might introduce some bias to the results. Nonetheless, these limitations have a minimal impact on the findings of this study. Based on this study, future research should include the following six aspects:

• When evaluating GWQ, it is recommended to use multiple methods, including the EWQI, the CPI, and the Nemerow index, and to promote the single parameter analysis method of the NI proposed in this study.
• It is encouraging to confirm the causal relationships between indicators and between the indicators and outcomes, ensuring that the associations identified through SHAP analysis are supported by robust evidence.
• Introduce and compare more models, including deep learning, reinforcement learning, and ensemble learning, to enhance the stability and accuracy of the results.
• Further promote the contribution of the PSR framework in spatial mapping and prediction for indicator selection to ensure the completeness of model construction.
• In addition to calculating spatial average probabilities and supplementing with additional groundwater samples, develop more methods to reduce spatial uncertainty to provide managers with more accurate mapping results.
• Further develop the application of EMLTs in groundwater management.

6. Conclusions

The spatial mapping and prediction of GWQ is essential for identifying pollution sources and informing comprehensive groundwater management strategies. However, this area has not yet been fully explored. The research gaps mainly include the inaccuracy of traditional spatial interpolation for spatial mapping, insufficient consideration of the geological environment and human activities in ML models, the limitation to single pollutants, and the lack of a systematic approach in the selection of indicators. By taking Guanzhong Plain as a case study, this study utilized the EWQI, the LightGBM model, the TPE optimization method, the PSR framework, and SHAP analysis for the spatial mapping and prediction of GWQ, aiming to address the aforementioned research gaps. Through analysis and discussion, we have made several interesting and important findings.

Firstly, according to the NI results for various parameters, NO₃⁻, Mg²⁺, SO₄²⁻, Na⁺ and Cr⁶⁺ should be prioritized for remediation. The skewed distribution of the EWQI indicates that the overall GWQ in the study area is generally good, but a few areas suffer from severely poor quality, which warrants urgent attention. Secondly, based on four randomly selected training and testing sets, although their model performances were high (with AUROC around 0.9), they exhibited spatial uncertainty, with the lowest spatial correlation being only 0.5365 (between selection 2 and 3). This issue is not limited to spatial GWQ mapping and prediction but also extends to other fields. The spatial averaging method and additional groundwater samples may be possible solutions for this issue, but further methods need to be explored. Thirdly, population and nighttime light are the most critical indicators, while the indicators of net recharge, LULC and degree of urbanization have the lowest importance. Combining SHAP values, we infer that economic development and population have both positive and negative impacts on GWQ, while point-source pollution is the main cause of the decline in GWQ in the study area. Additionally, we speculate that human expansion activities over the past decade have generally had a negative impact on GWQ.

Due to the limited research on the spatial mapping and prediction of GWQ, future studies should include six different aspects in this field, involving multi-method GWQ assessment; causal relationships between indicators and between the indicators and outcomes; the introduction and comparison of more spatial mapping and prediction models; the application of the PSR framework for indicator selection; the development of more methods to reduce spatial uncertainty; and the application of EMLTs in groundwater management. In this way, future research will support the development of the spatial mapping and prediction of GWQ from different perspectives, aiming to further assist groundwater managers in achieving sustainable groundwater management in the future.