Impact of Climate Change on the Spatio-Temporal Variation in Groundwater Storage in the Guangdong–Hong Kong–Macao Greater Bay Area

: The Guangdong–Hong Kong–Macao Greater Bay Area (GBA) is one of the world’s four major bay areas. Groundwater is indispensable in ensuring water supply for human production and living, as well as social and economic development. Studying the spatial–temporal variation in groundwater storage (GWS) and exploring the impact of climate change on GWS is of great signiﬁcance for water resource management in the GBA. In this work, we conducted a simulation using the Community Land Model version 5.0 (CLM5.0) and combined it with Gravity Recovery and Climate Experiment (GRACE) data to calculate GWS in the GBA. In addition, based on the multiple linear regression model, we quantitatively assessed the effects of different climate factors on the change in GWS in the GBA. Comparisons with groundwater wells, automatic weather stations, and satellite observations demonstrated reasonable results. Our results showed that precipitation and evapotranspiration are the main factors affecting the change in GWS in the GBA. Precipitation dominates GWS anomaly changes in areas where wetting and precipitation vary drastically, such as the northern part of Foshan. GWS is closely related to evapotranspiration, in which water and heat changes are signiﬁcant.


Introduction
The Guangdong-Hong Kong-Macao Greater Bay Area (GBA) is one of the world's four major bay areas. It is one of China's most economically and culturally developed regions and plays an important strategic role in its development. The sound development of the GBA requires the support of adequate water resources. Due to the developed internal river networks and broad drainage basins, the GBA is rich in water resources. However, the spatio-temporal distribution of precipitation is uneven; as it is mainly concentrated in the period between April and September [1]. There are more water resources in mountainous areas than plain areas, in eastern areas than western areas, and in southeastern coastal areas than northwestern inland areas [2]. Given economic and social development and population growth, the distribution of social-economic productivity does not match the allocation of water resources. Many cities in the GBA face severe water shortage, with the per capita water resources (the average amount of available water per person) below the threshold for severe water shortage (500 m 3 ) in populated areas [3], such as Hong Kong, Macao, Shenzhen, Zhuhai, and Guangzhou [4]. The massive discharge of waste water In addition, the influencing factors of GWS change in the GBA were analyzed based on a multiple linear regression method.

Study Area
The GBA is located in the south subtropical zone along the southeastern coast of Guangdong Province, China, extending from 21 • N to 25 • N and 111 • E to 116 • E. It covers an area of approximately 5.6 × 10 4 km 2 and includes nine municipalities directly administered by the governments of Guangzhou, Shenzhen, Foshan, Zhaoqing, Dongguan, Zhongshan, Huizhou, Jiangmen, and Zhuhai, as well as the special administrative regions of Hong Kong and Macao ( Figure 1). It is located in a tropical monsoon region and has an annual average temperature range of 21-23 • C. The average annual precipitation is 1300-2400 mm [27]. Under the influence of typhoons that occur in the East China Sea, rainfall increases from the northwest to the southeast. The annual relative humidity is 75-85% [2]. The hydrogeological characteristic of the GBA mainly comprises the Quaternary sequence and upper Mesozoic intrusion unit. The Quaternary sequence is characterized by anisotropic and heterogeneous groundwater flow [28]. Below the Quaternary sequence, the third series reconstruction layer sequence, karst layer, and cataclastic rocks constitute rich groundwater layers [9,11].
Sustainability 2023, 15, x FOR PEER REVIEW 3 of addition, the influencing factors of GWS change in the GBA were analyzed based on multiple linear regression method.

