Shenzhen is a coastal city located in southeast China, Guangdong Province, lying between 113°46′–114°37′E and 22°27′–22°52′N (Figure 1). The study region included 11 watersheds covering approximately 2000 km² area (Qiu et al., 2022). The region has a subtropical marine monsoon climate with an annual average precipitation of 1,830 mm. Precipitation is mainly concentrated from April to October, accounting for approximately 85% of the annual rainfall (Shenzhen Statistics Bureau, 2020). The study area terrain was high in the southeast and low in the northwest, mostly in low-hilly areas with gentle terraces. Groundwater types include quaternary loose deposits, pore water, bedrock fissure water, and karst water. The groundwater level was shallow and generally in the range of 1.5–4.5 m.

Groundwater quality data used in this study were obtained from the Water Authority of Shenzhen Municipality (http://swj.sz.gov.cn/). Water samples were collected from 30 monitoring wells in the study area from 2015 to 2021, covering the wet and dry seasons (684 samples; Supplementary Table S1). In this study, March and December represented dry season. June and September represented wet season. Locations of the sampling points are shown in Figure 1. The groundwater quality data included levels of pH, total dissolved solids (TDS), total hardness (TH), chloride (Cl⁻), fluoride (F⁻), ammonia-nitrogen (NH₃-N), nitrate (NO₃-N), nitrite (NO₂-N), sulphate (SO₄ ^2–), arsenic (As), total iron (Fe), total manganese (Mn), lead (Pb), cadmium (Cd), mercury (Hg), and dissolved oxygen (DO). Thirty monitoring wells were uniformly distributed throughout the study area.

Figure 2 shows the framework of the WQI-BMA model integrating four traditionally used WQI models based on the BMA method to assess groundwater quality. In this framework, the implementation of the WQI-BMA model involves two stages. First stage involved the construction of four traditionally used WQI models as members of the BMA. In the first stage, the groundwater quality parameters were determined. All measurable groundwater quality parameters were used for comprehensive groundwater quality assessment. Once the groundwater quality parameters were determined, a machine learning approach, called the extreme gradient boosting algorithm, was used to systematically assign weights based on the rank order centroid weight method. The sub-index functions are then calculated based on a uniform calculation rule. The final step was to calculate the aggregation functions based on the four WQI models, and the results were then used for stage two, BMA merging. In stage two, BMA merged the four WQI models from stage one to quantify their uncertainty and provide a final groundwater quality assessment. By combining WQI methods and Bayesian model, the WQI-BMA model can provide more accurate and reliable predictions of water quality. This is due to the fact that the WQI-BMA model combines the strengths of multiple WQI models. The Bayesian model can take into account the uncertainty of the WQI models, which can be used to generate more accurate predictions. Additionally, the WQI-BMA model is more robust and less sensitive to outliers than any single WQI model.

WQI models are widely used tools for assessing the water quality. They can analyze large amounts of temporally and spatially varying water quality data to inform decision-makers about the water quality. The construction of a WQI model typically involves four steps. First step is to collect and analyze the water quality data to determine the number of parameters needed to evaluate the water quality. In the second step, parameter weights are assigned based on the importance of the water quality. In the third step, these parameters are normalized to a 0–100 scale using sub-index functions. Finally, the normalized sub-indices and parameter weights are combined into an overall WQI value using a suitable aggregation function, and the water quality is inferred using a classification scheme that defines the water quality classes.

Parameter selection is the first step for the WQI model. The parameters are usually selected based on data availability, expert opinion, and the environmental impact of water quality (Gupta and Gupta, 2021). In this study, 16 parameters that commonly reflect the water quality were considered for the WQI models. These 16 parameters included the levels of pH, TDS, DO, TH, Cl⁻, F⁻, NH₃-N, NO₃-N, NO₂-N, SO₄ ^2–, and heavy metals (As, Fe, Mn, Pb, Cd, and Hg).

Parameter weight is a key component of WQI models. It is typically determined based on the environmental importance of the groundwater quality parameters, together with the recommended guidance value. Such parameter weight allocations are usually biased towards the subjectivity of experts and do not follow any systematic technique. In this study, an objective weighting approach based on machine learning (ML) (Uddin et al., 2022) was employed to assign weights to groundwater quality parameters. This approach differs from the subjective approach used in the literature (e.g., expert opinion) and aims to reduce the uncertainty caused by inappropriate weighting. Two processes are involved in the objective weighting approach: i) extreme gradient boosting (XGBoost) ranking of the parameter, and the ii) rank order centroid (ROC) weight method, which attributes weightings based on rank.

