A new AI paradigm for groundwater level prediction: integrating the AlphaEarth foundations model in the Leizhou Peninsula

Abstract

Accurate spatial prediction of groundwater level is crucial for sustainable resource management but remains challenging in coastal regions like the Leizhou Peninsula due to sparse monitoring wells and the complexity of integrating various environmental data. Traditional models often fail to capture nonlinear interactions or generalize effectively under these data-constrained conditions. To address these limitations, this study proposes a novel framework that integrates the AlphaEarth Foundation (AEF) large model with the Tabular Prior-Data Fitted Network (TabPFN) deep learning architecture. We utilize the shared embedding data of the AEF to fuse various earth observation data including optical imagery, radar, and climate simulations into unified 64-dimensional features. Using monitoring data from 2023, six spatial prediction models were systematically evaluated and optimized through Bayesian frameworks to map groundwater depth (GWD). The results demonstrate that the proposed hybrid model achieved superior performance with a coefficient of determination (R²) of 0.793 and a root mean square error (RMSE) of 2.079 m. In addition, ablation experiments, 95% confidence intervals, and prediction bias analysis reveal that it outperformed traditional methods such as XGBoost and Random Forest. Furthermore, global sensitivity and interpretability analyses revealed that AEF embeddings effectively encoded key topographical and climatic drivers. The generated 90 m resolution GWD map unveiled distinct spatial heterogeneity across the peninsula. However, given the inherent uncertainties of small sample sizes, this framework is best utilized as a cost-effective, macro-level preliminary screening tool to guide future monitoring network densification and regional water management, rather than an absolute substitute for dense field monitoring.

1 Introduction

Groundwater is one of the most important freshwater resources worldwide, accounting for approximately 30% of the total global freshwater reserves, and it plays an irreplaceable role in ensuring drinking water security, agricultural irrigation, and industrial production (Oki and Kanae, 2006; McCartney et al., 2022). In particular, under conditions of uneven spatiotemporal distribution of surface water resources and the increasing frequency of extreme climate events, groundwater has become a critical water source for sustaining socioeconomic development in many regions owing to its high stability and strong buffering capacity (Jódar et al., 2024). However, with rapid economic growth and accelerated urbanization, groundwater resources worldwide are increasingly threatened by overexploitation, continuous decline in groundwater level, and water quality deterioration, posing serious risks to regional ecological security and the sustainable use of water resources (Aeschbach-Hertig and Gleeson, 2012; Wang et al., 2022; Noori et al., 2023; Igwebuike et al., 2024; Jasechko et al., 2024; Kim et al., 2024). Consequently, accurately characterizing the spatiotemporal distribution of groundwater level is a fundamental prerequisite for rational groundwater allocation and scientific management.

Groundwater level monitoring and prediction are of great importance for different types of aquifers globally, with particular urgency in coastal regions. Coastal groundwater systems are influenced not only by natural factors such as precipitation, evapotranspiration, and topography, but also by multiple anthropogenic and marine-related disturbances, including seawater intrusion, intensive groundwater extraction, and land-use change (Chihi and Ben Cheikh Larbi, 2023; Richardson et al., 2024). As a result, groundwater level dynamics in coastal zones often exhibit pronounced nonlinearity and strong spatial heterogeneity. The Leizhou Peninsula, located along the southern coast of China, is characterized by relatively abundant groundwater resources and has long played a vital role in regional agricultural irrigation, domestic water supply, and industrial development (Chen et al., 2023). Nevertheless, due to large-scale historical groundwater exploitation, the vulnerability of the regional groundwater system has gradually emerged, manifested by intensified groundwater level fluctuations and the accumulation of ecological risks in certain areas. Therefore, conducting spatial prediction of groundwater level in the Leizhou Peninsula is not only of strong regional representativeness but also of substantial practical significance for groundwater management in China’s ecologically fragile coastal zones.

In recent years, with rapid advances in remote sensing technology, geographic information systems, and computational capacity, groundwater level prediction based on machine learning and deep learning has been widely applied (Zhou et al., 2018; Sun et al., 2022; Boo et al., 2024; Zhou and Zuo, 2025). Existing studies typically integrate multi-source environmental factors, such as digital elevation models (DEM), meteorological variables (e.g., precipitation, temperature, and evapotranspiration), land-use types, and vegetation indices. Data-driven models, including Random forest, Support vector machine, Gradient boosting trees, and Neural networks, were used to predict groundwater level, (Gharehbaghi et al., 2022; Tao et al., 2022; Li et al., 2023; Zhou and Xiao, 2024). These approaches have partially overcome the strict dependence of traditional physical models on parameters and boundary conditions and have demonstrated strong capabilities in capturing complex nonlinear relationships. However, several critical challenges remain. First, inconsistencies in spatial resolution, temporal scale, and data structure among multi-source heterogeneous datasets substantially reduce data integration efficiency. Second, in regions with sparse groundwater monitoring wells, model performance is highly sensitive to sample size, leading to limited predictive accuracy and generalization ability. Third, most data-driven models rely heavily on manual feature selection and construction, making it difficult to systematically capture the multi-scale coupled information embedded within groundwater systems. Collectively, these limitations severely constrain large-scale, high-resolution groundwater level prediction.

To tackle the challenges in data-scarce regions, various systematic methodological frameworks have emerged. For instance, Borzì, (2025) categorized the research on groundwater modeling in data-scarce regions into ten distinct emerging methodological approaches, which include physical groundwater models, machine learning techniques, geostatistical methods, remote sensing inversion techniques, kriging-based interpolation methods, and others. These studies have provided a clear methodological reference for relevant research. Geostatistical methods are the traditional mainstream techniques, whose strengths lie in enabling regional interpolation by leveraging the spatial autocorrelation of samples and quantifying prediction uncertainty. However, such methods are highly sensitive to the spatial uniformity of sample distribution. The interpolation accuracy will declines significantly when monitoring wells are extremely sparse or exhibit a clustered distribution (Hasan et al., 2021). The virtual sample generation method is often used for modeling groundwater pollutants in data-scarce regions. This method mitigates the limitations of small sample sizes through data augmentation, yet it is difficult to strictly ensure the physical consistency of the generated samples. This approach is prone to introducing spurious information, thereby leading to prediction biases (Pan et al., 2025). Machine learning techniques have been widely used in groundwater modeling due to their powerful nonlinear fitting ability (Rafik et al., 2023; Yang et al., 2025). Nevertheless, traditional machine learning methods generally depend on manual feature engineering, which requires researchers to manually integrate multi-source data including topography, climate and underlying surface conditions. This process is not only time-consuming and labor-intensive, but also likely to result in the omission of key information due to the subjectivity in feature design. Remote sensing inversion techniques and integrated modeling frameworks achieve large-scale spatial coverage by fusing satellite observation data, which effectively compensates for the inadequacy of ground monitoring data (Hasan et al., 2019). However, remote sensing data has limitations such as insufficient accuracy at the local scale and the necessity of validation with ground-based data.