Study Area
The GBA is located in the south subtropical zone along the southeastern coast Guangdong Province, China, extending from 21° N to 25° N and 111° E to 116° E. It cov an area of approximately 5.6 × 10 4 km 2 and includes nine municipalities directly admin tered by the governments of Guangzhou, Shenzhen, Foshan, Zhaoqing, Donggua Zhongshan, Huizhou, Jiangmen, and Zhuhai, as well as the special administrative regio of Hong Kong and Macao ( Figure 1). It is located in a tropical monsoon region and has annual average temperature range of 21-23 °C. The average annual precipitation is 130 2400 mm [27]. Under the influence of typhoons that occur in the East China Sea, rainf increases from the northwest to the southeast. The annual relative humidity is 75-85% [ The hydrogeological characteristic of the GBA mainly comprises the Quaternary sequen and upper Mesozoic intrusion unit. The Quaternary sequence is characterized by ani tropic and heterogeneous groundwater flow [28]. Below the Quaternary sequence, t third series reconstruction layer sequence, karst layer, and cataclastic rocks constitute ri groundwater layers [9,11].

Model
This study used CLM5.0 to simulate the hydrological components in the GBA in t period 2003-2016. CLM5.0 is the land surface model of the Community Earth Syste Model [29]. Based on previous versions, CLM5.0 incorporated new changes in photosy thesis and stomatal conductance, improved the CLM's snow physics, added MOSART river routing, described external inputs to the nitrogen cycle, and updated the CRO

Model
This study used CLM5.0 to simulate the hydrological components in the GBA in the period 2003-2016. CLM5.0 is the land surface model of the Community Earth System Model [29]. Based on previous versions, CLM5.0 incorporated new changes in photosynthesis and stomatal conductance, improved the CLM's snow physics, added MOSART for river routing, described external inputs to the nitrogen cycle, and updated the CROP option. [29][30][31]. In CLM5.0, the vertical soil moisture transport was calculated via the mass conservation equation and Darcy's law [29,32].

Meteorological Forcing Data
The Chinese Meteorological Forcing Dataset (CMFD) was adopted, which was developed by the Institute of Tibetan Plateau Research, the Chinese Academy of Sciences, for the study of land surface processes in China [33][34][35]. The CMFD was obtained from Princeton reanalysis data, National Aeronautics and Space Administration/Global Energy and Water Cycle Experiment-Surface Radiation Budget radiation data, Tropical Rainfall Measuring Mission precipitation data, and Global Land Data Assimilation System data, and it was combined with the observation data from nearly 700 meteorological stations managed by the China Meteorological Administration via ANU-Spline statistical interpolation. The time span of CMFD is 1979-2018, the time resolution is 3 h, and the spatial resolution is 0.1 • × 0.1 • . It contains seven near-surface meteorological element variables, including the air temperature at 2 m above the ground, surface pressure, specific humidity, the wind speed at 10 m above the ground, downward short-wave radiation, downward long-wave radiation, and the precipitation rate.

Soil Properties Data
The soil properties data were based on a gridded Chinese soil properties dataset developed by ShangGuan et al. (2013) [36]. Its spatial resolution is 30 arcs × 30 arcs, covering an area of approximately 1 km × 1 km, including eight soil layers. In the dataset, 9000 soil profiles from the second National Soil Survey (1979)(1980)(1981)(1982)(1983)(1984)(1985) of the China Soil Map (1:1,000,000) were used [37], and the spatial distribution of soil properties was deduced via the polygonal linkage method [38].

Plant Functional Types Data
The plant functional types data used in the model comprised the plant functional types map in China with a resolution of 1 km released by the Cold and Arid Regions Data Center [39]. These data were obtained by applying the climate rules for land cover and plant functional type conversion proposed by Bonan et al. (2002) to the Multi-Source Integrated Chinese Land Cover map [39][40][41][42][43]. The Multi-Source Integrated Chinese Land Cover data include Chinese land use, Chinese plant types, Chinese glacier distribution, Chinese marsh wetland, and land cover.

GRACE Satellite Data
GRACE is a twin satellite that was launched by NASA and the German Aerospace Centre in 2002 to monitor changes in the Earth's gravity [12]. The GRACE data mainly estimate the monthly mean global gravity field data based on the accurate K-band ranging measurements of the K-band ranging system [44]. These measurements were obtained from the release datasets of the Center for Space Research at the University of Texas [45,46]. The spatial resolution of the data was 0.25 • × 0.25 • , the temporal resolution was months, and the time span was 2003-2016.

In Situ Groundwater Table Depth Data
To validate the estimation of GWS changes, the monthly groundwater table depth data of 26 in situ observation stations were collected for the period 2005-2016 from the Geological Environment Monitoring Institute of the China Geological Survey. Outliers that deviated from the average by three standard deviations were removed, and the values were averaged from these 26 in situ observation stations. This analysis was performed because of the large number of missing observations (15 in situ observation stations with more than 45% missing data). It should be noted that since the 0.25 • resolution of the estimated GWSA can only reflect the average state within a grid cell, few observations were available and they were clustered together with similar hydrogeological information. Thus, we used the average value of observations to validate the estimated GWS with 0.25 • resolution. In addition, there are few fast recharge systems available in the locations of the groundwater table depth stations used in this study [47]; thus, it was reasonable to utilize the above processing method for groundwater table depth data. Most in situ observation stations were located in the northwestern area of the GBA, and there were a few in situ observation stations in the southeastern area. The locations of the sites are shown by black triangular dots in Figure 1 [48]. These data were obtained using a combination of active and passive microwave sensors that measured microwave radiation emitted or reflected from the surface. The spatial resolution of the data was 0.25 • , and the temporal resolution was 1 day. For the convenience of comparison and verification, the monthly average of the time dimension was obtained.

Calculation of GWSA
The TWS anomaly (TWSA), which was provided by GRACE, was the anomaly of TWS relative to the long-term (2004-2009) mean, and it was a vertically integrated hydrologic variable [49] that consisted of GWSA, SnWS anomaly (SnWSA), SMS anomaly (SMSA), CWS anomaly (CWSA), and SWS anomaly (SWSA). The filtering process used to smoothen the GRACE data weakened the geophysical signal, using a set of scaling factors for signal recovery [50,51], though this method ignored some of the effects of surface water [52]. In addition, most of the LSMs did not contain SWS [23]. Therefore, GWSA was calculated as follows [17,25]: where, SnWSA, CWSA, and SMSA were calculated using CLM5.0. All storage anomalies were those values relative to the mean value of the period 2004-2009 ( Figure 2).

Multiple Linear Regression
Multiple linear regression was used to evaluate the contribution of different meteorological elements to the GWS. Due to the problem of different units and scales between variables, each variable required standardization.
The equation for standard deviation normalization was as follows:

Multiple Linear Regression
Multiple linear regression was used to evaluate the contribution of different meteorological elements to the GWS. Due to the problem of different units and scales between variables, each variable required standardization.
The equation for standard deviation normalization was as follows: where x is the variable that needs to be standardized, x is the standardized variable, µ is the mean value of the variable, and σ is the standard deviation of the variable. After standardization, multiple linear regression was used to calculate the corresponding regression coefficient for each variable. The regression equation was as follows [53,54]: where, y is the monthly average value of GWS at a given grid cell, and x 1 , x 2 , and x 3 are the monthly average precipitation, evapotranspiration, and temperature of each grid cell, respectively. a 1 , a 2 , and a 3 are regression coefficients, and b is the residual term. The monthly average precipitation and the monthly average temperature of each grid cell were derived using CMFD forcing data, and the monthly average evapotranspiration was obtained via CLM5.0 simulation. The contribution of precipitation, temperature, and evapotranspiration to GWS was quantified via the following equation [54]: where a i is the contribution coefficient of each variable to GWS.

Mann-Kendall (MK) Test
This study used the MK trend test to determine whether the trend was significant [55][56][57]. The MK trend test could be used for mutation testing, significance analysis, and evaluating the monotonicity or variability in time series. Its advantages were that the data did not need to obey a specific distribution, allowed missing data, were not disturbed by a few outliers, and had strong applicability. Given a time series with the length of n, the null hypothesis of no trend assumed that the time series were independently distributed. The standardized test statistic was calculated as follows: where the test statistic S is calculated as follows: where a positive value of S indicates an increasing trend, while a negative value indicates a decreasing trend. The sgn(x) is a symbolic function, and x i and x j represent the values of sequence i and j, respectively. When n > 8 [55,56], S is approximately normally distributed, and variance Var(S) is calculated as follows: where q k is the number of data in the tied group, and p is the number of tied groups. The null hypothesis of no trend was rejected when the absolute value was larger than the theoretical value, albeit within the statistical significance level concerned (this figure was set to 0.05 in this study).

Root Mean Square Error (RMSE)
Following previous studies [58,59], several statistical indicators are commonly used to assess the accuracy of simulated SM: the RMSE, Pearson's correlation coefficient (R), and Taylor's skill score (TSS).
The RMSE was calculated as follows [60]: where N is the number of samples. D i A and D i P are the observed and predicted values, respectively. In this study, this equation could be used to verify the reliability of simulated SM.

Pearson's Correlation Coefficient (R)
The R was used to analyze the correlation between simulated SM and observed SM, as well as the correlation between the respective variables and dependent variables of the multiple linear regression. The R was calculated as follows: where n is the number of samples, X is the mean value of variable X, and Y is the mean value of variable Y. The value range of r XY is [-1,1]. When r XY > 0, X and Y variables were positively linearly correlated; when <0, X and Y variables were negatively linearly correlated.

Mean Squared Error (MSE)
The MSE was calculated as follows [61]: where n is the number of samples, and Y i andŶ i are the observed and predicted values, respectively. The smaller the MSE value, the more accurate the prediction result. In this study, this equation was used to verify the reliability of simulated SM.

TSS
TSS [62] was used to quantitatively evaluate model performance and calculated as follows [59,63]: where σ f is the ratio of standard deviation between simulations and observations, R is the correlation coefficient between simulations and observations, and R 0 is the maximum possible correlation coefficient. The value of TSS ranges from 0 to 1, with higher values reflecting better performance.

Experimental Design
The CMFD was used to drive the CLM5.0. The GBA surface data were obtained using the CLM5.0 tools, and the tools replaced the original soil properties and plant functional type data for each grid cell with high-resolution data (details in Sections 2.3.2 and 2.3.3 [36,39]) A cyclically spin-up run was conducted for 50 years to balance out the hydrological process before the formal simulation. Combined with the hydrological components from the CLM5.0 simulation and GRACE data, the GWSA from 2003 to 2016 relative to the average value from 2004 to 2009 were calculated in mm.

Spatio-Temporal Distribution of GWSA
Regarding the spatial distribution (Figure 3a), the GWSA of the GBA was negative in the central region, generally positive in the surrounding areas, and slightly negative in the southwest region, though the absolute value was relatively small. The trends of change ( Figure 3b) and spatial distribution were similar, showing a decreasing trend in the central region and an increasing trend in surrounding areas. The northwestern and northeastern areas are primarily mountainous, with high altitudes and fast groundwater flow rates, and are affected by short-term climate change; therefore, the changes here were more significant. In the central part of the GBA, increased precipitation results in increased evapotranspiration, and the groundwater shows a significant decreasing trend. Figure 3c shows the temporal variation in GWSA. Moreover, GWSA in the GBA showed a slowly increasing trend from 2003 to 2016, with a growth rate of 0.12 mm/month ( Figure 3c). Moreover, the uncertainty in the estimated GWS in the GBA mainly originated from the SM due to the small values of SnWSA and CWSA, according to Equation (1). The relative deviation in the simulated SM from the ESACCI data is approximately 14.66%, which directly indicates that the SM simulated in this study is reasonable and indirectly indicates the reasonableness of the estimated GWS.

Verification of SM
To verify the accuracy of the model simulation, the simulated SM was compared to the observed SM (Figure 4a). The observed SM was used for the top 100 cm (10, 20, 30, 40, 50, 60, 80, and 100 cm). The first nine soil layers (1, 4, 9, 16, 26, 40, 58, 80, and 106 cm) of the simulated SM were interpolated to the SM at the corresponding observed depth. The weighted averages of the observed and simulated SMs for the top 100 cm were then cal-

Verification of SM
To verify the accuracy of the model simulation, the simulated SM was compared to the observed SM (Figure 4a). The observed SM was used for the top 100 cm (10, 20,  30,40,50,60,80, and 100 cm). The first nine soil layers (1,4,9,16,26,40,58,80, and 106 cm) of the simulated SM were interpolated to the SM at the corresponding observed depth. The weighted averages of the observed and simulated SMs for the top 100 cm were then calculated and compared. The R and MSE between the observed and simulated results (see Figure 4a) indicated a reasonable simulated SM. The trend of the simulated SM was 0.03 mm 3 /mm 3 /decade, which was similar to the trend of observed SM (0.046 mm 3 /mm 3 /decade).  The SM simulated using CLM5.0 and ESA CCI SM v6.1 data were compared and verified via the spatial scale. The weighted average value of the first three soil layers (thickness: 2 cm, 4 cm, and 6 cm; weights: 0.2, 0.4, and 0.4) of the CLM5.0 simulation was compared to that of the ESA CCI. Figure 5 shows that the simulated SM was higher than the ESA CCI in central GBA, but lower in the northwestern mountainous area, as well as the northeastern and southwestern areas. The spatial distribution between the simulated SM and ESA CCI was consistent (Figure 5a,b). Figure 5c,d shows the R and RMSE of SM simulated via CLM5.0 and the SM simulated via ESA CCI, respectively. There was a good correlation between both methods (Table 1)  The SM simulated using CLM5.0 and ESA CCI SM v6.1 data were compared and verified via the spatial scale. The weighted average value of the first three soil layers (thickness: 2 cm, 4 cm, and 6 cm; weights: 0.2, 0.4, and 0.4) of the CLM5.0 simulation was compared to that of the ESA CCI. Figure 5 shows that the simulated SM was higher than the ESA CCI in central GBA, but lower in the northwestern mountainous area, as well as the northeastern and southwestern areas. The spatial distribution between the simulated SM and ESA CCI was consistent (Figure 5a,b). Figure 5c,d shows the R and RMSE of SM simulated via CLM5.0 and the SM simulated via ESA CCI, respectively. There was a good correlation between both methods (Table 1)

Verification of GWSA
The calculated GWSA was then compared to the observed values ( Figure 6). To validate the estimated GWSA, we processed the observed groundwater table depth data in the following way. Firstly, outliers that were more than three times the standard deviation of the average value were eliminated. We then calculated the average value of observations, and finally, for consistency with the estimated GWSA, the observation anomaly relative to the 2005-2009 mean value was calculated (considering the available time series of observed data). A detailed explanation of the reasonableness of the processing method is given in Section 2.3.5. Next, the observed groundwater table depth anomaly was converted to the groundwater level anomaly by multiplying it by −1. The observed station was at a single point, and the resolution of the calculated GWSA was 0.25°; some uncertainties may have occurred due to the mismatch between the different resolutions. Thus, the calculated GWSA was compared to the observed groundwater level anomaly from the GBA average value perspective. We noted that the groundwater level data obtained from the water wells needed to be multiplied by specific yields to estimate the GWS [19]. Due to unreliable specific yield data, GWSA was verified by comparing the annual variation and trends of observed and calculated data.

Verification of GWSA
The calculated GWSA was then compared to the observed values ( Figure 6). To validate the estimated GWSA, we processed the observed groundwater table depth data in the following way. Firstly, outliers that were more than three times the standard deviation of the average value were eliminated. We then calculated the average value of observations, and finally, for consistency with the estimated GWSA, the observation anomaly relative to the 2005-2009 mean value was calculated (considering the available time series of observed data). A detailed explanation of the reasonableness of the processing method is given in Section 2.3.5. Next, the observed groundwater table depth anomaly was converted to the groundwater level anomaly by multiplying it by −1. The observed station was at a single point, and the resolution of the calculated GWSA was 0.25 • ; some uncertainties may have occurred due to the mismatch between the different resolutions. Thus, the calculated GWSA was compared to the observed groundwater level anomaly from the GBA average value perspective. We noted that the groundwater level data obtained from the water wells needed to be multiplied by specific yields to estimate the GWS [19]. Due to unreliable specific yield data, GWSA was verified by comparing the annual variation and trends of observed and calculated data. There was an increasing trend of GWSA in the GBA over the period 2005-2016. consistency between the trend of calculated GWSA and the observed data indicates the combination of GRACE and CLM5.0 simulation can better reflect variations in GW However, differences occurred in 2009, which were possibly caused by the mismatch tween the calculated GWSA and the groundwater table depth observations. Since the tual area of the GBA is 56,000 km 2 , with strong heterogenic hydrological and geolog conditions, GRACE data with a coarse resolution cannot capture these heterogeneities original spatial resolution of GRACE data is 1° × 1° [64]). There was also a large ga spatial scale between the sites and GRACE data, and the spatial scale mismatch cau part of the difference [65]. Although GRACE data are filtered on the balance of accu and spatial resolution to reduce noise interference [66], signal loss also increases with creasing noise interference, meaning that smaller spatial scales will result in more err especially those associated with the leakage problem [67].