ML can quantify the importance of input variables based on their influence on model output variables (Uddin et al., 2022). Similarly, it can be used to quantify the importance of water quality parameters. In this study, 16 water quality parameters at 30 monitoring wells were used as inputs for the ML model, and the corresponding groundwater quality status estimated based on the Environmental Quality Standards for Surface Water (GB3838-2017) was used as the output. Then, an initial assessment was conducted using several ML methods according to 5-folds of the k-fold cross-validation, including decision tree (DT), random forest (RF), naive Bayes (NB), extremely randomized trees (ERT), and XGBoost algorithms. All ML methods in this study were applied using Python-based open-source scikit-learn codes.

In this study, an approach combining the XGboost and ROC weight methods was used to assign weights to provide a better allocation than the expert opinion method (Uddin et al., 2021; Uddin et al., 2022). XGBoost is a widely used ensemble ML method based on the gradient boosting algorithm (Zhu et al., 2022). Some studies have used XGBoost to extract key variables, including water quality (Naghibi et al., 2020), air pollution (Li et al., 2022), and wastewater treatment (Wang D. et al., 2022), for developing new models. XGBoost has lower prediction errors than other algorithms as its regularized boosting technique controls the complexity of the model and prevents overfitting by adding regularization terms to the loss function. XGBoost can rank features based on their relative importance using SHapley Additive exPlanations (SHAP) values. SHAP values are local interpretation methods that provide more detailed, additional, and individualized explanations and attributions than global interpretation techniques (Li, 2022). This study ranked the relative importance of the groundwater parameters based on XGBoost and SHAP methods for subsequent weight calculations.

Hyperparameters are important for determining the structure and prediction accuracy of the model. In this study, grid search optimization method (Abbaszadeh et al., 2022) was used to determine the hyperparameters in the XGBoost model to achieve high prediction accuracy. XGBoost performance was evaluated using the accuracy and R ² of the training and test datasets, where the input data were randomly divided into two groups: 80% of training dataset and 20% of testing dataset.

ROC weight method estimates the weights that minimize the maximum error of each weight by identifying the centroids of all possible consequences, while maintaining the rank order of objective importance. ROC weights of a set of n variables ranked from j = 1 to n were calculated using the following equation (Roszkowska, 2013): w j = 1 n ∑ i = j n 1 r i (1) where w _j is the jth weight and r _i is the rank position of w _j.

Sub-index converts the measured groundwater quality values into unitless values. Sub-index values (q _j . ranging from 0 to 100) were obtained using the Delphi technique (Bordalo et al., 2001), which was used to indicate whether the variable had a poor or good status. In this study, the linear interpolation rescaling approach by Uddin et al. (2022) was used to compute the sub-index for all groundwater variables using Eqs. (2) and (3): q j = S I u − S I l − S I u * C j S I T u − S I T l (2) q j = C j − S I T l S I T u − S I T l * S I u (3) where SI _l and SI _u are the lower and upper limits of the possible q _j values (0 and 100, respectively); SIT _l and SIT _u are the lower and upper threshold values, respectively; and C _j is the measured value of the water quality parameter. Threshold guidelines are listed in Supplementary Table S2.

Final step of the WQI model involves the aggregation function, which converts multiple water quality parameters into a single value to express the overall water quality status. Several studies have shown that different results from different WQI models are related to the aggregation function (Lumb et al., 2011; Juwana et al., 2016; Ma et al., 2020). In this study, we considered four common traditional weighted aggregation functions in the WQI models: NSF, WQM, SRDD and WJ (Table 1). The output of the final WQI model was recommended for the four water quality classes listed in Table 2, according to Uddin et al. (2022). A WQI value in the range of 0–29 indicates that the water quality is very poor, and the water cannot be used for any purpose. A WQI value in the range of 30–49 indicates that the water can only be used for a small number of specific cases, with a high the risk of water pollution. A WQI value in the range of 50–79 indicates that the water body has minor pollution, with some indicators meeting the water quality standards. A WQI value in the range of 80–100 indicates that the water quality is good, satisfying most of the water quality standards.

BMA is a statistical method for multi-model ensembles that has higher reliability and accuracy than the single-model method (Duan et al., 2007). BMA results are usually better than those of single models in terms of skill and reliability (Ashofteh et al., 2022). BMA has been used in many studies for groundwater vulnerability assessment (Gharekhani et al., 2022), mortality forecasting (Ashofteh et al., 2022), soil moisture estimation (Chen et al., 2020), merging of precipitation products (Yumnam et al., 2022), and streamflow simulation (Wang J. et al., 2022). BMS is a free R package for performing BMA using the open-source software R. In this study, the BMA model was applied using the BMS package in R version 4.1.2. For more information refers to BMS package is available at http://bms.zeugner.eu.