To address these challenges, the recently released AEF large model by Google provides a novel technical pathway for the unified representation and efficient utilization of multi-source Earth observation data. AEF embeds various Earth observation datasets into a unified 64-dimensional shared feature space. These datasets include optical remote sensing, synthetic aperture radar (SAR), DEM, and climate simulations, which enable deep spatial and semantic integration across different data sources (Brown et al., 2025; Tollefson, 2025). This framework effectively alleviates the long-standing issues of “data silos” and spatiotemporal misalignment commonly encountered in traditional multi-source data approaches. As these embedding features comprehensively reflect surface morphology, hydroclimatic conditions, and human activity intensity, they possess considerable potential to characterize groundwater level variations (Houriez et al., 2025). Therefore, constructing groundwater level prediction models based on AEF embedding features is expected to enable robust and efficient spatial prediction, particularly in data-sparse regions.

The AEF combined with Tabular Prior-Data Fitted Network (TabPFN) model proposed in this study serves as an integrated modeling method combination framework, integrating the advantages of machine learning technology and remote sensing technology. AEF fuses a 64-dimensional feature space and can automatically extract high-dimensional nonlinear features without the need for manual feature engineering. Based on a pre-training paradigm with the Transformer architecture, TabPFN learns general predictive patterns from massive synthetic tabular data. Compared with traditional machine learning methods, it saves a great deal of time devoted to hyperparameter tuning and demonstrates excellent predictive performance.

In terms of data quality, the AEF embedding features have undergone pre-training and standardization by the large foundation model, which can effectively reduce noise interference in multi-source data and improve data consistency, yet they are still constrained by the accuracy of the original remote sensing data. In terms of scale transferability, the AEF embeddings possess a unified spatial resolution. In terms of characterizing complex hydrogeological processes, the high-dimensional AEF features can capture information on topography, climate and other factors, but their capability to characterize subsurface processes such as deep aquifer structures and hydraulic conductivities is limited.

At present, AlphaEarth embedding data have demonstrated strong applicability and generalization ability in fields such as land use or land cover classification and urban air pollution prediction (Alvarez et al., 2025; Seydi, 2025; Wei et al., 2025). However, these approaches have not been systematically applied to groundwater research, especially for predicting groundwater level in coastal areas. Against this background, this study takes the Leizhou Peninsula as the study area and integrates measured groundwater level data from monitoring wells in 2023 with the 64-dimensional AEF embedding features. Multiple machine learning and deep learning models were developed, including Random Forest (RF), Extreme Gradient Boosting (XGBoost), Light Gradient Boosting Machine (LightGBM), Categorical Boosting (CatBoost), Support Vector Regression (SVR), and the TabPFN. Model hyperparameters were optimized using a combination of Bayesian optimization and Leave-One-Out Cross-Validation (LOOCV), and model performance was systematically evaluated using the coefficient of determination (R²), root mean square error (RMSE), and mean absolute error (MAE). In addition, Shapley Additive exPlanations (SHAP) were employed to interpret model behavior and elucidate the influence mechanisms of AEF embedding features on groundwater depth (GWD) prediction. Based on the best-performing model, a groundwater depth distribution map with a spatial resolution of 90 m was generated for the Leizhou Peninsula.

The objective of this study is to evaluate the applicability of the AEF for groundwater depth prediction in data-sparse coastal regions and to explore its potential in overcoming bottlenecks associated with multi-source data integration and feature construction. By providing a new technical paradigm for high-precision spatial prediction of groundwater depth, the findings offer a scientific basis for groundwater extraction control and refined management in the Leizhou Peninsula. Moreover, this framework can serve as a transferable reference for groundwater resource management and ecological risk assessment in other ecologically fragile coastal regions worldwide.

2 Study area and data

2.1 Study area

The Leizhou Peninsula is located at the southernmost part of mainland China, with a total area of approximately 12,470 km², making it the third largest peninsula in China (Figure 1a). It occupies a unique geographical position, bordered by the Lingnan Hills to the north, the South China Sea to the east, the Beibu Gulf to the west, and separated from Hainan Island by the Qiongzhou Strait to the south (Dong et al., 2020). The terrain is generally higher in the northern and southern distinct volcanic cones, transitioning to lower plains along the coast (Figure 1b). Tectonically, the Leizhou Peninsula is situated along the southern margin of the Eurasian Plate, and its geological evolution has been jointly influenced by the subduction of the Pacific Plate and the far-field tectonic stresses associated with the Indian Plate (Wang et al., 2024).

FIGURE 1

The study area lies within a low-latitude northern tropical monsoon climate zone, which is characterized by persistently high temperatures throughout the year. The long-term mean annual temperature ranges from 22 °C to 26 °C (Chai et al., 2025). Atmospheric moisture is primarily sourced from the western Pacific Ocean and the South China Sea; however, due to the generally low-relief topography and weak orographic uplift, moisture condensation efficiency is limited. As a result, the mean annual precipitation ranges from approximately 1200–1700 mm and exhibits pronounced spatial and temporal heterogeneity. Higher precipitation occurs in the northern and southern parts of the peninsula, followed by the eastern coastal zone and the southwestern inland areas, while the southwestern coastal region receives the lowest precipitation due to topographic effects (Zhang et al., 2021).

