Implementation of Groundwater Lateral Flow and Human Water Regulation in CAS‐FGOALS‐g3

Both groundwater lateral flow (GLF) and human water regulation (HWR) significantly impact hydrological processes, climate, and even socioeconomic sustainable development. Reasonably representing these processes in climate system models is vital for improving hydrological predication and climate modeling. In this study, schemes including GLF and HWR were implemented into the Flexible Global Ocean‐Atmosphere‐Land System model grid‐point version 3 (CAS‐FGOALS‐g3) to investigate the hydroclimatic effects of GLF and HWR. Three groups of simulations using CAS‐FGOALS‐g3 were conducted for the period from 1976 to 2010. Comparisons between the simulations and the observations show their good performance in reproducing the hydrological processes. Results show that soil moisture and latent heat flux increased when GLF was included in the western United States, northern Australia, and northern South America, along with a shallower water table depth. The largest increases in latent flux are located in regions without water and energy limitations. Increased summer precipitation occurred in the western United States due to the wetting and cooling effects of GLF. Latent heat flux significantly increases in three key regions of the world (central United States, north China plain, and northern India), caused by wetting surface soil due to irrigation. The atmosphere also responded to HWR, with cooling at the 850 hPa level over northern India and Pakistan. Decreased precipitation occurred in India because the upward movement was weaker as a result of HWR. GLF can replenish the groundwater depression cone caused by overexploitation, especially in thick aquifers.