Based on Bayes’ theorem and the total probability law, posterior distribution of the BMA predicted quantities Δ) for the given data D) is (Hoeting et al., 1999): p Δ D = ∑ k = 1 K p Δ M k , D p M k D (4) where p (Δ|D) is the probability of the prediction of data Δ) given the observed data D); p (Δ|M _k ,D) is the conditional probability of the predicted quantity given the observed data D) and given model (M _k ), which is the output of SRDD, NSF, WJ, and WQM models; and p (M _k |D) is the posterior probability of the model given the data D), which is the model weight. Posterior probability of the model (M _k ) is calculated using the following equation (Gharekhani et al., 2022): p M k D = p D M k p M k ∑ i = 1 K p D M i p M i (5) where, p D M k = ∫ p D θ k , M k p θ k M k d θ k (6)

Here, p (D|M _k ) is the integrated likelihood of model M _k ; θ _k is the vector of parameters of model M _k ; p (θ _k |M _k ) is the prior density of θ _k under model M _k ; p (D|θ _k , M _k ) is the likelihood, and p (M _k ) is the prior model probability for the model M _k .

Statistical analysis results of the groundwater quality parameters in the study area are listed in Table 3. The parameters include many variables, such as the minimum value (Min), maximum value (Max), mean, standard deviation (SD), and coefficient of variation (CV). NH₃-N concentration values ranged from 0 to 56.40 mg/L with a mean value of 1.40 mg/L, which exceeds the limit of 0.5 mg/L set by the World Health Organization. pH varied from 5.14 to 8.60 with a mean value of 6.84, indicating that most groundwater within the study area was slightly acidic. Mean concentrations of As, Fe, and Mn were 0.012, 0.49, and 0.20 mg/L, respectively. Mean concentrations of the other parameters in the groundwater were less than the thresholds proposed in the WHO and Environmental Quality Standards for Surface Water (GB3838-2002) guidelines (Supplementary Table S2). CV for heavy metals ranged from 618.94% to 140.70% in the order As > Pb > Fe > Cd > Mn > Hg, indicating a large difference in the distribution of heavy metals in groundwater. SD of the groundwater quality data indicated high variability; therefore, variation in the mean value of the data was also high (Elbeltagi et al., 2022).

Supplementary Table S3 lists the findings of an initial assessment conducted using several ML methods (DT, RF, NB, ERT, and XGBoost) according to 5-folds of the k-fold cross-validation. Accuracy is a metric for evaluating ML classification problems and indicates the ratio of all correctly classified samples to the total number of samples. Among the 4 ML methods, the XGBoost model exhibited the highest prediction accuracy (83%). Prediction accuracies of DT (82%), RF (82%), ERT (77%), and NB (65%) were slightly lower than that of XGBoost. Thus, the XGBoost model was selected to rank the parameters according to the relative importance of the selection of important groundwater quality parameters. Table 4 lists the final adjusted hyperparameter values for the XGBoost model obtained using the grid-search optimization method. Combination of the hyperparameters in Table 4 resulted in the highest R ² and accuracy in the predicted and observed values. A maximum tree depth of seven and a learning rate of 0.1 were used to increase the maximum effect of interactions between variables. The training and testing accuracies of the XGBoost model were 1 and 0.95, respectively, and the training and testing R ² were 0.99 and 0.92, respectively. Results of the selection of important groundwater quality parameters using the XGBoost model integrated with SHAP values are shown in Supplementary Figure S1. At the study site, the most important factors influencing the water quality status were NH₃-N and Mn (Supplementary Figure S1). The order of importance of the groundwater parameters was NH₃-N > Mn > pH > NO₃-N > TDS > Pb > Fe > As > Hg > DO > SO₄ ^2– > Cd > Cl⁻ > F⁻ > TH > NO₂-N.

Table 5 lists the weight values of the groundwater quality parameters calculated using the ROC method. Weight values followed the same order as the rankings of the groundwater quality parameters, with NH₃-N, Mn, and pH exhibiting the highest weight values and F⁻, TH, and NO₂-N exhibiting the lowest weight values.