The groundwater system of the Leizhou Peninsula can be broadly classified into four main types: volcanic rock pore–fracture water, unconsolidated porous water, carbonate rock fracture–karst water, and bedrock fracture water. Among these, the widely distributed Quaternary basalt pore–fracture aquifer constitutes the most representative groundwater system in the region. Shallow groundwater is primarily recharged by atmospheric precipitation infiltration and, driven by the topographic gradient from the volcanic plateau to the coastal plain, flows radially from the central area toward the surrounding coastal zones or gradually follows the terrain toward the shoreline (Lu et al., 2015; Chen et al., 2023). Shallow groundwater vertically recharges deeper aquifers through overflow processes or primary basaltic fracture pathways and is ultimately discharged into the ocean in coastal areas in the form of spring groups or submarine groundwater discharge (Wang et al., 2024).

A total of 87 groundwater depth samples were collected across the region to ensure spatial representativeness. As shown in Figure 1c, the groundwater depth in the study area exhibits significant variability, ranging from approximately 1.2 m–17.8 m. The dataset shows a mean depth of 9.51 m and a median of 10.64 m, with a standard deviation of 4.14 m, indicating a relatively dispersed distribution with a slight deviation from a perfect normal distribution.

2.2 Data sources

The dependent variable in this study is the measured groundwater depth in the Leizhou Peninsula, derived from a dataset containing 87 groundwater monitoring samples collected in 2023. To better understand the hydrogeological characteristics of the study area, we conducted a statistical analysis of the sampling data (Figure 1c). As illustrated in the distribution histogram, the GWD exhibits significant spatial variability, ranging from a minimum of 1.17 m to a maximum of 17.8 m. The dataset shows a mean depth of 9.507 m and a median of 10.640 m, with a standard deviation of 4.136 m, indicating a relatively high dispersion in the water table levels across different monitoring sites. Although the data approximates an approximate normal distribution trend, a high frequency of samples is concentrated in the 10–13 m range, reflecting the heterogeneity of the local aquifer systems.

The independent variables consist of 64-dimensional satellite embedding features derived from AEF, a cutting-edge geospatial foundation model developed by Google DeepMind (Table 1). These features are not direct physical measurements; instead, they are condensed latent feature vectors generated by a self-supervised deep learning model. Trained on a massive corpus of Earth observation data, the AlphaEarth model integrates large-scale multimodal data, including optical imagery (Sentinel-2 and Landsat), SAR data (Sentinel-1), as well as topographic and climatic information (Brown et al., 2025). Unlike traditional remote sensing indices, AEF embeddings are designed to learn universal representations of the Earth’s surface, effectively capturing complex nonlinear characteristics and semantic relationships related to land surface cover, terrain morphology, and hydro-meteorological conditions even in data-sparse scenarios.

3 Methodology

Given the limited sample size (N = 87), the strategy was designed as follows: (1) The 87 groundwater depth samples were randomly partitioned into a training set (69 samples) and an independent test set (18 samples) at an 8:2 ratio. (2) The results of LOOCV were used as the target evaluation index of Bayesian hyperparameter optimization using only the training set data. The independent test set did not participate in any part of this process at all. (3) Based on the global optimal hyperparameters obtained by hyperparameter optimization, all training set samples were used to complete the final training of the model. (4) The final trained model was applied to the independent test set for evaluation, and the performance indicators such as R², RMSE, MAE were calculated to complete the final evaluation of the generalization ability of the model. (5) Based on the final trained model, the AEF feature grid with global resolution of 90 m in Leizhou Peninsula was deduced, and the spatial distribution map of regional groundwater depth and the prediction uncertainty distribution map were generated (Figure 2).

FIGURE 2

3.1 Predictive model construction

3.1.1 TabPFN

Tabular Prior-Data Fitted Network (TabPFN) is a tabular data prediction model based on the concept of prior learning (Hollmann et al., 2023). TabPFN adopts a Transformer architecture, in which both training samples and prediction samples are jointly encoded as sequence inputs. By learning a general prior distribution of tabular data during training, TabPFN can maintain relatively stable predictive performance even when the available sample size is limited.

3.1.2 Random Forest

Random Forest is an ensemble learning method based on the bagging strategy. It achieves robust estimation of the target variable by constructing multiple mutually independent regression trees and aggregating their predictions (Breiman, 2001). In groundwater level prediction, this mechanism allows Random Forest to effectively capture nonlinear relationships between groundwater level and environmental factors in high-dimensional feature space, while reducing variance and improving prediction robustness.

3.1.3 XGBoost

XGBoost is an ensemble learning model based on the gradient boosting framework. Its iterative optimization process enables XGBoost to continuously correct prediction residuals, thereby accurately characterizing the response of groundwater level to multi-source environmental features (Chen and Guestrin, 2016).

3.1.4 LightGBM

LightGBM is also based on the gradient boosting decision tree framework, and its basic prediction formulation is consistent with that of XGBoost. However, in terms of algorithmic implementation, LightGBM employs histogram-based feature binning and a leaf-wise tree growth strategy (Ke et al., 2017). This strategy prioritizes splitting the leaf node that contributes most to the reduction of the loss function, allowing lower training error to be achieved under the same tree depth.

3.1.5 CatBoost

CatBoost (Categorical Boosting) is an ensemble learning model based on Gradient Boosting Decision Trees (GBDT). Its core feature lies in the improvement of the gradient estimation process, which effectively mitigates the prediction bias problem introduced by target variable involvement in feature encoding in traditional gradient boosting methods (Prokhorenkova et al., 2018).

3.1.6 Support Vector Regression

Support Vector Regression (SVR) predicts groundwater level by constructing an optimal regression function in feature space. The optimization goal of SVR is to control the prediction error within a given tolerance range while ensuring model smoothness (El Bilali et al., 2023). This philosophy implies that SVR does not pursue an exact fit for all samples but focuses on characterizing the overall trend.

3.2 Model evaluation metrics

To quantitatively evaluate and compare the predictive performance of different spatial prediction models for groundwater depth, three widely used statistical metrics were adopted in this study: the coefficient of determination (), RMSE, and MAE. These metrics assess model performance from complementary perspectives, including explanatory power, overall error magnitude, and average deviation, and have been extensively applied in groundwater modeling and spatial prediction studies.