Spatio-Temporal Variation in Precipitation, Temperature, and Evapotranspiration
Precipitation and evapotranspiration affect the availability of water infiltration in subsurface as recharge, and temperature is associated with evapotranspiration [68] deepen understanding of GWSA changes in the GBA, changes in climate factors were vestigated. Figure 7   There was an increasing trend of GWSA in the GBA over the period 2005-2016. The consistency between the trend of calculated GWSA and the observed data indicates that the combination of GRACE and CLM5.0 simulation can better reflect variations in GWSA. However, differences occurred in 2009, which were possibly caused by the mismatch between the calculated GWSA and the groundwater table depth observations. Since the actual area of the GBA is 56,000 km 2 , with strong heterogenic hydrological and geological conditions, GRACE data with a coarse resolution cannot capture these heterogeneities (the original spatial resolution of GRACE data is 1 • × 1 • [64]). There was also a large gap in spatial scale between the sites and GRACE data, and the spatial scale mismatch caused part of the difference [65]. Although GRACE data are filtered on the balance of accuracy and spatial resolution to reduce noise interference [66], signal loss also increases with decreasing noise interference, meaning that smaller spatial scales will result in more errors, especially those associated with the leakage problem [67].

Spatio-Temporal Variation in Precipitation, Temperature, and Evapotranspiration
Precipitation and evapotranspiration affect the availability of water infiltration in the subsurface as recharge, and temperature is associated with evapotranspiration [68]. To deepen understanding of GWSA changes in the GBA, changes in climate factors were investigated. Figure 7  Regarding the spatial distribution (Figure 8a,c,e), precipitation was higher in the east and lower in the west, as well as higher on the southeast coast than in the northwest inland area, which is similar to the findings of Cao et al. (2021) [2]. The temperature in the Pearl River estuary is higher than that in other areas, especially in places such as Guangzhou, where the temperature differs significantly from the surrounding areas, and the average summer temperature can reach 31 • C. The spatial distribution of evapotranspiration was not apparent, though high-value areas were primarily located in the regions with greater precipitation and higher temperatures. Regarding the spatial distribution (Figure 8a,c,e), precipitation was higher in the east and lower in the west, as well as higher on the southeast coast than in the northwest inland area, which is similar to the findings of Cao et al. (2021) [2]. The temperature in the Pearl River estuary is higher than that in other areas, especially in places such as Guangzhou, where the temperature differs significantly from the surrounding areas, and the average summer temperature can reach 31 °C. The spatial distribution of evapotranspiration was not apparent, though high-value areas were primarily located in the regions with greater precipitation and higher temperatures. Figure 8b shows that precipitation had an increasing trend over a wide area in the GBA, especially in the central region, with a maximum increase of 0.51 mm/month. Only the precipitation in the eastern part of the GBA had a decreasing trend of −0.14 mm/month. Temperature variation (Figure 8e) Figure 8b shows that precipitation had an increasing trend over a wide area in the GBA, especially in the central region, with a maximum increase of 0.51 mm/month. Only the precipitation in the eastern part of the GBA had a decreasing trend of −0.14 mm/month. Temperature variation (Figure 8e) was generally slight. However, there was a cooling trend in the central region of the GBA of approximately 0.8 • C per decade, while the temperature typically increased in other regions, with an increasing trend of up to 0.83 • C per decade in the southwestern GBA. Evapotranspiration (Figure 8f) increased in the central region of the GBA at 0.17 W/m 2 per month, while it decreased in the northwestern and southeastern regions at −0.04 W/m 2 per month. From the time series perspective, there was a slight increase in precipitation and evapotranspiration, as well as some fluctuations in temperature, but no significant trends, which is similar to a previous study's results [27]. a slight increase in precipitation and evapotranspiration, as well as some fluctuations in temperature, but no significant trends, which is similar to a previous study's results [27].