In this study, four WQI models based on different weighted aggregation functions were used to calculate WQI values and infer the groundwater status. Figure 3 shows the range of the WQI values calculated for each function within the study area. Overall, there were significant differences in the WQI values obtained using different WQI models (Figure 3). WQM provided the smallest distribution of scores and the highest mean score. Conversely, WJ provided the largest distribution of scores and the lowest mean score. The mean score for NSF was lower than that for WQM but greater than that for SRDD. For the NSF and WQM models, the groundwater quality evaluation grade was higher than that for the SRDD and WJ models. In addition, most of the groundwater quality in the study area was above the grade “poor,” indicating that the groundwater quality in the study area was in the marginal to good level. In this study, all WQI models used the same groundwater quality parameters and sub-index functions, except the aggregation function, to calculate the WQI values. Therefore, the differences were caused by the different aggregation functions.

Figure 4 shows the groundwater quality classification results for WQI-BMA value distribution predicted using BMA. Groundwater quality status can be categorized into four types based on WQI values: poor (0–29), marginal (29–49), fair (50–79), and good (80–100) (Table 2). WQI-BMA values ranged from 35.01 to 98.45 based on BMA prediction. According to the groundwater quality classification of WQI values (Table 2), the groundwater quality status was classified into three categories: marginal, fair, and good. Based on the WQI-BMA model, 17.6% of the groundwater fell under the marginal category, 45.7% fell under the fair category, and 36.7% fell under the good groundwater quality category.

Figure 5 shows the percentages of groundwater quality classification in the study site over the dry and wet seasons from 2015 to 2020. Percentage of groundwater quality classification refers to the percentage of wells classified as marginal, fair, or good. During the period 2015–2020, the fairness category had the highest percentage of all categories, ranging from 43.3% to 83%, with an annual average of 54.7%. In contrast, the marginal category had the lowest percentage, ranging from 3.3% to 40%, with an annual average of 16.7%. The percentage of good categories was between fair and marginal, ranging from 13.3% to 43.3%, with an annual average of 28.6%. Regarding seasonal variation, there were differences in the groundwater quality status over the dry and wet seasons during the study period (2015–2020). For example, the seasonal variation in the groundwater quality in the fair category was small, ≤7%, except in 2015. However, groundwater quality in the marginal and good categories varied the most in 2020 at 20% and 13%, respectively.

Geographic information system-based spatiotemporal analysis was used to determine the spatiotemporal distribution of the WQI map in the study area from 2015 to 2020 (Figure 6). Spatially, most sites fell under the fair-to-good category, with a few scattered areas falling under the marginal category, indicating that groundwater quality of the study area has been well maintained over the past few years. Specifically, the best areas for groundwater quality classification are the eastern parts of the city, such as Dapeng and Yantian. This area is mostly mountainous and hilly, with abundant forest resources. Due to the lack of heavy industry, the maintenance of groundwater quality is better. There are a total of 12 monitoring wells in the region. In 2015, there were two wells fell under the good category, nine wells fell under the fair category, and one well fell under the marginal category. With the development of smart water management, groundwater quality has gradually improved by 2019. There were seven wells fell under the good category, four wells fell under the fair category, and one well fell under the marginal category. The western region of the study area is an industrial agglomeration, the central part is the city center area with dense population, and the quality of groundwater is worse than that of the eastern region. Subsequent monitoring of groundwater quality should be strengthened in the central and western regions, especially in the coastal areas of the west. None of the samples in the study area fell under the poor and inappropriate for drinking categories. These results indicate that during the scarcity of surface water, groundwater can be used as a recharge source. Overall, the spatial distribution of groundwater quality in the study area showed that most of the groundwater in the study area was above the fair category and suitable for most applications.

Temporally, the groundwater quality in the study area showed seasonal fluctuations from 2015 to 2020. Relatively, the number of wells in the marginal groundwater quality category first increased, then decreased, and subsequently increased over time. Notably, 2020 had the highest number of sites in the marginal category. The number of wells in the good groundwater quality category first increased and then stabilized, followed by a decrease over time. The number of wells in the good groundwater quality category increased significantly compared to those in 2015. In contrast, the fair groundwater quality category remained at a relatively stable level. Overall classification of groundwater in the study area was mainly above the marginal category and suitable for domestic and agricultural use. Basic data for the model input were groundwater monitoring data. For emerging pollutants, the monitoring data of the relevant pollutants can be collected and introduced into the model for evaluation; however, currently, there is no threshold standard for the limit value of emerging pollutants in groundwater. Therefore, threshold standards need to be determined in the future along with relevant toxicological studies.