To quantify the uncertainty of model performance under small-sample conditions, we employed a Bootstrap resampling approach (1000 iterations with replacement). For each resample, a model was constructed using the identical hyperparameters and training procedures as the baseline model. Performance metrics were then evaluated on a fixed, independent test set. Finally, the mean and 95% confidence intervals (CI) for each indicator were calculated based on the 1000 simulation results to quantify the range of performance fluctuations.

In order to analyze the influence of the spatial autocorrelation between the test set and the training set on the model, the study used the spatial leave-one-region-out cross-validation, which divided the sample of the study area into four geographical partitions, and retained one partition each time as the validation set, and the other three partitions as the training set. Average performance metrics were computed to assess the spatial generalization ability of the model.

3.3 Interpretability analysis

To interpret the influence mechanisms of AEF features on groundwater level predictions and to enhance the transparency of machine learning models, SHAP were employed in this study. It is a unified, game-theoretic framework that quantifies the contribution of each input feature to an individual model prediction by fairly distributing the prediction output among all features. SHAP analysis was applied to the best-performing groundwater depth prediction model to quantify the relative importance of AEF features and to reveal their nonlinear effects on groundwater level variations. Global feature importance was assessed by aggregating the absolute SHAP values across all samples, while local explanations were used to investigate how specific AEF features influence groundwater level predictions at individual locations.

4 Results

4.1 Predictive model performance comparison

The 95% confidence intervals of the evaluation metrics calculated by the bootstrap method show that the TabPFN model performs most robustly, followed by XGBoost (Figure 3). TabPFN has the highest bootstrap mean estimate of R² (0.76), and the lower limit of the 95% confidence interval [0.55, 0.90] is even higher than the mean of the traditional Random Forest model (0.44). In terms of error indicators, the RMSE ([1.22, 3.27]) and MAE ([0.90, 2.44]) of TabPFN not only have the lowest mean values, but also have a relatively narrow span of the confidence intervals.

Figure 4 further visualizes the agreement between observed and predicted groundwater depth using scatter plots with marginal density histograms. Consistent with the statistical metrics, TabPFN (Figure 4a) exhibited the most robust generalization capability. Its regression fit line aligns almost perfectly with the 1:1 diagonal, achieving an of 0.793 and a Pearson correlation () of 0.894. The marginal histograms also indicate that the distribution of predictions by TabPFN closely matches the observed data, effectively preserving the data’s variance.

FIGURE 3

FIGURE 4

XGBoost (Figure 4b) ranked second (), maintaining a relatively tight cluster around the identity line, though slight deviations appeared in the lower value range. In comparison, CatBoost (Figure 4c) and SVR (Figure 4d) showed moderate performance. of 0.601 and 0.560 respectively. While these models captured the general trend, their scatter plots reveal a systematic bias: a tendency to overestimate low GWD values and underestimate high GWD values, indicated by the fit lines having slopes flatter than the 1:1 line.

Conversely, LightGBM (Figure 4e) and Random Forest (Figure 4f) displayed significant dispersion and limited predictive accuracy on this dataset, with values dropping below 0.50 ( and , respectively). As observed in their marginal histograms, the predicted distributions are narrower than the observed ones, clustering heavily between 6 m and 12 m. This suggests that these models suffered from a “smoothing effect”, failing to adequately capture extreme hydrological conditions.

The study quantified the bias of the predictions. The results demonstrate that all models exhibit a consistent trend of overestimating shallow groundwater depth and underestimating deep groundwater depth. Among all models, the TabPFN model outperforms others in controlling the magnitude of extreme-value prediction bias (Figure 5).

FIGURE 5

4.2 Cross ablation experiments

This study designed cross ablation experiments involving full-combination tests between AEF features and traditional manually engineered features, as well as between the TabPFN architecture and traditional ensemble learning models. The experimental results demonstrate that while keeping the underlying model unchanged, models utilizing AEF features outperform those using traditional manually engineered features significantly. This indicates that AEF, as a pre-trained foundation model-based feature extractor, can effectively integrate multi-source Earth observation data. With the same traditional feature input, the performance of TabPFN is superior to that of other models. The most significant performance improvement is achieved when AEF features are combined with the TabPFN architecture (Figure 6).

FIGURE 6

4.3 Spatial distribution of groundwater depth

Figure 7 illustrates the frequency distribution histograms and Kernel Density Estimation (KDE) curves of GWD predictions derived from six machine learning models (TabPFN, XGBoost, CatBoost, SVR, LightGBM, and Random Forest). In general, most models exhibited quasi-normal and unimodal distributions, such as TabPFN, XGBoost, CatBoost and Random Forest. Their predicted means were primarily clustered between 8.76 m and 9.31 m, with medians consistently surpassing the means, indicating a slight negative skewness. Notably, the TabPFN model predicted the deepest average depth (9.31 m) and demonstrated a relatively broad distributional span. In contrast, the SVR model displayed the most distinct characteristics, characterized by a pronounced leptokurtic distribution. Its predictions were highly concentrated around a lower mean of 7.88 m with the narrowest range, suggesting a tendency toward value fixation. While LightGBM generally aligned with the trends of other ensemble models, it exhibited minor fluctuations in the lower value range.

FIGURE 7

From a spatial perspective, the predicted GWD exhibits significant heterogeneity characterized by a distinct land-sea gradient (Figure 8a). In general, the groundwater depth tends to be deeper in the inland hinterlands and elevated terrains, contrasting sharply with the shallower level found along the coastal belts and low-lying plains. This pattern shows a strong consistency with the actual hydrogeological conditions of the Leizhou Peninsula.

FIGURE 8

Specifically, the regions shaded in dark blue represent areas with substantial GWD approaching 15 m, which are primarily aggregated in the northern interior of the peninsula and select inland volcanic highlands in the central and southern districts. Conversely, the light blue regions signify shallow GWD ranging from roughly 0–2 m. These shallow zones are widely dispersed across the jagged coastal margins, river deltas, and central depressions. The occurrence of shallow groundwater in these locations is largely driven by the prevailing flat topography and their immediate proximity to surface water bodies and marine boundaries.

