Groundwater quality evolution across China

Domain: earth_science
Published: 2025-03-14
Editorial outcome: major_revision

Abstract
China is facing a severe groundwater quality crisis amid economic development and climate change, yet the extent and trajectory of this crisis remain largely unknown. Here we developed a machine-learning model, incorporating natural and social-economic factors, to construct annual probabilistic maps of poor groundwater quality (PGQ, i.e., Class V based on the Chinese groundwater quality standard) across China from 1980 to 2100. Alarmingly, our findings indicate a concerning escalation in PGQ area ratio, rising from 17.3% in 1980 to 30.1% in 2000, and surging to 40.8% by 2020, adversely affecting 6.8%, 17.5%, and 36.0% of the Chinese population, respectively. The predominant drivers of this degradation were identified as agricultural discharge (contributing to 10.7% growth in PGQ area ratio), followed by groundwater exploitation (5.6%), industrial discharge (5.3%), domestic discharge (1.7%), climate change (0.5%), and land use change (-0.3%). By 2050, the PGQ area ratio could range from 37.9% to 48.3% under different socio-economic and climate scenarios. Our study highlights the urgent need for effective water resources management and conservation measures to mitigate the deteriorating trend of groundwater quality and address the challenges posed by socio-economic development and climate change, thereby safeguarding water security for China and the global community.

Introduction
Groundwater provides almost half of all drinking water worldwide, particularly in rural areas1. Yet, it faces substantial threats from socio-economic development and climate change2,3. Serving as a crucial source of drinking water for over 400 cities across China, groundwater is especially vital in the northern regions where it accounts for two-thirds of the drinking water, half of the industrial water, and one-third of the irrigation water4,5. As the world’s largest industrial powerhouse and a major agricultural producer, China has experienced considerable groundwater quality problems over recent decades6. Extensive pollution discharges from agricultural and industrial activities, coupled with large-scale over-exploitation, have led to widespread groundwater pollution across the nation7,8. In 2020, a report by China’s Ministry of Ecology and Environment, indicated that 33.7% of the tested 10,242 sites (dominated by shallow groundwater) had “marginal” groundwater quality, while another 43.6% were classified as “poor”9.

In response to the groundwater quality crisis, nations including China have developed environmental monitoring networks aimed at assessing and mitigating the negative impacts of contamination and over-exploitation on this vital resource10. These networks are important for the early detection of pollutants, comprehending groundwater flow dynamics, and formulating strategies to increase groundwater sustainability11. However, a key limitation of these networks is their sparse spatial distribution of monitoring sites, which fail to provide comprehensive national coverage, often both costly and time-consuming12. This deficiency is particularly acute in emerging economies like China, where the low density of monitoring sites inadequately captures the spatial variability of groundwater quality and pollution incidents, compromising the accuracy and completeness of monitoring efforts13. Furthermore, considerable challenges persist in comprehensively quantifying the driving factors behind groundwater quality deterioration and in projecting future changes.

To address the lack of comprehensive data for groundwater quality assessment, researchers have increasingly utilized machine learning (ML) to enhance the mapping of contaminants like arsenic and fluoride, which often stem from natural geological sources14,15. Pioneering works by Amini et al.16 and Rodríguez-Lado et al.17 have incorporated environmental factors into ML models to create predictive maps. These efforts are notable for their ability to interpret complex environmental data—including climate, soil, geology, and topography—into useful predictive tools for assessing groundwater quality18–21. However, these studies primarily focused on static environmental factors as predictors, which may not sufficiently capture the dynamic and evolving impacts of human activities and climate change on groundwater systems.

In this study, we employed extensive groundwater survey data, along with natural and socio-economic variables across China, to develop an ML model for predicting groundwater quality. Utilizing this model, we mapped the annual probabilities of poor groundwater quality (PGQ, i.e., Class V based on the Chinese groundwater quality standard22; Table S1) from 1980 to 2100, under various future scenarios of socio-economic development and climate change. This approach allowed us to illuminate evolving patterns of groundwater quality and assess the profound influence of human activities and climate change on China’s groundwater quality. Our analysis offers crucial insights and guidance for crafting strategies to address the environmental threats that China’s groundwater resources are facing. We underscore the importance of adapting to changes in development patterns and climate conditions to effectively mitigate these challenges.

Results
Random forest modeling
In our study, we extracted groundwater quality data from 1977 published surveys, culminating in a dataset. This dataset incorporated geospatial information, including location and land-use type, as well as temporal data represented by the sampling year. The spatial distribution and temporal variation of the survey data are presented in Fig. S1. Recognizing the need for a robust and unbiased model, we systematically segregated these 1977 groundwater quality surveys into two distinct datasets: 90% for training and 10% for validation. To counterbalance any potential issues arising from the relatively small sample size, we utilized a strategic data augmentation technique. This involved the generation of modified duplicates of the existing data, thereby enhancing our training dataset. The augmentation of the validation dataset was conducted independently of the training dataset, ensuring no overlap and maintaining the integrity of our model’s validation process.

To establish a quantitative relationship between predictors and groundwater quality (PGQ or non-PGQ), we employed a non-parametric supervised ML technique known as random forest (RF). To develop a parsimonious classification model, an initial selection of 51 potential predictor variables underwent the elimination of redundant variables (Table S2). This process of variable pruning, which was guided by evaluating the collinearity and independence (Fig. S2), ensured that only those variables that did not compromise the prediction accuracy of subsequent RF models were retained (“Model predictor selection” in Methods). Addressing multicollinearity in these parameters is crucial, as it helps prevent adverse effects such as overfitting, reduced interpretability, and ultimately, compromised predictive performance of the model23. The finally retained 25 predictors (indicated in bold in Table S2) includes factors related to soil properties, geographical and hydrogeological conditions, climate change, groundwater exploitation, pollution discharge, and land-use change. The out-of-bag error of each predictor, reflecting its impact on the model’s predictive accuracy24, is used to assess the importance of each predictor (Fig. S3). Five most notable predictors are identified: Depth-based groundwater type, air temperature, aridity index, precipitation, and sand content of soil (100–200 cm).

The performance of the RF model on the validation set (10% of the data, which was randomly selected while maintaining the relative distribution of PGQ and non-PGQ) is summarized in the confusion matrix in Table S3. Despite a prevalence of PGQ of only 43% in the dataset, the model performs well in predicting both PGQ (sensitivity: 0.82) and non-PGQ (specificity: 0.87) at a probability cutoff of 0.50. The accuracy is correspondingly high at 0.85. Likewise, the model’s area under the receiver operating characteristic curve (AUC), which considers the full range of possible cutoffs, has a very high value of 0.88 on the validation set (Fig. S4). Additionally, the model achieves substantial agreement beyond chance in its classifications with Cohen’s kappa coefficient of 0.69. While slightly below the ideal, this value falls within 0.61–0.80, which Landis and Koch defined as substantial agreement25. This slightly lower kappa value reflects the inherent complexity and variability of the environmental data and processes analyzed in this study. The concordance between our model’s predictions and the PGQ maps presented in “2021 Annual Report of China’s Groundwater Monitoring Project”26, attest to the model’s accuracy and reliability (Fig. S5).

China’s groundwater quality during 1980–2020
Utilizing the developed model, we drew annual maps of the probability of PGQ across China from 1980 to 2020, with a detailed resolution of 1 km. In 1980, regions with PGQ probability exceeding 0.5 were primarily located in Southwest China, Northwest China, and parts of Northeast China, showing diverse levels of PGQ probability (Fig. 1a). As the 21st century unfolded, the situation evolved, with a substantial increase in PGQ probability observed in North China and Central China (Fig. 1b). These regions, which previously had relatively good groundwater quality, experienced a notable degradation by 2020 (Fig. 1c). The probability of being classified as having PGQ exceeds 0.9 in up to 2.1% of the national area, particularly across most regions of North China, signaling a notable deterioration in groundwater quality.

In 1980, only about 17.3% of the national area was affected by PGQ (Fig. 2a). By 2000, this percentage had slightly increased to 22.2%. However, there was a substantial rise to a peak of 41.6% in 2019, followed by a slight decrease to 40.8% in 2020.

Figure 3a illustrates the spatial distribution of groundwater quality category changes in China from 1980 to 2020, detailing regions of deterioration (non-PGQ changed to PGQ), improvement (PGQ changed to non-PGQ), and stability (unchanged). Across the country, 25.3% of the national area experienced deterioration in groundwater quality, 1.8% showed improvement, and 72.9% remained unchanged. Figs. S6a and S6b show that the area with deteriorated groundwater quality expanded from 6.8% during 1980–2000 to 20.9% during 2000–2020, suggesting an accelerated deterioration over recent decades. Between 1980 and 2000, regions with deteriorated groundwater quality were primarily scattered across Central China, North China, and coastal regions (Fig. S6a). In the last two decades, this situation spreads to encompass more regions, including North China, Northeast China, Northwest China, and Southeast China (Fig. S6b).

Compared to the expansion of PGQ area, the populations affected have undergone a more pronounced change. In the last century, most of the regions with PGQ were not densely populated, with the majority having a population density less than 100 people per km² (Fig. 1d). However, as we entered the new century, there has been a noticeable increase in the population density within the newly emerged PGQ regions (Fig. 1e). Approximately one-third of these PGQ regions now have a population density exceeding 100 people per km² (Fig. 1f). From 1980 to 2020, the area ratio suffering from PGQ has doubled, while the proportion of affected population experienced a fourfold escalation. This demographic change highlights an increase from affecting 6.8% of the population in 1980, to 17.5% in 2000, peaking at 38.3% in 2018 before slightly decreasing to 36.0% in 2020 (Fig. 2b).

By analyzing drivers of groundwater quality changes from 1980 to 2020 (“Analysis of drivers” in Methods), we have obtained the distribution of dominant factors (Fig. S7a) and their relative contributions (Fig. S7b). This analysis showed a heterogeneous geographical pattern, illustrating the spatial heterogeneity of the driver impacts.

Pollution-dominant regions encompassed the largest area, covering 16.2% of the national area (Fig. 3b). Among the pollution sources, agricultural discharge emerged as the predominant contributor, affecting 9.6% of the national area, mainly located in North China, Northeast China, and Central China. Industrial discharge as the dominant driver was concentrated primarily in Northeast China and Central China, affecting 4.4% of national area. In contrast, regions predominantly influenced by domestic discharge were concentrated in Southeast China, comprising 2.2% of the national area. Additionally, regions where groundwater exploitation was the main driver constituted 4.6% of national area, notably impacting Northeast China and Northwest China. The distribution of regions dominated by land use was fragmented, with the smallest area ratio of 0.2%.

During 1980–2020, the drivers of groundwater quality change in China have evolved spatially and temporally. Notably, in regions experiencing groundwater deterioration, agricultural discharge and climate change emerged as the primary drivers from 1980 to 2000, affecting 2.7% and 1.8% of the national area, respectively (Fig. S6c). The period from 2000 to 2020 saw a change in dominant drivers, with agricultural discharge and groundwater exploitation taking precedence, controlling 7.6% and 5.7% of the national area, respectively (Fig. S6d). This period witnessed a marked escalation in the cumulative impact of various environmental factors on the change of PGQ area ratio. Specifically, agricultural discharge experienced a dramatic increase in its contribution to the PGQ area ratio, climbing from 0.04% in 1985 to 10.7% in 2020, while industrial discharge rose from 0.09% to 5.3% over the same period (Fig. 3c). Similarly, groundwater exploitation marked a notable increase, escalating from 0.07% to 5.6%. Domestic discharge, though initially negligible, also notably impacted PGQ area ratio, contributing an increase of 1.7% by 2020 from 0.1% in 1985. Concurrently, climate change demonstrated a fluctuating yet persistent moderate influence, starting from a 0.1% contribution in 1985, peaking at 0.9% in 2000, and stabilizing at 0.5% in 2020. In contrast, land use changes had a consistently minimal and predominantly positive impact on groundwater quality dynamics.

Post-2020 projections of groundwater quality changes
We analyzed the temporal evolution of groundwater quality from 2021 to 2100 under four distinct scenarios, enabling an assessment of long-term influence of these key drivers on groundwater quality (Fig. 4). In Scenario I, we maintained the intensity of all influencing factors at their 2020 levels. To further investigate the responses of groundwater quality to future socio-economic development and climate change, we used the Shared Socio-economic Pathways (SSP) scenarios, i.e., SSP1-1.9 (Scenario II), SSP2-4.5 (Scenario III), and SSP5-8.5 (Scenario IV), from the Sixth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC AR6)27 (“Scenario settings” in Methods). Specifically, Fig. 4a-d shows the PGQ probability in 2050 under four scenarios (I, II, III, and IV), while Fig. 4e-h depicts the affected populations, respectively. Additionally, Fig. S8 explores changes in PGQ probability up to 2050, and Fig. S9 displays the changes in groundwater quality categories between 2020 and 2050 under four scenarios, whereas Fig. S10 extends to 2100, detailing the distribution of PGQ probability as well.

Scenario I projects a near-constant rate of deterioration, with the PGQ area ratio projected to rise from 40.8% in 2020 to 42.9% in 2050, and slightly increase further to 43.1% in 2100 (affecting ~35.3% of the population) (Fig. 2). Key regions of concern include Southeast China, Southwest China, and Central China by 2050 (Fig. S9a), particularly in South China where deterioration is expected to expand by 2100 (Fig. S10i).

In contrast, Scenario II projects substantial improvements, driven by proactive environmental policies and reduced emissions. This optimistic projection is evidenced by the PGQ area ratio’s decline to 37.9% in 2050, with expectations for a further reduction to 33.2% in 2100, illustrating a sustained positive trend (Fig. 2). The affected population is expected to decrease from 31.1% in 2050 to 21.1% in 2100. This improvement is predominantly observed in Northeast China and North China in 2050 (Fig. S9b), with additional progress anticipated in South China in 2100 (Fig. S10j).

Scenario III presents a concerning outlook, revealing a troubling trend with the PGQ area ratio rising from 47.2% in 2050 to 49.1% in 2100 (Fig. 2). Approximately 39.1% of the population could be exposed to PGQ in 2050, further increasing to 40.5% in 2100. This scenario reflects a situation of medium pollution and persistent emissions, leading to widespread deterioration, particularly in Central China and Southwest China in 2050 (Fig. S9c). Moreover, the projection indicates an expansion of the affected regions, extending groundwater quality deterioration to Northwest China and South China in 2100 (Fig. S10k).

Scenario IV, which represents the highest emission trajectory, projects an even more alarming increase in PGQ area ratio (48.3% in 2050 and 50.8% in 2100) (Fig. 2). Concurrently, approximately 40.7% of the population could be exposed to PGQ in 2050, with this proportion slightly decreasing to about 39.5% in 2100. The spatial distribution of PGQ expands notably, observed in the more industrialized and populated regions like Northeast and Central China in 2050 (Fig. S9d), and later spreads to additional regions in Northwest China and South China in 2100 (Fig. S10l).

Across all four scenarios, relative to the spatial distribution pattern observed in 2020, the probability of encountering PGQ in northern regions, especially in North China and Northeast China, is expected to slightly decrease by 2050 (Fig. S8) and 2100 (Figs. S10e-h), though the overall groundwater quality is projected to remain poor (Figs. S9 and S10i-l).

Our analysis reveals a complex interplay of factors (i.e., climate change, groundwater exploitation, land use and pollution discharge) driving changes in China’s groundwater quality during 2020–2100 (“Analysis of drivers” in Methods). The spatial distribution of dominant factors affecting PGQ probability in China during 2020-2050 is detailed in Fig. S11, specifically focusing on areas experiencing deterioration (Fig. 5a, b, d, e) or improvement (Fig. 5c). Temporal trends of these drivers across China are presented in Fig. S12.

Pollution discharge plays a critical role, especially in regions like Central China and Southwest China, where industrial, agricultural, and domestic activities are prevalent, with Scenario IV controlling the largest area, with 3.0% of the national area (Fig. 5e). In comparison, pollution affects 2.2% of the national area under Scenario III (Figs. 5d), and 2.6% in Scenario I (Fig. 5a). During 2020–2050, pollution is the leading contributor to the increase in PGQ area ratio, particularly noted in Scenario IV with an increase of 4.0% (Fig. 5f). This trend not only persists but also intensifies by 2100, with pollution contributing an additional increase of 4.2% in PGQ area ratio relative to 2020, reflecting the escalating effects of pollution under high-emission scenarios.

Land use, whose changes are driven by agricultural expansion, urbanization, and other human activities, predominantly impacts groundwater quality in Central China and Eastern China, particularly in regions such as the North China Plain, the middle and lower reaches of the Yangtze River, and the Northeast Plain, where intensive agricultural and urban development are prevalent. Under Scenario IV, land use dominates groundwater quality change in 5.6% of the national area, while Scenario III sees a controlling area ratio of 5.5% (Fig. 5d, e). Figure 5f reveals that land use is one of the key contributors to changes in PGQ area ratio, with Scenario IV leading to a 1.5% increase, and Scenario III showing a 1.2% increase by 2100. Notably, Scenario II stands out with a substantial decrease in land-use-related impacts, reducing the PGQ area ratio by 0.9% in 2050 and further to 1.8% in 2100, indicating a potential positive effect of controlled changes in land use on groundwater quality in certain scenarios.

Groundwater exploitation notably affects regions with high water demand, such as Northeast China, Northwest China, and North China. Under all scenarios, it contributes to groundwater quality deterioration in up to 0.3% of the national area. However, Scenario II stands out for its positive effects, improving groundwater quality across ~4.0% of the national area (Fig. 5c). From 2020 to 2050, reductions in groundwater exploitation, particularly under Scenario II, result in a slight improvement in PGQ area ratio by 0.2% (Fig. 5f). In 2100, continued efforts to reduce groundwater exploitation under Scenario II are projected to further reduce the PGQ area ratio by up to 0.6%, underscoring the importance of sustainable groundwater management.

Climate change exerts a moderate yet notable influence on groundwater quality across all scenarios, especially in Northwest China, Central China, and Southeast China. In Scenario I, climate change affects 1.2% of the national area (Fig. 5a), increasing to 3.1% in Scenario IV (Fig. 5e). The impact is similarly substantial in Scenario III, where it influences 3.1% of the national area. In 2050, climate change has a modest impact on the PGQ area ratio, subtly altering it by less than 1.0% across various scenarios (Fig. 5f). However, in 2100, its effect becomes slightly more pronounced, potentially increasing the PGQ area ratio by up to 1.8% in Scenario IV, reflecting a steady, though gradual, worsening of groundwater quality driven by this factor.

Methods
Data collection and pre-processing
The groundwater quality data were obtained from published articles by searching in China National Knowledge Infrastructure (CNKI) and Web of Science. Initially, a total of 31,535 studies comprising 25,200 articles in Chinese and 6335 articles in English were retrieved using the following search formula: “Groundwater quality” (Topic) and “China” (Topic). Given the heterogeneity of the data, the retrieved publications were further filtered according to the following rules: (i) at least five groundwater quality indicators were determined; (ii) sampling locations were provided at prefecture level; (iii) there were more than five sampling points and sufficient statistical information (including mean at least); (iv) land-use type was provided; (v) samples were taken between 1990 and 2020. To effectively capture the contemporary groundwater quality status, data predating 1990 were deliberately omitted due to its limited representativeness and inadequacy in reflecting the prevailing conditions. Finally, 753 articles were identified to establish our dataset, which contains a total of 1977 surveys, as shown in Fig. S1.