Contributions of Climatic Factors to GWSA
GWSA, precipitation, temperature, and evapotranspiration were standardized for each grid cell. A multiple linear regression model was established that used the standardized GWSA as the dependent variable and precipitation, temperature, and evapotranspiration as the independent variables. The contributions of each climatic factor to GWS were calculated based on Equation (1). Figure 9 shows that the average contribution of evapotranspiration was 58.53%, while average contributions of precipitation and temperature were 26.83% and 14.63%, respectively. The GWS of the GBA was mainly controlled via precipitation and evapotranspiration, while temperature had a relatively slight effect. The precipitation-controlled regions were located in the wetting areas of the GBA, while evapotranspiration-controlled regions were located in the northwestern mountainous area and the central region of the GBA.

Contributions of Climatic Factors to GWSA
GWSA, precipitation, temperature, and evapotranspiration were standardized for each grid cell. A multiple linear regression model was established that used the standardized GWSA as the dependent variable and precipitation, temperature, and evapotranspiration as the independent variables. The contributions of each climatic factor to GWS were calculated based on Equation (1). Figure 9 shows that the average contribution of evapotranspiration was 58.53%, while average contributions of precipitation and temperature were 26.83% and 14.63%, respectively. The GWS of the GBA was mainly controlled via precipitation and evapotranspiration, while temperature had a relatively slight effect. The precipitationcontrolled regions were located in the wetting areas of the GBA, while evapotranspirationcontrolled regions were located in the northwestern mountainous area and the central region of the GBA.
The northwest region of the GBA is mountainous and composed of basement units. The permeability decreases rapidly with increasing depth, and the precipitation is low. The precipitation mainly replenishes the river with surface runoff, and the remaining water enters the soil, which is greatly affected by evapotranspiration. Therefore, evapotranspiration has a significant effect on GWS. In the southwestern region of the GBA, precipitation was high, and temperature increased significantly; therefore, GWSA in the southwest of GBA was closely related to evapotranspiration and temperature. In northern Foshan and western Guangzhou, precipitation increased significantly; hence, precipitation in this region played a dominant role in GWS changes. In the central region of the GBA, precipitation and temperature increased significantly, and evapotranspiration also changed drastically under sufficient water-heat recharge. The significant decrease in GWS in this region was closely related to evapotranspiration. In the southeastern coastal area of the GBA, precipitation was high, and GWSA was not limited by moisture; therefore, the influence of evapotranspiration here was more pronounced. Overall, precipitation dominated GWSA changes in areas where wetting and precipitation varied drastically, such as northern Foshan. The GWS was closely related to evapotranspiration, where water and heat changes were significant. The northwest region of the GBA is mountainous and composed of basement units. The permeability decreases rapidly with increasing depth, and the precipitation is low. The precipitation mainly replenishes the river with surface runoff, and the remaining water enters the soil, which is greatly affected by evapotranspiration. Therefore, evapotranspiration has a significant effect on GWS. In the southwestern region of the GBA, precipitation was high, and temperature increased significantly; therefore, GWSA in the southwest of GBA was closely related to evapotranspiration and temperature. In northern Foshan and western Guangzhou, precipitation increased significantly; hence, precipitation in this region played a dominant role in GWS changes. In the central region of the GBA, precipitation and temperature increased significantly, and evapotranspiration also changed drastically under sufficient water-heat recharge. The significant decrease in GWS in this region was closely related to evapotranspiration. In the southeastern coastal area of the GBA, precipitation was high, and GWSA was not limited by moisture; therefore, the influence of evapotranspiration here was more pronounced. Overall, precipitation dominated GWSA changes in areas where wetting and precipitation varied drastically, such as northern Foshan. The GWS was closely related to evapotranspiration, where water and heat changes were significant.