To further validate the spatial reliability of the model, observed monitoring data were overlaid onto the prediction map. As visualized in Figure 8a and the zoomed-in inset of the Zhanjiang downtown area (Figure 8b), there is a high degree of spatial agreement between the predicted background and the observation points. The color of the observation dots generally matches the underlying prediction layer; for instance, red dots (indicating deep GWD) are predominantly located within dark blue zones, while yellow dots (shallow GWD) align well with light blue areas.

This spatial consistency is quantitatively confirmed in Figure 8c, which compares the observed GWD with the corresponding predicted grid values. The scatter plot reveals a strong positive correlation (), indicating that the model successfully captures the spatial variance of the groundwater table across the study area. Although the extreme high values are moderately underpredicted, the overall trend demonstrates the model’s robustness in mapping regional groundwater distribution.

The uncertainty map of GWD predictions demonstrates that uncertainty is maintained at a low level across the vast majority of the region. Conversely, zones characterized by high uncertainty are primarily distributed along the periphery of the study area or in regions with sparse or absent monitoring well coverage. This map not only explicitly delineates the reliability bounds of the model predictions but also provides a scientific foundation for the future optimization and densification of the groundwater monitoring network.

We conducted global and local Moran’s I of the model prediction residuals. As shown in Figure 9b, the global Moran’s I index was 0.044 (z-score = 0.787, p-value = 0.431), which was not statistically significant. This suggests the prediction residuals were spatially random with no significant spatial autocorrelation. Similarly, the local Moran’s I analysis (Figure 9a) showed no significant spatial clustering of residuals at most monitoring wells. These results confirm that the model sufficiently captured the environmental covariates driving the spatial variability of groundwater depth.

FIGURE 9

A decrease in R² was observed using spatial cross-validation (Figure 9c). The main reason is that spatial cross-validation reduces the spatial similarity between training and test samples, requiring the model to make predictions under more challenging conditions. This realistically reflects the influence of spatial heterogeneity on groundwater depth prediction. Such a decrease revealed the potential constraints of the uneven distribution of monitoring wells on prediction accuracy.

Although the rigorous statistical examinations corroborated the robustness of the present model framework within the study area, it should be noted that the model is affected by the spatial autocorrelation of sample distribution, given the limited dataset of only 87 monitoring wells. Accordingly, the generalizability of the conclusions drawn in this study may be limited when extrapolating to other extensive regions with more complex and highly heterogeneous hydrogeological conditions. Future research is recommended to increase the sample size, integrate multi-source remote sensing data, couple with physically based mechanistic models, or employ transfer learning techniques to further alleviate the limitations induced by small sample sizes.

4.4 SHAP interpretation results

Figure 10a integrates feature importance ranking with SHAP value distribution to visually demonstrate the global contribution and non-linear influence of various driving factors on the XGBoost model’s predictions of GWD. Based on the mean absolute SHAP values, the top three most critical features are identified as A51, A39, and A54, respectively. This indicates that these variables play a dominant role in the model’s ability to capture the variation patterns of GWD. Notably, the importance of A51 is the most significant, with its mean SHAP value far exceeding that of the other features.

FIGURE 10

In the beeswarm plot, red points represent high feature values while blue points represent low values, facilitating the analysis of the positive or negative impact of features on GWD. As the primary dominant factor, A51 exhibits a clear positive influence. The red dots (high values) are predominantly distributed along the positive SHAP axis (to the right of the zero axis), whereas the blue dots (low values) are located on the negative axis. This suggests that higher A51 values significantly increase the predicted GWD, while lower A51 values tend to result in shallower depths.

Conversely, A39 shows a distinct negative influence. The blue dots representing low values are concentrated in the positive SHAP region, and the rightmost points exceed 0.5. By contrast, the red dots for high values cluster near zero and in the negative region. This implies that lower A39 values act as a significant driving force for increased GWD.

To validate the hydrogeological and physical significance of the AEF embedding features, this study performed Pearson correlation analysis between the top ten features identified by SHAP and several key environmental variables, including the DEM, precipitation, Normalized Difference Vegetation Index (NDVI), and land use types (Figure 10b). The results strongly corroborate the physical interpretability of the AEF features. Specifically, feature A51, which exhibited high importance in the SHAP analysis, showed a significant correlation with the actual DEM (r = −0.47, p < 0.01), while A32 was significantly positively correlated with the DEM (r = 0.51, p < 0.01). It should be noted that during high-dimensional feature extraction via the AEF, the mapping direction in the latent space is arbitrary. Therefore, these significant correlations indicate that A51 and A32 effectively capture the spatial characteristics of topographic elevation within the study area. Similarly, feature A54 exhibited a significant positive correlation with precipitation (r = 0.41, p < 0.01), confirming its role in characterizing climatic conditions. A05 was significantly positively correlated with NDVI (r = 0.51, p < 0.01). These quantitative assessments demonstrate that AEF embeddings are not physically uninterpretable “black-box” variables; rather, they successfully encode and compress complex real-world environmental drivers into representative latent features.

4.5 Sobol sensitivity analysis

After the SHAP analysis, the variance-based Sobol global sensitivity analysis (Figure 11) was performed. Features A51 and A39 exhibit the highest total sensitivity index. More importantly, the vast majority of the variance contribution of these two features originates from first-order effects. This further corroborates our conclusions in the correlation analysis. Sobol analysis quantifies the interaction effects between AEF embedded features. For example, for features A13 and A46, their interaction effects account for a considerable proportion of their total sensitivities; However, for features A07 and A37, the interaction effects clearly surpass their own first-order effects. This suggests that the influence of these characteristics on groundwater depth does not act in isolation but depends on the synergy of other variables.

FIGURE 11

5 Discussion

5.1 Advantages of AEF and TabPFN in data-scarce regions