The groundwater quality data were classified into two categories: non-PGQ and PGQ, based on the criteria outlined in the Groundwater Quality Standards of the People’s Republic of China (GB/T 14848-2017)22. This classification process involves the following four steps:
Step 1: Collect all relevant groundwater quality indicators from the dataset, ensuring they contain more than five indicators, as specified by established testing standards.
Step 2: Assess each indicator individually by classifying it as Class V or not, based on Table S1.
Step 3: Determine the overall groundwater quality with the worst indicator.
Step 4: Categorize the groundwater quality into two types: non-PGQ (Class I to IV) and PGQ (Class V).

In our analysis, we meticulously built the relationship between predictors and groundwater quality category (PGQ or non-PGQ). As listed in Table S2, the predictors included two classes of steady-state factors (soil properties and geographic and hydrogeologic characteristics) and four classes of time-variant factors (i.e., climate change, pollution discharge, groundwater exploitation, and land-use change).

Groundwater quality is influenced by various environmental and chemical processes, including redox reactions, contaminant transport, and natural filtration mechanisms51. Given China’s vast and diverse landscape, along with notable regional variations in data availability and quality, we prioritized accessible and quantifiable predictors to ensure broader applicability and consistency in groundwater quality predictions across different regions. Accordingly, we selected predictors that not only capture redox-related influences indirectly but also reflect other critical mechanisms affecting groundwater quality, providing a comprehensive foundation for reliable assessment. To effectively capture these complex processes, we selected predictors such as soil properties (e.g., porosity, permeability, soil type, and soil composition across different soil layers) to reflect contaminant adsorption and biological degradation, geographic and hydrogeologic characteristics (e.g., topographic wetness index, water-table ratio, and groundwater type) to indicate contaminant accumulation and natural dilution, and climatic factors (e.g., temperature, precipitation, and aridity) to represent groundwater recharge, evaporation, and pollutant mobility17,19,30. Additionally, pollution discharge (agricultural, industrial, and domestic sources) can substantially affect groundwater quality by increasing contaminant loads, while land-use change, can disrupt the natural infiltration and recharge processes52. Notably, groundwater exploitation (e.g., groundwater supply and groundwater exploitation rate) emerges as a crucial anthropogenic perturbation factor53,54. Excessive groundwater extraction lowers water tables, concentrates contaminants, and, in coastal areas, can trigger seawater intrusion, increasing salinity35,55,56. It can also mobilize naturally occurring contaminants like arsenic and fluoride, further degrading water quality30.

To classify groundwater quality at the prefecture level, we implemented a series of data preprocessing steps. First, regional data were converted into 1 km² grid cells based on regional land-use types, with all cells within a prefecture assigned the corresponding groundwater quality data to ensure consistency. The groundwater quality data were then binarized, assigning a value of zero to non-PGQ and one to PGQ. The purpose of this approach was twofold: (i) to prioritize the fundamental health aspect of safety (non-PGQ) or unsafety (PGQ) of groundwater for drinking; and (ii) to address variations in precision resulting from diverse analysis methods employed in different data sources. Next, natural factors in the predictors were extracted at a 1 km² resolution, assigning each grid cell a unique value representing local environmental conditions. Socio-economic factors were uniformly applied to all grid cells within a province, based on the sampling year. The dataset on land-use types covers specific benchmark years: 1980, 1990, 1995, 2000, 2005, 2010, 2015, 2020, and every five years between 2020 and 2100. For other years, land-use data are assumed to be identical to the adjacent prior benchmark year. For instance, land-use data in 1989 are taken to be the same as those in 1980. For other dynamic factors, their cumulative average values from 1980 to the target year were calculated as the model inputs. This accounts for the lag in pollutant migration and effectively captures the long-term impacts on groundwater quality. Fertilizer consumption and pesticide consumption were averagely allocated per km2 of agricultural land, with non-agricultural land assigned zero. Industrial wastewater discharge, industrial solid waste discharge, domestic wastewater discharge, domestic solid waste discharge, and groundwater supply were allocated per km2 based on population density for each province. Per capita data for these variables were calculated by dividing the total value by the population of each province, then distributed per km2 based on population density. Groundwater exploitation rates were uniformly distributed across each province, with regions having a population density below 1 person/km² assigned a rate of zero. Other factors were allocated per km2 based on the resolution of raster data or vector data. Finally, groundwater quality categories for each year and region were matched with the corresponding natural and socio-economic variables to create the training and validation datasets for the RF model.

Model development and evaluation
In this study, the RF model was constructed to predict the spatio-temporal dynamic of groundwater quality from natural and social-economic factors. The RF model was chosen for its ability to handle large datasets with numerous predictor variables while efficiently capturing complex interactions among them24. The RF model uses the bootstrap resampling method to build multiple decision trees. For the construction of each tree, samples are independently selected; however, the distributions for all trees in the forest are the same, which guarantees the robustness of the model. In addition, the RF model offers key advantages such as preventing overfitting by introducing randomness at each decision node, which enhances generalization. By selecting random feature subsets for each tree, RF effectively addresses the challenges of high dimensionality in the data. It also reduces experimental noise and improves prediction accuracy through an ensemble approach, averaging predictions from multiple decision trees. Furthermore, RF provides inherent feature importance measures, helping identify key variables, and handles missing data efficiently, either through imputation or splitting based on available features. With minimal parameter tuning required compared to other machine learning algorithms, RF presents itself as an ideal model for this study14,57,58.

We began by optimizing the hyperparameters of the RF model, with a particular focus on the number of decision trees. By systematically testing tree counts in increments of 10, from 10 to 200, we identified 50 trees as the optimal balance between accuracy and computational efficiency (Fig. S21). Other parameters, such as minimum leaf size, were kept at their default settings, as these have proven effective in similar applications.

To further validate the model’s performance and safeguard against overfitting, we conducted 10-fold cross-validation. The dataset was randomly partitioned into 10 subsets, with 9 used for training and 1 for validation in each iteration. This process was repeated across all subsets, ensuring comprehensive evaluations. The averaged performance metrics (Table S4 and Fig. S22) confirmed the model’s ability to generalize effectively, without overfitting.

After confirming the model’s robustness through cross-validation, we proceeded to retrain it using 90% of the data (1779 surveys, comprising 67,315 points) for training and 10% (198 surveys, comprising 15,335 points) for validation, each preserving the proportion of PGQ and non-PGQ of the full dataset. This approach maximizes the model’s learning capacity by providing a larger training set while preserving an independent validation set to assess performance on unseen data. This final validation confirmed that the model maintained strong predictive accuracy, further enhancing its generalization ability.

To precisely and comprehensively establish the relationship between model predictors and groundwater quality, 51 predictor variables were initially included. We conducted Pearson correlation analysis to reduce redundancy in the predictor variables. Initially, all original variables were included in the RF model training. Among the list of candidate variables, one is opted out if its Pearson correlation coefficient23 with any other variables was larger than 0.7. However, certain variables were retained despite exceeding this correlation threshold. This decision was made given that excluding these variables notably reduced the model’s predictive performance, demonstrating their essential contributions to accurately capturing the system’s dynamics. This step-by-step approach ensures that the final set of variables used in modeling or further analysis is less likely to include redundant information, thus improving the parsimony and interpretability of the model. Finally, 25 predictors were reserved to predict the groundwater quality in China. The data sources and detailed information about these predictors were described in the Supplementary Information.

The final RF model was trained with 50 trees using the selected 25 predictors as inputs to classify each 1 km2 grid as being PGQ or not. Each decision tree in the forest casts a ‘vote’ for a class label based on the input data, with the final class label determined by the majority vote across all trees. The model also calculates class-label probabilities based on the proportion of trees voting for each class. The performance of the RF model is systematically validated with the following metrics.

Prevalence (Prev) quantifies the proportion of actual positives in the dataset and is calculated as follows59:
Prev = (TP+FN) / (TP+FN+TN+FP) (1)

where TP (true positives) are samples correctly identified as PGQ; FN (false negatives) are samples incorrectly identified as non-PGQ; TN (true negatives) are samples correctly identified as non-PGQ; FP (false positives) are samples incorrectly identified as PGQ.

General statistical indices, such as accuracy (Acc), sensitivity (Sen), precision (Prec) and specificity (Spec), were employed to quantify the bias of the RF model, with higher values indicating better model effectiveness60. The definitions of TP, TN, FP, and FN shed light on the seemingly subtle distinctions among these statistics, which play a crucial role in assessing the model’s performance and its ability to minimize false positives and negatives59:

Acc = (TP+TN) / (TP+FN+TN+FP) (2)

Sen = TP / (TP+FN) (3)

Prec = TP / (TP+FP) (4)

Spec = TN / (TN+FP) (5)

The performance of the trained RF model was then assessed using the area under the receiver operating characteristic (ROC) curve (AUC). The ROC curve is primarily used to evaluate the performance of a classification model by showing the trade-off between the sensitivity and specificity across different thresholds61. AUC is a single measure that summarizes the overall performance of RF for classification. The higher the AUC, the better the model is at distinguishing between PGQ and non-PGQ. A perfect classifier would have an AUC of 1, while a completely random one would have an AUC of 0.5.

Additionally, Cohen’s kappa coefficient (K), which measures the agreement between the observed and predicted classifications corrected for chance, is used to further assess the model’s accuracy59. Kappa coefficient of 0.6–0.8 is generally considered good, and values above 0.8 are deemed excellent, calculated by:

K = (Po - Pe) / (1 - Pe) (6)

where Po is the observed agreement (Acc), and Pe is the expected agreement by chance:

Pe = ((TP+FP) * (TP+FN) + (FN+TN) * (FP+TN)) / (TP+FN+TN+FP)^2 (7)

Apart from evaluating the predictive abilities of the RF model, the spatial patterns of the estimates were validated by comparing them with nationally recognized maps, ensuring their consistency and reliability.

To assess predictor importance, we utilized the error of out-of-bag (OOB) samples among trees, a common metric employed in RF models. This approach involves randomly permuting the values of a specific predictor variable and calculating the resulting change in OOB error. OOB error refers to the prediction error computed on samples that were not included in the bootstrap sample used to train each decision tree in the RF model. If permuting a predictor leads to a considerable increase in OOB error, it suggests that the variable plays a critical role in prediction accuracy. This approach offers a reliable means of identifying the important features for model prediction.

To further assess the sensitivity of each influencing factor on the model’s predictions62, we conducted a sensitivity analysis by systematically increasing or decreasing each continuous predictor by 10%. This approach allowed us to observe how these perturbations influenced the predicted average probability of PGQ over the historical period from 1980 to 2020. By comparing these modified predictions to the baseline scenario, where no changes to the predictors were applied, we were able to quantify the degree to which each factor contributed to variability in PGQ probability, providing insights into which predictors most notably affect groundwater quality predictions under different conditions.

Following the identification of key predictors, partial dependence plots (PDPs) were generated to visualize each predictor’s marginal effect on the probability of PGQ. The PDPs were generated by holding the focal predictor constant across its empirical range while integrating out the effects of all other predictors, thereby isolating the unique contribution of the predictor to the response63. This was achieved by averaging the model predictions over the distribution of the data, providing a clear visual representation of the dependency between the predictors and the probability of PGQ.

In addition, feature interaction charts (FICs) were constructed to explore the combined effect of two predictors on the predicted probability of PGQ63. It takes into account the range of values for the two concerned predictors. However, instead of holding one variable constant, FICs vary both variables simultaneously over their empirical ranges. The model’s predictions are then averaged over the joint distribution of the two predictors, and the combined effect is visualized in a two-dimensional color-coded map. This approach highlights how the interaction between two predictors influences the model output, thus providing a deeper understanding of the relationships within the data.

PGQ prediction
The RF model was used to reconstruct and predict the spatiotemporal distribution of groundwater quality from 1980 to 2100, enabling us to infer long-term trends in populations potentially exposed to PGQ. We apply a probability cutoff (Prob_cut) of 0.5 to identify the populations at risk in areas exceeding this threshold. Then the affected population living in each grid (1 km2) was calculated by multiplying the grid’s population (Pop) by its probability of PGQ (Prob_poor). The calculation of the potentially affected population (Pop_affect) is summarized in the following equations:

Pop_affect = Pop * Prob_poor, Prob_poor > Prob_cut (8)
Pop_affect = 0, Prob_poor <= Prob_cut

The dynamics of PGQ probability among different periods were controlled by the four categories of time-variant drivers (i.e., climate change, pollution discharge, groundwater exploitation, and land-use change). To quantify the contribution of each driver on PGQ, we designed factorial experiments between 1980 and 2020. In the experiment, we just kept one specific driver category (e.g., climate) constant at their initial states in 1980, while allowing the other driver categories vary with time as they are (i.e., the time series in the real world). For instance, the Prob_Climate means that all drivers in the climate category remained constant at values in 1980 but with real (time-variant) features of drivers in pollution and groundwater categories.

By comparing the simulations between experiments, we could isolate the impact of a target driver category and quantify its contribution. To implement the comparison, we defined the differences between simulation driven by all variables in period y (Prob_All,y) and that driven by partial time-variant driver categories (Prob_Vi,y) as the actual contribution of the target variable i (Vi) and period y (Con_Vi,y):

Con_Vi,y = Prob_All,y - Prob_Vi,y (9)

The relative contribution (%) of variable i (i=1, 2, 3,…, n) and period y (y=1, 2, 3,…, k) (RelCon_Vi,y) was defined as:

RelCon_Vi,y = |Con_Vi,y| / sum(|Con_Vi,y|) * 100% (10)

To quantify the contribution of each driver category to changes in PGQ area ratio on the national scale, we define the differences between the area ratio of PGQ driven by all variables in period y (AR_All,y) and that driven by partial time-variant driver categories (AR_Vi,y) as the primitive contribution of the target variable i (Vi) and period y (ConAR_Vi,y):

ConAR_Vi,y = AR_All,y - AR_Vi,y (11)

The differences between the area ratio of PGQ driven by all variables in period y (AR_All,y) and that driven by all variables in the baseline year (AR_All,BY) are defined as the actual growth:

ActAR_All,y = AR_All,y - AR_All,BY (12)

Then, the relative contributions of each driver category to changes in the PGQ area ratio are calculated as:

RelAR_Vi,y = (ConAR_Vi,y / sum(ConAR_Vi,y)) * ActAR_All,y (13)

The grid annual mean relative contributions of the four driver categories in groundwater were used to visualize the spatial patterns of the influential intensity of different driver categories on a national scale. We defined the single dominant driver category of each pixel as the variable with a maximum relative contribution. As pollution discharges are composed of agricultural, industrial, and domestic pollution, the dominant driver category of pollution discharge for each pixel is defined as the secondary driver category with the largest relative contributions among the three secondary driver categories.

To further elucidate the temporal dynamics of groundwater quality, the differences in the projected groundwater quality categories for the years 2020, 2000, and 1980 were calculated, generating maps of groundwater quality changes for the periods 1980–2000, 2000–2020, and 1980–2020. Regions of groundwater quality deterioration (non-PGQ to PGQ) during these periods were extracted as masks. This approach yielded detailed maps of the driving factors behind deteriorated groundwater quality regions across the three distinct periods.

In future scenarios, the analysis of dominant driver focuses exclusively on four primary categories: climate change, groundwater exploitation, land use and pollution discharges, without individually detailing secondary pollution categories. We kept one specific driver category (e.g., climate) constant at its states in 2020, while allowing the other driver categories to vary with time.

We evaluated the temporal evolution of groundwater quality from 2021 to 2100 in four different scenarios, which allowed us to assess the long-term changes in groundwater quality by considering the continuity of the underlying factors over the specified time frame. To evaluate the impacts of climate change on groundwater quality, we utilized the SSP scenarios, namely SSP1-1.9, SSP2-4.5, and SSP5-8.5, as outlined in the IPCC AR664. These scenarios provide a framework for understanding and projecting different socio-economic and climate conditions, allowing us to examine the potential responses of groundwater quality to various development and climate change trajectories.

In Scenario I (Baseline), we assumed a consistent development intensity for all drivers. Climate, groundwater exploitation, and pollution discharges were maintained at the same levels as the cumulative average values in 2020. Similarly, land-use type in Scenario I was also assumed to remain unchanged from 2020. In contrast, Scenarios II-IV incorporated projections to account for future variability. Specifically, climate data from one CMIP6 model—EC-Earth3 were used as model input to project future PGQ probabilities65–67. Additionally, land-use data in these scenarios were derived from the SSP-RCP global 1km land-use simulation dataset (2020-2100) for SSP1-1.9, SSP2-4.5, and SSP5-8.5 scenarios, updated every five years68. To maintain consistency in representation, these land-use data were aligned with those used for the historical period from 1980 to 2020.

In the three scenarios, annual change rates of −3%, 1%, and 3% were assigned as gradients for groundwater exploitation and pollution discharge, respectively. These rates were determined based on historical trends69,70 and future projections71,72 related to industrial activities, agricultural expansion, and environmental policies. The low rate of −3% reflects a proactive approach toward reducing resource use and pollution discharge, aligning with the sustainability-focused SSP1-1.9 scenario73. In contrast, the medium rate of 1% corresponds to a moderate development pathway, as seen in SSP2-4.5. This scenario embodies a balanced approach, allowing socio-economic development to progress while placing emphasis on environmental protection74. The high rate of 3% aligns with the SSP5-8.5 scenario, characterized by rapid economic growth and minimal environmental regulation75. The scenarios thus represent a spectrum of future development and industrial expansion. We acknowledge that future development may exhibit more intricate patterns, and our current study may not fully capture this complexity and uncertainty. Then, these rates were applied to the corresponding values from 2021 to 2100, and the cumulative average was calculated as the model inputs for each year. Please refer to Table S5 for detailed information on the specific scenario assumptions.

To evaluate the impact of uncertainties in climate change, groundwater exploitation, and pollution charge on model predictions, we applied Monte Carlo (MC) simulations to quantify the variability in groundwater quality predictions due to input uncertainties76. This analysis was conducted for the period 1980–2100, with data processed at 5-year intervals for 1980–2015, annually from 2016 to 2020, and at 5-year intervals for future projections from 2025 to 2100 under four scenarios.

In each scenario-year combination, MC simulations with 50 members were run to introduce variability in key natural and socio-economic predictors. Perturbations were applied to these dominant factors, such as climate, groundwater exploitation and pollution charge, reflecting their real-world uncertainties. For instance, temperature varied within a ±2°C range, while predictors like precipitation, fertilizer, and pesticide use were perturbed by up to ±20%.