Model Performance and Uncertainty
Combined GRACE data and CLM5.0 simulation in this work were successfully applied to GWSA estimation in the GBA. Many factors will affect the accuracy of the estimation. The initial resolution of GRACE was 1° × 1°, and the area of each grid cell was approximately 10,000 km 2 [45], while the area of the GBA was 56,000 km 2 . The coarse resolution of GRACE could lead to uncertainties for a small region. The leakage problem of GRACE could also influence the accuracy of TWSA. During GRACE data processing, mass changes in the study area were affected by mass changes outside of the area, an effect

Model Performance and Uncertainty
Combined GRACE data and CLM5.0 simulation in this work were successfully applied to GWSA estimation in the GBA. Many factors will affect the accuracy of the estimation. The initial resolution of GRACE was 1 • × 1 • , and the area of each grid cell was approximately 10,000 km 2 [45], while the area of the GBA was 56,000 km 2 . The coarse resolution of GRACE could lead to uncertainties for a small region. The leakage problem of GRACE could also influence the accuracy of TWSA. During GRACE data processing, mass changes in the study area were affected by mass changes outside of the area, an effect known as the leakage problem, which caused deviation in the retrieved TWSA [70,71].
In addition, the accuracy of meteorological data, land cover data, and model structure could lead to some bias in the simulation of hydrological components (e.g., SM) [72]. Since CLM5.0 is a land component used for climate simulations and does not fully consider geological information, there may be biases in the simulation of deep soils. Furthermore, the interpolation process introduces uncertainties during data processing because of the different resolutions between CLM5.0 simulations and GRACE.