To address the challenges of sparse groundwater monitoring wells and the spatiotemporal scale mismatch among multi-source data in the Leizhou Peninsula, this study innovatively applied the AEF for feature engineering, coupled with the TabPFN model for spatial prediction. The results demonstrate that the TabPFN model, utilizing AEF embeddings, achieved superior performance on the independent testing set (, RMSE = 2.079 m), significantly outperforming traditional ensemble learning models such as XGBoost and Random Forest. This advantage is quantitatively supported by cross-ablation experiments (Figure 6), which demonstrated that integrating AEF features consistently outperformed manually engineered features across different algorithms. Furthermore, the TabPFN architecture maximized the utility of these embeddings.

This superiority is primarily attributed to the powerful representational capability of the AEF foundation model. Traditional methods rely on limited, manually extracted features, which often struggle to capture the complex non-linear relationships between surface environments and groundwater systems. In contrast, AEF utilizes a 64-dimensional shared embedding space to generate a unified representation of optical imagery, SAR, DEM, and climate data. This approach effectively resolves the issue of “data silos” and uncovers high-dimensional non-linear features latent within multi-source Earth observation data, thereby enabling the model to extract richer information even from limited training samples.

Crucially, the architectural advantage of TabPFN for small-sample data played a pivotal role in enhancing prediction performance. Given the scarcity of groundwater monitoring points, traditional gradient boosting tree models, particularly LightGBM and Random Forest, exhibited clear signs of overfitting and poor generalization in data-sparse zones. In contrast, TabPFN is built upon a Transformer-based architecture utilizing In-Context Learning (ICL) and a Prior-Data Fitting (PFN) mechanism. Consequently, it demonstrates superior robustness and generalization capabilities on small datasets without requiring extensive iterative training or hyperparameter tuning, thereby significantly reducing computational overhead and improving overall modeling efficiency.

5.2 Interpretability analysis of driving factors

Through SHAP analysis, this study elucidated the driving influence of AEF features on GWD. The results indicate that features A51 and A39 are the most critical factors influencing GWD in the Leizhou Peninsula. Inferred from the spatial distribution patterns and SHAP dependence plots, A51 (the primary dominant factor) exhibits a significant positive influence, which is highly consistent with the controlling effect of topographic elevation on GWD. The northern platform of the peninsula is characterized by higher terrain and greater GWD, where high values of A51 coincide with deep groundwater level. This indicates that the AEF model has effectively encoded geomorphological features from multi-modal data. In contrast, A39 exhibits a negative influence, with low values corresponding to greater depths. This is likely associated with regional climatic conditions or hydrological recharge characteristics (e.g., precipitation or infiltration rates). Ultimately, these findings confirm that the embedding features generated by AEF have successfully captured the topographic, climatic, and hydrogeological characteristics controlling groundwater dynamics through self-supervised learning.

5.3 Implications for regional management

The 90 m resolution spatial distribution map, generated based on the optimal model, validates the reliability of the model for regional extrapolation across the Leizhou Peninsula. For this ecologically fragile zone, this high-precision map holds significant managerial implications for implementing differentiated groundwater protection strategies. In coastal and estuarine delta areas with shallow GWD, vigilance is required regarding the risk of seawater intrusion induced by overexploitation. Conversely, in inland platform areas with greater depths, where water extraction costs for agricultural irrigation are higher, it is necessary to optimize groundwater extraction management and promote water-saving irrigation technologies. Future monitoring stations should be prioritized in transitional zones where the GWD gradient is steep to capture dynamic changes more effectively. The framework proposed in this study offers a novel, low-cost, and highly efficient pathway for refined groundwater management in data-scarce regions.

Furthermore, mitigating systematic bias at extreme values carries profound practical implications for regional water engineering and ecological risk management. In deep inland zones, severe underestimation of GWD can lead to inadequate drilling depth planning for agricultural wells, causing significant shortfalls in cost estimates. Conversely, in shallow coastal zones, overestimating GWD may create a false impression of abundant subsurface buffer capacity, which compromises flood early warning capabilities and masks the true vulnerability to seawater intrusion. While the significant reduction of prediction bias by the TabPFN model corroborates its practical utility, managers must remain aware that residual uncertainties persist due to small-sample constraints. Therefore, the generated high-resolution maps are best utilized as a macro-level decision-support reference to prioritize areas for future monitoring network densification, rather than for microscale engineering design.

5.4 Limitations and prospects

Despite the demonstrated potential of the AEF model in groundwater spatial prediction, this study acknowledges certain limitations. While global sensitivity analyses identified the dominant contributions of specific AEF embeddings, their exact physical mechanisms remain partially implicit. This study relies on spatial snapshot data from 2023 and does not account for temporal dynamics driven by seasonal precipitation or anthropogenic extraction. Without the validation of multi-temporal and long-term series data, the current results cannot support the long-term management decision of regional groundwater. Furthermore, the small sample size limits the model’s spatial generalizability, and the uneven distribution of sampling points also exacerbates this limitation. This can be clearly observed from the performance degradation that occurred during spatial cross-validation.

In future research, we will integrate long-term continuous monitoring data from the Leizhou Peninsula, expand the sample size across the entire study domain, enhance the spatial uniformity of sample distribution, and improve overall sample representativeness. Coupling these AI-driven features with measured hydrogeological parameters and physically based mechanistic models will further alleviate the constraints. These efforts are intended to further strengthen the model’s robustness and generalizability, thereby addressing the current limitations associated with data scarcity and spatial sampling bias.

6 Conclusion

To address the challenges of scarce groundwater monitoring data and the heterogeneity of multi-source data in the coastal ecological zone of the Leizhou Peninsula, this study proposed a GWD spatial prediction framework based on the AEF embeddings and the TabPFN. AEF facilitates the seamless fusion of multi-source heterogeneous data. This effectively overcomes the limitations of traditional manual feature engineering, which is often tedious and lacking generalization capability, providing a novel data representation paradigm for groundwater research in data-scarce regions.

The TabPFN model demonstrated optimal performance under small-sample conditions. Even with sparse monitoring well data, TabPFN, combined with AEF features, achieved high-precision spatial prediction. It achieved the highest predictive accuracy on the independent test set (R² = 0.793, RMSE = 2.079 m). This validated the superiority of this methodological approach in scenarios characterized by small sample sizes and high-dimensional data.