Introduction
Groundwater interacts with soil water, rivers, lakes, and oceans and plays an important role in the global water cycle. Groundwater dynamics affect soil moisture and hence water balance and energy fluxes at the land surface (Liang et al., 2003;Niu et al., 2007;York et al., 2002). Groundwater receives surplus water from soil and also sustains base flow in humid regions (Arnold et al., 2000;Fan et al., 2007). For arid and semi-arid regions, groundwater increases water content at root zones and receives river seepage, sustaining the channel ecosystem (Chen & Hu, 2004;Scanlon et al., 2012;Zeng et al., 2016a). The river water replenishment program that has occurred since 2000 for the Tarim River basin was carried out to recharge the surrounding aquifer through groundwater lateral flow (GLF), to maintain vegetation growth, and therefore to improve the river ecological environment Zhu et al., 2019). GLF also influences latent heat flux and transpiration partitioning through groundwater redistribution (Maxwell & Condon, 2016), which affects climate at local and possibly larger scales (Maxwell et al., 2007;Maxwell & Kollet, 2008). Hence, GLF performs a pivotal part in the hydrological process and energy flux simulation in climate models.
Human water regulation (HWR), such as the water used for farmland irrigation, industrial production, and life, can affect the hydrological process and climate. Groundwater constitutes about 30% of the global fresh water and plays a critical role in satisfying the basic needs of agriculture, domestic, and industrial water use (Niu et al., 2007;Wada et al., 2014). Groundwater was exploited continuously, and water has been consumed rapidly to sustain growing food demand and an increasing standard of life (Haddeland et al., 2006;Wada et al., 2014;Wisser et al., 2010). With the increase of water demand, many regions in the world, such as the north China plain, the central and western United States, and the Middle East, face the problem of water resource shortages (Rodell et al., 2009;Shi et al., 2013;. Excessive groundwater extraction will reduce the recharge of groundwater to rivers, intensifying hydrological drought (the occurrence of anomalously low streamflow) (Wada et al., 2013) and even reducing the area of lakes (Coe & Foley, 2001). Plant water uptake through roots can be reduced and the stability of the ecosystem can be damaged by the influence of excessive groundwater extraction on soil moisture (Zeng, Xie, Yu, et al., 2016). In addition, groundwater extraction for irrigation can enhance the local latent heat flux, affecting the atmospheric circulation Zou et al., 2014). The changed groundwater level induced by excessive groundwater extraction also changes the GLF. There is evidence that GLF recharges the groundwater depletion at a maximum rate of 40%, especially in plain regions (Zeng, Xie, Yu, et al., 2016a). This evidence indicates that the anthropogenic activities and lateral hydrological processes are inseparable. Therefore, it is necessary to reasonably represent the GLF and HWR in climate system models and promote the balance between human water resource management and the stability of the ecosystem.
Land surface models (LSMs) and climate models are effective tools for explaining how the refined hydrological structure affects soil moisture, land-atmosphere flux, and regional climate (Fan et al., 2019;Maxwell et al., 2007;Maxwell & Kollet, 2008;Miguez-Macho & Fan, 2012a, 2012b. As the lower boundary in the climate system models, current LSMs ignore the GLF, which is driven by topography. However, groundwater affects convection, advection, and precipitation by affecting the water-heat flux between land and atmosphere (Chen & Hu, 2004;Haddeland et al., 2006;. In recent years, a number of studies discuss incorporating GLF into models (Fan et al., 2013;Maxwell & Condon, 2016;. De Graaf et al. (2015) constructed an equilibrium groundwater table map with a high-resolution global-scale groundwater model using the modular hydrologic model MODFLOW. Zeng et al. (2018) incorporated a GLF module into the Community Land Model 4.5 (CLM 4.5). These studies do not consider the effects of human water use, which ignore the influence of groundwater extraction on groundwater horizontal movement. Some studies began to consider human water management schemes in LSMs and climate models. Ozdogan et al. (2010) considered irrigation and crop type in an LSM and evaluated the effects of irrigation on evapotranspiration and sensible heat flux. Döll et al. (2012) used the Water Global Assessment and Prognosis model to show that water withdrawal affected water storage variations. Pokhrel et al. (2012) incorporated the anthropogenic water regulation modules into an LSM and found that irrigation can cause a maximum increase in the latent heat flux averaged over June-August. Zou et al. (2015) used the Regional Climate Model to quantify the wetting and cooling effects of groundwater exploitation. These works either focus on irrigation water, use non-comprehensive water management plans, or do not incorporate more realistic groundwater flow processes, making it difficult to fully quantify the degree of human influence on water resources and the resulting hydrological and climatic effects.
To better understand the Earth system responses to external forcing changes and promote model development, the Coupling Model Working Group of the World Climate Research Programme proposed the Coupled Model Intercomparison Project (CMIP). As one of the participated models in CMIP6 (Eyring et al., 2016), Flexible Global Ocean-Atmosphere-Land System model grid-point version 3 (CAS-FGOASL-g3) is developed by the State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics, Institute of Atmospheric Physics, Chinese Academy of Sciences (LASG/ IAP), and considered to be a crucial modeling tool in climate science research (Zhou et al., 2018). The objective of this paper is to implement GLF and HWR in CAS-FGOALS-g3 for investigating the hydroclimatic effects of GLF and HWR and quantifying the relationship between GLF and HWR on a global scale. The model development and experimental design are presented in section 2. In section 3, we show the validation and evaluation. Finally, some concluding remarks are given in section 4.

Groundwater Lateral Flow
In this study, only the unconfined groundwater system was considered. The GLF model is based on Darcy's law and the Dupuit approximation (Bear, 1972), and it is expressed as where R is the net groundwater lateral discharge rate (when the value is greater than zero) or recharge rate (when the value is less than zero) per second per unit area (L/T); x and y are distances in the longitude and latitude directions (L), respectively; T is the transmissivity (integration of hydraulic conductivity over flow depth) (L 2 /T); and h is the groundwater table head (L). Each grid cell in the model exchanges water with its surrounding eight grids (Figure 1), and the numerical formulation of GLF for each grid cell is expressed as where i and j represent the row and column number of the model grid, respectively. R i,j is the net lateral discharge rate from the surrounding eight grid cells (L/T). When the groundwater table depth is deeper than the value of depth to bedrock (obtained from Shangguan et al., 2017), the grid cell can only receive water from the neighboring cells. n is the index of the eight neighboring grid cells. T i,j is the transmissivity of the grid (L 2 /T), w i,j is the width of the flow cross section of the grid cell (L), h i,j is the groundwater table head of the grid cell (L), h n is the groundwater table head of the number n neighbor of the grid cell (L), S i,j is the grid cell area (L 2 ), and l n is the center-to-center distance between the grid cell and its neighbor (L). The calculations of T i,j and w i,j refer to Zeng et al. (2018). After calculating R i,j , the groundwater table depth d′(L) and the unconfined aquifer water storage W′(L) are written as where Δt is the time step (T), s is the aquifer-specific yield, and d and W are the groundwater table depth (L) and unconfined aquifer water storage (L) before the calculation of GLF.
We used a similar scheme as Fan et al. (2013) to represent the conversion of groundwater and surface water. When the groundwater level is higher than the ground surface, that is, the groundwater converges on the surface to form a ponded area, and the excess water is routed to a river network to make the groundwater table depth equal to zero. Currently, the model does not describe the process of rivers recharging groundwater. Additional experiments were conducted (in the supporting information) to prove that this omission has a small effect on the global groundwater table depth patterns ( Figure S1).

Human Water Regulation
HWR includes surface water extraction and groundwater extraction. In an earlier study, it was coupled with the CLM 4.5 (Zeng et al., 2017;Zou et al., 2015). The HWR scheme is described in details in . Only the groundwater extraction scheme is considered in the present study. The groundwater extraction is expressed as where W′ and W ′′ are the water storage of aquifer before and after groundwater extraction (L), respectively, and q g is the groundwater extraction rate (L/T). After groundwater extraction, the groundwater table depth is written as where d′ and d ′′ are the groundwater table depth before and after groundwater extraction (L), respectively.
HWR was classified into three components: farmland irrigation, industrial, and domestic uses. The water for irrigation was visualized as effective rain into the ground surface, and the industrial and domestic water use have two sinks; the waste water is treated as runoff and the net water consumption is treated as evaporation Zou et al., 2014).

Incorporating GLF and HWR Into CAS-FGOALS-g3
The  is the Grid-Point Atmospheric Model of LASG/IAP version 3 (GAMIL3) (Tang et al., 2019), which uses a finite difference dynamical core with a horizontal resolution of about 2°and model top at 2.194 hPa (Li, Lin et al., 2013;Wang et al., 2004). The oceanic component is LICOM3, which uses a tripolar grid and a new dynamic framework . The sea ice component is CICE4 for CAS-FGOALS-g3 (Liu, 2010;Wang et al., 2009). The land component of CAS-FGOALS-g3 is the LSM developed by the Chinese Academy of Sciences (CAS-LSM) and contains GLF, HWR, and the frost and thaw fronts (FTFs) . CAS-LSM utilizes a subgrid hierarchy with multiple columns in a grid cell ( Figure 1b) and discretizes the soil profile with 15 layers, which are similar to CLM4.5. The river component of CAS-LSM was the River Transport Model (RTM), which used a linear transport scheme to route water from each grid to its downstream neighboring grid based on the steepest slope (Branstetter & Erickson, 2003).
Previous studies have shown that a reasonable running GLF model requires high model resolution, and the model grid length is preferably chosen to be around 1 km Zeng, Xie, Yu, et al., 2016). However, the calculation of LSMs with similar resolutions at a global scale can be very large, not to mention coupled with the climate system model. Therefore, the approach adopted in this study was to run a land surface process model at 2°resolution, run the GLF model at 1-km resolution, and use the conversion scheme in Zeng et al. (2018).
The conversion scheme was based on the subgrid structure of CAS-LSM ( Figure 1b). The soil water and groundwater exchange flux q recharge calculated by the CAS-LSM is input to the GLF as the upper boundary condition, and then, the GLF calculates the lateral flow and transmits the GLF flux q lateral to the CAS-LSM to update the water table. From the land model to the GLF downscaling, using the subgrid structure of the land model, the q recharge calculated by each plant functional type (PFT) in the 2°grid is assigned to the corresponding PFT in the 1-km grid. Bilinear interpolation was used on all GLF 1-km grids to achieve a smooth q recharge . From GLF to CAS-LSM upscaling, the model accumulates the lateral flux of all GLF 1-km grids in each 2°grid and divides by the number of grids to find the lateral flux q lateral of the 2°grid. In addition, a parallel approach was used to speed up the calculation, and a detailed description of the parallel running scheme can be found in Zeng et al. (2018).
It is important to point out that GLF can output groundwater table depth at a high resolution of 1 km, and the HWR is based on the LSM CAS-LSM with a 2°resolution. Therefore, the effect of HWR on the 1-km grid also needs to be considered, that is, the groundwater extraction rate q g in Equation 5 is downscaled from the 2°grid to the 1-km grid. If we still used the above downscaling scheme, it will lead to the extraction effects occurring in many non-extraction areas (such as forests and lakes) on the GLF model grid, which is inconsistent with the real situation. Since nearly 80% of the water abstraction is used for agricultural irrigation, here we approximately think that groundwater extraction of the GLF model grid occurs in farmland. Through the above schemes, we incorporated GLF and HWR into CAS-LSM and CAS-FGOALS-g3. The schematic diagram describing the coupled model is shown in Figure 1a.
In order to calculate the GLF and HWR modules, the input data needed are shown in Table 1. The depth to bedrock was used to control the GLF and the groundwater distribution (Shangguan et al., 2017). The global near-surface permeability and elevation were used to calculate Equation 2, obtained from Gleeson et al. (2011Gleeson et al. ( , 2014 and International Steering Committee for Global Mapping (Miliaresis & Argialas, 1999. The global land cover data were used for the resolution conversion and obtained from the Global Land Cover by National Mapping Organizations (Friedl et al., 2002;Tateishi et al., 2011). The initial value of groundwater table depth was obtained from Fan et al. (2013) and used for the initial condition before the model spin-up. The irrigation data were provided by , which was based on the Food and Agriculture Organization (FAO) of the United Nations and the Global Map of Irrigation Areas version 5.0 (GMIA5; Siebert et al., 2005). The GMIA5 data set provided the irrigated area and the area irrigated with groundwater for each grid. HWR from industrial and domestic sectors was estimated over the 2°global grid using the proportion of water use for each component provided by FAO.

Experimental Design
To quantify the impacts of GLF and HWR on land hydrological variables and energy fluxes, as well as to better understand the mechanism of climate feedback, three offline simulations from 1976 to 2010 were conducted. All of the three offline simulations have the same parameterizations, except the different modules activated. As shown in Table 2, the first simulation (CTL-off) activated the FTFs module only ("off" means offline model configuration). The second simulation (GC-off) only used the HWR and FTF modules. The last simulation (LGC-off) activated the GLF, HWR, and FTFs. In this paper, CTL only, including FTFs, was considered as the baseline. Thus, the GC-off runs do not include the GLF, whereas the LGC-off runs do, and the LGC-off-minus-GC-off difference would provide a measure of the impact of the GLF. The spatial resolution of the offline was set to 2°× 2°, and the GLF module was set to a resolution of~1 km at the equator. The time steps used in the LSM and the GLF module were 1,800 s and 6 hr, respectively (to reduce computation time). In addition, to reach a quasi-steady state before the formal simulation, we "spin up" the offline model for 100 years. The equilibrium groundwater table depth pattern (Fan et al., 2013) was the initial condition at 1-km resolution in the spin-up. The sea level was assumed to be 0 m, and lake water levels were set to the elevation of land surfaces as the boundary condition. The CRUNCEP forcing data set was used to run the offline model (Harris et al., 2014), which has a 6-hr time step at 0.5°horizontal spatial resolution.
Three groups of online simulations were conducted using CAS-FGOALS-g3: CTL, LC, and LGC. Note that only the atmospheric (GAMIL) and land (CAS-LSM) components were integrated in these experiments, and the monthly observations of sea surface temperatures (SSTs) and sea ice concentration from the HadISST data set (Rayner et al., 2003) were used for prescribing the ocean surface conditions. All simulations were conducted for a 41-year duration from 1976 to 2010 (using the period of 1970-1975 for spin-up) with a horizontal spacing of 2°× 2°and a time step of 10 min for both the GAMIL and CAS-LSM. To reduce the internal noise and enhance the forced signal caused by the HWR and GLF, the ensemble averages are used here. Each group of experiments contained five ensemble simulations using different initial conditions and then averaged the results over the ensemble for evaluation (Koster et al., 2002(Koster et al., , 2006). In the current study, five individual simulations, typically differing only in their initial atmospheric and land surface conditions, comprised each ensemble (CTL, GC, or LGC). To set the five sets of different "restart files," a 500-year 1850 control was run by the CAS-FGOALS-g3 and then continued with a 50-year 1850 control run, and a separation time of 10 years between the restart conditions was used (Koster et al., 2006). Simulations with corresponding ensemble members in the three sets of simulations use the same restart files. For instance, CTL1, GC1, and LGC1 shared the same initial conditions. Ensemble averages of the simulations for the years from 1976 to 2010 were analyzed to study the effects of HWR and GLF.

Validation Data
In order to evaluate the capacity of the model to simulate groundwater variation, we used 58 sites worldwide (7 in China, 11 in the United States, 13    Zeng et al. (2018). In addition to China, the rest of the continuous observation time covers more than 20 years. For each station, the nearest modeled groundwater table depth from the station is compared with the observations. To eliminate the errors associated with scale mismatch between simulated groundwater table depth and observations, we normalized the groundwater depth for both values. After standardization, the systematic error of resolution mismatch is also eliminated, while the groundwater dynamics (i.e., the variability of the groundwater table depth time series) are retained. As shown in Table 1, other observation data applied in this study include the groundwater table depth data for~1.6 million stations around the world, offered by Fan et al. (2013); the satellite-based soil moisture product of the European Space Agency's Climate Change Initiative (ESA-CCI) Gruber et al., 2017;Gruber et al., 2019); the in situ soil moisture data for 140 stations from the International Soil Moisture Network (ISMN) (Dorigo et al., 2011;Dorigo et al., 2013;Robock et al., 2015); the terrestrial water storage derived from the GRACE satellite data (Rodell et al., 2009;Tapley et al., 2004;Tiwari et al., 2009); and the latent heat flux product from the FLUXNET-MTE, which is generated by using a multitree ensemble (MTE) method and FLUXNET tower flux observations (Jung et al., 2019;Tramontana et al., 2016).

Validation 3.1.1. Soil Moisture
Soil moisture is a key participant in land-atmosphere interactions. To better analyze the climate effects later, we first validate the reasonability of simulated soil moisture. Figure 2 compares the climatology spatial distribution of the simulated top 10-cm soil moisture with the ESA CCI product between 1979 and 2010. The top 10-cm soil moisture is a weighted average of the first four soil layer thicknesses (the thicknesses are 1.75, 2.76, 4.55, and 7.5 cm, and the weights are 0.175, 0.276, 0.455, and 0.094, respectively). Figures 2a  and 2b show that the simulated results agree with the ESA CCI in most areas (the absolute value of the difference is less than 0.05 mm 3 mm −3 ). Figure 2c shows the differences between simulations and ESA CCI data. ESA-CCI had lower soil moisture compared with the simulated results from south India and Southeast Asia. This deviation is partly from the precipitation deviation of atmospheric simulations. Results from LGC were also compared with in situ soil moisture observations to evaluate the model capability in soil moisture variation simulations. Excluding sites with a lot of missing data during the analysis period, we used a total of 140 sites (48 in United States, 40 in China, 8 in Australia, 12 in France, 30 in Russia, and 2 in Niger networks). Considering the coarse resolution of the modeled results and the closeness of the sites, we compared the regional averages of the soil moisture time series. Figures S4-S9 Figure 3 shows the temporal correlation coefficients between the annual time series from the LGC simulation and observations for each site in different regions. At the same time, Figure 3 also shows the root mean square error (RMSE) of each site (different colors represent the magnitude of RMSE). As seen from Figure 3, with the exception of Africa, the correlation coefficients of most sites were significantly positive. In China and Peru in particular, almost all sites passed the 95% significance test, demonstrating the ability of the model to simulate the annual variation of groundwater. The performance of the model is better in China, Peru, and France (smaller RMSE and significant correlation). While in Portugal and some African countries, the correlation coefficients are rather low, which may be related to the errors caused by atmospheric simulation and measurements.

Journal of Geophysical Research: Atmospheres
The modeled time-averaged groundwater table depth was also compared with observations from~1.6 million stations worldwide (Fan et al., 2013) to validate the model. Figure 4 shows the groundwater table depth of the LGC simulation minus the observations in Africa, Asia, the United States, Europe, Canada, Australia, and South America. Generally, the simulated water table depth is shallower than the observations except in three severe groundwater extraction regions (central United States, north China plain, and northern India), a possible explanation is insufficient HWR. Figure 4f indicates that the simulated water table depth is deeper than observed in western Australia, which is most likely due to the contribution of the poorly defined boundary condition (here, the sea level was assumed to be 0 m). The above comparison illustrates the ability of the model to represent reasonable groundwater changes and spatial patterns and also emphasizes the importance of boundary conditions for groundwater simulation.

Latent Heat Flux
Changes in groundwater will affect latent heat flux through interaction with the soil, so we further validate latent heat flux. Figure 5 shows the climatology spatial distribution of the latent heat flux from the CTL minus products from FLUXNET-MTE, GC minus CTL, and LGC minus CTL. As Figure 5 shows, CTL generally underestimates latent heat flux relative to the FLUXNET-MTE data set, in particular over South America and southern Africa. However, Figure 5a shows that the strong latent heat flux indicated by FLUXNET-MTE in northern India and north China plain was not accurately simulated in the CTL. In contrast, after the HWR module was included, the underestimation was improved (the difference between the GC simulation minus CTL becomes positive). The contrasting accuracy of latent heat flux simulations in CTL, GC, and LGC scenarios occurred because, according to HWR, areas in northern India are equipped with intensive farmland irrigation, and the irrigated water would have significantly enhanced the latent heat flux. Moreover, the underestimation of latent heat flux simulation in South America and Africa was alleviated when considering GLF and HWR (Figure 5c). This comparison stressed the importance of including GLF and HWR when studying the water-energy cycle of regions and confidently validated the ability of the CAS-FGOALS-g3 in hydrological variables.

Terrestrial Water Storage Change
As shown in Equations 4 and 5, GLF and HWR can directly affect the terrestrial water storage. Then, we compared the simulated terrestrial water storage anomalies (TWSAs) with GRACE satellite data. Figure 6a shows the patterns of the temporal correlation coefficients between the monthly TWSA from the LGC simulations  Figure 6a show that the model is able to capture the TWS seasonality for most regions in the world, with the exception of the cold regions at high latitudes (e.g., Iceland and Russia). These deficiencies are most likely because the glacial process has not been represented sufficiently well. Also, the changes of hydrological parameters (e.g., the transmissivity) caused by frozen conditions are not considered, which may cause errors. Figures 6b-6d compare the trends in simulated TWSA with those from GRACE, estimated as the average change in centimeters per month over the period In general, the LGC simulation can capture major groundwater-depleted regions (e.g., the north China plain, Pakistan, India, central United States, and parts of the Middle East). For these regions, the changes of TWSA for the LGC simulation are larger than the GRACE. This is mainly due to the imperfect HWR, as well as the fact that the spatial averaging used in the GRACE may cause deviation. Moreover, there are some clear differences over southern Africa, western Australia, and some parts of Europe. A possible explanation is errors in atmospheric simulation, which are not very accurate in these regions. Finally, the trend from the LGC simulation was −0.0314 cm month −1 . This value was much closer to the actual trend of −0.0252 cm month −1 indicated by the GRACE data than to −0.0184 cm month −1 predicted by the GC simulation. The trend from the CTL run was 0.0095 cm month −1 . The deviations between the TWSA trend predicted by the CTL simulation and the observed TWSA trend indicated the rationality of GLF and HWR being considered in the climate system model.

Hydroclimate Response to GLF
The previous section demonstrated reasonability in the hydrological variable simulation when GLF and HWR are considered in the model. In this section, the effects of GLF on land surface hydrology and climate simulations are investigated at the global scale. To gain the direct changes in land surface variables due to GLF and HWR regardless of the atmospheric feedbacks and to better explain the mechanism of climate feedback, three offline simulation results were used to distinguish the hydrologic response. The global distributions of differences for groundwater table depth, surface soil moisture, and latent heat flux over 1976-2010 between LGC-off and GC-off are presented in Figure 7. Figure 7a shows that GLF enhances the spatial variability of the groundwater table depth globally. The groundwater level has risen about 1-5 m in eastern North America and eastern South America, while in most parts of Europe and Australia, groundwater table depth is deepened to varying degrees after considering GLF. The largest increases in groundwater table depth are in the low-latitude regions of the Northern Hemisphere, especially in the Sahara region (Figure 7a). This is mainly due to the fact that the arid region receives less recharge, and the GLF is more pronounced. In most parts of the world, the differences of groundwater table depth passed the 95% confidence level in the t test (stippled areas). Changes in groundwater table depth cause corresponding changes in soil moisture. Compared to the GC-off simulation, global latent heat flux increases by~2% in the LGC-off simulation with the largest difference in the regions located between 15°S and 30°N.
From the above results, we found that the latent heat flux change is basically consistent with the change of groundwater level, but the individual areas, such as the rainforests of South America, do not have the consistency expected. We further evaluated the difference of latent heat flux for the Köppen-Geiger climate types (Kottek et al., 2006). A brief description of the climate classification is as follows: the first letter refers to the climate types: tropical (A), arid (B), temperate (C), and cold (D). The second letter indicates the precipitation conditions: rainforest (f), monsoon (m), and savannah (s) in tropical and desert (W) and steppe (S) in arid, dry summer (s), dry winter (w), and without dry season (f) in temperate and cold climates. The last letter refers to hot (h) and cold (k) in arid and hot climates (a), warm summer (b), cold summer (c), and very cold summer (d) in temperate and cold climates. Moreover, we used the climate regions defined by Feng and Zhang (2015), where the first climate letter labeled "arid" was the arid regions, the second letter "f" was the humid regions, and the other regions were the transitional regions. As Figure 7d shows, the largest increases of latent flux are located in the Aw, Cs, and Cw subregions, which are characterized by a tropical or temperate climate and a shallow water table depth. The difference of latent heat flux in the BW and BS subregions, except BWk, which are characterized by hot deserts and steppe and accompanied by deep water table depth, is relatively small. Small changes of latent heat flux occur in cold regions (such as Ds and Dw subregions), which is energy limited. In humid regions (such as Cf and Df subregions), the change in moisture caused by GLF is much smaller than the climatic effect, so the latent heat changes little. These results indicate that LGC-off-GC-off) of latent heat flux for 28 Köppen-Geiger climate types. For each type, the color of the bar represents the water table depth. The dots indicate the difference passed the 95% confidence level in the Student's t test (the null hypothesis is that there is no difference between the results from two simulations, and the significance level is 0.05).

10.1029/2019JD032289
Journal of Geophysical Research: Atmospheres the largest increases of latent heat flux due to GLF are located in regions without water and energy limitations.
To examine the climate response to GLF, we compared the LGC and GC runs on an annual bias. The 850-hPa temperature in the LGC run is cooler than the GC run by about 0.1°C (with p value < 0.05) when ensemble averaged over the northern part of the Middle East and western America. Note that we will not focus on the Middle East later due to the lack of interaction between lake and groundwater in the model and the uncertainty that may exist. The cooling effect of GLF is mainly due to the increase in latent heat flux caused by the reduction of the groundwater table depth. Figure 8b indicates that the precipitation in summer is increased in western America. To understand the precipitation response to GLF, we analyzed the 850-hPa wind field and accumulated the horizontal vapor flux in summer (Figures 8c and 8d). GLF decreased the groundwater table depth, increased soil moisture and evapotranspiration, and transported water from the West Coast to the inland regions of the western United States (Figure 8c) through water vapor advection, resulting in an increase in summer precipitation within the western United States, about 0.1 mm d −1 . Moreover, we mapped the LGC-minus-GC differences for June-August (JJA) horizontal vapor flux integrated from the 850 to 200 hPa level; the result is shown in Figure 8d. The accumulated horizontal vapor flux is calculated as where F is the flux of vertically accumulated horizontal vapor (kg m −1 s −1 ); z1 and z2 are the heights of the 850 and 200 hPa pressure surfaces (m), respectively; q is the specific humidity in each layer (kg m −3 ); v is the wind vector in each layer (m s −1 ); dz is the thickness of each layer (m); and Δ represents the difference between the LGC and GC ensemble simulations. As shown in Figure 8d, there is indeed a horizontal flow of water vapor from the coast to the inland caused by GLF in western America, and the summer precipitation was subsequently increased.

Hydroclimate Response to HWR
Similarly, we analyze the hydroclimate effects of HWR. Figure 9 shows hydrological variable changes due to HWR. Compared to the effect of GLF, HWR effects are of similar magnitude, but HWR has large differences in magnitude around the world. Figure 9a shows that a significantly dropping groundwater table appears in the central United States, parts of the north China plain, and northern India and Pakistan due to the long-term exploitation. Global average annual mean groundwater table depth increases less than 1 m and increases 2-6 m in severe extraction regions (Figure 9d) because of high water demand. The shaded areas in Figures 9a-9c are at the 95% confidence level based on the t test. The differences in water table depth are correlated with the water demand. Differences in surface soil moisture are shown in Figure 9b, which has a similar spatial pattern with the changes of the water table depth. Areas that are extensively irrigated (northern India and China) would be expected to have higher soil water contents (Figure 9d) than areas with less irrigation (central United States). This is mainly because the irrigation water use was added to the topsoil as precipitation. Note that the magnitude of soil moisture change is not large, but some areas have passed significant tests. In addition, HWR accelerates or decelerates soil moisture trends under climate change especially in regions with high water demand for irrigation, production, and human consumption (Wang et al., 2019). The effect of HWR on global latent heat flux is slightly smaller than GLF, as latent heat flux mainly increases over irrigated areas. Globally, the increases of latent heat flux are almost negligible (less than 0.2 W m −2 ), while in the main extraction regions, the effects of HWR cannot be ignored (about 3 W m −2 ). The changes of latent heat flux correspond to the changes of surface soil moisture (Figure 9c). The increase of latent heat flux is mainly due to the increase of soil moisture as a result of irrigation. The dots indicate the difference passed the 95% confidence level in the Student's t test (the null hypothesis is that there is no difference between the results from two simulations, and the significance level is 0.05).

Journal of Geophysical Research: Atmospheres
demonstrating that the cooling effect was caused by groundwater irrigation. Therefore, HWR, especially irrigation, increases the surface soil moisture and evapotranspiration and leads to an increase for latent heat flux. The HWR process contributes to changes in the water and energy budgets, resulting in a decreasing effect on the surface temperature. Note that the cooled area in Russia is not the area experiencing severe groundwater extraction, which may be because the prevailing wind of India in the irrigation seasons takes the cooling air farther north and causes the reduction in surface temperature (Atlas et al., 1996;Young, 1999;. Such a climate response could not be simulated by an offline model (Döll et al., 2012), which stresses the necessity of using a global-scale land-atmosphere coupling model to study the climate effects of groundwater extraction. Moreover, we find an increase in the surface temperature in the Tibetan Plateau in eastern China bordering northern India (more than 0.4°C). The enhanced temperature may be a result of the increase in the air density and the local geopotential height induced by the irrigation of India (Zeng, Xie, & Zou, 2016b).
We next examine the JJA precipitation response to HWR by calculating the precipitation changes between the GC and CTL runs (Figure 10b). The results show that the simulated JJA precipitation increased over the groundwater extraction regions, such as parts of Europe, India, and eastern China. Over the groundwater extraction regions, the precipitation increased by about 0.1 mm d −1 . The increased precipitation of the north China plain is consistent with the strong summer land-atmosphere coupling interaction due to increases in local soil moisture after HWR (Dirmeyer et al., 2013;. Groundwater extraction can affect precipitation by changing the surface soil moisture, while the opposite effects of groundwater extraction on precipitation occurred in India. We further compare the changes in the JJA wind field and horizontal vapor water flux. Figures 10c and 10d show that the opposite wind anomaly direction appears at 850 and 200 hPa, a wind anomaly from the Indian subcontinent to the Arabian Sea at 850 hPa and from the sea to continent at 200 hPa. The results indicate that a thermal circulation was produced by HWR between the Arabian Sea and Indian subcontinent. Although the water vapor is flowing from the Arabian Sea to the northern Indian subcontinent (Figure 10e), the thermal circulation may reduce the effects of water vapor as most of the water vapor was located in the lower atmosphere.
To further understand how HWR affects the precipitation in India, we examine the changes in the vertical profile of the air temperature, vertical wind, and specific humidity ( Figure 11). Figure 11a shows that HWR causes a cooling effect, especially in the near-surface layer (about 0.23°C). The vertical wind anomaly was shown in Figure 11b (down is positive), which is calculated as where w z1 is the vertical wind (m s −1 ) at z1 (m), p is the height of the pressure surface (hPa), and u and v are the meridional wind and latitudinal wind (m s −1 ), respectively. The average horizontal divergence takes the arithmetic mean of the horizontal divergence at z1 and z2, and the horizontal divergence is calculated according to the central difference method.

Journal of Geophysical Research: Atmospheres
As shown in Figure 11c, the cooling effect in India may enhance the atmospheric stability due to the anomalous downdraft, which may reduce the upward water vapor transport, and the moisture decreased. The moisture increased in the lower troposphere, which was mainly induced by the increased surface soil moisture. The wetter and cooling effects of HWR have the opposite effect on precipitation, the wetter effect increases the amount of precipitation and effective energy of convection, while the cooling effect causes an abnormality in the downdraft, enhances atmospheric stability, reduces upward water vapor transport, and has a negative effect on precipitation. Note that the magnitude of the increase of water vapor is very small (0.2 g/kg) and almost negligible. The results reveal that the cooling effect of HWR in India, which is related to the anomalous downdraft, played a crucial role in the decrease in JJA precipitation. Our results on these regional-scale precipitation responses are consistent with previous studies (DeAngelis et al., 2010;Lo & Famiglietti, 2013;Zhang et al., 2017;Zou et al., 2015).

The Relationship Between GLF and HWR
As mentioned previously, excessive groundwater extraction will cause a lower local groundwater level than the surrounding area, which will cause a groundwater depression cone. Under the effect of the hydraulic gradient, the groundwater in the surrounding area will replenish the extraction area, slowing down the groundwater level decrease. Therefore, we further study the offset effect of GLF on the groundwater depression cone on a global scale. We additionally conducted another simulation called LC-off. The model configuration of LC-off was the same as CTL-off except that the GLF module was activated. Based on the difference between LGC-off and LC-off, we can evaluate the relationship between the change of GLF and groundwater extraction. Considering that this replenishment effect occurs locally, we use 1 km results for evaluation. On the whole, as the amount of groundwater extraction increases, the amount of GLF also increases ( Figure 12a). Especially in areas where the groundwater extraction exceeds 300 mm year −1 , the increase of GLF is obvious. However, the relationship between the change of GLF and the amount of groundwater extraction is not a significant linear. Even in overexploited areas, the change of GLF can be small. Zeng, Xie, Yu, et al. (2016) indicated that the terrain slope is a key factor determining the GLF offset in the extraction area. At the same time, a steep (large terrain slope) area generally has a thin aquifer, and the recharge from surrounding aquifers will be small due to low water transmissivity. Thus, we further explore the

10.1029/2019JD032289
Journal of Geophysical Research: Atmospheres relationship between the change of GLF and aquifer thickness. As shown in Figure 12a, a significant relationship appears only when groundwater extraction exceeds 300 mm year −1 . With the increase of the aquifer thickness, the replenishment effect of GLF also increases. When the aquifer thickness is greater than 7 m, the replenishment effect of GLF is basically stable. To better combine the replenishment effects of GLF with water management, we further explore the relationship between offset rate (ratio of the change of GLF to the amount of groundwater extraction) and aquifer thickness (Figure 12b). We find that the offset rate increases with aquifer thickness increases and up to 40%. Thus, a relatively high offset rate can be obtained by extraction in thick aquifers, which can be used as a reference to sustainable development of water resource.

Conclusions and Discussion
In the present study, schemes including GLF and HWR were incorporated into climate system model CAS-FGOALS-g3. Three groups of global simulations for the years 1976-2010 using the developed model were then conducted to quantify the effects of GLF and HWR on the hydroclimate. Comparing the simulations with the soil moisture from ISMN and ESA CCI, the in situ observational groundwater table depth from throughout the world, the latent heat flux from FLUXNET-MTE, and the terrestrial water storage from GRACE demonstrated that the coupling model is able to reliably represent hydrological variables. Results illustrate the hydroclimate response to GLF and HWR and the relationship between GLF and HWR. GLF can essentially enhance the spatial variability of groundwater table depth and change the soil moisture, evaporation, and latent heat flux. In the western United States, GLF decreased the groundwater table depth, increased soil moisture and evapotranspiration, and transported water from the West Coast to the inland regions through water vapor advection, resulting in an increase in summer precipitation. The HWR process contributes to changes in the water and energy budgets, resulting in a decreasing effect on surface temperature, especially a significant cooling effect in northern India. The north China plain experienced increased precipitation due to the wetter surface soil, while precipitation in northern India decreased, which is related to the cooling effect of HWR in India, may lead to the anomalous downdraft, enhance the atmospheric stability, and reduce the upward water vapor transport. The groundwater depression cone area caused by overexploitation can be replenished by the surrounding grid under the effect of the hydraulic gradient. With the increase of the aquifer thickness, the replenishment effect of GLF also increases. These results suggest that HWR and GLF, which are generally simplified or excluded in climate system models, may provide missing parts for modeling water and energy fluxes and improving water management.
Our study highlights the importance of considering GLF and HWR in climate system models. Both HWR and GLF are significantly impacting hydrology and climate through changing the water and energy cycles. The parameterization schemes can be expected to improve with the increasing abundance of global groundwater information. This coupled model system may provide a more comprehensive platform for water resource management. However, the limitations in the current study should be noted. Besides the uncertainties within the CAS-FGOALS-g3 (Zhou et al., 2018), the scheme GLF is highly simplified and is difficult to accurately express due to the complex subsurface properties. The permeability, land cover, and other data used for the initial condition of the GLF module also induced the uncertainties. In addition, groundwater-surface water interactions have not been fully considered. Results of additional experiments (in the supporting information) indicated that the impact of river-groundwater interactions requires more detailed research in typical riparian areas. For example, the groundwater table depth changed significantly in the Amazon River basin while the impact of river on the groundwater was small in the Nile River basin ( Figure S3). However, since the groundwater changed due to river-groundwater interactions mainly occurred in the riparian areas ( Figures S2 and S3), and the effects of the river-groundwater interactions on the global groundwater table depth pattern was not significantly ( Figure S1), the simplified scheme used in this study was reasonable. Future studies should focus on those aspects mentioned above.