Limitations
However, it should be noted that there are still some limitations. The impact of human activities needs to be considered when estimating the GWSA for smaller areas. According to the Pearl River Basin Water Resources Bulletin, the groundwater supply of several major cities in the GBA only accounts for 0.616% of the total water supply, which has a negligible impact on the GWSA. Considering that the primary purpose of this study was to estimate the changes in GWSA at a large scale, with the results showing that our approach is reasonable, the role of human activities was not further explored, though it is worth investigating in the future.
Future work could focus on the sustainable development and management of groundwater resources in the GBA; for example, the potential for groundwater extraction and changes in groundwater quality could be considered. Groundwater extraction potential is influenced by geological conditions, land use, topographic density, drainage density, slope, soil type, rainfall, elevation, and groundwater table fluctuations [73]. The hydrogeologic structure includes the hydro-lithologic characteristics of aquifers, which affect the variation in the groundwater table [74]. The suitability of groundwater for extraction is also related to groundwater quality, with pollutants from population growth and industrialization increasing the pressure on surface and groundwater [75]. Future studies could focus on strategies to address these aspects.

Conclusions
In this study, GRACE data were combined with the CLM5.0 simulation to calculate GWSA in the GBA, and the reliability of simulated SM and calculated GWSA was verified. Based on the multiple linear regression model, the influence of different climate factors on GWS changes in the GBA were quantitatively evaluated. The following findings were obtained.
The simulated SM using CLM5.0 was highly correlated with the observed data and could reflect the spatio-temporal changes in SM. A combination of GRACE data and CLM5.0 land surface simulation used to calculate GWSA in the GBA had a certain reliability. GWSA in the GBA showed a slowly increasing trend from 2003 to 2016, with a growth rate of 0.12 mm/month. Precipitation and evapotranspiration were the main factors affecting GWS changes in the GBA. Precipitation dominated GWSA changes in areas where wetting and precipitation varied drastically, such as northern Foshan and western Guangzhou. GWS was closely related to evapotranspiration in regions with significant water and heat changes.