Furthermore, the robustness tests and global sensitivity analyses confirmed that the AEF embeddings physically encoded the terrain and climate driving factors, such as A51 and A54. It is notable that the AEF-TabPFN framework performed better than other models in mitigating the systematic overestimation and underestimation biases of groundwater depth, providing more cautious guidance for regional flood risk warning and irrigation cost planning. While limitations of small-sample sizes and spatial generalization persist, the proposed framework provides a highly cost-effective, macro-level preliminary tool.

In summary, this study validates the application potential of the deep learning paradigm fusing AEF and TabPFN in the field of hydrogeology. This approach not only provides a scientific basis for the sustainable management of groundwater resources in the Leizhou Peninsula but also offers a valuable reference for other data-scarce, ecologically fragile coastal zones worldwide.

References

• 1

  Aeschbach-HertigW.GleesonT. (2012). Regional strategies for the accelerating global problem of groundwater depletion. Nat. Geosci.5, 853–861. 10.1038/ngeo1617

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AlvarezC. I.Ulloa VacaC. A.Echeverria LlumipantaN. A. (2025). Machine learning for urban air quality prediction using google AlphaEarth foundations satellite embeddings: a case study of quito, Ecuador. Remote Sens.17, 3472. 10.3390/rs17203472

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  BooK. B. W.El-ShafieA.OthmanF.KhanMd. M. H.BirimaA. H.AhmedA. N. (2024). Groundwater level forecasting with machine learning models: a review. Water Res.252, 121249. 10.1016/j.watres.2024.121249

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  BorzìI. (2025). Modeling groundwater resources in data-scarce regions for sustainable management: methodologies and limits. Hydrology12, 11. 10.3390/hydrology12010011

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  BreimanL. (2001). Random forests. Mach. Learn.45, 5–32. 10.1023/A:1010933404324

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  BrownC. F.KazmierskiM. R.PasquarellaV. J.RucklidgeW. J.SamsikovaM.ZhangC.et al (2025). AlphaEarth foundations: an embedding field model for accurate and efficient global mapping from sparse label data. 10.48550/arXiv.2507.22291

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  ChaiX.LiuB.GuoF.ChenY.LinY.LiY.et al (2025). Machine learning-based remote sensing inversion and spatiotemporal characterization of Chl-a concentration in the Leizhou Peninsula coastal waters. J. Mar. Sci. Eng.13, 1787. 10.3390/jmse13091787

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  ChenT.GuestrinC. (2016). “XGBoost: a scalable tree boosting system,” in Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD ’16) (New York, NY: Association for Computing Machinery), 785–794. 10.1145/2939672.2939785

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  ChenY.SunY.ShiH.ZhaoH.GaoH.PanG.et al (2023). Environmental capacity and fluxes of land-sourced pollutants around the Leizhou Peninsula in the summer. Front. Mar. Sci.10, 1280753. 10.3389/fmars.2023.1280753

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  ChihiH.Ben Cheikh LarbiI. (2023). Spatial prediction of groundwater levels using machine learning and geostatistical models: a case study of coastal faulted aquifer systems in southeastern Tunisia. Hydrogeol. J.31, 1387–1404. 10.1007/s10040-023-02686-7

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  DongL.XuW.QianF.FanZ.WangJ.FuC.et al (2020). Evaluation of vertical groundwater flux from borehole temperatures: a case study in the Leizhou Peninsula, South China. Environ. Earth Sci.79, 171. 10.1007/s12665-020-08918-3

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  El BilaliA.AbdeslamT.AyoubN.LamaneH.EzzaouiniM. A.ElbeltagiA. (2023). An interpretable machine learning approach based on DNN, SVR, Extra Tree, and XGBoost models for predicting daily pan evaporation. J. Environ. Manage.327, 116890. 10.1016/j.jenvman.2022.116890

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  GharehbaghiA.GhasemlouniaR.AhmadiF.AlbajiM. (2022). Groundwater level prediction with meteorologically sensitive gated recurrent unit (GRU) neural networks. J. Hydrology612, 128262. 10.1016/j.jhydrol.2022.128262

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  HasanE.TarhuleA.HongY.IiiB. M. (2019). Assessment of physical water scarcity in Africa using GRACE and TRMM satellite data. Remote Sens.11, 904. 10.3390/rs11080904

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  HasanK.PaulS.ChyT. J.AntipovaA. (2021). Analysis of groundwater table variability and trend using ordinary kriging: the case study of Sylhet, Bangladesh. Appl. Water Sci.11, 120. 10.1007/s13201-021-01454-w

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  HollmannN.MüllerS.EggenspergerK.HutterF. (2023). TabPFN: a transformer that solves small tabular classification problems in a second. 10.48550/arXiv.2207.01848

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  HouriezL.PilarskiS.VahediB.AhmadalipourA.ScullyT. H.AflittoN.et al (2025). Scalable geospatial data generation using AlphaEarth foundations model. 10.48550/arXiv.2508.11739

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  IgwebuikeN.AjayiM.OkolieC.KanyerereT.HalihanT. (2024). Application of machine learning and deep learning for predicting groundwater levels in the West Coast aquifer system, South Africa. Earth Sci. Inf.18, 6. 10.1007/s12145-024-01623-w

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  JasechkoS.SeyboldH.PerroneD.FanY.ShamsudduhaM.TaylorR. G.et al (2024). Rapid groundwater decline and some cases of recovery in aquifers globally. Nature625, 715–721. 10.1038/s41586-023-06879-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  JódarJ.UrrutiaJ.HerreraC.CustodioE.Martos-RosilloS.LambánL. J. (2024). The catastrophic effects of groundwater intensive exploitation and Megadrought on aquifers in central Chile: global change impact projections in water resources based on groundwater balance modeling. Sci. Total Environ.914, 169651. 10.1016/j.scitotenv.2023.169651

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  KeG.MengQ.FinleyT.WangT.ChenW.MaW.et al (2017). “LightGBM: a highly efficient gradient boosting decision tree,” in (La jolla: neural information processing systems (nips)). GuyonI.LuxburgU. V.BengioS.WallachH.FergusR.VishwanathanS.et alAdv. Neural Inf. Process. Syst. 10.5555/3294996.3295074

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  KimJ.MaY.MaxwellR. M. (2024). Integrating groundwater pumping data with regression-enhanced random forest models to improve groundwater monitoring and management in a coastal region. Front. Water6, 1509945. 10.3389/frwa.2024.1509945

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  LiW.FinsaM. M.LaskeyK. B.HouserP.Douglas-BateR.LiW.et al (2023). Groundwater level prediction with machine learning to support sustainable irrigation in water scarcity regions. Water15, 3473. 10.3390/w15193473

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  LuY.TangC.ChenJ.ChenJ. (2015). Groundwater recharge and hydrogeochemical evolution in Leizhou Peninsula, China. J. Chem.2015, 427579. 10.1155/2015/427579

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  McCartneyM. P.RexW.YuW.UhlenbrookS.GnechtenR. (2022). Change in global freshwater storage. Colombo, Sri Lanka: International Water Management Institute IWMI. Available online at: https://cgspace.cgiar.org/items/8512e846-1e94-4f7d-a2c0-1d0dad3017a8 (Accessed September 28, 2025).

  • Google Scholar
  • View reference in article