The results of these simulations can be used to evaluate the range of possible groundwater quality outcomes based on uncertainties in input data. For each year and scenario, we calculated the mean and standard deviation of the predicted probabilities based on the ensemble outcomes, providing a robust estimate of both expectation and variability in groundwater quality predictions.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper addresses a critical environmental issue, exploring the groundwater quality crisis across China using extensive data and Random Forest (RF) modeling technique. The study is well-motivated and holds significant implications for environmental management and policy-making. The manuscript is generally well-structured, with a clear introduction, and comprehensive results and discussion sections. However, there are several areas where the manuscript could be strengthened to enhance its clarity, robustness, and impact. Specifically, the methodology section needs more detail on the hyperparameter tuning process, data preprocessing steps, and the rationale for choosing the RF model. Additionally, the results section would benefit from a more detailed interpretation of the key predictors and a thorough uncertainty assessment in the scenario analysis. Finally, the discussion should provide a more nuanced analysis of the interactions between various drivers of groundwater quality and offer specific policy recommendations based on the study's findings. Addressing these areas will significantly improve the manuscript's scientific rigor and practical relevance.\nMajor concerns:\n1) The current methodology provides limited details on the hyperparameter tuning process for the Random Forest (RF) model. Specifically, it lacks a thorough explanation of how the number of trees, maximum depth, and other critical parameters were optimized. Elaborate on the hyperparameter tuning process by describing the specific methods used, such as grid search or randomized search. Include a discussion on the range of values tested for each parameter, the criteria for selecting the optimal values, and any cross-validation techniques employed to ensure robust model performance. For instance, detailing the grid search process might include explaining the range of the number of trees tested (e.g., 50 to 500 in increments of 50), the range of maximum depth (e.g., 5 to 30), and the evaluation metric used to select the best combination (e.g., highest AUC score on validation sets). Additionally, explain how overfitting was addressed, such as by using cross-validation techniques and evaluating the model on a separate validation set to ensure it generalizes well to unseen data.\nDiscuss the rationale for selecting the RF model over other machine learning algorithms, emphasizing its strengths in handling large datasets with numerous predictor variables and its ability to capture complex interactions between predictors without requiring extensive parameter tuning. Moreover, provide details on any data preprocessing steps undertaken, such as normalization or standardization of the predictor variables, to ensure that they are on a comparable scale and to enhance the model's performance. Explain the rationale for these preprocessing steps and how they contribute to the overall robustness of the model. Providing a comprehensive explanation of the hyperparameter tuning process, including how overfitting was checked, the rationale for selecting the RF model, and the data preprocessing steps, will enhance the transparency and reproducibility of the study. It will demonstrate the rigor of the model development process and increase confidence in the robustness and reliability of the model predictions, ultimately strengthening the scientific contribution of the research.\n2) While the methodology mentions using out-of-bag (OOB) error to assess the importance of predictor variables, it lacks a detailed analysis of how these predictors contribute to the model and their implications for groundwater quality. Conduct a more detailed analysis of the top predictors identified by the RF model. For each top predictor, discuss its relative importance and provide insights into how it influences groundwater quality. For instance, if soil type and industrial wastewater discharge are among the top predictors, explain their potential mechanisms of influence on groundwater quality. Additionally, include visual aids such as partial dependence plots to illustrate the relationship between these predictors and the predicted groundwater quality outcomes. A detailed analysis of predictor variable importance will enhance the interpretability of the model and provide valuable insights into the factors driving groundwater quality, thereby strengthening the scientific contribution of the study and informing policymakers and stakeholders about the critical variables that need to be managed to improve groundwater quality.\n3) The methodology describes the use of factorial experiments to isolate the impacts of specific driver categories under different scenarios, but the implementation details of these scenarios and the handling of uncertainties in future projections are insufficiently described. Expand on the scenario analysis methodology by providing a detailed description of how each scenario is constructed, including the data sources, assumptions, and adjustments made for future projections. Incorporate a discussion on how uncertainties in climate, exploitation, and pollution data are addressed, such as using sensitivity analysis or Monte Carlo simulations to quantify the impact of these uncertainties on model predictions. Present the results of these uncertainty assessments to highlight the range of possible outcomes and their implications. A comprehensive and transparent scenario analysis, coupled with a thorough uncertainty assessment, will significantly strengthen the study's credibility by providing a clearer understanding of the potential future states of groundwater quality under various conditions. This will aid in the development of robust and informed management strategies, enhancing the study's value to policymakers by offering a nuanced understanding of the risks and uncertainties involved.\n4) The discussion section effectively identifies major drivers of groundwater quality deterioration, including pollution discharge and groundwater exploitation. However, it should provide a more nuanced analysis of the interactions between these drivers and climatic factors. Addressing the limitations of the study, such as the assumption of constant land-use types in future scenarios and the exclusion of certain predictor variables due to high correlation, is crucial. Moreover, the policy implications of the findings should be explored in greater depth, offering specific recommendations for mitigation strategies to address the identified drivers of groundwater quality degradation. Highlighting the socio-economic implications and feasibility of these strategies in different regions of China will enhance the manuscript's practical relevance for policymakers and stakeholders involved in groundwater management. By integrating these improvements, the results and discussion sections will present a more comprehensive and impactful analysis, strengthening the scientific rigor and applicability.", "reviewer_2": "The study of groundwater quality evolution is of significant importance and relevance, given its pivotal role in sustaining freshwater resources and maintaining ecosystem health. The application of machine learning, particularly the random forest model, to analyze and predict the spatiotemporal dynamics of groundwater quality across entire China is innovative and represents a valuable contribution to the field. The manuscript is well organized and easy to follow. However, there are several areas that require clarification and additional detail to enhance the robustness and transparency of this study. Below are some comments and suggestions aimed at improving the manuscript:\n1. Definition of Poor Groundwater Quality (PGQ): The study does not clearly distinguish between PGQ caused by human activities and that resulting from geogenic origins, two distinct categories influenced by different driving forces. Moreover, different groundwater pollutants have different sources and influencing factors. For example, nitrogen pollution from non-point sources is significantly affected by surface input and migration through the unsaturated zone, yet the machine learning model does not account for the impact of the unsaturated zone. Similarly, heavy metals are largely influenced by redox conditions, especially in cases where naturally poor-quality groundwater is more affected by hydrogeological conditions than by surface processes. Consequently, the machine learning model may overemphasize the role of surface soil, potentially overlooking other critical factors.\n2. Hydrogeological Factors: While this study developed a machine learning model to construct annual probabilistic maps of poor groundwater quality across China using 25 natural and social-economic factors, few of these are explicitly related to hydrogeological processes. Groundwater flow and solute transport are clearly significant for predicting groundwater quality. As an example, the authors mention a factor called 'Aquifer-media-based groundwater type' in the supporting information, but its establishment is not clearly explained. Does this refer to distinctions such as confined vs. unconfined, porous media vs fractured media, or other characteristics?\n3. Chemical Processes: As noted in comment #1, chemical processes is critical in determining groundwater quality. For example, nitrate (NO3), one of the indicators considered in this study for groundwater quality, is redox-sensitive and undergoes denitrification under reducing conditions. Nitrate concentrations in groundwater are therefore strongly dependent on redox conditions, regardless of nitrogen input at the land surface. How redox conditions are considered in assessing the groundwater quality is an area that warrants further exploration.\n4. Temporal Variations: This study aims to investigate the temporal changes in Chinese groundwater quality; however, only social-economic factors were considered temporally dynamic, while soil and environmental factors were treated steady-state. It’s noteworthy that environmental conditions vary over time, especially the depth to the water table, which in turn affects contaminant transport from the land surface to the underlying aquifer system. Considering these factors as steady state may underestimate their impacts, and this should be addressed in the paper.\n5. Data Set Clarity: The manuscript mentions the number of surveys derived from a certain number of articles but does not specify the exact number of data points within the groundwater quality dataset. For clarity, could you please provide the total number of data points and detail how they are divided into training and validation sets? This information is crucial for understanding the scale of the analysis and the potential for model overfitting.\n6. Spatial Partitioning Methodology (L.365-372): The method of spatial partitioning into 1KM*1KM grids requires further elaboration. Specifically, it is unclear how the presence of both PGQ and non-PGQ within a single grid is resolved. Additionally, please elucidate how point-scale driver variables are expanded to fit the 1KM*1KM grid format. A clear explanation of these methodologies will greatly aid in understanding the spatial analysis performed.\n7. Model Parameters and Tuning: The manuscript would benefit from a detailed description of the random forest model parameters, including the settings for each parameter and the methods used for tuning. If parameter optimization was conducted, please describe the process and the rationale behind the chosen methods and settings.\n8. Uncertainty Analysis: An uncertainty analysis of the model results would be beneficial. Providing an assessment of the model’s uncertainty will give readers a better understanding of the reliability of the predictions.\n9. Data Processing Details: It would be beneficial to include a section on data preprocessing. Specifically, readers would benefit from knowing whether the input driver variable data were standardized or normalized before being fed into the model. Standardization is a common practice in machine learning to ensure that each feature contributes proportionately to the final prediction.\n10. Rationale for future scenario settings (L.497-514): The setting of future scenarios is crucial for the interpretation of model results. Given that variables such as fertilizer and pesticide usage appear to have a significant impact on the model, as indicated by Fig. S3 and the results section, the rationale behind the high (3%), medium (1%), and low (-3%) annual change rates should be discussed. This discussion could include the basis for these specific rates and how they align with projected socio-economic and environmental trends.\nAddressing these points will significantly improve the manuscript and provide readers with a clearer understanding of your research methodology and findings.", "reviewer_3": "Based on a lot of collected data, the manuscript used the random forest model to construct the annual probabilistic map of poor groundwater quality across China from 1980 to 2100. Although the method is not new, the idea they use it to model the groundwater quality is good. The new insights from the study are helpful for the policy-makers to take reasonable measures to protect China’s groundwater. However, there are some issues that need further addressing:\n1. Title of the paper didn’t fully reflet the contents, I think it should changed.\n2. Climate includes temperature, precipitation, evaporation etc., they usually have different impacts to groundwater quality. For example, precipitation is usually beneficial for groundwater quality, while evaporation is detrimental to groundwater quality. What does climate mean in your paper?\n3. Land use changes have different types, for example, from forest to farmland, or from farmland to construction land, etc. What do you mean in your study?\n4. How the groundwater exploitation affects its quality, the links between these should be clarified.\n5. Whether it is agricultural or industrial pollution, there is usually a lag time for it to migrate to groundwater. The time can be several years or tens of years. How do you treat this in you study?\n6. This study should consider the influence of geological background on water quality.\n7. Is a kappa coefficient of 0.69 a bit low? Generally, we believe that a kappa coefficient higher than 0.8 is best.\n8. In table S1, ug/L should be μg/L.\n9. In the assessment of groundwater quality in lines 335-336, does the use of a single factor assessment method make the results too pessimistic? Please explain.\n10. Since there is a lot of missing data before 1990, why not start the study directly from 1990?\n11. The ground water Standard they listed in Table S1 is 2017, but the reference they cited is 1993.\n12. Line 337, when we evaluate the water quality, we usually think that the water above Class III is beyond the standard. What is the basis for choosing Class V water as the boundary in this manuscript?\n13. The spatial pattern of future land use is affected by both human activities and climate change. In your article, you assume that the land use pattern of different scenarios is in 2020. Is this assumption too arbitrary? Different scenarios of future land use are also available in the Qinghai-Tibet Plateau database. Please explain why you don't use them?\n14. The research conclusion should be summarized."}, "author_response": "This paper addresses a critical environmental issue, exploring the groundwater quality crisis across China\nusing extensive data and Random Forest (RF) modeling technique. The study is well-motivated and holds\nsignificant implications for environmental management and policy-making. The manuscript is generally\nwell-structured, with a clear introduction, and comprehensive results and discussion sections. However,\nthere are several areas where the manuscript could be strengthened to enhance its clarity, robustness, and\nimpact. Specifically, the methodology section needs more detail on the hyperparameter tuning process,\ndata preprocessing steps, and the rationale for choosing the RF model. Additionally, the results section\nwould benefit from a more detailed interpretation of the key predictors and a thorough uncertainty\nassessment in the scenario analysis. Finally, the discussion should provide a more nuanced analysis of\nthe interactions between various drivers of groundwater quality and offer specific policy\nrecommendations based on the study’s findings. Addressing these areas will significantly improve the\nmanuscript’s scientific rigor and practical relevance.\nMajor concerns:\n1) The current methodology provides limited details on the hyperparameter tuning process for the\nRandom Forest (RF) model. Specifically, it lacks a thorough explanation of how the number of\ntrees, maximum depth, and other critical parameters were optimized. Elaborate on the\nhyperparameter tuning process by describing the specific methods used, such as grid search or\nrandomized search. Include a discussion on the range of values tested for each parameter, the\ncriteria for selecting the optimal values, and any cross-validation techniques employed to ensure\nrobust model performance. For instance, detailing the grid search process might include explaining\nthe range of the number of trees tested (e.g., 50 to 500 in increments of 50), the range of maximum\ndepth (e.g., 5 to 30), and the evaluation metric used to select the best combination (e.g., highest\nAUC score on validation sets).\nResponse: Thank you for your advice. In our modeling approach with random forest (RF), we focused\nprimarily on optimizing the number of trees. We tested different values, increasing the number of trees\nin increments of 10, ranging from 10 to 200. Through this process, we identified 50 trees as the optimal\nnumber, delivering the highest model accuracy (Fig. S21). To ensure the evaluation’s stability and\nreliability, we utilized 10-fold cross-validation, which allowed us to evaluate performance across\ndifferent data subsets, providing a robust assessment and helping to mitigate the risk of overfitting.\nOther hyperparameters of the RF model, such as the minimum leaf size for classification trees (set to 1),\nwere left at their default settings as provided by the MATLAB Statistics and Machine Learning Toolbox.\nThese defaults are widely regarded as effective for many applications, and in some cases, the toolbox\nenforces fixed values to prevent alteration. Given this, we did not find it necessary to adjust these\nparameters further. Additionally, since we did not use a deep RF model, tree depth was not treated as a\ntunable hyperparameter.\nMoreover, performance of the RF model is systematically evaluated with various metrics, please refer to\nour response to your subsequent comment on “how overfitting was addressed”.\nTo enhance clarity, we have incorporated the following texts in “Splitting the training and validation sets”\nsection of “Methods”:\n“We began by optimizing the hyperparameters of the RF model, with a particular focus on the\nnumber of decision trees. By systematically testing tree counts in increments of 10, from 10 to\n200, we identified 50 trees as the optimal balance between accuracy and computational\nefficiency (Fig. S21). Other parameters, such as minimum leaf size, were kept at their default\nsettings, as these have proven effective in similar applications.”\n0.79\n0.78\n0.77\n0.76\n0.75\n0.74\n0.73\n0 20 40 60 80 100 120 140 160 180 200\nFig. S21| Validation accuracy (Acc) of the RF model with different numbers of decision trees.”\nAdditionally, explain how overfitting was addressed, such as by using cross-validation techniques\nand evaluating the model on a separate validation set to ensure it generalizes well to unseen data.\nResponse: Thank you for your suggestion. RF is inherently resistant to overfitting, as it builds an\nensemble of decision trees, each trained on random subsets of the data with randomly selected features\nat each split. To further ensure generalization, we applied systematic cross-validation to assess the\nmodel’s performance and tested the trained RF model on an independent validation set (Table S4 and\nFig. S22). This approach confirms the model’s robustness on unseen data, minimizing the risk of\noverfitting and ensuring reliable performance beyond the training set.\nIn the revised manuscript, we have provided a detailed explanation of how overfitting was addressed in\n“Splitting the training and validation sets” section of “Methods”:\n“To further validate the model’s performance and safeguard against overfitting, we conducted\n10-fold cross-validation. The dataset was randomly partitioned into 10 subsets, with 9 used\nfor training and 1 for validation in each iteration. This process was repeated across all subsets,\nensuring comprehensive evaluations. The averaged performance metrics (Table S4 and Fig.\nS22) confirmed the model’s ability to generalize effectively, without overfitting.\nAfter confirming the model’s robustness through cross-validation, we proceeded to retrain it\nusing 90% of the data (1,779 surveys, comprising 67,315 points) for training and 10% (198\nsurveys, comprising 15,335 points) for validation, each preserving the proportion of PGQ and\nccA\nNumber of decision trees\nnon-PGQ of the full dataset. This approach maximizes the model’s learning capacity by\nproviding a larger training set while preserving an independent validation set to assess\nperformance on unseen data. This final validation confirmed that the model maintained strong\npredictive accuracy, further enhancing its generalization ability.”\nTable S4. Statistics summarizing the results of applying the RF model to the validation set by 10-fold\ncross-validation at a probability cutoff of 0.50. The definitions of the statistics can be found in the\nMethods section.\nStatistics Value\nPrevalence 0.44 (0.40 ~ 0.48)\nAccuracy 0.77 (0.72 ~ 0.83)\nSensitivity 0.70 (0.60 ~ 0.79)\nPrecision 0.77 (0.72 ~ 0.83)\nSpecificity 0.83 (0.77 ~ 0.89)\nAUC 0.85 (0.79 ~ 0.89)\nCohen’s kappa coefficient 0.54 (0.43 ~ 0.64)\nNote: The value is presented as the mean (min ~ max).\nFig. S22| The ROC and the AUC of the RF model based on the validation set by 10-fold cross-validation\n(a-j).\nRPT\nAUC = 0.79\n1\n0.9\n0.8\n0.7\n0.6\n0.5\n0.4\n0.3\n0.2\n0.1\n0\n0 0.2 0.4 0.6 0.8 1\nFPR\nRPT\n0.9\n0.8\n0.7\n0.6\n0.5\n0.4\n0.3\n0.2\n0.1\n0\n0 0.2 0.4 0.6 0.8 1\nFPR\nRPT\nRPT\n0.9\n0.8\n0.7\n0.6\n0.5\n0.4\n0.3\n0.2\n0.1\n0\n0 0.2 0.4 0.6 0.8 1\nFPR\nRPT\na b\nAUC = 0.86\n1\n0.9\n0.8\n0.7\n0.6\n0.5\n0.4\n0.3\n0.2\n0.1\n0\n0 0.2 0.4 0.6 0.8 1\nc FPR d\ne f\ng h\ni j\nDiscuss the rationale for selecting the RF model over other machine learning algorithms,\nemphasizing its strengths in handling large datasets with numerous predictor variables and its\nability to capture complex interactions between predictors without requiring extensive parameter\ntuning.\nResponse: Thank you for this suggestion. In the revised manuscript, we have provided a detailed\nexplanation in “Splitting the training and validation sets” section of “Methods”, outlining the rationale\nbehind selecting the RF model over other machine learning algorithms:\n“In this study, the RF model was constructed to predict the spatio-temporal dynamic of\ngroundwater quality from natural and socio-economic factors. The RF model was chosen for\nits ability to handle large datasets with numerous predictor variables while efficiently\ncapturing complex interactions among them24. The RF model uses the bootstrap resampling\nmethod to build multiple decision trees. For the construction of each tree, samples are\nindependently selected; however, the distributions for all trees in the forest are the same, which\nguarantees the robustness of the model. In addition, the RF model offers key advantages such\nas preventing overfitting by introducing randomness at each decision node, which enhances\ngeneralization. By selecting random feature subsets for each tree, RF effectively addresses the\nchallenges of high-dimensionality in the data. It also reduces experimental noise and improves\nprediction accuracy through an ensemble approach, averaging predictions from multiple\ndecision trees. Furthermore, RF provides inherent feature importance measures, helping\nidentify key variables, and handles missing data efficiently, either through imputation or\nsplitting based on available features. With minimal parameter tuning required compared to\nother machine learning algorithms, RF presents itself as an ideal model for this study14,57,58.”\nReference:\nmodeling-A comprehensive review. Water Research233, 119745 (2023).\nNature625, 715-721 (2024).\ninspired algorithms for mapping groundwater nitrate concentration in a coastal multi-layer aquifer\nsystem. Journal of Cleaner Production343, 130900 (2022).\nMoreover, provide details on any data preprocessing steps undertaken, such as normalization or\nstandardization of the predictor variables, to ensure that they are on a comparable scale and to\nenhance the model’s performance. Explain the rationale for these preprocessing steps and how\nthey contribute to the overall robustness of the model. Providing a comprehensive explanation of\nthe hyperparameter tuning process, including how overfitting was checked, the rationale for\nselecting the RF model, and the data preprocessing steps, will enhance the transparency and\nreproducibility of the study. It will demonstrate the rigor of the model development process and\nincrease confidence in the robustness and reliability of the model predictions, ultimately\nstrengthening the scientific contribution of the research.\nResponse: Thank you for your suggestions. In general, we do not need to normalize or scale the\npredictors for an RF model. RF is based on decision trees, which splits the data based on feature\nthresholds. Since the splits are determined independently for each feature, the model is not sensitive to\nthe magnitude or distribution of the predictors. In the revised manuscript, we have included a detailed\nexplanation of the data preprocessing workflow, clearly outlining each step from raw data collection to\nthe generation of training samples:\n“To classify groundwater quality at the prefecture level, we implemented a series of data\npreprocessing steps. First, regional data were converted into 1 km² grid cells based on\nregional land-use types, with all cells within a prefecture assigned the corresponding\ngroundwater quality data to ensure consistency. The groundwater quality data were then\nbinarized, assigning a value of zero to non-PGQ and one to PGQ. The purpose of this\napproach was twofold: (i) to prioritize the fundamental health aspect of safety (non-PGQ) or\nunsafety (PGQ) of groundwater for drinking; and (ii) to address variations in precision\nresulting from diverse analysis methods employed in different data sources. Next, natural\nfactors in the predictors were extracted at a 1 km² resolution, assigning each grid cell a unique\nvalue representing local environmental conditions. Socio-economic factors were uniformly\napplied to all grid cells within a province, based on the sampling year. The dataset on land-\nuse types covers specific benchmark years: 1980, 1990, 1995, 2000, 2005, 2010, 2015, 2020,\nand every five years between 2020 and 2100. For other years, land-use data are assumed to\nbe identical to the adjacent prior benchmark year. For instance, land-use data in 1989 are\ntaken to be the same as those in 1980. For other dynamic factors, their cumulative average\nvalues from 1980 to the target year were calculated as the model inputs. This accounts for the\nlag in pollutant migration and effectively captures the long-term impacts on groundwater\nquality. Fertilizer consumption and pesticide consumption were averagely allocated per km2\nof agricultural land, with non-agricultural land assigned zero. Industrial wastewater\ndischarge, industrial solid waste discharge, domestic wastewater discharge, domestic solid\nwaste discharge, and groundwater supply were allocated per km2 based on population density\nfor each province. Per capita data for these variables were calculated by dividing the total\nvalue by the population of each province, then distributed per km2 based on population density.\nGroundwater exploitation rates were uniformly distributed across each province, with regions\nhaving a population density below 1 person/km² assigned a rate of zero. Other factors were\nallocated per km2 based on the resolution of raster data or vector data. Finally, groundwater\nquality categories for each year and region were matched with the corresponding natural and\nsocio-economic variables to create the training and validation datasets for the RF model.”\nFor explanation of other hyperparameter tuning process, including how overfitting was checked, and the\nrationale for selecting the RF model, please refer to our previous responses.\n2) While the methodology mentions using out-of-bag (OOB) error to assess the importance of\npredictor variables, it lacks a detailed analysis of how these predictors contribute to the model and\ntheir implications for groundwater quality. Conduct a more detailed analysis of the top predictors\nidentified by the RF model. For each top predictor, discuss its relative importance and provide\ninsights into how it influences groundwater quality. For instance, if soil type and industrial\nwastewater discharge are among the top predictors, explain their potential mechanisms of\ninfluence on groundwater quality. Additionally, include visual aids such as partial dependence plots\nto illustrate the relationship between these predictors and the predicted groundwater quality\noutcomes. A detailed analysis of predictor variable importance will enhance the interpretability of\nthe model and provide valuable insights into the factors driving groundwater quality, thereby\nstrengthening the scientific contribution of the study and informing policymakers and stakeholders\nabout the critical variables that need to be managed to improve groundwater quality.\nResponse: Thank you for your valuable suggestions. We have expanded our analysis of the top predictors\nidentified by the RF model using OOB error. These predictors include depth-based groundwater type,\nair temperature, aridity index, precipitation, and soil sand content (100-200 cm). We also provide detailed\nexplanations of how each predictor influences groundwater quality, supported by partial dependence\nplots to enhance interpretability. The revised manuscript includes the following texts to “Discussion”:\n“In this study, our model ranks the factors that significantly influence the prediction accuracy\nof groundwater quality, identifying depth-based groundwater type as the most critical\npredictor (Fig. S3). Sensitivity analysis, involving a ±10% change in each variable,\nhighlighted permeability as the most sensitive factor (Fig. S13). In addition, we analyzed the\ninterrelationships among the predictors (Fig. S2), as well as the partial dependence plots and\nfeature interaction charts (Figs. S14 and S15), demonstrating the complex interactions among\nthe predictors. Shallow groundwater is more vulnerable to surface contaminants, while deeper\ngroundwater is less impacted27. Regions with lower precipitation, higher temperatures, and\nhigher groundwater exploitation rates showed higher PGQ probabilities (Fig. S15). Climate-\ndriven changes directly affect recharge rates, influence groundwater replenishment, while\nindirectly alter groundwater usage patterns28. In arid and semi-arid regions, reduced rainfall\nleads to over-extraction and higher pollutant concentrations due to decreased dilution,\nwhereas humid regions see increased mobilization and dilution of pollutants29-30. As climate\nchange continues to reshape global temperature and precipitation patterns, these dynamics\nare expected to become increasingly complex. Additionally, soil sand content (100-200 cm)\nsignificantly affects recharge and contamination, with higher content increasing susceptibility\nto surface contaminants (Fig. S14f)31. It is essential to note that this analysis is based on a\nnationwide perspective, and the most critical predictors vary spatially, as shown in Fig. S7a.\nThese insights highlight the importance of informing policymakers and stakeholders about\nregion-specific factors that must be managed to improve groundwater quality.”\nAnd we have detailed the steps in “Assessment of predictor importance” section of “Methods”:\n“To further assess the sensitivity of each influencing factor on the model’s predictions62, we\nconducted a sensitivity analysis by systematically increasing or decreasing each continuous\npredictor by 10%. This approach allowed us to observe how these perturbations influenced\nthe predicted average probability of PGQ over the historical period from 1980 to 2020. By\ncomparing these modified predictions to the baseline scenario, where no changes to the\npredictors were applied, we were able to quantify the degree to which each factor contributed\nto variability in PGQ probability, providing insights into which predictors most significantly\naffect groundwater quality predictions under different conditions.\nFollowing the identification of key predictors, partial dependence plots (PDPs) were\ngenerated to visualize each predictor’s marginal effect on the probability of PGQ. The PDPs\nwere generated by holding the focal predictor constant across its empirical range while\nintegrating out the effects of all other predictors, thereby isolating the unique contribution of\nthe predictor to the response63. This was achieved by averaging the model predictions over the\ndistribution of the data, providing a clear visual representation of the dependency between the\npredictors and the probability of PGQ.\nIn addition, feature interaction charts (FICs) were constructed to explore the combined effect\nof two predictors on the predicted probability of PGQ63. It takes into account the range of\nvalues for the two concerned predictors. However, instead of holding one variable constant,\nFICs vary both variables simultaneously over their empirical ranges. The model’s predictions\nare then averaged over the joint distribution of the two predictors, and the combined effect is\nvisualized in a two-dimensional color-coded map. This approach highlights how the\ninteraction between two predictors influences the model output, thus providing a deeper\nunderstanding of the relationships within the data.”\nReference:\nIllinois. Vol. 532 (Illinois State Geological Survey, 1984).\nHydrology518, 250-266 (2014).\nEnvironmental Science and Technology51, 2895-2933 (2021).\ngroundwater in China: A review. Science of the Total Environment643, 967-993 (2018).\nJournal of Hydrology332, 187-205 (2007).\nparameters. Reliability Engineering & System Safety 93, 1563-1573 (2008).\na 1 b1\n2\n0.8 6\n12\n15 0.6 17\n19\n0.4 20 22\n24\n0.2 25\n26\n0 27\n28\n29\n-0.2 30\n35\n-0.4 36 38\n39\n-0.6 41\n42\n43\n-0.8 44\n51\n-1 52\nFig. S2| Pearson correlation analysis of predictors (a) before and (b) after variable pruning. The\nvariable pruning is carried out to remove redundant variables to prevent overfitting. The numbers in the\nx and y axis indicate different variables: 1. Porosity; 2. Permeability; 3. Parent material; 4. Sand content\nof soil (0-20cm); 5. Silt content of soil (0-20cm); 6. Clay content of soil (0-20cm); 7. pH (0-5cm); 8. Bulk\ndensity (BD) (0-5cm); 9. Cation exchange capacity (CEC) (0-5cm); 10. Soil and sedimentary deposit\n1 2 3 4 5 6 7 8 9 01 11 21 31 41 51 61 7811 91 02 12 22 32 42 52 62 72 82 92 03 13 23 33 43 53 63 73 83 93 04 14 24 34 44 54 64 74 84 94 05 15 25 1 2 6 21 51 71 91 02 22 42 52 62 72 82 92 03 53 63 83 93 14 24 34 44 15 25\n123 1\n1111 0123\n456789\n0\n0\n.\n.\n6\n8\n111\n111\n456\n789 0.4\n222 012\n0.2\n23 2222 4567\n0\n2233333333 8901234567\n-\n-\n0\n0\n.\n.\n4\n2\n3344444444\n44555\n8901234567\n89012\n-\n-\n-\n1\n0\n0\n.\n.\n8\n6\nthickness; 11. Soil organic carbon (SOC) (0-100cm); 12. Sand content of soil (100-200cm); 13. Silt\ncontent of soil (100-200cm); 14. Clay content of soil (100-200cm); 15. pH (100-200cm); 16. BD (100-\n200cm); 17. CEC (100-200cm); 18. Coarse fragment volume (CFVO) (100-200cm); 19. Soil type; 20.\nSeawater intrusion; 21. Digital elevation model (DEM); 22. Normalized difference vegetation index\n(NDVI); 23. Slope; 24. Topographic wetness index (TWI); 25. Water-table ratio (WTR); 26. Aquifer-\nmedia-based groundwater type; 27. Depth-based groundwater type; 28. Precipitation; 29. Temperature;\n30. Aridity; 31. Actual evapotranspiration (ET); 32. Priestley-Taylor alpha coefficient; 33. Humidity; 34.\nParticulate matter 10 micrometers or less in diameter (PM ); 35. Fertilizer consumption; 36. Pesticide\n10\nconsumption; 37. Coal consumption; 38. Industrial wastewater discharge; 39. Industrial solid waste\ndischarge; 40. Industrial gross domestic product (GDP); 41. Domestic wastewater discharge; 42.\nDomestic solid waste discharge; 43. Groundwater supply; 44. Groundwater exploitation rate; 45. Crude\noil consumption; 46. Natural gas consumption; 47. Non-ferrous metal mining; 48. GDP per capita; 49.\nEnvironmental investment; 50. Total wastewater discharge; 51. Land-use type; and 52. Groundwater\nquality category.\nFig. S3| Importance of the predictors in the final random forest (RF) model based on out-of-bag error\nanalysis.\na b\n0.05\n0.00\n-0.05\n-0.10\n-0.15\n-0.20\n-0.25\n1 2 3 4 5 6 7 8 9 1011121314151617181920\nPredictor\nFig. S13| Boxplots showing the difference in average national PGQ probability from 1980 to 2020\ncompared to the baseline scenario under increased (+10%) (a) and decreased (-10%) (b) conditions for\nvarious factors. The numbers in the x-axis represent the following factors: 1. Porosity; 2. Permeability;\n3. Clay content of soil (0-20 cm); 4. CEC (100-200 cm); 5. pH (100-200 cm); 6. Sand content of soil\n(100-200 cm); 7. NDVI; 8. TWI; 9. WTR; 10. Aridity; 11. Temperature; 12. Precipitation; 13.\nGroundwater supply; 14. Groundwater exploitation rate; 15. Fertilizer consumption; 16. Pesticide\nconsumption; 17. Industrial wastewater discharge; 18. Industrial solid waste discharge; 19. Domestic\nwastewater discharge; and 20. Domestic solid waste discharge.\nFig. S15| Feature interaction charts of the RF model. Subfigure (a) shows the combined effects of\nprecipitation and temperature, and (b) of precipitation and groundwater exploitation rate, on the\npredicted probability of PGQ. Both variables are evaluated across their full empirical ranges and are\nallowed to vary simultaneously. Model predictions are averaged over the joint distribution of these\nvariables. The resulting interactions are displayed on a two-dimensional, color-coded map, where the\ngradient from blue to red visually represents the shift from lower to higher probabilities of PGQ.\nQGP\nfo\nytilibaborp\nni\negnahC\n25th Percentile ~ 75th Percentile 1.5 IQR range Median Mean\n0.05\n0.00\n-0.05\n-0.10\n-0.15\n-0.20\n-0.25\n1 2 3 4 5 6 7 8 9 1011121314151617181920\nPredictor\nQGP\nfo\nytilibaborp\nni\negnahC\n25th Percentile ~ 75th Percentile 1.5 IQR range Median Mean\na b\na b c d\ne f g h\ni j k l\nm n o p\nq r s t\nFig. S14| Partial dependence plots for PGQ probability across multiple features. Subfigures (a-t) depict\nthe partial dependence of PGQ probability on different factors: (a) porosity, (b) permeability, (c) clay\ncontent of soil (0-20cm), (d) CEC (100-200cm), (e) pH (100-200cm), (f) sand content of soil (100-200cm),\n(g) NDVI, (h) TWI, (i) WTR, (j) aridity, (k) temperature, (l) precipitation, (m) groundwater supply, (n)\ngroundwater exploitation rate, (o) fertilizer consumption, (p) pesticide consumption, (q) industrial\nwastewater discharge, (r) industrial solid waste discharge, (s) domestic wastewater discharge, and (t)\ndomestic solid waste discharge, respectively. Each subfigure shows the relationship between the feature\nand the predicted PGQ probability, indicating how changes in each factor influence PGQ probability.\n3) The methodology describes the use of factorial experiments to isolate the impacts of specific\ndriver categories under different scenarios, but the implementation details of these scenarios and\nthe handling of uncertainties in future projections are insufficiently described. Expand on the\nscenario analysis methodology by providing a detailed description of how each scenario is\nconstructed, including the data sources, assumptions, and adjustments made for future projections.\nIncorporate a discussion on how uncertainties in climate, exploitation, and pollution data are\naddressed, such as using sensitivity analysis or Monte Carlo simulations to quantify the impact of\nthese uncertainties on model predictions. Present the results of these uncertainty assessments to\nhighlight the range of possible outcomes and their implications. A comprehensive and transparent\nscenario analysis, coupled with a thorough uncertainty assessment, will significantly strengthen\nthe study’s credibility by providing a clearer understanding of the potential future states of\ngroundwater quality under various conditions. This will aid in the development of robust and\ninformed management strategies, enhancing the study’s value to policymakers by offering a\nnuanced understanding of the risks and uncertainties involved.\nResponse: Thank you for your suggestions.\n(1) In the revised manuscript, we have provided a more detailed explanation of the scenario construction\nprocess in “Scenario settings” section of “Methods”, outlining the data sources, underlying assumptions,\nand adjustments made for future projections:\n“In Scenario I (Baseline), we assumed a consistent development intensity for all drivers.\nClimate, groundwater exploitation, and pollution discharges were maintained at the same\nlevels as the cumulative average values in 2020. Similarly, land-use type in Scenario I was\nalso assumed to remain unchanged from 2020. In contrast, Scenarios II-IV incorporated\nprojections to account for future variability. Specifically, climate data from one CMIP6\nmodel——EC-Earth3 were used as model input to project future PGQ probabilities65-67.\nAdditionally, land-use data in these scenarios were derived from the SSP-RCP global 1km\nland-use simulation dataset (2020-2100) for SSP1-1.9, SSP2-4.5, and SSP5-8.5 scenarios,\nupdated every five years68. To maintain consistency in representation, these land-use data were\naligned with those used for the historical period from 1980 to 2020.\nIn the three scenarios, annual change rates of -3%, 1%, and 3% were assigned as gradients\nfor groundwater exploitation and pollution discharge, respectively. These rates were\ndetermined based on historical trends69,70 and future projections71,72 related to industrial\nactivities, agricultural expansion, and environmental policies. The low rate of -3% reflects a\nproactive approach toward reducing resource use and pollution discharge, aligning with the\nsustainability-focused SSP1-1.9 scenario73. In contrast, the medium rate of 1% corresponds\nto a moderate development pathway, as seen in SSP2-4.5. This scenario embodies a balanced\napproach, allowing socio-economic development to progress while placing emphasis on\nenvironmental protection74. The high rate of 3% aligns with the SSP5-8.5 scenario,\ncharacterized by rapid economic growth and minimal environmental regulation75. This\nscenario anticipates high energy demands, with significant reliance on fossil fuels and\nindustrial expansion. We acknowledge that future development may exhibit more intricate\npatterns, and our current study may not fully capture this complexity and uncertainty. Then,\nthese rates were applied to the corresponding values from 2021 to 2100, and the cumulative\naverage was calculated as the model inputs in each year. Please refer to Table S5 for detailed\ninformation on the specific scenario assumptions.”\nReference:\n(2022).\n2100. National Tibetan Plateau / Third Pole Environment Data Center (2022).\n2100. National Tibetan Plateau / Third Pole Environment Data Center (2022).\nChina. (http://www.geodata.cn).\nPress, 1985-2021).\nStatistics Press, 1981-2021).\ngas emissions implications: An overview. Global Environmental Change42, 153-168 (2017).\nworld futures in the 21st century. Global Environmental Change42, 169-180 (2017).\ngreen growth paradigm. Global Environmental Change42, 237-250 (2017).\nof-the-road scenario for the 21st century. Global Environmental Change42, 251-267 (2017).\nfor the 21st century. Global Environmental Change42, 297-315 (2017).\nTable S5. Scenario settings.\nScenari Climat Descriptio Descriptio Land Descriptio Descriptio\nExploitation Pollution\no e n n use n n\nBaselin Baselin\nBaseline Baseline\ne e land-\nexploitation pollution\nⅠ climate Historical Historical use Historical Historical\nscenario scenario\nscenario scenario\n(EB) (PB)\n(CB) (LB)\nClimate\nClimate change\nExploitation Pollution\nchange Low land- Low\nⅡ SSP1-1.9 reduction SSP1-1.9 reduction\nscenario （-3%） use （-3%）\nscenario (E1) scenario (P1)\n(C1) scenario\n(L1)\nClimate\nClimate Exploitation change Pollution\nchange intensificatio Medium land- intensificatio Medium\nⅢ SSP2-4.5 SSP2-4.5\nscenario n scenario （+1%） use n scenario （+1%）\n(C2) (E2) scenario (P2)\n(L2)\nClimate\nClimate Exploitation change Pollution\nchange intensificatio High land- intensificatio High\nⅣ SSP5-8.5 SSP5-8.5\nscenario n scenario （+3%） use n scenario （+3%）\n(C3) (E3) scenario (P3)\n(L3)\n(2) To evaluate the impact of uncertainties in climate, groundwater exploitation, and pollution data on\nmodel predictions, we conducted Monte Carlo simulations using n=50 sets of randomized input data\nrepresenting various uncertain future scenarios. This approach enables us to quantify output uncertainty\nand highlight the range of possible outcomes and their implications. By incorporating these uncertainty\nassessments, we offer a clearer understanding of potential future groundwater quality under different\nconditions, providing policymakers with valuable insights into associated risks.\nWe have detailed the steps and methods used for uncertainty analysis in the “Uncertainty analysis”\nsection of “Methods”:\n“Uncertainty analysis\nTo evaluate the impact of uncertainties in climate change, groundwater exploitation, and\npollution charge on model predictions, we applied Monte Carlo (MC) simulations to quantify\nthe variability in groundwater quality predictions due to input uncertainties76. This analysis\nwas conducted for the period 1980-2100, with data processed at 5-year intervals for 1980-\n2015, annually from 2016 to 2020, and at 5-year intervals for future projections from 2025 to\n2100 under four scenarios.\nIn each scenario-year combination, MC simulations with 50 members were run to introduce\nvariability in key natural and socio-economic predictors. Perturbations were applied to these\ndominant factors, such as climate, groundwater exploitation and pollution charge, reflecting\ntheir real-world uncertainties. For instance, temperature varied within a ±2°C range, while\npredictors like precipitation, fertilizer, and pesticide use were perturbed by up to ±20%.\nThe results of these simulations can be used to evaluate the range of possible groundwater\nquality outcomes based on uncertainties in input data. For each year and scenario, we\ncalculated the mean and standard deviation of the predicted probabilities based on the\nensemble outcomes, providing a robust estimate of both expectation and variability in\ngroundwater quality predictions.”\nIn addition, we have included an uncertainty analysis in “Discussion”:\n“Although visual comparisons and statistical evaluations suggest that the model estimates are\ngenerally reasonable, notable uncertainties remain. Monte Carlo simulations were used to\nassess uncertainties in PGQ probability predictions, providing insights into spatial and\ntemporal variability. As shown in Fig. S18(d-f) and Fig. S19(e-h), uncertainties are moderate\noverall but higher in North and Northeast China from 1980 to 2100, likely due to socio-\neconomic variability and complex groundwater interactions. Under Scenario IV,\nuncertainties become more pronounced (Fig. S20), likely driven by increased emissions and groundwater\nexploitation. This highlights the need for adaptive groundwater management strategies,\nparticularly in high-emission scenarios, to mitigate adverse impacts and enhance prediction\nreliability in vulnerable regions.”\nReference:\nEditorial Note: The maps in Fig S18 of this Peer Review File are from Standard Map Service of the\nMinistry of Natural Resources, People’s Republic of China [Map Review Number: GS(2019)1822].\na\nd\nb e\nc f\nMean value of the probability of PGQ Standard deviation of the probability of PGQ\nFig. S18| Spatial distribution of mean and standard deviation of PGQ probability from 1980 to 2020\nbased on Monte Carlo simulations. Subfigures (a-c) show the spatial distribution of the mean of PGQ\nprobability for the years 1980, 2000, and 2020, respectively, and the insets in (a-c) indicate the area\nratios of different PGQ probability levels; Subfigures (d-f) present the standard deviation of PGQ\nprobability, representing the prediction uncertainties for the same years, and the insets in (d-f) indicate\nthe area ratios of different uncertainty levels of PGQ probability predictions.\n0891\n0002\nMinistry of Natural Resources, People’s Republic of China [Map Review Number: GS(2019)1822]\na\ne\nb f\nc g\nd h\nMean value of the probability of PGQ Standard deviation of the probability of PGQ\nFig. S19| Spatial distribution of mean and standard deviation of PGQ probability under four scenarios\nin 2050 based on Monte Carlo simulations. Subfigures (a-d) show the mean of PGQ probability under\nthe four scenarios in 2050, and the insets in (a-d) indicate the area ratios of different PGQ probability\nlevels; Subfigures (e-h) present the uncertainty (standard deviation) of PGQ probability under the four\nscenarios in 2050, and the insets in (e-h) indicate the area ratios of different uncertainty levels in PGQ\nprobability predictions.\nIoiranecS\nIIoiranecS\nIIIoiranecS\nVIoiranecS\na b c\nd e\nOutlier\n1.5 IQR range\n75th Percentile\nMean\nMedian\n25th Percentile\n1.5 IQR range\nOutlier\nFig. S20| Temporal variation of PGQ probability from 1980 to 2100 across different scenarios based on\nMonte Carlo simulations. (a) Probability distribution from 1980 to 2020. (b-e) Probability distributions\nunder Scenarios I, II, III, and IV from 2025 to 2100, respectively.\n4) The discussion section effectively identifies major drivers of groundwater quality deterioration,\nincluding pollution discharge and groundwater exploitation. However, it should provide a more\nnuanced analysis of the interactions between these drivers and climatic factors.\nResponse: Thank you for your feedback. In the revised manuscript, we have provided a more nuanced\nanalysis of the interactions among major drivers in “Discussion”:\n“In this study, our model ranks the factors that significantly influence the prediction accuracy\nof groundwater quality, identifying depth-based groundwater type as the most critical\npredictor (Fig. S3). Sensitivity analysis, involving a ±10% change in each variable,\nhighlighted permeability as the most sensitive factor (Fig. S13). In addition, we analyzed the\ninterrelationships among the predictors (Fig. S2), as well as the partial dependence plots and\nfeature interaction charts (Fig. S14 and S15), demonstrating the complex interactions among\nthe predictors. Shallow groundwater is more vulnerable to surface contaminants, while deeper\ngroundwater is less impacted27. Regions with lower precipitation, higher temperatures, and\nhigher groundwater exploitation rates showed higher PGQ probabilities (Fig. S15). Climate-\ndriven changes directly affect recharge rates, influence groundwater replenishment, while\nindirectly alter groundwater usage patterns28. In arid and semi-arid regions, reduced rainfall\nleads to over-extraction and higher pollutant concentrations due to decreased dilution,\nwhereas humid regions see increased mobilization and dilution of pollutants29-30. As climate\nchange continues to reshape global temperature and precipitation patterns, these dynamics\nare expected to become increasingly complex. Additionally, soil sand content (100-200 cm)\nsignificantly affects recharge and contamination, with higher content increasing susceptibility\nto surface contaminants (Fig. S14f)31. It is essential to note that this analysis is based on a\nnationwide perspective, and the most critical predictors vary spatially, as shown in Fig. S7a.\nThese insights highlight the importance of informing policymakers and stakeholders about\nregion-specific factors that must be managed to improve groundwater quality.”\nReference:\nIllinois. Vol. 532 (Illinois State Geological Survey, 1984).\nHydrology518, 250-266 (2014).\nEnvironmental Science and Technology51, 2895-2933 (2021).\ngroundwater in China: A review. Science of the Total Environment643, 967-993 (2018).\nJournal of Hydrology332, 187-205 (2007).\na 1 b1\n2\n0.8 6\n12\n15 0.6 17\n19\n0.4 20\n22\n24\n0.2 25\n26\n0 27\n28\n29\n-0.2 30\n35\n-0.4 36 38\n39\n-0.6 41\n42\n43\n-0.8 44\n51\n-1 52\nFig. S2| Pearson correlation analysis of predictors (a) before and (b) after variable pruning. The\nvariable pruning is carried out to remove redundant variables to prevent overfitting. The numbers in the\nx and y axis indicate different variables: 1. Porosity; 2. Permeability; 3. Parent material; 4. Sand content\nof soil (0-20cm); 5. Silt content of soil (0-20cm); 6. Clay content of soil (0-20cm); 7. pH (0-5cm); 8. Bulk\ndensity (BD) (0-5cm); 9. Cation exchange capacity (CEC) (0-5cm); 10. Soil and sedimentary deposit\nthickness; 11. Soil organic carbon (SOC) (0-100cm); 12. Sand content of soil (100-200cm); 13. Silt\ncontent of soil (100-200cm); 14. Clay content of soil (100-200cm); 15. pH (100-200cm); 16. BD (100-\n200cm); 17. CEC (100-200cm); 18. Coarse fragment volume (CFVO) (100-200cm); 19. Soil type; 20.\nSeawater intrusion; 21. Digital elevation model (DEM); 22. Normalized difference vegetation index\n(NDVI); 23. Slope; 24. Topographic wetness index (TWI); 25. Water-table ratio (WTR); 26. Aquifer-\nmedia-based groundwater type; 27. Depth-based groundwater type; 28. Precipitation; 29. Temperature;\n30. Aridity; 31. Actual evapotranspiration (ET); 32. Priestley-Taylor alpha coefficient; 33. Humidity; 34.\nParticulate matter 10 micrometers or less in diameter (PM ); 35. Fertilizer consumption; 36. Pesticide\n10\nconsumption; 37. Coal consumption; 38. Industrial wastewater discharge; 39. Industrial solid waste\ndischarge; 40. Industrial gross domestic product (GDP); 41. Domestic wastewater discharge; 42.\nDomestic solid waste discharge; 43. Groundwater supply; 44. Groundwater exploitation rate; 45. Crude\noil consumption; 46. Natural gas consumption; 47. Non-ferrous metal mining; 48. GDP per capita; 49.\nEnvironmental investment; 50. Total wastewater discharge; 51. Land-use type; and 52. Groundwater\nquality category.\n1 2 3 4 5 6 7 8 9 01 11 21 31 41 51 61 7811 91 02 12 22 32 42 52 62 72 82 92 03 13 23 33 43 53 63 73 83 93 04 14 24 34 44 54 64 74 84 94 05 15 25 1 2 6 21 51 71 91 02 22 42 52 62 72 82 92 03 53 63 83 93 14 24 34 44 15 25\n1 1\n1 111111111 0 123456789\n23456789\n0\n0\n0\n.\n.\n.\n4\n6\n8\n2222222222333 0123456789012\n-\n0\n0\n0\n.2\n.2\n33333334 34567890\n-0.4\n444444444 123456789\n-\n-\n0\n0\n.\n.\n8\n6\n555 012\n-1\nFig. S3| Importance of the predictors in the final random forest (RF) model based on out-of-bag error\nanalysis.\na b\n0.05\n0.00\n-0.05\n-0.10\n-0.15\n-0.20\n-0.25\n1 2 3 4 5 6 7 8 9 1011121314151617181920\nPredictor\nFig. S13| Boxplots showing the difference in average national PGQ probability from 1980 to 2020\ncompared to the baseline scenario under increased (+10%) (a) and decreased (-10%) (b) conditions for\nvarious factors. The numbers in the x-axis represent the following factors: 1. Porosity; 2. Permeability;\n3. Clay content of soil (0-20 cm); 4. CEC (100-200 cm); 5. pH (100-200 cm); 6. Sand content of soil\n(100-200 cm); 7. NDVI; 8. TWI; 9. WTR; 10. Aridity; 11. Temperature; 12. Precipitation; 13.\nGroundwater supply; 14. Groundwater exploitation rate; 15. Fertilizer consumption; 16. Pesticide\nconsumption; 17. Industrial wastewater discharge; 18. Industrial solid waste discharge; 19. Domestic\nwastewater discharge; and 20. Domestic solid waste discharge.\nQGP\nfo\nytilibaborp\nni\negnahC\n25th Percentile ~ 75th Percentile 1.5 IQR range Median Mean\n0.05\n0.00\n-0.05\n-0.10\n-0.15\n-0.20\n-0.25\n1 2 3 4 5 6 7 8 9 1011121314151617181920\nPredictor\nQGP\nfo\nytilibaborp\nni\negnahC\n25th Percentile ~ 75th Percentile 1.5 IQR range Median Mean\na b c d\ne f g h\ni j k l\nm n o p\nq r s t\nFig. S14| Partial dependence plots for PGQ probability across multiple features. Subfigures (a-t) depict\nthe partial dependence of PGQ probability on different factors: (a) porosity, (b) permeability, (c) clay\ncontent of soil (0-20cm), (d) CEC (100-200cm), (e) pH (100-200cm), (f) sand content of soil (100-200cm),\n(g) NDVI, (h) TWI, (i) WTR, (j) aridity, (k) temperature, (l) precipitation, (m) groundwater supply, (n)\ngroundwater exploitation rate, (o) fertilizer consumption, (p) pesticide consumption, (q) industrial\nwastewater discharge, (r) industrial solid waste discharge, (s) domestic wastewater discharge, and (t)\ndomestic solid waste discharge, respectively. Each subfigure shows the relationship between the feature\nand the predicted PGQ probability, indicating how changes in each factor influence PGQ probability.\na b\nFig. S15| Feature interaction charts of the RF model. Subfigure (a) shows the combined effects of\nprecipitation and temperature, and (b) of precipitation and groundwater exploitation rate, on the\npredicted probability of PGQ. Both variables are evaluated across their full empirical ranges and are\nallowed to vary simultaneously. Model predictions are averaged over the joint distribution of these\nvariables. The resulting interactions are displayed on a two-dimensional, color-coded map, where the\ngradient from blue to red visually represents the shift from lower to higher probabilities of PGQ.\nAddressing the limitations of the study, such as the assumption of constant land-use types in future\nscenarios and the exclusion of certain predictor variables due to high correlation, is crucial.\nResponse: Thank you for your valuable suggestion.\n(1) We have enhanced our model to account for dynamic land use changes from 2020 to 2100 across\nmultiple scenarios. Additionally, we have conducted an analysis of uncertainties associated with these\nfuture projections. For further details, please refer to our response to your third comment.\n(2) To improve the parsimony of the model, we excluded certain predictor variables that exhibited high\ncorrelation, a measure implemented to avoid multicollinearity. However, we acknowledge this as a\nlimitation, as it may exclude potentially relevant factors. In the revised manuscript, we have added some\ndiscussion about the impact of excluding some predictors and the potential implications for interpreting\nour results:\n“The use of regional averages instead of point-specific values, variability in sampling methods,\nand the exclusion of highly correlated predictor variables—though necessary for model\nparsimony—may have led to an incomplete representation of complex interactions, potentially\nunderestimating the contributions of correlated factors.”\nMoreover, the policy implications of the findings should be explored in greater depth, offering\nspecific recommendations for mitigation strategies to address the identified drivers of groundwater\nquality degradation. Highlighting the socio-economic implications and feasibility of these strategies\nin different regions of China will enhance the manuscript’s practical relevance for policymakers\nand stakeholders involved in groundwater management. By integrating these improvements, the\nresults and discussion sections will present a more comprehensive and impactful analysis,\nstrengthening the study’s scientific rigor and applicability.\nResponse: Thank you for your valuable suggestions. In the revised manuscript, we have expanded the\ndiscussion to include specific mitigation strategies concerning the key drivers of groundwater quality\ndegradation. Additionally, we have emphasized the socio-economic implications and feasibility of these\nstrategies across various regions of China, taking into account the unique environmental and economic\nconditions in different regions. The following texts have been added to “Discussion”:\n“In response to diverse development and climate change trajectories, China urgently needs\nproactive and flexible strategies to manage groundwater resources and control industrial\npollution, with a focus on enforcing strict compliance and discharge standards for\nmanufacturing enterprises, while prioritizing sustainable development goals and integrating\nclimate resilience tactics to ensure the long-term viability of these assets, particularly in\nvulnerable areas43-45. In North China, particularly the North China Plain, reducing\ngroundwater over-exploitation through stricter regulations and promoting water-saving\nirrigation is crucial46. Additionally, incentivizing sustainable agricultural practices, such as\ncrop rotation and controlled fertilizer application, can also mitigate groundwater\ncontamination47. In Northwest China, policies should target reducing industrial water use and\nenhancing recharge programs, while encouraging water recycling and less water-intensive\ntechnologies48. In Southern China, particularly in regions prone to extreme weather events,\neffective flood control and improved land-use planning are needed to reduce contaminant\ntransport during floods, alongside stronger waste management regulations to minimize\ncontamination risks49,50.”\nReference:\nNature Climate Change3, 621-626 (2013).\n(2023).\nGlobal threats, opportunities and realising the potential of groundwater. Science of the Total\nEnvironment811, 152471 (2022).\ntypical alluvial-pluvial fans, North China Plain. Environmental Science and Pollution Research, 1-18\n(2022).\nof widespread adoption of potential technology solutions to reduce water use and pollution:\nApplication to China’s textile industry. Environmental Impact Assessment Review79, 106293 (2019).\ndownstream flooding based on a land use change analysis and the concept of ecosystem-based\ndisaster risk reduction. Journal of Environmental Management287, 112341 (2021).\n(World Health Organization. Regional Office for Europe, 2011).\nReviewer #2 (Remarks to the Author):\nThe study of groundwater quality evolution is of significant importance and relevance, given its pivotal\nrole in sustaining freshwater resources and maintaining ecosystem health. The application of machine\nlearning, particularly the random forest model, to analyze and predict the spatiotemporal dynamics of\ngroundwater quality across entire China is innovative and represents a valuable contribution to the field.\nThe manuscript is well organized and easy to follow. However, there are several areas that require\nclarification and additional detail to enhance the robustness and transparency of this study. Below are\nsome comments and suggestions aimed at improving the manuscript:\n1. Definition of Poor Groundwater Quality (PGQ): The study does not clearly distinguish between\nPGQ caused by human activities and that resulting from geogenic origins, two distinct categories\ninfluenced by different driving forces. Moreover, different groundwater pollutants have different\nsources and influencing factors. For example, nitrogen pollution from non-point sources is\nsignificantly affected by surface input and migration through the unsaturated zone, yet the\nmachine learning model does not account for the impact of the unsaturated zone. Similarly, heavy\nmetals are largely influenced by redox conditions, especially in cases where naturally poor-quality\ngroundwater is more affected by hydrogeological conditions than by surface processes.\nConsequently, the machine learning model may overemphasize the role of surface soil, potentially\noverlooking other critical factors.\nResponse: Thank you for your valuable feedback on the definition of PGQ and the potential limitations\nof our approach.\n(1) Our study defines PGQ based on the latest Chinese groundwater quality standards, classifying\ngroundwater as either poor or not poor for drinking purposes. Since the original groundwater quality data\ncollected from the literature lacks explicit information about contamination sources, we were unable to\nclearly differentiate between PGQ caused by human activities and geogenic factors when training the\nmodel. However, we examined the spatial distribution of the dominant drivers of groundwater quality\nchanges across China, as shown in Fig. S7(a). This analysis provides a rough distinction between PGQ\ndriven by anthropogenic activities and that arising from natural processes.\nIn the revised manuscript, we have added the following text to “Discussion” to acknowledge the\nlimitations of our approach:\n“Furthermore, the lack of distinction between geogenic and anthropogenic sources of PGQ\nand depth differentiation among phreatic, confined, and unspecified groundwater samples\ncould affect the accuracy of predictions.”\nEditorial Note: The map in Fig S7a of this Peer Review File is from Standard Map Service of the\nMinistry of Natural Resources, People’s Republic of China [Map Review Number: GS(2019)1822].\na b\nc d\nFig. S7| Dominant drivers of groundwater quality changes and their relative contributions from 1980 to\n2020: (a) Dominant factors influencing the changes in PGQ probability over China from 1980 to 2020;\n(b) The relative contributions (RelCon) of each driver category from 1980 to 2100 over China; (c) The\narea ratios of regions with different dominant drivers over China; (d) The ratios of population affected\nby different dominant drivers over China. Here, ‘C’, ‘E’, ‘L’, ‘A’, ‘D’ and ‘I’ represent\nclimate, exploitation, land use, agricultural discharge, domestic discharge, and industrial discharge,\nrespectively.\n(2) Additionally, we recognize that different pollutants come from distinct sources and are influenced by\nvarious factors, with hydrogeological processes being just as important as surface processes. To address\nthis, our machine learning model integrates subsurface factors alongside surface soil characteristics. For\nexample, we incorporate soil properties at depths of 100-200 cm and account for various groundwater\ntypes based on aquifer media, such as pore water in loose sediments, karst water in carbonate rocks, and\nfissure water in igneous and metamorphic rocks (Table S2). This approach ensures that deeper\nhydrogeological influences are adequately considered without overemphasizing surface soil factors.\nIn the revised manuscript, we have clarified the principle of predictor collection in the “Predictor\ncollection and pre-processing” section of “Methods”:\n“In our analysis, we meticulously built the relationship between predictors and groundwater\nquality category (PGQ or non-PGQ). As listed in Table S2, the predictors included two classes\nof steady-state factors (soil properties and geographic and hydrogeologic characteristics) and\nfour classes of time-variant factors (i.e., climate change, pollution discharge, groundwater\nexploitation, and land-use change).\nGroundwater quality is influenced by various environmental and chemical processes,\nincluding redox reactions, contaminant transport, and natural filtration mechanisms51. Given\nChina’s vast and diverse landscape, along with significant regional variations in data\navailability and quality, we prioritized accessible and quantifiable predictors to ensure\nbroader applicability and consistency in groundwater quality predictions across different\nregions. Accordingly, we selected predictors that not only capture redox-related influences\nindirectly but also reflect other critical mechanisms affecting groundwater quality, providing\na comprehensive foundation for reliable assessment. To effectively capture these complex\nprocesses, we selected predictors such as soil properties (e.g., porosity, permeability, soil type,\nand soil composition across different soil layers) to reflect contaminant adsorption and\nbiological degradation, geographic and hydrogeologic characteristics (e.g., topographic\nwetness index, water-table ratio, and groundwater type) to indicate contaminant accumulation\nand natural dilution, and climatic factors (e.g., temperature, precipitation, and aridity) to\nrepresent groundwater recharge, evaporation, and pollutant mobility17,19,29. Additionally,\npollution discharge (agricultural, industrial, and domestic sources) can significantly affect\ngroundwater quality by increasing contaminant loads, while land-use change, can disrupt the\nnatural infiltration and recharge processes52. Notably, groundwater exploitation (e.g.,\ngroundwater supply and groundwater exploitation rate) emerges as a crucial anthropogenic\nperturbation factor53,54. Excessive groundwater extraction lowers water tables, concentrates\ncontaminants, and, in coastal areas, can trigger seawater intrusion, increasing salinity34,55,56.\nIt can also mobilize naturally occurring contaminants like arsenic and fluoride, further\ndegrading water quality29.”\nReference:\n866-868 (2013).\nEnvironmental Science and Technology51, 2895-2933 (2021).\nreaches deeper depths in heavily pumped aquifer systems. Nature Communications13, 5263 (2022).\nJournal of Biomedical Research27, 3855-3867 (2024).\n279-310 (Springer Netherlands, 2013).\nScience & Technology57, 5125-5136 (2023).\n7, 413-422 (2024).\nNature625, 715-721 (2024).\nin the eastern coastal of Laizhou Bay, China. Environmental Geology54, 1763-1770 (2008).\n“Table S2. Predictors for the groundwater qualityclassification model. The 25 variables ultimately\nused in the final prediction model are indicated in bold. (excerpt)\nResolution\nPredictor Period\nTemporal Spatial\nSoil\nPorosity7 steady-state Polygon /\nPermeability7 steady-state Polygon /\nSoil type8 steady-state 1:4 000 000 /\nParent material8 steady-state 1:1 000 000 /\nSand content of soil (0-20cm)9 steady-state 1 km /\nSilt content of soil (0-20cm)9 steady-state 1 km /\nClay content of soil (0-20cm)9 steady-state 1 km /\npH (0-5cm)9 steady-state 1 km /\nBD (0-5cm)9 steady-state 1 km /\nCEC (0-5cm)9 steady-state 1 km /\nSoil and sedimentary deposit thickness9 steady-state 1 km /\nSOC (0-100cm)9 steady-state 1 km /\nSand content of soil (100-200cm)9 steady-state 1 km /\nSilt content of soil (100-200cm)9 steady-state 1 km /\nClay content of soil (100-200cm)9 steady-state 1 km /\npH (100-200cm)9 steady-state 1 km /\nBD (100-200cm)9 steady-state 1 km /\nCEC (100-200cm)9 steady-state 1 km /\nCFVO (100-200cm)9 steady-state 1 km /\nGeography and hydrogeology\ncoastal buffer zones\nSeawater intrusion steady-state /\n(distance≤ 10 km)\nDEM10 steady-state 90 m /\nNDVI10 steady-state 500 m /\nSlope11 steady-state 1 km /\nTWI11 steady-state 1 km /\nWTR12 steady-state 1 km /\nAquifer-media-based groundwater type13 steady-state 1:4 000 000 /\nDepth-based groundwater type steady-state as recorded /\nNote: The population density data was not utilized directly as a predictive variable. Instead, it was\nused to derive the spatial distribution of groundwater exploitation and pollution discharge, which\nwe only have data at provincial scale. Specifically, we downscaled these data by assuming they are\nproportional to population density within each province. Aquifer-media-based groundwater types\ninclude loose sediment pore water, carbonate rock fissure cavern water, clastic rock void fissure\nwater, magmatic rock metamorphic rock fissure water, and frozen layer water. Groundwater can\nalso be classified by depth into phreatic, confined, and unspecified types. Additionally, land-use\ntypes encompass agricultural land, forest, grassland, water body, built-up land, and unused land.”\n2. Hydrogeological Factors: While this study developed a machine learning model to construct\nannual probabilistic maps of poor groundwater quality across China using 25 natural and social-\neconomic factors, few of these are explicitly related to hydrogeological processes. Groundwater\nflow and solute transport are clearly significant for predicting groundwater quality. As an example,\nthe authors mention a factor called ‘Aquifer-media-based groundwater type’ in the supporting\ninformation, but its establishment is not clearly explained. Does this refer to distinctions such as\nconfined vs. unconfined, porous media vs fractured media, or other characteristics?\nResponse: Thank you for your feedback. In our model, the variable “Aquifer-Media-Based Groundwater\nType” classifies aquifers into several types, including pore water in loose sediments, fissure cavern water\nin carbonate rocks, void fissure water in clastic rocks, fissure water in magmatic and metamorphic rocks,\nand frozen layer water. We also account for the depth to the water table, distinguishing among unconfined,\nconfined, and unspecified groundwater types. These categories reflect the intrinsic characteristics of the\naquifer media and groundwater, which significantly influence water flow and solute transport. These\ncategories are detailed in the revised Table S2.\n“Table S2. Predictors for the groundwater quality classification model. The 25 variables ultimately used\nin the final prediction model are indicated in bold.(excerpt)\nResolution\nGeography and hydrogeology Period\nTemporal Spatial\ncoastal buffer zones\nSeawater intrusion steady-state /\n(distance≤ 10 km)\nDEM10 steady-state 90 m /\nNDVI10 steady-state 500 m /\nSlope11 steady-state 1 km /\nTWI11 steady-state 1 km /\nWTR12 steady-state 1 km /\nAquifer-media-based groundwater type13 steady-state 1:4 000 000 /\nDepth-based groundwater type steady-state as recorded /\nNote: The population density data was not utilized directly as a predictive variable. Instead, it was used\nto derive the spatial distribution of groundwater exploitation and pollution discharge, which we only\nhave data at provincial scale. Specifically, we downscaled these data by assuming they are proportional\nto population density within each province. Aquifer-media-based groundwater types include loose\nsediment pore water, carbonate rock fissure cavern water, clastic rock void fissure water, magmatic rock\nmetamorphic rock fissure water, and frozen layer water. Groundwater can also be classified by depth into\nphreatic, confined, and unspecified types. Additionally, land-use types encompass agricultural land,\nforest, grassland, water body, built-up land, and unused land.”\n3. Chemical Processes: As noted in comment #1, chemical processes are critical in determining\ngroundwater quality. For example, nitrate (NO3), one of the indicators considered in this study for\ngroundwater quality, is redox-sensitive and undergoes denitrification under reducing conditions.\nNitrate concentrations in groundwater are therefore strongly dependent on redox conditions,\nregardless of nitrogen input at the land surface. How redox conditions are considered in assessing\nthe groundwater quality is an area that warrants further exploration.\nResponse: Thank you for your insightful comment. In response to your observation regarding the role\nof chemical processes, particularly redox conditions, in determining groundwater quality, we would like\nto clarify the following:\n(1) We fully acknowledge the significant role redox conditions play in shaping groundwater quality,\nespecially for redox-sensitive indicators like nitrate. While our model does not explicitly include a redox\npotential variable, certain predictors—such as climate conditions, soil properties, and groundwater\ntypes—indirectly capture variations that influence denitrification and other redox-related processes.\n(2) We recognize that redox potential is particularly crucial for pollutants like heavy metals, which are\nsensitive to these conditions. However, many groundwater quality indicators are less impacted by redox\npotential. Our model prioritizes factors that affect a broad range of indicators, as outlined in the latest\nChinese groundwater quality standards, rather than focusing exclusively on contaminants sensitive to\nredox changes. Additionally, the limited availability and consistency of redox data across China led us\nto focus on more accessible and quantifiable factors.\n(3) We emphasize that our study operates on a large national scale with a spatial resolution of 1 km²,\nwhich restricts the ability to capture fine-scale biochemical processes, such as site-specific denitrification.\nSuch processes require high-resolution data for accurate representation, which is beyond the scope of\nthis study. Therefore, we relied on indirect predictors that can capture these processes at a broader scale.\nIn the revised manuscript, we have added the following text to clarify these points:\nDiscussion:\n“Lastly, our model prioritizes factors influencing overall groundwater quality rather than\nspecific indicators, potentially differing from more targeted assessments.”\nPredictor collection and pre-processing:\n“In our analysis, we meticulously built the relationship between predictors and groundwater\nquality category (PGQ or non-PGQ). As listed in Table S2, the predictors included two classes\nof steady-state factors (soil properties and geographic and hydrogeologic characteristics) and\nfour classes of time-variant factors (i.e., climate change, pollution discharge, groundwater\nexploitation, and land-use change).\nGroundwater quality is influenced by various environmental and chemical processes,\nincluding redox reactions, contaminant transport, and natural filtration mechanisms51. Given\nChina’s vast and diverse landscape, along with significant regional variations in data\navailability and quality, we prioritized accessible and quantifiable predictors to ensure\nbroader applicability and consistency in groundwater quality predictions across different\nregions. Accordingly, we selected predictors that not only capture redox-related influences\nindirectly but also reflect other critical mechanisms affecting groundwater quality, providing\na comprehensive foundation for reliable assessment. To effectively capture these complex\nprocesses, we selected predictors such as soil properties (e.g., porosity, permeability, soil type,\nand soil composition across different soil layers) to reflect contaminant adsorption and\nbiological degradation, geographic and hydrogeologic characteristics (e.g., topographic\nwetness index, water-table ratio, and groundwater type) to indicate contaminant accumulation\nand natural dilution, and climatic factors (e.g., temperature, precipitation, and aridity) to\nrepresent groundwater recharge, evaporation, and pollutant mobility17,19,29. Additionally,\npollution discharge (agricultural, industrial, and domestic sources) can significantly affect\ngroundwater quality by increasing contaminant loads, while land-use change, can disrupt the\nnatural infiltration and recharge processes52. Notably, groundwater exploitation (e.g.,\ngroundwater supply and groundwater exploitation rate) emerges as a crucial anthropogenic\nperturbation factor53,54. Excessive groundwater extraction lowers water tables, concentrates\ncontaminants, and, in coastal areas, can trigger seawater intrusion, increasing salinity34,55,56.\nIt can also mobilize naturally occurring contaminants like arsenic and fluoride, further\ndegrading water quality29.”\nReference:\n866-868 (2013).\nEnvironmental Science and Technology51, 2895-2933 (2021).\nreaches deeper depths in heavily pumped aquifer systems. Nature Communications13, 5263 (2022).\nJournal of Biomedical Research27, 3855-3867 (2024).\n279-310 (Springer Netherlands, 2013).\nScience & Technology57, 5125-5136 (2023).\n7, 413-422 (2024).\nNature625, 715-721 (2024).\nin the eastern coastal of Laizhou Bay, China. Environmental Geology54, 1763-1770 (2008).\n4. Temporal Variations: This study aims to investigate the temporal changes in Chinese\ngroundwater quality; however, only social-economic factors were considered temporally dynamic,\nwhile soil and environmental factors were treated steady-state. It’s noteworthy that environmental\nconditions vary over time, especially the depth to the water table, which in turn affects contaminant\ntransport from the land surface to the underlying aquifer system. Considering these factors as\nsteady state may underestimate their impacts, and this should be addressed in the paper.\nResponse: Thank you for the valuable suggestion. We would like to address your concern as follows:\n(1) Due to limitations in the availability and validity of dynamic data for soil and hydrogeological\nvariables, we were unable to incorporate temporal changes for these factors directly. The lack of\nconsistent, high-resolution temporal data across the study period constrained our ability to model their\nfull dynamism. Therefore, we used available steady-state data to represent these factors within the model.\n(2) To account for the influence of fluctuations in soil and environmental factors on groundwater quality\nthrough changes in contaminant transport, we conducted a sensitivity analysis by introducing a ±10%\nperturbation to these predictors. This approach helps us understand the impact of potential temporal\nvariability of predictors within the constraints of available data.\nWe have included the following text in “Discussion” to address this consideration and to acknowledge\nthe limitations of our approach:\n“In this study, our model ranks the factors that significantly influence the prediction accuracy\nof groundwater quality, identifying depth-based groundwater type as the most critical\npredictor (Fig. S3). Sensitivity analysis, involving a ±10% change in each variable,\nhighlighted permeability as the most sensitive factor (Fig. S13).”\n“Additionally, assuming steady-state conditions for environmental variables like soil\nproperties and water table depth may overlook key temporal variations, limiting the capture\nof dynamic contaminant transport over time.”\nIn addition, we have detailed the steps and methods used for the sensitivity analysis in the\n“Assessment of predictor importance” section of “Methods”:\n“To further assess the sensitivity of each influencing factor on the model’s predictions62, we\nconducted a sensitivity analysis by systematically increasing or decreasing each continuous\npredictor by 10%. This approach allowed us to observe how these perturbations influenced\nthe predicted average probability of PGQ over the historical period from 1980 to 2020. By\ncomparing these modified predictions to the baseline scenario, where no changes to the\npredictors were applied, we were able to quantify the degree to which each factor contributed\nto variability in PGQ probability, providing insights into which predictors most significantly\naffect groundwater quality predictions under different conditions.\nReference:\nparameters. Reliability Engineering & System Safety 93, 1563-1573 (2008).\na b\n0.05\n0.00\n-0.05\n-0.10\n-0.15\n-0.20\n-0.25\n1 2 3 4 5 6 7 8 9 1011121314151617181920\nPredictor\nFig. S13| Boxplots showing the difference in average national PGQ probability from 1980 to 2020\ncompared to the baseline scenario under increased (+10%) (a) and decreased (-10%) (b) conditions for\nvarious factors. The x-axis represents the following factors from 1 to 20: (1) porosity, (2) permeability,\n(3) clay content of soil (0-20 cm), (4) CEC (100-200 cm), (5) pH (100-200 cm), (6) sand content of soil\nQGP\nfo\nytilibaborp\nni\negnahC\n25th Percentile ~ 75th Percentile 1.5 IQR range Median Mean\n0.05\n0.00\n-0.05\n-0.10\n-0.15\n-0.20\n-0.25\n1 2 3 4 5 6 7 8 9 1011121314151617181920\nPredictor\nQGP\nfo\nytilibaborp\nni\negnahC\n25th Percentile ~ 75th Percentile 1.5 IQR range Median Mean\n(100-200 cm), (7) NDVI, (8) TWI, (9) WTR, (10) aridity, (11) temperature, (12) precipitation, (13)\ngroundwater supply, (14) groundwater exploitation rate, (15) fertilizer consumption, (16) pesticide\nconsumption, (17) industrial wastewater discharge, (18) industrial solid waste discharge, (19) domestic\nwastewater discharge, (20) domestic solid waste discharge.\n5. Data Set Clarity: The manuscript mentions the number of surveys derived from a certain\nnumber of articles but does not specify the exact number of data points within the groundwater\nquality dataset. For clarity, could you please provide the total number of data points and detail\nhow they are divided into training and validation sets? This information is crucial for\nunderstanding the scale of the analysis and the potential for model overfitting.\nResponse: Thank you for your suggestion. The groundwater quality dataset contains 82,650 data points\ncollected from 1,977 surveys. We allocated 90% of the surveys (1,779 surveys, comprising 67,315 points)\nfor training the RF model, and the remaining 10% (198 surveys, comprising 15,335 points) for validation.\nThis division ensures a balance between providing ample training data and maintaining a robust\nvalidation process, thereby strengthening the credibility and reliability of our model’s predictions.\nAlthough RF models are inherently resistant to overfitting, we also applied cross-validation to rigorously\nassess model performance. Overfitting was not observed during these evaluations. In the revised\nmanuscript, we have clarified this with the following text:\n“To further validate the model’s performance and safeguard against overfitting, we conducted\n10-fold cross-validation. The dataset was randomly partitioned into 10 subsets, with 9 used\nfor training and 1 for validation in each iteration. This process was repeated across all subsets,\nensuring comprehensive evaluations. The averaged performance metrics (Table S4 and Fig.\nS22) confirmed the model’s ability to generalize effectively, without overfitting.\nAfter confirming the model’s robustness through cross-validation, we proceeded to retrain it\nusing 90% of the data (1,779 surveys, comprising 67,315 points) for training and 10% (198\nsurveys, comprising 15,335 points) for validation, each preserving the proportion of PGQ and\nnon-PGQ of the full dataset. This approach maximizes the model’s learning capacity by\nproviding a larger training set while preserving an independent validation set to assess\nperformance on unseen data. This final validation confirmed that the model maintained strong\npredictive accuracy, further enhancing its generalization ability.”\n6. Spatial Partitioning Methodology (L.365-372): The method of spatial partitioning into\n1KM*1KM grids requires further elaboration. Specifically, it is unclear how the presence of both\nPGQ and non-PGQ within a single grid is resolved. Additionally, please elucidate how point-scale\ndriver variables are expanded to fit the 1KM*1KM grid format. A clear explanation of these\nmethodologies will greatly aid in understanding the spatial analysis performed.\nResponse:Thank you for your insightful comments. We have thoroughly reviewed our original dataset\nand confirmed that no instances exist where a single 1 km² grid contains both PGQ and non-PGQ\nclassifications. In the event such a scenario arises, the grid would be classified as PGQ, as the presence\nof groundwater degradation within the grid takes precedence in determining the classification.\nIn the revised manuscript, we have clearly explained the process of partitioning both the predictors and\ntarget variables into the 1 km² grids:\n“To classify groundwater quality at the prefecture level, we implemented a series of data\npreprocessing steps. First, regional data were converted into 1 km² grid cells based on\nregional land-use types, with all cells within a prefecture assigned the corresponding\ngroundwater quality data to ensure consistency. The groundwater quality data were then\nbinarized, assigning a value of zero to non-PGQ and one to PGQ. The purpose of this\napproach was twofold: (i) to prioritize the fundamental health aspect of safety (non-PGQ) or\nunsafety (PGQ) of groundwater for drinking; and (ii) to address variations in precision\nresulting from diverse analysis methods employed in different data sources. Next, natural\nfactors in the predictors were extracted at a 1 km² resolution, assigning each grid cell a unique\nvalue representing local environmental conditions. Socio-economic factors were uniformly\napplied to all grid cells within a province, based on the sampling year. The dataset on land-\nuse types covers specific benchmark years: 1980, 1990, 1995, 2000, 2005, 2010, 2015, 2020,\nand every five years between 2020 and 2100. For other years, land-use data are assumed to\nbe identical to the adjacent prior benchmark year. For instance, land-use data in 1989 are\ntaken to be the same as those in 1980. For other dynamic factors, their cumulative average\nvalues from 1980 to the target year were calculated as the model inputs. This accounts for the\nlag in pollutant migration and effectively captures the long-term impacts on groundwater\nquality. Fertilizer consumption and pesticide consumption were averagely allocated per km2\nof agricultural land, with non-agricultural land assigned zero. Industrial wastewater\ndischarge, industrial solid waste discharge, domestic wastewater discharge, domestic solid\nwaste discharge, and groundwater supply were allocated per km2 based on population density\nfor each province. Per capita data for these variables were calculated by dividing the total\nvalue by the population of each province, then distributed per km2 based on population density.\nGroundwater exploitation rates were uniformly distributed across each province, with regions\nhaving a population density below 1 person/km² assigned a rate of zero. Other factors were\nallocated per km2 based on the resolution of raster data or vector data. Finally, groundwater\nquality categories for each year and region were matched with the corresponding natural and\nsocio-economic variables to create the training and validation datasets for the RF model.”\n7. Model Parameters and Tuning: The manuscript would benefit from a detailed description of the\nrandom forest model parameters, including the settings for each parameter and the methods used\nfor tuning. If parameter optimization was conducted, please describe the process and the rationale\nbehind the chosen methods and settings.\nResponse: Thank you for your comment. In our modeling approach, we focused primarily on optimizing\nthe number of trees. We tested different values, increasing the number of trees in increments of 10,\nranging from 10 to 200. Through this process, we identified 50 trees as the optimal number, delivering\nthe highest model accuracy. To ensure the evaluation’s stability and reliability, we utilized 10-fold cross-\nvalidation, which allowed us to evaluate performance across different data subsets, providing a robust\nassessment and helping to mitigate the risk of overfitting.\nOther hyperparameters of the RF model, such as the minimum leaf size for classification trees (set to 1),\nwere left at their default settings as provided by the MATLAB Statistics and Machine Learning Toolbox.\nThese defaults are widely regarded as effective for many applications, and in some cases, the toolbox\nenforces fixed values to prevent alteration. Given this, we did not find it necessary to adjust these\nparameters further. Additionally, since we did not use a deep RF model, tree depth was not treated as a\ntunable hyperparameter.\nTo enhance clarity, we have incorporated the following texts in the “Splitting the training and validation\nsets” section of “Methods”:\n“We began by optimizing the hyperparameters of the RF model, with a particular focus on the\nnumber of decision trees. By systematically testing tree counts in increments of 10, from 10 to\n200, we identified 50 trees as the optimal balance between accuracy and computational\nefficiency (Fig. S21). Other parameters, such as minimum leaf size, were kept at their default\nsettings, as these have proven effective in similar applications.”\n0.79\n0.78\n0.77\n0.76\n0.75\n0.74\n0.73\n0 20 40 60 80 100 120 140 160 180 200\nFig. S21| Validation accuracy (Acc) of the RF model with different numbers of decision trees.”\n8. Uncertainty Analysis: An uncertainty analysis of the model results would be beneficial.\nProviding an assessment of the model’s uncertainty will give readers a better understanding of the\nreliability of the predictions.\nResponse: To evaluate the uncertainties of our model results, we conducted Monte Carlo simulations\nusing n=50 sets of randomized input data representing various uncertain future scenarios. This approach\nenables us to quantify output uncertainty and highlight the range of possible outcomes and their\nimplications. By incorporating these uncertainty assessments, we offer a clearer understanding of\npotential future groundwater quality under different conditions, providing policymakers with valuable\ninsights into associated risks.\nWe have detailed the steps and methods used for uncertainty analysis in “Methods”:\n“Uncertainty analysis\nTo evaluate the impact of uncertainties in climate change, groundwater exploitation, and\npollution charge on model predictions, we applied Monte Carlo (MC) simulations to quantify\nthe variability in groundwater quality predictions due to input uncertainties76. This analysis\nwas conducted for the period 1980-2100, with data processed at 5-year intervals for 1980-\n2015, annually from 2016 to 2020, and at 5-year intervals for future projections from 2025 to\n2100 under four scenarios.\nIn each scenario-year combination, MC simulations with 50 members were run to introduce\nvariability in key natural and socio-economic predictors. Perturbations were applied to these\ndominant factors, such as climate, groundwater exploitation and pollution charge, reflecting\ntheir real-world uncertainties. For instance, temperature varied within a ±2°C range, while\npredictors like precipitation, fertilizer, and pesticide use were perturbed by up to ±20%.\nThe results of these simulations can be used to evaluate the range of possible groundwater\nquality outcomes based on uncertainties in input data. For each year and scenario, we\ncalculated the mean and standard deviation of the predicted probabilities based on the\nensemble outcomes, providing a robust estimate of both expectation and variability in\ngroundwater quality predictions.”\nIn addition, we have included an uncertainty analysis in “Discussion”:\n“Although visual comparisons and statistical evaluations suggest that the model estimates are\ngenerally reasonable, notable uncertainties remain. Monte Carlo simulations were applied to\nassess uncertainties in PGQ probability predictions, providing insights into spatial and\ntemporal variability. As shown in Figs. S18(d-f) and S19(e-h), uncertainties are moderate\noverall but higher in North and Northeast China from 1980 to 2100, likely due to complex\ninteractions between socio-economic factors and groundwater systems. Under Scenario IV,\nuncertainties become more pronounced (Fig. S20), likely driven by increased discharges and\ngroundwater exploitation. This highlights the need for adaptive groundwater management and\npollution control strategies, particularly in high-emission scenarios, to mitigate adverse\nimpacts and enhance prediction reliability in vulnerable regions.”\nReference:\nEditorial Note: The maps in Fig S18 of this Peer Review File are from Standard Map Service of the\nMinistry of Natural Resources, People’s Republic of China [Map Review Number: GS(2019)1822].\na\nd\nb e\nc f\nMean value of the probability of PGQ Standard deviation of the probability of PGQ\nFig. S18| Spatial distribution of mean and standard deviation of PGQ probability from 1980 to 2020\nbased on Monte Carlo simulations. Subfigures (a-c) show the spatial distribution of the mean of PGQ\nprobability for the years 1980, 2000, and 2020, respectively, and the insets in (a-c) indicate the area\nratios of different PGQ probability levels; Subfigures (d-f) present the standard deviation of PGQ\nprobability, representing the prediction uncertainties for the same years, and the insets in (d-f) indicate\nthe area ratios of different uncertainty of PGQ probability levels.\n0891\n0002\nMinistry of Natural Resources, People’s Republic of China [Map Review Number: GS(2019)1822].\na\ne\nb f\nc g\nd h\nMean value of the probability of PGQ Standard deviation of the probability of PGQ\nFig. S19| Spatial distribution of mean and standard deviation of PGQ probability under four scenarios\nin 2050 based on Monte Carlo simulations. Subfigures (a-d) show the mean of PGQ probability under\nthe four scenarios in 2050, and the insets in (a-d) indicate the area ratios of different PGQ probability\nlevels; Subfigures (e-h) present the uncertainty (standard deviation) of PGQ probability under the four\nscenarios in 2050, and the insets in (e-h) indicate the area ratios of different uncertainty of PGQ\nprobability levels.\nIoiranecS\nIIoiranecS\nIIIoiranecS\nVIoiranecS\na b c\nd e\nOutlier\n1.5 IQR range\n75th Percentile\nMean\nMedian\n25th Percentile\n1.5 IQR range\nOutlier\nFig. S20| Temporal variation of PGQ probability from 1980 to 2100 across different scenarios based on\nMonte Carlo simulations. (a) Probability distribution from 1980 to 2020. (b-e) Probability distributions\nunder Scenarios I, II, III, and IV from 2025 to 2100, respectively.\n9. Data Processing Details: It would be beneficial to include a section on data preprocessing.\nSpecifically, readers would benefit from knowing whether the input driver variable data were\nstandardized or normalized before being fed into the model. Standardization is a common practice\nin machine learning to ensure that each feature contributes proportionately to the final prediction.\nResponse: Thank you for your valuable suggestions. Typically, it is unnecessary to normalize or scale\nthe predictors when using a RF model. This is because RF relies on decision trees, which divide the data\nbased on specific feature thresholds. Each feature is treated independently during the splitting process,\nmeaning the model’s performance is not affected by the scale or distribution of the predictors. In the\nrevised manuscript, we have provided a comprehensive explanation of the data preprocessing steps,\nclearly detailing the workflow from raw data collection through to the generation of training samples.\n“To classify groundwater quality at the prefecture level, we implemented a series of data\npreprocessing steps. First, regional data were converted into 1 km² grid cells based on\nregional land-use types, with all cells within a prefecture assigned the corresponding\ngroundwater quality data to ensure consistency. The groundwater quality data were then\nbinarized, assigning a value of zero to non-PGQ and one to PGQ. The purpose of this\napproach was twofold: (i) to prioritize the fundamental health aspect of safety (non-PGQ) or\nunsafety (PGQ) of groundwater for drinking; and (ii) to address variations in precision\nresulting from diverse analysis methods employed in different data sources. Next, natural\nfactors in the predictors were extracted at a 1 km² resolution, assigning each grid cell a unique\nvalue representing local environmental conditions. Socio-economic factors were uniformly\napplied to all grid cells within a province, based on the sampling year. The dataset on land-\nuse types covers specific benchmark years: 1980, 1990, 1995, 2000, 2005, 2010, 2015, 2020,\nand every five years between 2020 and 2100. For other years, land-use data are assumed to\nbe identical to the adjacent prior benchmark year. For instance, land-use data in 1989 are\ntaken to be the same as those in 1980. For other dynamic factors, their cumulative average\nvalues from 1980 to the target year were calculated as the model inputs. This accounts for the\nlag in pollutant migration and effectively captures the long-term impacts on groundwater\nquality. Fertilizer consumption and pesticide consumption were averagely allocated per km2\nof agricultural land, with non-agricultural land assigned zero. Industrial wastewater\ndischarge, industrial solid waste discharge, domestic wastewater discharge, domestic solid\nwaste discharge, and groundwater supply were allocated per km2 based on population density\nfor each province. Per capita data for these variables were calculated by dividing the total\nvalue by the population of each province, then distributed per km2 based on population density.\nGroundwater exploitation rates were uniformly distributed across each province, with regions\nhaving a population density below 1 person/km² assigned a rate of zero. Other factors were\nallocated per km2 based on the resolution of raster data or vector data. Finally, groundwater\nquality categories for each year and region were matched with the corresponding natural and\nsocio-economic variables to create the training and validation datasets for the RF model.”\n10. Rationale for future scenario settings (L.497-514): The setting of future scenarios is crucial for\nthe interpretation of model results. Given that variables such as fertilizer and pesticide usage\nappear to have a significant impact on the model, as indicated by Fig. S3 and the results section,\nthe rationale behind the high (3%), medium (1%), and low (-3%) annual change rates should be\ndiscussed. This discussion could include the basis for these specific rates and how they align with\nprojected socio-economic and environmental trends.\nResponse:Thank you for your suggestion. The high (3%), medium (1%), and low (-3%) annual change\nrates for groundwater exploitation and pollution discharge were established based on historical trends\nand the Shared Socioeconomic Pathways (SSPs) framework. These rates were selected to strike a balance\nbetween past trends and future projections, taking into account key factors such as industrial activities,\nagricultural expansion, and environmental policies, as emphasized in relevant studies and SSP\nassessments. In the revised manuscript, we have discussed the basis and limitations of our approach:\nIn the three scenarios, annual change rates of -3%, 1%, and 3% were assigned as gradients\nfor groundwater exploitation and pollution discharge, respectively. These rates were\ndetermined based on historical trends69,70 and future projections71,72 related to industrial\nactivities, agricultural expansion, and environmental policies. The low rate of -3% reflects a\nproactive approach toward reducing resource use and pollution discharge, aligning with the\nsustainability-focused SSP1-1.9 scenario73. In contrast, the medium rate of 1% corresponds\nto a moderate development pathway, as seen in SSP2-4.5. This scenario embodies a balanced\napproach, allowing socio-economic development to progress while placing emphasis on\nenvironmental protection74. The high rate of 3% aligns with the SSP5-8.5 scenario,\ncharacterized by rapid economic growth and minimal environmental regulation75. This\nscenario anticipates high energy demands, with significant reliance on fossil fuels and\nindustrial expansion. We acknowledge that future development may exhibit more intricate\npatterns, and our current study may not fully capture this complexity and uncertainty. Then,\nthese rates were applied to the corresponding values from 2021 to 2100, and the cumulative\naverage was calculated as the model inputs in each year. Please refer to Table S5 for detailed\ninformation on the specific scenario assumptions.”\nReference:\nPress, 1985-2021).\nStatistics Press, 1981-2021).\ngas emissions implications: An overview. Global Environmental Change42, 153-168 (2017).\nworld futures in the 21st century. Global Environmental Change42, 169-180 (2017).\ngreen growth paradigm. Global Environmental Change42, 237-250 (2017).\nof-the-road scenario for the 21st century. Global Environmental Change42, 251-267 (2017).\nfor the 21st century. Global Environmental Change42, 297-315 (2017).\nAddressing these points will significantly improve the manuscript and provide readers with a clearer\nunderstanding of your research methodology and findings.\nReviewer #3 (Remarks to the Author):\nBased on a lot of collected data, the manuscript used the random forest model to construct the annual\nprobabilistic map of poor groundwater quality across China from 1980 to 2100. Although the method is\nnot new, the idea they use it to model the groundwater quality is good. The new insights from the study\nare helpful for the policy-makers to take reasonable measures to protect China’s groundwater. However,\nthere are some issues that need further addressing:\n1. Title of the paper didn’t fully reflect the contents; I think it should change.\nResponse:We appreciate your feedback on the title of the paper. To more accurately capture the scope\nand content of the study, we propose revising the title to: “Groundwater Quality Evolution Across\nChina”.\n2. Climate includes temperature, precipitation, evaporation etc., they usually have different\nimpacts to groundwater quality. For example, precipitation is usually beneficial for groundwater\nquality, while evaporation is detrimental to groundwater quality. What does climate mean in your\npaper?\nResponse: Thank you for your feedback. As outlined in Table S2, the term “climate” encompasses\nseveral key variables, including temperature, precipitation, aridity, actual evapotranspiration (ET),\nthe\nPriestley-Taylor alpha coefficient, humidity, and particulate matter 10 micrometers or less in diameter\n(PM ).\n10\n“Table S2 Predictors for the groundwater quality classification model. The 25 variables ultimately used\nin the final prediction model are indicated in bold.(excerpt)”\nResolution\nClimate Period\nTemporal Spatial\nTemperature9,10 yearly 1 km 1980-2100\nPrecipitation11,12 yearly 1 km 1980-2100\nActual evapotranspiration (ET)13 steady-state 1 km /\nAridity (Potential ET / Precipitation)14 steady-state 1 km /\nPriestley‐Taylor alpha coefficient (Actual ET / Potential ET)13 steady-state 1 km /\nHumidity3 steady-state 1 km /\nParticulate matter 10 micrometers or less in diameter (PM )3 steady-state 1 km /\nprocesses such as recharge, contaminant dilution, and pollutant concentration. Moreover, they interact\nwith a range of natural and social factors. In the revised manuscript, we have added a detailed discussion\non the influence of these various factors on groundwater quality:\n“In this study, our model ranks the factors that significantly influence the prediction accuracy of\ngroundwater quality, identifying depth-based groundwater type as the most critical predictor (Fig. S3).\nSensitivity analysis, involving a ±10% change in each variable, highlighted permeability as the most\nsensitive factor (Fig. S13). In addition, we analyzed the interrelationships among the predictors (Fig.\nS2), as well as the partial dependence plots and feature interaction charts (Figs. S14 and S15),\ndemonstrating the complex interactions among the predictors. Shallow groundwater is more vulnerable\nto surface contaminants, while deeper groundwater is less impacted27. Regions with lower precipitation,\nhigher temperatures, and higher groundwater exploitation rates showed higher PGQ probabilities (Fig.\nS15). Climate-driven changes directly affect recharge rates, influence groundwater replenishment, while\nindirectly alter groundwater usage patterns28. In arid and semi-arid regions, reduced rainfall leads to\nover-extraction and higher pollutant concentrations due to decreased dilution, whereas humid regions\nsee increased mobilization and dilution of pollutants29-30. As climate change continues to reshape global\ntemperature and precipitation patterns, these dynamics are expected to become increasingly complex.\nAdditionally, soil sand content (100-200 cm) significantly affects recharge and contamination, with\nhigher content increasing susceptibility to surface contaminants (Fig. S14f)31. It is essential to note that\nthis analysis is based on a nationwide perspective, and the most critical predictors vary spatially, as\nshown in Fig. S7a. These insights highlight the importance of informing policymakers and stakeholders\nabout region-specific factors that must be managed to improve groundwater quality.”\nReference:\nIllinois. Vol. 532 (Illinois State Geological Survey, 1984).\nHydrology518, 250-266 (2014).\nEnvironmental Science and Technology51, 2895-2933 (2021).\ngroundwater in China: A review. Science of the Total Environment643, 967-993 (2018).\nJournal of Hydrology332, 187-205 (2007).\na 1 b1\n2\n0.8 6\n12\n15 0.6 17\n19\n0.4 20\n22\n24\n0.2 25\n26\n0 27\n28\n29\n-0.2 30\n35\n-0.4 36 38\n39\n-0.6 41\n42\n43\n-0.8 44\n51\n-1 52\nFig. S2| Pearson correlation analysis of predictors (a) before and (b) after variable pruning. The\nvariable pruning is carried out to remove redundant variables to prevent overfitting. The numbers in the\nx and y axis indicate different variables: 1. Porosity; 2. Permeability; 3. Parent material; 4. Sand content\nof soil (0-20cm); 5. Silt content of soil (0-20cm); 6. Clay content of soil (0-20cm); 7. pH (0-5cm); 8. Bulk\ndensity (BD) (0-5cm); 9. Cation exchange capacity (CEC) (0-5cm); 10. Soil and sedimentary deposit\nthickness; 11. Soil organic carbon (SOC) (0-100cm); 12. Sand content of soil (100-200cm); 13. Silt\ncontent of soil (100-200cm); 14. Clay content of soil (100-200cm); 15. pH (100-200cm); 16. BD (100-\n200cm); 17. CEC (100-200cm); 18. Coarse fragment volume (CFVO) (100-200cm); 19. Soil type; 20.\nSeawater intrusion; 21. Digital elevation model (DEM); 22. Normalized difference vegetation index\n1 2 3 4 5 6 7 8 9 01 11 21 31 41 51 61 7811 91 02 12 22 32 42 52 62 72 82 92 03 13 23 33 43 53 63 73 83 93 04 14 24 34 44 54 64 74 84 94 05 15 25 1 2 6 21 51 71 91 02 22 42 52 62 72 82 92 03 53 63 83 93 14 24 34 44 15 25\n1 1\n1 111111111 0 123456789\n23456789\n0\n0\n0\n.\n.\n.\n4\n6\n8\n2222222222333 0123456789012\n-\n0\n0\n0\n.2\n.2\n33333334 34567890\n-0.4\n444444444 123456789\n-\n-\n0\n0\n.\n.\n8\n6\n555 012\n-1\n(NDVI); 23. Slope; 24. Topographic wetness index (TWI); 25. Water-table ratio (WTR); 26. Aquifer-\nmedia-based groundwater type; 27. Depth-based groundwater type; 28. Precipitation; 29. Temperature;\n30. Aridity; 31. Actual evapotranspiration (ET); 32. Priestley-Taylor alpha coefficient; 33. Humidity; 34.\nParticulate matter 10 micrometers or less in diameter (PM ); 35. Fertilizer consumption; 36. Pesticide\n10\nconsumption; 37. Coal consumption; 38. Industrial wastewater discharge; 39. Industrial solid waste\ndischarge; 40. Industrial gross domestic product (GDP); 41. Domestic wastewater discharge; 42.\nDomestic solid waste discharge; 43. Groundwater supply; 44. Groundwater exploitation rate; 45. Crude\noil consumption; 46. Natural gas consumption; 47. Non-ferrous metal mining; 48. GDP per capita; 49.\nEnvironmental investment; 50. Total wastewater discharge; 51. Land-use type; and 52. Groundwater\nquality category.\nFig. S3| Importance of the predictors in the final random forest (RF) model based on out-of-bag error\nanalysis.\na b\n0.05\n0.00\n-0.05\n-0.10\n-0.15\n-0.20\n-0.25\n1 2 3 4 5 6 7 8 9 1011121314151617181920\nPredictor\nFig. S13| Boxplots showing the difference in average national PGQ probability from 1980 to 2020\ncompared to the baseline scenario under increased (+10%) (a) and decreased (-10%) (b) conditions for\nQGP\nfo\nytilibaborp\nni\negnahC\n25th Percentile ~ 75th Percentile 1.5 IQR range Median Mean\n0.05\n0.00\n-0.05\n-0.10\n-0.15\n-0.20\n-0.25\n1 2 3 4 5 6 7 8 9 1011121314151617181920\nPredictor\nQGP\nfo\nytilibaborp\nni\negnahC\n25th Percentile ~ 75th Percentile 1.5 IQR range Median Mean\nvarious factors. The x-axis represents the following factors from 1 to 20: (1) porosity, (2) permeability,\n(3) clay content of soil (0-20 cm), (4) CEC (100-200 cm), (5) pH (100-200 cm), (6) sand content of soil\n(100-200 cm), (7) NDVI, (8) TWI, (9) WTR, (10) aridity, (11) temperature, (12) precipitation, (13)\ngroundwater supply, (14) groundwater exploitation rate, (15) fertilizer consumption, (16) pesticide\nconsumption, (17) industrial wastewater discharge, (18) industrial solid waste discharge, (19) domestic\nwastewater discharge, (20) domestic solid waste discharge.\na b c d\ne f g h\ni j k l\nm n o p\nq r s t\nFig. S14| Partial dependence plots for PGQ probability across multiple features. Subfigures (a-t) depict\nthe partial dependence of PGQ probability on different factors: (a) porosity, (b) permeability, (c) clay\ncontent of soil (0-20cm), (d) CEC (100-200cm), (e) pH (100-200cm), (f) sand content of soil (100-200cm),\n(g) NDVI, (h) TWI, (i) WTR, (j) aridity, (k) temperature, (l) precipitation, (m) groundwater supply, (n)\ngroundwater exploitation rate, (o) fertilizer consumption, (p) pesticide consumption, (q) industrial\nwastewater discharge, (r) industrial solid waste discharge, (s) domestic wastewater discharge, and (t)\ndomestic solid waste discharge. Each subfigure shows the relationship between the feature and the\npredicted PGQ probability, indicating how changes in each factor influence groundwater quality\nprobabilities.\na b\nFig. S15| Feature interaction charts of the RF model. Subfigure (a) shows the combined effects of\nprecipitation and temperature, and (b) of precipitation and groundwater exploitation rate, on the\npredicted probability of PGQ. Both variables are evaluated across their full empirical ranges and are\nallowed to vary simultaneously. Model predictions are averaged over the joint distribution of these\nvariables. The resulting interactions are displayed on a two-dimensional, color-coded map, where the\ngradient from blue to red visually represents the shift from lower to higher probabilities of PGQ.\n3. Land use changes have different types, for example, from forest to agricultural land, or from\nagricultural land to construction land, etc. What do you mean in your study?\nResponse:Thank you for your comment. In our study, we utilize land use types—agricultural land, forest,\ngrassland, water body, built-up land, and unused land—at various times and locations as predictors (Table\nS2). Land use changes are identified by comparing these land use types across different years.\n“Table S2 Predictors for the groundwater quality classification model. The 25 variables ultimately used\nin the final prediction model are indicated in bold.(excerpt)\nResolution\nLand use Period\nTemporal Spatial\n1980, 1990, 1995, 2000, 2005, 2010, 2015,\nLand-use type23,24 yearly 1 km 2020, and every five years between 2020\nand 2100\nNote: The population density data was not utilized directly as a predictive variable. Instead, it was used\nto derive the spatial distribution of groundwater exploitation and pollution discharge, which we only\nhave data at provincial scale. Specifically, we downscaled these data by assuming they are proportional\nto population density within each province. Aquifer-media-based groundwater types include loose\nsediment pore water, carbonate rock fissure cavern water, clastic rock void fissure water, magmatic rock\nmetamorphic rock fissure water, and frozen layer water. Groundwater can also be classified by depth into\nphreatic, confined, and unspecified types. Additionally, land-use types encompass agricultural land,\nforest, grassland, water body, built-up land, and unused land.”\n4. How the groundwater exploitation affects its quality, the links between these should be clarified.\nResponse: Thank you for your comment. In the revised manuscript, we have addressed this point with\nthe following explanation:\n“Notably, groundwater exploitation (e.g., groundwater supply and groundwater exploitation\nrate) emerges as a crucial anthropogenic perturbation factor53,54. Excessive groundwater\nextraction lowers water tables, concentrates contaminants, and, in coastal areas, can trigger\nseawater intrusion, increasing salinity34,55,56. It can also mobilize naturally occurring\ncontaminants like arsenic and fluoride, further degrading water quality29.”\nReference:\nEnvironmental Science and Technology51, 2895-2933 (2021).\nreaches deeper depths in heavily pumped aquifer systems. Nature Communications13, 5263 (2022).\nScience & Technology57, 5125-5136 (2023).\n7, 413-422 (2024).\n625, 715-721 (2024).\neastern coastal of Laizhou Bay, China. Environmental Geology54, 1763-1770 (2008).\n5. Whether it is agricultural or industrial pollution, there is usually a lag time for it to migrate to\ngroundwater. The time can be several years or tens of years. How do you treat this in you study?\nResponse:Thank you for your insightful comment. In our study, we calculate cumulative average values\nof dynamic pollution factors starting from 1980, leveraging long-term historical data to account for the\ndelayed migration of pollutants into groundwater. This approach enables our model to effectively capture\nchanges in groundwater quality over timescales of several years to decades.\nIn the revised manuscript, we have clarified this process:\n“For other dynamic factors, their cumulative average values from 1980 to the target year were\ncalculated as the model inputs. This accounts for the lag in pollutant migration and effectively\ncaptures the long-term impacts on groundwater quality.”\n6. This study should consider the influence of geological background on water quality.\nResponse: Thank you for highlighting the importance of geological background in groundwater quality\nassessment. As reflected in Table S2, geological factors have been incorporated into our modeling\napproach, allowing us to account for the diverse influences of geological background on groundwater\nquality across different regions.\n“Table S2 Predictors for the groundwater quality classification model. The 25 variables ultimately used\nin the final prediction model are indicated in bold.(excerpt)\nResolution\nGeography and hydrogeology Period\nTemporal Spatial\ncoastal buffer zones\nSeawater intrusion steady-state /\n(distance≤ 10 km)\nDEM10 steady-state 90 m /\nNDVI10 steady-state 500 m /\nSlope11 steady-state 1 km /\nTWI11 steady-state 1 km /\nWTR12 steady-state 1 km /\nAquifer-media-based groundwater type13 steady-state 1:4 000 000 /\nDepth-based groundwater type steady-state as recorded /\nNote: The population density data was not utilized directly as a predictive variable. Instead, it was used\nto derive the spatial distribution of groundwater exploitation and pollution discharge, which we only\nhave data at provincial scale. Specifically, we downscaled these data by assuming they are proportional\nto population density within each province. Aquifer-media-based groundwater types include loose\nsediment pore water, carbonate rock fissure cavern water, clastic rock void fissure water, magmatic rock\nmetamorphic rock fissure water, and frozen layer water. Groundwater can also be classified by depth into\nphreatic, confined, and unspecified types. Additionally, land-use types encompass agricultural land,\nforest, grassland, water body, built-up land, and unused land.”\n7. Is a kappa coefficient of 0.69 a bit low? Generally, we believe that a kappa coefficient higher\nthan 0.8 is best.\nResponse: Thank you for your insightful comment. While the kappa coefficient of 0.69 in our analysis\nfalls short of the ideal 0.80, it still reflects substantial agreement according to the standards set by\nLandis and Koch (1977)25, who classified kappa values between 0.61 and 0.80 as indicating substantial\nagreement. This is particularly significant given the inherent complexity and variability of environmental\ndata, which can affect classification consistency across different samples or contexts. As such, achieving\na kappa value above 0.8 is challenging in our study. In the revised manuscript, we have included the\nfollowing clarification:\n“Additionally, the model achieves substantial agreement beyond chance in its classifications\nwith Cohen’s kappa coefficient of 0.69. While slightly below the ideal, this value falls within\n0.61-0.80, which Landis and Koch (1977) defined as substantial agreement25. This slightly\nlower kappa value reflects the inherent complexity and variability of the environmental data\nand processes analyzed in this study.”\nReference:\nof majority agreement among multiple observers. Biometrics, 363-374 (1977).\n8. In table S1, ug/L should be μg/L.\nResponse: Thank you for pointing out this typo. We have corrected the notation from “ug/L” to “μg/L”\nthroughout the manuscript and the supporting information.\n9. In the assessment of groundwater quality in lines 335-336, does the use of a single factor\nassessment method make the results too pessimistic? Please explain.\nResponse: Thank you for your feedback. In our study, we adhere to the latest Chinese groundwater\nquality standard (GB/T 14848-2017), which classifies overall groundwater quality based on the worst\nindicator in the dataset. This approach closely aligns with the barrel theory, which asserts that the overall\nperformance or effectiveness of a system is limited by its weakest component, similar to how a barrel’s\ncapacity is determined by the height of its shortest stave. We acknowledge that this approach may lead\nto more stringent assessments; however, it is crucial to ensure that all potential risks to groundwater\nquality are thoroughly considered. Furthermore, by aligning with established standards, our results can\nprovide more targeted information for government agencies.\nIn the revised manuscript, we have added the following clarification to “Discussion”:\n“Our single-factor assessment focuses on the worst indicator, which may lead to stricter\nevaluations but offers a cautious approach to groundwater management.”\n10. Since there is a lot of missing data before 1990, why not start the study directly from 1990?\nResponse: Thank you for your suggestion. It is indeed true that there is a significant amount of missing\ndata prior to 1990. Similar to our predictions for future groundwater quality, this study also aims to\nreconstruct groundwater quality for the years preceding the range of our collected data. We have reliable\npredictor variables available for the decade from 1980 to 1990, which makes this reconstruction feasible.\nBy including this earlier period, we can capture the historical evolution of groundwater conditions,\nproviding a more comprehensive assessment of water management policies over an extended timeline\nand enhancing our understanding of long-term shifts in groundwater quality.\nIn the revised manuscript, we have added the following clarification in the “Calculation of the Affected\nPopulation” section of “Methods”:\n“The RF model was used to reconstruct and predict the spatiotemporal distribution of\ngroundwater quality from 1980 to 2100, enabling us to infer long-term trends in populations\npotentially exposed to PGQ. We apply a probability cutoff (Prob ) of 0.5 to identify the\ncut\npopulations at risk in areas exceeding this threshold. Then affected population living in each\ngrid (1 km2) was calculated by multiplying the grid’s population (Pop) by its probability of\nPGQ (Prob ).”\npoor\n11. The ground water Standard they listed in Table S1 is 2017, but the reference they cited is 1993.\nResponse: Thank you for pointing out this error. We have corrected the mistake to accurately reflect the\n2017 standards. We appreciate your attention to detail.\nReference:\nSAC. Standard for Groundwater Quality (GB/T 14848-2017). (General Administration of Quality\nSupervision, Inspection and Quarantine of the People’s Republic of China; Standardization\nAdministration of China, Beijing, China, 2017).\n12. Line 337, when we evaluate the water quality, we usually think that the water above Class III\nis beyond the standard. What is the basis for choosing Class V water as the boundary in this\nmanuscript?\nResponse: Thank you for your comment. We classify groundwater quality as poor or not poor using\nClass V as the threshold, determining whether the water is suitable for drinking, with or without treatment.\nAccording to the latest groundwater quality standards (GB/T 14848-2017), water in Classes I to IV is\nconsidered fit for domestic drinking purposes, possibly with some treatment, whereas Class V water is\ndeemed unsuitable for such use. We recognize that applying Class III as a stricter boundary would\nhighlight more critical groundwater quality challenges. In the revised manuscript, we have added the\nfollowing clarification:\n“Moreover, the use of Class V as a threshold for PGQ is conservative, and stricter thresholds\nlike Class IV could be considered.”\n13. The spatial pattern of future land use is affected by both human activities and climate change.\nIn your article, you assume that the land use pattern of different scenarios is in 2020. Is this\nassumption too arbitrary? Different scenarios of future land use are also available in the Qinghai-\nTibet Plateau database. Please explain why you don’t use them?\nResponse: Thank you for your valuable suggestion. In the revised manuscript, we have enhanced our\nmodels to incorporate dynamic land use changes from 2020 to 2100 across various scenarios. This\nimprovement allows for more reasonable projections of groundwater quality changes in response to shifts\nin land use patterns driven by human activities and climate change. To evaluate the impact of land use\nchange on groundwater quality, we have intentionally retained the 2020 baseline for Scenario Ⅰ, which\nreflects our original results without accounting for future land use changes. This approach provides a\nclear contrast between static and dynamic land use scenarios, enhancing both the interpretability and\npractical relevance of our findings.\nIn the revised manuscript, we have included the following updates in the “Scenario settings” section and\nthe “Analysis of drivers” section of “Methods”,and revised the results for future scenarios accordingly:\nScenario settings:\n“In Scenario I (Baseline), we assumed a consistent development intensity for all drivers.\nClimate, groundwater exploitation, and pollution discharges were maintained at the same\nlevels as the cumulative average values in 2020. Similarly, land-use type in Scenario I was\nalso assumed to remain unchanged from 2020. In contrast, Scenarios II-IV incorporated\nprojections to account for future variability. Specifically, climate data from one CMIP6\nmodel——EC-Earth3 were used as model input to project future PGQ probabilities65-67.\nAdditionally, land-use data in these scenarios were derived from the SSP-RCP global 1km\nland-use simulation dataset (2020-2100) for SSP1-1.9, SSP2-4.5, and SSP5-8.5 scenarios,\nupdated every five years68. To maintain consistency in representation, these land-use data were\naligned with those used for the historical period from 1980 to 2020.\nAnalysis of drivers:\n“In future scenarios, the analysis of dominant driver focuses exclusively on four primary\ncategories: climate change, groundwater exploitation, land use and pollution discharges,\nwithout individually detailing secondary pollution categories. We kept one specific driver\ncategory (e.g., climate) constant at its states in 2020, while allowing the other driver\ncategories to vary with time.”\nReference:\n(2022).\n2100. National Tibetan Plateau / Third Pole Environment Data Center (2022).\n2100. National Tibetan Plateau / Third Pole Environment Data Center (2022).\nChina. (http://www.geodata.cn).\n14. The research conclusion should be summarized.\nResponse: Thank you for your suggestion. As the journal does not support including a “conclusion”\nsection in the article, we have incorporated a summary of the key findings and implications of our study\nwithin the “Discussion” section, as follows:\n“Despite the inherent limitations, our approach offers distinct advantages over traditional\nmethods, which are often labor-intensive and costly. By integrating regional data from\nprevious research with accessible natural and socio-economic predictors, our methodology\npresents a comprehensive analytical framework for assessing groundwater qualities. This\nstrategy provides detailed patterns of groundwater quality evolution across China from 1980\nto 2100, revealing that pollution discharge, groundwater over-exploitation, and land-use\nchange are significant contributors to its degradation, particularly in high-risk regions\naffecting large populations. Furthermore, our projections indicate that proactive\nenvironmental policies can substantially improve groundwater quality under diverse future\nscenarios. This approach is especially valuable for developing countries, where limited\ngroundwater data and PGQ are common. By offering a cost-effective solution, our model not\nonly fills critical data gaps but also provides actionable insights into resource management\nand strategic planning. Although uncertainties remain—especially regarding socio-economic\nmodels and emissions paths—the findings underscore the need for adaptive, region-specific\ngroundwater management practices to mitigate socio-economic and climate-related impacts.\nFuture research can build on this model by refining data resolution and examining specific\npollutant pathways to enhance predictive reliability and support long-term, sustainable\ngroundwater policies.”"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I am satisfied with the revised version of the manuscript, as the authors have thoroughly addressed all my comments and suggestions. The revisions have significantly improved the clarity and quality of the work. Therefore, I recommend the manuscript for publication.", "reviewer_2": "The authors have done a commendable job throughly addressing my comments from the first round of review. I am satisfied with their detailed responses and the corresponding revisions made to the manuscript. I recommend the paper for publication.", "reviewer_3": "The authors comprehensively revised the manuscript. I have no further suggestions."}, "author_response": "Reviewer #1 (Remarks to the Author):\nI am satisfied with the revised version of the manuscript, as the authors have thoroughly addressed all\nmy comments and suggestions. The revisions have significantly improved the clarity and quality of the\nwork. Therefore, I recommend the manuscript for publication.\nResponse: Thank you for your positive feedback and for recognizing the efforts we have made to\naddress your comments and suggestions. We are pleased to hear that the revisions have improved the\nclarity and quality of the manuscript. We appreciate your recommendation for publication.\nReviewer #2 (Remarks to the Author):\nThe authors have done a commendable job thoroughly addressing my comments from the first round of\nreview. I am satisfied with their detailed responses and the corresponding revisions made to the\nmanuscript. I recommend the paper for publication.\nResponse: We are grateful to you for your constructive comments during the first round of review and\nfor acknowledging the revisions made to the manuscript. We are glad that our responses and revisions\nmeet your expectations. Thank you for recommending our paper for publication.\nReviewer #3 (Remarks to the Author):\nThe authors comprehensively revised the manuscript. I have no further suggestions.\nResponse: Thank you for your thorough review and for confirming that the manuscript has been\ncomprehensively revised. We appreciate your time and effort in evaluating our work and are pleased to\nhear that no further suggestions are needed."}]

Decision letter