• 26

  NooriR.MaghrebiM.JessenS.BateniS. M.HeggyE.JavadiS.et al (2023). Decline in Iran’s groundwater recharge. Nat. Commun.14, 6674. 10.1038/s41467-023-42411-2

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  OkiT.KanaeS. (2006). Global hydrological cycles and world water resources. Science313, 1068–1072. 10.1126/science.1128845

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 28

  PanZ.GuoZ.ChenK.LuW.ZhengC. (2025). A deep adaptive bidirectional generative adversarial neural network (Bi-GAN) for groundwater contamination source estimation. J. Hydrology653, 132753. 10.1016/j.jhydrol.2025.132753

  • CrossRef
  • Google Scholar

• 29

  ProkhorenkovaL.GusevG.VorobevA.DorogushA. V.GulinA. (2018). “CatBoost: unbiased boosting with categorical features,” in Advances in neural information Processing systems 31 nips 2018. La Jolla: neural Information Processing Systems nips.

  • Google Scholar

• 30

  RafikA.Ait BrahimY.AmazirhA.OuaraniM.BargamB.OuatikiH.et al (2023). Groundwater level forecasting in a data-scarce region through remote sensing data downscaling, hydrological modeling, and machine learning: a case study from Morocco. J. Hydrology Regional Stud.50, 101569. 10.1016/j.ejrh.2023.101569

  • CrossRef
  • Google Scholar

• 31

  RichardsonC. M.DavisK. L.Ruiz-GonzálezC.GuimondJ. A.MichaelH. A.PaldorA.et al (2024). The impacts of climate change on coastal groundwater. Nat. Rev. Earth Environ.5, 100–119. 10.1038/s43017-023-00500-2

  • CrossRef
  • Google Scholar

• 32

  SeydiS. T. (2025). Deep learning-based burned area mapping using Bi-Temporal siamese networks and AlphaEarth foundation datasets. 10.48550/arXiv.2509.07852

  • CrossRef
  • Google Scholar

• 33

  SunJ.HuL.LiD.SunK.YangZ. (2022). Data-driven models for accurate groundwater level prediction and their practical significance in groundwater management. J. Hydrology608, 127630. 10.1016/j.jhydrol.2022.127630

  • CrossRef
  • Google Scholar

• 34

  TaoH.HameedM. M.MarhoonH. A.Zounemat-KermaniM.HeddamS.KimS.et al (2022). Groundwater level prediction using machine learning models: a comprehensive review. Neurocomputing489, 271–308. 10.1016/j.neucom.2022.03.014

  • CrossRef
  • Google Scholar

• 35

  TollefsonJ. (2025). Google AI model mines trillions of images to create maps of Earth ‘at any place and time. Nature644, 313–314. 10.1038/d41586-025-02412-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36

  WangL.JiaB.XieZ.WangB.LiuS.LiR.et al (2022). Impact of groundwater extraction on hydrological process over the beijing-tianjin-hebei region, China. J. Hydrology609, 127689. 10.1016/j.jhydrol.2022.127689

  • CrossRef
  • Google Scholar

• 37

  WangY.CaoX.YeS.ZhouY.ZhangK.ChenB.et al (2024). Spatial distribution characteristics of groundwater quality and hydrochemical genesis of iron, manganese, and pH exceeding standards in the Leizhou Peninsula. Environ. Monit. China40, 183–197.

  • Google Scholar

• 38

  WeiY.ZhongR.YangY.WeiY.ZhongR.YangY. (2025). Groundwater fluoride prediction for sustainable water management: a comparative evaluation of machine learning approaches enhanced by satellite embeddings. Sustainability17, 8505. 10.3390/su17188505

  • CrossRef
  • Google Scholar

• 39

  YangJ.DaiH.LiuH.YeM.JiaoT.LiuZ.et al (2025). Advancing groundwater vulnerability assessment to nitrate contamination: a comprehensive evaluation of index-based, statistical, and machine learning approaches with sensitivity analysis. J. Hydrology663, 134189. 10.1016/j.jhydrol.2025.134189

  • CrossRef
  • Google Scholar

• 40

  ZhangP.DaiP.ZhangJ.LiJ.ZhaoH.SongZ. (2021). Spatiotemporal variation, speciation, and transport flux of TDP in Leizhou Peninsula coastal waters, South China Sea. Mar. Pollut. Bull.167, 112284. 10.1016/j.marpolbul.2021.112284

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  ZhouY.XiaoF. (2024). A glimpse into the new advances in Earth science research based on artificial intelligence and big data. Geosci. Front.31 (4), 1–6.

  • Google Scholar

• 42

  ZhouY.ZuoR. (2025). Application of big data mining, machine learning and artificial intelligence in ore deposits (Basel, Switzerland: MDPI). 10.3390/books978-3-7258-5464-6

  • CrossRef
  • Google Scholar

• 43

  ZhouY.ZhangL.ZhangO.WangJ. (2018). Big data mining and machine learning in Earth sciences. Guangzhou: Sun Yat-sen University Press, 1–269p.

  • Google Scholar