# Influence of hillslope flow on hydroclimatic evolution under climate change

Pedro Arboleda-Obando, Agnès Ducharne, Frederique Cheruy, Anne Jost, Josefine Ghattas, Jeanne Colin, Camille Nous

► To cite this version:

Pedro Arboleda-Obando, Agnès Ducharne, Frederique Cheruy, Anne Jost, Josefine Ghattas, et al. Influence of hillslope flow on hydroclimatic evolution under climate change. Earth's Future, 2022, 10, pp.e2021EF002613. 10.1029/2021ef002613. hal-03797413

HAL Id: hal-03797413

https://hal.science/hal-03797413

Submitted on 4 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Earth's Future

## RESEARCH ARTICLE

10.1029/2021EF002613

**Key Points:**

* Hillslope flow sustains higher soil moisture, enhancing evapotranspiration and precipitation rates and cooling down the air temperature

* Hillslope flow slightly attenuates global warming

* At regional scale, hillslope flow attenuates climate change trends of all hydrological variables except for evapotranspiration increases

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

P. F. Arboleda Obando,
pedro.arboleda_obando@upmc.fr

**Citation:**

Arboleda Obando, P. F., Ducharne, A., Cheruy, F., Jost, A., Ghattas, J., Colin, J., & Nous, C. (2022). Influence of hillslope flow on hydroclimatic evolution under climate change. *Earth's Future, 10*, e2021EF002613. https://doi.org/10.1029/2021EF002613

Received 16 DEC 2021

Accepted 10 AUG 2022

© 2022. The Authors. Earth's Future published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

## Influence of Hillslope Flow on Hydroclimatic Evolution Under Climate Change

Pedro Felipe Arboleda Obando¹,² ID, Agnès Ducharne¹,² ID, Frédérique Cheruy²,³ ID, Anne Jost¹,² ID, Josefine Ghattas², Jeanne Colin⁴, and Camille Nous¹,²

¹Laboratoire METIS (UMR 7619, Sorbonne Université, CNRS, EPHE), Paris, France, ²Institut Pierre Simon Laplace (FR 636, Sorbonne Université, CNRS), Paris, France, ³Laboratoire de Météorologie Dynamique (UMR 8539, Sorbonne Université, CNRS), Paris, France, ⁴Centre National de Recherches Météorologiques (UMR 3589, Météo-France, CNRS, Université Fédérale de Toulouse), Toulouse, France

**Abstract** We analyzed the influence of hillslope flow on projections of climate change by comparing two transient climate simulations with the IPSL climate model between 1980 and 2100. Hillslope flow induces a reorganization and increment of soil moisture (+10%), which increases evapotranspiration (+4%) and precipitation (+1%) and decreases total runoff (-3%) and air temperature (-0.1 °C) on an annual average over land for 1980–2010 when compared to simulation not representing hillslope flow. These changes in land/ atmosphere fluxes are not homogenous and depend on regional climate and surface conditions. Hillslope flow also influences climate change projections. On average over land, it amplifies the positive trend of soil moisture (+23%), evapotranspiration (+50%), and precipitation (+7%) and slightly attenuates global warming (-1%), especially for daily maximum air temperature. The role of hillslope flow in supporting surface/atmosphere fluxes is more evident at a regional scale. Where precipitation is projected to decrease, hillslope flow is shown to attenuate the related declines in evapotranspiration, precipitation, and total runoff, regardless of aridity conditions and mean air temperature. Where precipitation is projected to increase, hillslope flow amplifies evapotranspiration enhancement but attenuates the increase in precipitation and total runoff. Warming is generally attenuated, especially in semiarid and cold areas, and humid and warm/temperate regions, but the signal is weak. These results demonstrate the role of hillslope flow in enhancing water and energy fluxes between the surface and the atmosphere. They also suggest that including hillslope flow in climate models would weaken the projected intensification of hydrological extreme events.

**Plain Language Summary** We analyze how the flows caused by topography, called hillslope flow, affect the evolution of climate using simulations from a climate model. Results show that hillslope flow increases soil moisture in the valleys. More soil moisture enhances evapotranspiration and precipitation and decreases total runoff and air temperature for the period 1980–2010. But the increase in water exchanges between land and atmosphere is not homogenous. In the future, hillslope flow amplifies positive trends of climate change for soil moisture, evapotranspiration and precipitation, while global warming is minimally slowed. The role of hillslope flow in sustaining exchanges of water and energy in the future is most evident in the regions. Where precipitation decreases in the future, evapotranspiration and precipitation declines are less intense when hillslope flow is included, and that is the case for the decline in total runoff as well. In regions where precipitation increases in the future, evapotranspiration increases faster, but precipitation and runoff increases slower. Effect of hillslope flow on warming is weak, but in general, the air warms up more slowly, especially in both semiarid/cold regions, and humid and warm/temperate regions. The results highlight the role of hillslope flows in increasing water exchange between the surface and atmosphere

## 1. Introduction

There is strong evidence that soil moisture plays a critical role in the evolution of climate (Seneviratne et al., 2010). For instance, at the seasonal scale, the GLACE project (Koster et al., 2006) revealed that soil moisture anomalies in the Northern Hemisphere, the Sahel and equatorial Africa can heavily affect seasonal values of precipitation and temperature (Koster, 2004). In Europe, soil moisture plays a role in the northward propagation of Mediterranean drought and in the occurrence of extreme summertime temperatures (Quesada et al., 2012; Zampieri et al., 2009); preexisting dry soil conditions are a necessary prerequisite for the occurrence of heat waves (Horton
et al., 2016; Perkins, 2015). In the United States, irrigation-induced changes in soil moisture could help explain part of the modeling bias in summer temperature and precipitation (Al-Yaari et al., 2019; Barlage et al., 2021).

Further efforts in the GLACE-CMIP5 project analyzed the feedbacks of climate change and land-atmosphere coupling. Early results revealed that the projected decrease in soil moisture would result in a stronger decrease in evapotranspiration and precipitation and a higher increase in temperature at the end of the 21st century (Seneviratne et al., 2013). Other results pointed to feedback between decreasing soil moisture and land surface temperature, relative humidity, and precipitation, thus explaining an increase in aridity under climate change (Berg et al., 2016). By the end of the 21st century, in some regions, decreasing trends in soil moisture due to climate change can also lead to further increases in the intensity, frequency, and duration of temperature extremes and in dry precipitation extremes (Lorenz et al., 2016).

For land-atmosphere coupling, knowledge gaps exist regarding the effect of topography heterogeneity on energy and water fluxes between land and atmosphere. The topographic gradient induces a hillslope flow of surface water and groundwater between the upland area and the lowland valley (Fan et al., 2019), which can create a wetter riparian wetland (Fan & Miguez-Macho, 2011). The role of groundwater is especially important for hillslope flow, because groundwater is the largest continental reservoir (Gleeson et al., 2016), and its slow flow influences a longer soil moisture memory in areas with shallow water tables (Cuthbert et al., 2019; Gleeson et al., 2011; Martinez et al., 2016a, 2016b; Martínez-De La Torre & Miguez-Macho, 2019). In water-limited regions, wetter soil induced by groundwater increases evapotranspiration (ET) (Fan et al., 2019; Maxwell & Condon, 2016). Higher ET directly affects surface water and energy budgets and can lead to increased precipitation as response from the atmosphere (Lo & Famiglietti, 2011; Wang et al., 2018).

Modeling is a good option for disentangling the link between hillslope flow and land-atmosphere fluxes. Coupling a land surface model (LSM) to a general circulation model (GCM) allows us to explore the evolution of the land and atmospheric components at the global (Wang et al., 2018) and regional scales (Anyah et al., 2008; Campoy et al., 2013; Fan et al., 2007). In LSMs, the representation of groundwater storage and its contribution to river discharge is now commonly included, with multiple approaches from a simple reservoir to physically based representations with water table dynamics (Gleeson et al., 2021). If implemented at a high enough spatial resolution, the latter approaches allow one to explicitly simulate hillslope flow and its influence on soil moisture heterogeneities, but they remain very rare in climate models, and only in regional ones (Anyah et al., 2008; Furusho-Percot et al., 2019). In coarse resolution LSMs classically coupled to atmospheric models, hillslope processes remain a challenge (Clark et al., 2015). Most attempts make use of the TOPMODEL formalism (Band et al., 1993; Beven & Kirby, 1979), either in a diagnostic mode (Gedney & Cox, 2003), or with a full coupling of subgrid water table depth distribution with the ones of soil moisture, runoff, and ET (Ducharne et al., 2000; Koster et al., 2000; Walko et al., 2000). These models, however, may overestimate the moistening of soils by capillary fluxes from the water table, as it is assumed to always be shallow in TOPMODEL (Beven et al., 2021; Gascoin et al., 2009). Interactions with deeper groundwater systems can also be described in climate models, usually via 1D capillary fluxes to the soil, which overlooks the modulations by small scale (subgrid) topography. The most comprehensive description of groundwater-soil moisture interactions in a global climate model is presently offered by the ISBA-CTRIP LSM (Decharne et al., 2019)), combining 2D horizontal groundwater flow between grid cells with vertical capillary rise in the fraction of each grid cell with the lowest elevation.

Therefore, the feedback of hillslope flow on climate change projections is not clear. Climate change, especially warming, is expected to have an impact on groundwater recharge and storage (Markovich et al., 2016; Smerdon, 2017; Wu et al., 2020). Such changes in groundwater storage and hillslope flow may support higher ET rates in transition zones between wet and dry climates for some time, but not indefinitely (Condon et al., 2020). However, the evidence provided is regional, and the modeling efforts overlook dynamical interactions with the atmosphere. Because land-surface coupling has an impact on the projection of climate change, the reorganization of soil moisture induced by hillslope flow may play a role in projections of some hydroclimatic variables on a global scale.

Here, we present evidence of the influence of hillslope flow on global climate change projections of hydroclimatic variables. We contrast the results of two transient land-atmosphere simulations for the period 1979–2100, one of which uses the default land surface representation of the ORCHIDEE LSM (Section 2.1.1), while the other includes a novel subgrid hillslope flow parameterization (Section 2.1.2). As described in Section 2.2, we use LMDZOR, the coupled land-atmosphere component of the IPSL-CM6 climate model from the Institut Pierre
Simon Laplace (Boucher et al., 2020; Cheruy et al., 2020). First, we focus the analysis on the historical period (1980–2010) to explore the effects of hillslope flow on yearly and seasonal average values in terms of both sensitivity and realism against observations (Section 3.1). We then explore the effect of hillslope flow on climate change trends, which we call hillslope flow modulation. To calculate the hillslope flow modulation, we first estimate the long-term trend for each simulation during the 21st century; then the trend of the difference of the two simulations at the global scale in Section 3.2; and finally, the spatial distribution of the hillslope flow modulation in Section 3.3. To further analyze hillslope flow modulation at the regional scale, we use a simple climate classification inspired by the Köppen-Geiger and Thornthwaite classifications (Beck et al., 2018; Feddema, 2005) in Section 3.4. The latter allows for defining which regions are prone to hillslope flow modulation and what is the type of modulation in these areas. In Section 4, we discuss the main limitations of our results, and we close in Section 5 with the main conclusions and perspectives.

## 2. Materials and Methods

### 2.1. Description of the Land Surface Model

#### 2.1.1. ORCHIDEE-REF

ORCHIDEE (Organizing Carbon and Hydrology in Dynamic Ecosystems) is a process-based model that describes the fluxes of mass, momentum, and heat between the surface and the atmosphere (Krinner et al., 2005). The version used here as a reference (called ORCHIDEE-REF in the following) corresponds to version 2.0, included in the IPSL-CM6 climate model for CMIP6 simulations and described in Cheruy et al. (2020), Boucher et al. (2020), and Tafasca et al. (2020). Here, we summarize the main characteristics of the model. In Section 2.1.2, we summarize the hillslope flow parametrization.

In each grid cell, vegetation is represented by a mosaic of up to 15 plant functional types (PFTs), including bare soil. Each PFT is characterized by a set of parameters (Boucher et al., 2020; Mizuochi et al., 2021), and fractions are described by the LUHv2 data set (Lurton et al., 2020). Plant phenology is controlled by the STOMATE module, which couples photosynthesis and the carbon cycle and computes the evolution of the leaf area index (LAI) (Krinner et al., 2005). It means that CO₂ influences plant growth and phenology. PFTs are grouped into three soil columns according to their physiological behavior: high vegetation (forest, eight PFTs), low vegetation (grasses and crops, six PFTs), and bare soil. A separate water budget is calculated independently for each soil column, in order to prevent forest PFTs from depriving the other PFTs of soil moisture but within each soil column, the uptake of water for transpiration considers the root distribution of the corresponding PFT (de Rosnay et al., 2002). In contrast, the energy balance is calculated for the whole grid cell (Boucher et al., 2020).

Evapotranspiration is represented by a classical bulk aerodynamic approach with four subfluxes: snow sublimation, interception loss, bare soil evaporation, and transpiration. The first two proceed at potential rates and originate from the snow-covered area and from the PFT fraction effectively covered by foliage. Bare soil evaporation is limited by upward water diffusion. Transpiration originates from the PFT fractions effectively covered by foliage, with no intercepted water, and is controlled by stomatal resistance, which depends on soil moisture and vegetation parameters.

Vertical soil water flow is represented by a 1-D Richards equation (Figure 1a) coupled to a mass balance (de Rosnay et al., 2002), and soil depth is set to 2 m, here discretized into 22 layers to finely model lower layers implicated in drainage (Campoy et al., 2013; Wang et al., 2018). Surface infiltration is limited by the hydraulic conductivity of surface layers and is represented by a sharp wetting front based on the Green and Ampt model (D'Orgeval et al., 2008; Tafasca et al., 2020). The resulting increase in top soil moisture is delivered to the Richards redistribution scheme as a boundary condition, while the lower one consists in free drainage equal to the hydraulic conductivity of the deepest node. Lateral fluxes between cells are neglected (Campoy et al., 2013; de Rosnay et al., 2002).

Soil is assumed to be homogeneous inside every grid cell and represented by the dominant USDA soil texture, as taken from the Zobler (1986) map. Soil parameters are a function of the soil texture, following Carsel and Parrish (1988), while the unsaturated values of hydraulic conductivity and diffusivity depend on soil moisture using the Van Genuchten-Mualem model (Mualem, 1976; van Genuchten, 1980). Even if the soil texture is assumed to be uniform inside every grid cell, the saturated hydraulic conductivity decreases with depth, following D'Orgeval et al. (2008).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematization of the main hydrological fluxes in ORCHIDEE-REF (a) and ORCHIDEE-HSL (b), blue arrows represent surface and subsurface flows, orange arrows represent fluxes to the atmosphere. Map of lowland fraction from Tootchi et al. (2019) interpolated to ORganizing Carbon and Hydrology in Dynamic EcosystEmS (ORCHIDEE) resolution (c), corresponding to 2.5 × 1.3° grid size.</figcaption></figure>

The routing scheme (Figure 1a) transfers surface runoff and drainage from land to the ocean through a cascade of linear reservoirs (Guimberteau et al., 2012; Ngo-Duc et al., 2007). Each grid cell is split into subbasins according to a flow direction map based on the work of Vörösmarty et al. (2000), enhanced over the polar regions by Oki et al. (1999) with a resolution of 0.5°. Every subbasin includes three reservoirs (Figure 1a), corresponding to groundwater, overland storage, and stream storage, with decreasing residence times (Ngo-Duc et al., 2007). The groundwater reservoir collects drainage from the soil column, while the overland reservoir collects surface runoff. There is no feedback from these two reservoirs on soil moisture. In particular, there is no capillary rise from the groundwater reservoir, so that groundwater cannot influence the atmosphere through the soil column. These two reservoirs are internal to each subbasin, and they both feed the stream reservoir, which also collects streamflow from the upstream subbasins, hence contributing to large-scale routing across subbasins and grid cells.

Eventually, the surface energy and water budget are computed at the same time step as the atmospheric model, that is, 15 min. We impose the same time step on the routing scheme for consistency with ORCHIDEE-HSL. In contrast, the carbon and plant phenology processes in STOMATE are solved with a daily time step.
### 2.1.2. ORCHIDEE-HSL

ORCHIDEE-HSL describes the effect of hillslope flow along topography at the subgrid scale. To this end, it introduces a “lowland” fraction, which buffers the flow between the upland fraction and the river system (Figure 1b) and has its own water budget, while the remaining “upland” fraction behaves like the soil of the entire grid cell in the reference version. Surface runoff and drainage replenish overland and groundwater reservoirs, which now feed the lowland fraction instead of the stream reservoir, but if there is no lowland fraction within the grid cell, overland and groundwater flows pass directly to the river, as in ORCHIDEE-REF. The lowland fraction can thus be seen as a topographically driven riparian wetland (Fan & Miguez-Macho, 2011), fed by the convergence of both surface water and groundwater, with a potentially higher evapotranspiration than the upland fraction. It is worth noting that this approach is equivalent to the use of the representative hillslope concept (Fan et al., 2019; Swenson et al., 2019), and is suitable for large-scale modeling efforts with redistribution at subgrid scale, but not for high-resolution simulations, for which intercell flow (2D or even 3D) is necessary (Felfelani et al., 2021; Markovich et al., 2016).

In contrast to many LSMs describing topographically driven wetlands with a variable area, such as in TOPMODEL (Band et al., 1993; Beven & Kirby, 1979), we suppose here for simplicity that the lowland fraction remains constant in a grid cell, and only its soil moisture changes with time. In this framework, the lowland fraction is prescribed from a 500-m resolution global-scale wetland map recently designed for this purpose (Tootchi et al., 2019). It combines open-water and inundation imagery and high-resolution groundwater modeling (Fan et al., 2013) and is interpolated to the ORCHIDEE resolution according to Section 2.2 (Figure 1c).

Even if some vegetation types may be favored/prevented in lowland areas, this is overlooked due to the lack of guiding rules to do otherwise (Fan et al., 2019). Therefore, the lowland fraction has the same PFT composition as the upland fraction (and as the full grid cell in the reference configuration). A given PFT, however, undergoes weaker water stress in the wetter lowland fraction, which therefore produces higher transpiration and a higher LAI. Bare soil evaporation is also enhanced by the higher soil moisture in the lowland fraction.

In this fraction, the overland flow from the upland fraction is added to throughfall and snowmelt. Depending on soil moisture and hydraulic conductivity, the resulting amount of water can either infiltrate into the soil or produce surface runoff, which directly flows to the river. The flow from the groundwater reservoir is injected at the bottom of the soil column, which is supposed to be impermeable to vertical drainage, so that the deep layers can gradually saturate, thus forming a water table that drains horizontally to the stream reservoir as baseflow. The water table is practically defined as the top of the uppermost saturated soil layer when starting from the impermeable soil bottom. If the whole soil column becomes saturated, the excess water will add to surface runoff and flow into the river (Tootchi, 2019). The soil depth is kept at 2 m in the lowland fraction, and it is discretized into 22 layers (with increasing height from 1 mm in the top layer to 12.5 mm in the eighth layer and all layers below) to accurately simulate the water table depth and the overlying soil moisture gradients and resulting water fluxes according to the Richards equation (Campoy et al., 2013). This discretization is also imposed for consistency in the upland fraction and in the full grid cell of ORCHIDEE-REF.

Baseflow $Q_{base}$ (mm/s) to the streams comes from the lowland fraction. Following Darcy's law, it originates from the saturated layers below the water table, and is given by a solution of the long-term linearized Boussinesq equation at the catchment scale (Brutsaert, 2005; Tootchi, 2019) in equivalent water depth:

$$Q_{base} = \sum_{i=wtl}^{22} q_i = E_F * \frac{\pi^2}{4} * \frac{\Delta h}{B} * \frac{1}{B} * \sum_{i=wtl}^{22} \Delta z_i * K_i \quad (1)$$

where $E_F$ [-] is the exchange factor multiplied by $\pi^2/4$, which accounts for variations of hydraulic conductivity along horizontal direction, anisotropy and riverbed clogging, based on water table shape (Brutsaert, 2005). The term $\Delta h/B$ corresponds to the mean gradient along hillslopes in the lowland fraction. It depends on $\Delta h$ [L], the height of water table above the bottom of the soil column, that is, 2 m depth, and on $B$ [L], the mean aquifer breadth from the streams to the divides. This term is equal to $-1/(2*\delta)$, where $\delta$ [L⁻¹] is the drainage density in the lowland fraction, assumed here to be the same as in the entire grid cell by lack of specific information. The term $\Delta z_i * K_i$ corresponds to the layer transmissivity, dependent on the layer thickness $\Delta z_i$ [L] and the layer hydraulic conductivity $K_i$ [L/T]. Finally $wtl$ is the layer that contains the water table.
Both higher values of the exchange $E_F$ factor and of the drainage density $\delta$ increase the baseflow rate: the former because it retains water in the lower zone for a shorter period, and the latter because it decreases the width of the aquifer, which facilitates the flow toward the river. In our simulation, $E_F$ is set equal to 1, and $\delta$ is set to $0.535 \times 10^{-3} \text{ m}^{-1}$, the average over the Seine river basin based on the global $\delta$ map of Schneider et al. (2017a).

## 2.2. Coupled Simulation and Experiment Setup

The ORCHIDEE LSM is coupled to the LMDZ6A atmospheric model (Hourdin et al., 2020), as embedded in the IPSL-CM6A-LR climate model (Boucher et al., 2020). This land-atmosphere coupled model is often referred to as LMDZOR (Cheruy et al., 2020). Here, the resolution of LMDZOR is defined by 144 × 143 points in longitude and latitude, which correspond to $2.5 \times 1.3^\circ$, respectively, and by 79 vertical levels. The sole difference in our simulations compared to the IPSL-CM6A-LR setup is that land cover is kept constant and representative of the early 21st century conditions (year 2000) and permafrost is turned off. It should be noted that these simulations exclude irrigation.

Two coupled simulations using LMDZOR were run for the period 1979–2100: the REF simulation, with no lowland representation, and the HSL simulation, which includes the subgrid representation of hillslope flow through a lowland fraction (Section 2.1.2). Spin up was conducted for each simulation to reach storage equilibrium. The SST/SIC forcing data sets consist of bias-corrected values from a fully coupled (land-ocean-atmosphere) simulation by the CNRM-CM6-1 (Voldoire et al., 2019) climate model against observed data. The simulation was performed for ScenarioMIP (Tebaldi et al., 2021) under historical and SSP5-8.5 radiative forcing. The radiative forcing was set differently for historical (1979–2014) and future (2015–2100) periods: historical corresponds to observed values, while future corresponds to the shared socioeconomic pathway (SSP) SSP5-8.5 (O'Neill et al., 2016)), selected to obtain a strong climate change signal.

# 3. Results

## 3.1. The Impact of Hillslope Flow on Historical Climate

In this section, yearly and seasonal mean values for the period 1980–2010 were compared for both simulations. The statistical significance of the mean difference between both simulations was assessed at each pixel with Student's t test at the confidence level of 5%.

Figure 2 shows the difference between REF and HSL simulations for annual averages covering the historical period (1980–2010). Hillslope flow increases soil moisture (SM) almost everywhere (Figure 2a), except in arid and semiarid areas where there are fewer lowland fractions and changes are not detected by the statistical test. Evapotranspiration (ET) is also enhanced (Figure 2b), but the increase is limited to moisture-limited regions (Seneviratne et al., 2010). Precipitation (P) increases, but the increment is weak and clustered in a few hotspots (Figure 2c). In particular, the inclusion of hillslope flow does not change the zonal distribution of P (Figure S1a in Supporting Information S1) or the positive P bias against observed data (for yearly and seasonal values), especially in the tropics. This bias impacts ET values (Figure S1b in Supporting Information S1), which tend to be excessive in the tropics as well. In areas with a positive P bias, hillslope flow does not enhance ET, as these regions are not water-limited but energy-limited, so there is no coupling between SM and ET (Seneviratne et al., 2010).

Total runoff (R) decreases in most zones (Figure 2d) as a result of ET enhancement, except in areas where P increases (R will increase if the P rise is greater than the ET rise). Mean air temperature (Tas, Figure 2e) cools down as a result of increased ET, with the exception of Alaska, where it increases, most likely linked to internal variability of the model. The cooling down signal is weak, -0.09 °C over land, and the correlation between ET enhancement and Tas cooling does not necessarily fit well (for instance, in the northwestern USA, the Siberian Far East, or the Caucasus). The average daily maximum temperature (Tasmax) shows a stronger cooling in magnitude over a wider area than Tas (Figure 2f). The impact on the Tasmax is linked to the cooling effect of ET and to changes in the soil thermal inertia (Cheruy et al., 2017). The control that soil moisture exerts on thermal inertia may explain the greater impact on the Tasmax compared to Tas, as well as an increase in the average daily minimum temperature (Tasmin, Figure S3 in Supporting Information S1). The increase in Tasmin as a result of wetter soil is due to a reduction in both surface temperature variability and nocturnal soil cooling.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Change in mean yearly values between HSL and REF for the period 1980–2010: (a) SM, (b) ET, (c) P, (d) R, (e) Tas, (f) Tasmax. Statistical significance of the mean differences is tested at each point with a Student's test ($p = 0.05$). The areas with insignificant changes are left gray.</figcaption></figure>

Figures 3a–3d show the differences between both simulations for boreal summer. During JJA, the ET increase is stronger than for yearly values (Figure 3a), with higher values in both the Northern and Southern Hemispheres. The statistical test does not detect any change in areas such as India (where there is an increment in yearly means), and there is a decrease in some tropical areas (e.g., the northern Amazonian basin) owing to a reduction in downwelling radiation at the surface, which is consistent with the results from Wang et al. (2018). The P increment is also stronger during JJA (Figure 3b), especially in the Sahelian band, even if P decreases in small areas (such as in the Democratic Republic of Congo). R decreases (Figure 3c) following increases in ET, except in areas where the P rise counteracts the ET rise (e.g., in the transition areas between the Sahel and the wetter southern areas) or where ET decreases (north of South America). Tasmax cooling is stronger during boreal summer as well (Figure 3d) and fits well with areas where ET is enhanced, with the exception of the Siberian Far East and North China, where cooling down and ET enhancement do not match.

Comparisons of monthly climatological values are shown for Southern Hudson Bay (SHB, Figure 3e) in Canada and for Sahelian Sudanese Band (SSB, Figure 3f) in Africa. Both regions correspond to different latitudes, seasonality and climates, but both show an increase in ET and air cooling by the inclusion of hillslope flow at yearly scales. Monthly values with detectable differences according to the *t* test and 5% confidence are marked by a shaded column. In both regions, there is an increase in ET during the wetter months (July in SHB, August in SSB), but in SSB, the test detects a change during the driest months, indicating that hillslope flow keeps sustaining ET. For both regions, increases in ET do not translate into a change in the shape but rather into a change in the ranges of the ET curve. In SHB, changes in P occur during the wettest months, while in SSB, the statistical test detects changes in June but not in August when the peak takes place, even if the peak value slightly increases. These results are consistent with (Lo & Famiglietti, 2011; Wang et al., 2018), who state that a shallow water table affects the intensity and extent of rain belts and enhances local convection but does not affect the main P patterns. For SHB, the increase in P may be due to positive feedback between ET and P from increased convection, while in SSB, the changes may be due to enhanced Hadley circulation and poleward movement of the rain belts (which also explains local decreases in P during boreal winter, see Figure S2 in Supporting Information S1).
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

**Figure 3.** Same as Figure 2 but for (a) ET, (b) P, (c) R, (d) Tasmax, and boreal summer (JJA) mean values. Regional monthly multiyear values of ET, P, R, Tas, and Tasmax for Southern Hudson Bay (SHB, e) and Sahelian Sudanese Band (SSB, f) for 1980–2010. Statistical significance of the mean differences was tested for each month with Student's *t* test (*p* = 0.05). Months with significant changes are shaded in blue. In the Tasmax and Tas caption, the top blue and red lines correspond to Tasmax.

R shows larger changes between the two regions. In SHB, the R peak is higher and occurs earlier due to hillslope flow. In SSB, R decreases during the dry season and increases during the wet months, but the statistical test only detects changes during the dry season and not changes during the peak runoff months. In both regions, the decrease in R is linked to ET enhancement, which leaves less water to runoff, but for the case of SHB, the decrease is partially counteracted by an increase in P. Changes in R peak timing are explained by the characteristics of the lowland fraction. Since the lowland fraction does not have drainage, it can become saturated by groundwater flow and infiltration more easily. If the lowland becomes saturated, all excess water will then become surface runoff
and will flow directly to the river. This means that the lowland fraction may produce surface runoff sooner, with a larger peak. For the air temperature, cooling is mostly observed, but it does not occur with the same intensity all the time. Sometimes, this decrease occurs in conjunction with increases in ET (although higher ET does not necessarily impact temperature, e.g., during the hottest months in SHB). In others, it is not due to changes in ET but rather to changes in thermal inertia due to increases in SM (e.g., in February in SHB or March in SSB).

The above changes in the regional hydrology are consistent with the results obtained by Tootchi (2019) based on offline simulations in the Seine River basin (Northern France, temperate humid climate). This author also illustrates the sensitivity of the hillslope parametrization to soil depth and the exchange factor $E_F$ (Equation 1), and to the extent of the lowland fraction: when the latter increases from 5% to 75% in the Seine River basin, annual mean SM and ET increase by 21% and 6%, respectively, while annual mean R decrease by 24%.

## 3.2. Global Influence of Hillslope Flow on Climate Change Projections

### 3.2.1. Global Trends

In this section, and hereafter, long-term trends (1980–2100) were computed using Sen's Kendall slope estimator, a nonparametric method, robust enough against outliers (Burn & Hag Elnur, 2002; Perkins & Alexander, 2013; Sen, 1968). The slope was computed for each simulation and also for the difference between them to isolate the effect of hillslope flow on the climate change trend. The statistical significance of the trends was tested with the Mann-Kendall tau test at the confidence level of 5%.

Figure 4 shows the land average time series between 1980 and 2100 for both simulations and the difference between HSL and REF (Diff. in the figure). Linear long-term trends are calculated for the three time series. REF simulation includes the effect of climate change, HSL includes the effect of climate change and hillslope flow, and Diff. isolates the effect of hillslope flow on the projection. The spatial distribution of climate change (calculated as the difference between averages in future, 2070–2100, and in historical, 1980–2010) is shown in Figures S4 and S5 in Supporting Information S1 for REF. The simulation demonstrates the main characteristics of climate change: global warming and increment of the water cycle (more P, ET, and R), but with regional increases of aridity (drier soil, less P, ET, and R), for example, in Southern Europe and the Mediterranean area (with a stronger and wider increase of aridity during boreal summer).

The REF simulation shows strong warming for the Tasmax (Figure 4a) as a result of climate change (+6.2 °C per century), but when considering climate change and hillslope flow (HSL simulation), warming is slightly less important (+6 °C per century). The trend of the difference (Diff.) between HSL and REF is negative because hillslope flow slows down warming and is detected by the statistical test. The Diff. trend indicates the magnitude and direction of the hillslope flow impact in the long term. Additionally, the trend of the difference (–0.08 °C per century) does not exactly match the differences between HSL and REF trends (–0.2 °C) as a result of numerical errors in the slope calculation induced by nonlinearity, but these errors are generally small.

For ET (Figure 4b), the rate is enhanced by climate change (+0.02 mm/day/century), and it increases faster if the hillslope is included (+50.5% faster). For the case of yearly P (Figure 4c), climate change increases P at the end of the century (+0.35 mm/day/century), and hillslope flow amplifies the increment of P, that is, P trend increases faster (+6.9% faster). During the boreal summer, ET decreases (–0.03 mm/day/century) due to climate change (Figure 4d), but hillslope flow opposes this trend by enhancing ET and attenuating the reduction (–39.2%).

### 3.2.2. Hillslope Flow Modulation

With respect to the climate change trend (REF simulation), the isolated effect of hillslope flow on long-term trends (HSL-REF time series, Diff. slope in Figure 4) can adequately reflect the influence of hillslope flow on the evolution of a variable. This influence is referenced here as hillslope flow modulation and is defined in Table 1b. A single color is associated with each of the possible combinations of the hillslope flow effect and climate change trend, and we use this color matrix in further analysis. When the trend increases/decreases more rapidly, that is, REF and Diff. slopes have the same sign, we speak of amplification. When the trend increases/decreases more slowly, that is, the REF and Diff. slopes have different signs, we speak of attenuation. Finally, if the REF and HSL trends have opposite signs, indicating that the hillslope flow counteracts climate change, we speak of an inversion. A significant modulation means that the HSL-REF trend is significant according to the Mann-Kendall tau test at 5%.
<figure><img src="image_11.png" /><figcaption>**Figure 4.** Global land average for both simulations, and difference between them, for the period 1980–2100, for daily maximum temperature (a), evapotranspiration (b), and precipitation (c), at yearly scale, and evapotranspiration (d) for boreal summer. The statistical significance of the trend was tested with a Mann-Kendall test (*p* = 0.05). Significant trends are written in bold with an *. Diff. corresponds to the trend of the difference between HSL and REF.</figcaption></figure>
<table><caption>**Table 1**
Trends of Global Land Average Values for 1980–2100 and Hillslope Flow Modulation in % (a) and Matrix of Hillslope Flow Modulation With Respect to Climate Change (C.C.) Trend and Corresponding Color (b)</caption><thead><tr><th>Season</th><th>Variables</th><th>SM<br>mm/Cent.</th><th>ET<br>mm/day/Cent.</th><th>P<br>mm/day/Cent.</th><th>R<br>mm/day/Cent.</th></tr></thead><tbody><tr><td rowspan="2">Yearly</td><td>C.C. Trend</td><td>8.0</td><td>0.02</td><td>0.3</td><td>0.3</td></tr><tr><td>Modulation</td><td>23.6%</td><td>50.5%</td><td>6.9%</td><td>3.0%</td></tr><tr><td rowspan="2">DJF</td><td>C.C. trend</td><td>14.1</td><td>0.05</td><td>0.5</td><td>0.4</td></tr><tr><td>Modulation</td><td>2.6%</td><td>-2.0%</td><td>1.0%</td><td>8.8%</td></tr><tr><td rowspan="2">MAM</td><td>C.C. trend</td><td>15.6</td><td>0.05</td><td>0.3</td><td>0.4</td></tr><tr><td>Modulation</td><td>-0.07%</td><td>41.1%</td><td>2.1%</td><td>-11.2%</td></tr><tr><td rowspan="2">JJA</td><td>C.C. trend</td><td>-1.5</td><td>-0.03</td><td>0.1</td><td>0.1</td></tr><tr><td>Modulation</td><td>-279.3%</td><td>-39.2%</td><td>11.5%</td><td>0.2%</td></tr><tr><td rowspan="2">SON</td><td>C.C. trend</td><td>3.0</td><td>0.02</td><td>0.4</td><td>0.4</td></tr><tr><td>Modulation</td><td>112.0%</td><td>91.5%</td><td>10.1%</td><td>13.2%</td></tr><tr><td rowspan="2">Season</td><td>Variables</td><td>Tasmax<br>°C/Cent.</td><td>Tas<br>°C/Cent.</td><td>Tasmin<br>°C/Cent.</td><td>SHF<br>W/m<sup>2</sup>/Cent.</td></tr><tr><td>Yearly</td><td>C.C. trend</td><td>6.2</td><td>6.4</td><td>6.9</td><td>4.2</td></tr><tr><td rowspan="2">DJF</td><td>C.C. trend</td><td>-1.3%</td><td>-1.2%</td><td>-1.1%</td><td>-4.5%</td></tr><tr><td>Modulation</td><td>5.9</td><td>6.6</td><td>7.4</td><td>2.0</td></tr><tr><td rowspan="2">MAM</td><td>C.C. trend</td><td>5.4</td><td>5.7</td><td>6.2</td><td>3.8</td></tr><tr><td>Modulation</td><td>-3.9%</td><td>-3.6%</td><td>-3.0%</td><td>-11.4%</td></tr><tr><td rowspan="2">JJA</td><td>C.C. trend</td><td>6.6</td><td>6.6</td><td>6.7</td><td>6.6</td></tr><tr><td>Modulation</td><td>-0.9%</td><td>-0.8%</td><td>-0.8%</td><td>-2.0%</td></tr><tr><td rowspan="2">SON</td><td>C.C. trend</td><td>6.6</td><td>6.8</td><td>7.2</td><td>4.3</td></tr><tr><td>Modulation</td><td>-0.02%</td><td>0.2%</td><td>0.5%</td><td>-6.0%</td></tr></tbody></table><table><thead><tr><th></th><th>Negative CC</th><th>Positive CC</th></tr></thead><tbody><tr><td>Amplification</td><td>Negative HSL-REF</td><td>Positive HSL-REF</td></tr><tr><td>Attenuation</td><td>Positive HSL-REF</td><td>Negative HSL-REF</td></tr><tr><td>Inversion</td><td>Positive HSL-REF</td><td>Negative HSL-REF</td></tr></tbody></table>

Note. The statistical significance of the trend was tested with a Mann-Kendall test (p = 0.05). Significant trends are written in italic bold.
We show the modulation of hillslope flow in climate change trends in Table 1a for global land average values at yearly and seasonal scales. Each variable presents the slope of climate change and the color of the hillslope flow modulation. The percentage corresponds to the ratio of the Diff. slope (HSL-REF time series) to the climate change slope (REF simulation). It indicates the magnitude and direction of hillslope flow modulation compared to climate change impact. For the case of SM, hillslope flow amplifies a positive trend for yearly values (+24%), as well as for ET (+51%) and P (+7%), while the R increase is attenuated. Sensible heat flux (SHF) rise is attenuated. For air temperature, we observe a small reduction in warming for the Tasmax, Tas, and Tasmin but it is only detected by the statistical test for the Tasmax (−1.3%).

These modulations may change at the seasonal scale, for example, during boreal summer (JJA) for SM and ET. The enhanced water storage in the lowland fraction at all times leads to amplified SM and ET increases in boreal winter and on an annual average, and helps counteract the decrease of these variables that is induced by climate change in boreal summer. Such an increased soil moisture memory, by means of water redistribution at the landscape scale, increases the resilience of water resources and land surface fluxes to enhanced droughts caused by climate change. Another example of modulation inversion at the seasonal scale compared to the annual scale is found for R, Tas, and Tasmin. These variables all exhibit positive long-term trends whichever the season, which are usually attenuated but in boreal autumn (SON). For Tas and Tasmin these seasonal amplifications are small and insignificant, unlike the one of R which continues in boreal winter, and is probably related to the higher soil moisture in the lowland fraction, favoring runoff in the rainy season, as noted in the historical period (Section 3.1).

### 3.3. Spatial Distribution of Hillslope Flow Modulation on Climate Change Projections

Many areas show an amplification of ET increase due to hillslope flow (red color, Figure 5a), particularly in the Northern Hemisphere, especially in areas with high lowland fraction (e.g., in Southern Hudson Bay and Central Siberia, detected by the statistical test). Other areas show a reduction in ET decline (green color, like Eastern Europe), although their detection by statistical tests is less frequent. P presents a similar modulation distribution (mostly amplification of P increases or attenuation of P declines, Figure 5b), but in this case, the hillslope flow modulation is undetected in most cases by the statistical test. Geographically, the type of ET and P modulation seems to agree with each other. For instance, attenuation of a decline (in green) is seen for both variables in Southern Europe and California. Amplification of an increase (in red) also matches (in Northern Australia and the Amazonian basin).

For yearly values of R (Figure 5c), attenuation/amplification of an increase (purple and red color, respectively) are the most frequent modulations. The modulation of hillslope flow on R appears to depend on both the magnitude and direction of ET and P modulations. For example, in North Africa, a slower reduction in R (green color) matches with a slower reduction in ET (which should accelerate the R reduction) and a slower reduction in P (which should slow down the R reduction if it is greater than changes in ET). For the case of yearly Tasmax (Figure 5d), we observe mostly a slowing down of warming (purple color) rarely detected by the test. The statistical test also detects an acceleration of warming (red color) (for instance, in the USA South Atlantic (SA USA), Figure S8 in Supporting Information S1). In general, slower warming (purple) occurs in areas where ET increases faster (red color) or decreases slower (green color). Faster warming seems to match slower ET increases (purple color), such as in Southeast Asia. However, in other regions, this agreement does not match.

During the boreal summer, the distribution of hillslope flow modulation shifts for ET (Figure 6a), and attenuation of ET declines becomes more common (green color, e.g., in Eastern Europe). The inversion of the negative trend of ET (i.e., hillslope flow counteracts the signal of climate change from negative to positive, brown color) is also important in the Northern Hemisphere. For P (Figure 6b), changes in the modulation distribution are similar to ET but again are mostly undetected by the statistical test. For R (Figure 6c), attenuation of increases (purple) is still the most common modulation, but the attenuation of R decline (in green) and the inversion of positive trends to negative (in magenta) occur more often during JJA and are more likely to be detected by the statistical test. The type of modulation for R appears to remain dependent on ET and P changes, but changes in the timing and intensity of the R peak due to the presence of the lowland fraction may also affect seasonal trends. The modulation distribution during boreal winter is similar to those of yearly means for R, with the exception of Northern Russia (Figure S6 in Supporting Information S1).

For the Tasmax (Figure 6d), hillslope flow still slows down warming (purple color) during the boreal summer, but the statistical test detects this modulation almost exclusively in the Northern Hemisphere. For JJA means,
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 5.** Spatial distribution of hillslope flow modulation in yearly mean values for period 1980–2100, ET (a), P (b), R (c), and Tasmax (d). The areas with insignificant changes in hillslope flow modulation have transparent colors. Donuts show the distribution of hillslope flow modulation (outer circle) and the corresponding partitioning of the significant and insignificant modulation (inner circle). We use a Mann-Kendall test (*p* = 0.05). The table below the maps define the modulation types and corresponding colors. Regional annual average values of ET, P, and Tas and Tasmax in Southern Hudson Bay, SHB (e) and Western Europe and the Mediterranean, WEM (f) for 1980–2100. Statistical significance of the modulation was tested with a Mann-Kendall test (*p* = 0.05), and significant trends are shown by a * in the vertical titles. In the Tasmax and Tas panel, the top blue and red lines correspond to Tasmax.</figcaption></figure>

attenuation of ET declines or amplification of ET increases (green and red colors, respectively) seem to fit well with areas where warming slows down (in Europe), while amplification of ET declines and attenuation of ET increases (blue and purple colors, respectively) agree with warming acceleration. Sensible heat flux also depicts similar modulations in some regions (Figure S7 in Supporting Information S1). During winter, slower warming is still the main hillslope flow modulation, but the statistical test mostly detects warming acceleration, probably linked to long-term trends of thermal inertia, following modulations on SM (Figures S6 and S7 in Supporting Information S1).
<figure><img src="image_14.png" /><figcaption>**Figure 6.** Same as Figure 5 in boreal summer (JJA) mean values for 1980–2100. In panel b, areas in white have very low JJA precipitation over the entire period, which prevents extracting a trend.</figcaption></figure>

In SHB (Figure 5e), ET increases more rapidly due to hillslope flow (modulation detected by the statistical test). This increase in ET does not cause a change in the P projection (although P remains higher throughout the period when hillslope flow is included). The air temperature projection does not change either (although the HSL simulation continues to show lower values). In the case of Western Europe and Mediterranean, WEM (Figure 5f), ET decline is slower (but the statistical test does not detect the modulation). This attenuation of ET decline does not agree with an attenuation of P decline, but it slows down air temperature warming (especially for the Tasmax, whose modulation is detected by the statistical test).
<figure><img src="image_15.png" /><figcaption>**Figure 7.** Aridity index (a) and mean temperature (b) from REF simulation for the period 1980–2010. Simple climatic classification from both variables (c). Areas with negative P trends (P−) for the period 1980–2100 are stippled. The piechart shows the distribution of climate classes, separated among positive and negative P trends under climate change (identified by + and −). Climate classes correspond to arid (a), semiarid and warm (SA-w), humid and warm (H-w), semiarid and temperate (SA-t), humid and temperate (H-t), semiarid and cold (SA-c), humid and cold (H-c), and subpolar (SP). Refer to the text for definition of each class.</figcaption></figure>

During the boreal summer in SHB (Figure 6e), the ET trend inverts (from negative to positive) due to hillslope flow (the modulation is detected by the statistical test). The change in the ET trend does not attenuate the P reduction and does not slow down warming, but HSL simulation consistently presents higher values for P and lower values for Tas and Tasmax. Surprisingly, we detected an acceleration of warming in SHB during boreal summer, despite the increase in ET, but the hillslope flow modulation in this case is small and undetected by the statistical test. In WEM (Figure 6f), ET decline is slower (modulation detected by the statistical test). The attenuation does not attenuate P decline, but it slows down warming, and modulation is detected for both the Tas and Tasmax.

**3.4. Influence of Hillslope Flow Modulation on Regional Climate Projections**

**3.4.1. Climate Classification for Regional Analysis**

To further analyze whether climate characteristics have an impact on hillslope flow modulation, we separated regions based on a simple climate classification. This classification is based on the use of a moisture factor and a thermal factor, following Feddema (2005). The aridity index (AI) and mean annual air temperature (MAT) are chosen as moisture and thermal factors, respectively, since they are continuous variables and are simple to calculate from the outputs of climate models. We do not consider a seasonality factor or the use of other variables, as other classifications, for the sake of simplicity. The objective is to reproduce the main features of other climate classifications, such as the Köppen-Geiger classification (Beck et al., 2018) and revised Thornthwaite classification (Feddema, 2005).

Both the aridity index and the mean air temperature are calculated for the reference period 1980–2010 from the REF simulation (Figures 7a and 7b). The aridity index AI is computed following the World Atlas of Desertification (Cherlet et al., 2018), that is, $AI = P/PET$, P is the annual mean precipitation, and PET is the annual mean potential evapotranspiration. To compute the potential ET from the model's outputs, we follow Greve et al. (2019)
and Milly and Dunne (2016), which calculate PET = 0.8R<sub>N</sub>/λ, with R<sub>N</sub> being the net radiation at the surface and λ being the latent heat of vaporization.

We classify AI into three groups: arid (A), semiarid (SA) and humid (H), for AI < 0.35, 0.35 < AI < 1.125 or AI > 1.125, respectively. MAT is classified into four groups: subpolar (SP), cold (c), temperate (t), and warm (w) for *Tas* < -5, -5 < *Tas* < 7, 7 < *Tas* < 20, and *Tas* > 20°C, respectively. We group all the grid cells from the subpolar group into a single class, SP, regardless of the value of AI. For the remaining cells, class is the result of merging AI and MAT groups, except for arid Group A, which is considered a single class. The resulting classification is shown in Figure 7c. The inclusion of the effect of climate change on the water cycle is important, as there is evidence that it can take different trajectories even within the same climatic region, as in West Africa (Gaetani et al., 2020). To include the long-term effects of climate change on aridity, we divide the classes according to their long-term trend in P for the period 1980–2100. P+ is added to the class code if the trend is positive or P− if it is negative (P− is marked as * in Figure 7c) even when the P trend is not detected by the test, except for subpolar (SP) and humid cold (H-c) because negative P trends rarely occur in these areas. Finally, to isolate the effect of proximity to the ocean, we consider an additional class called “Coastal,” which corresponds to grid cells with a fraction of land less than 50%.

The final classification includes 9 classes where the P trend is not considered and 15 classes where it is consid-
ered (Figure 8). The climate classification captures the main features of the Köppen-Geiger classification, as
calculated by Beck et al. (2018), and the revised Thornthwaite classification by Feddema (2005), but there are
biases and mismatches linked to the positive P bias detected in the tropics and midlatitudes. For instance, the
climate classification does not capture the arid region of Central Australia well due to a wet bias, while the semi-
arid areas of India and the Sahel are too large. The use of HSL simulation to calculate the classification produces
a similar distribution, with small regional differences (in Australia and northern high latitudes, Figure S9 in
Supporting Information S1), because the inclusion of hillslope flow does not drastically alter the distribution of
P and ET biases.

3.4.2. Hillslope Flow Modulation at Regional Scale

For ET and P (Figures 8a–8d), we observe that hillslope flow attenuates the reduction of yearly means due to
climate change for all classes with P- (P reduces due to climate change). For ET, attenuation of yearly decline
ranges from 6% to 67%, and for P, it ranges from 16% to 48%. During summer, we observe the same modulation
with the exception of arid areas (A P-), which amplify the reduction of both ET and P, and humid warm areas
(H-w P-), where ET increase is enhanced. Some of these modulations are detected by the statistical test, and
they match well with attenuation of SM decline and attenuation of SHF increase (Figure S10 in Supporting
Information S1).

When P increases due to climate change, hillslope flow mostly enhances ET increment (except for humid temper-
ate, H-t P+, where reduction is attenuated, but absolute trend values are small). In the case of P, hillslope flow
mostly attenuates the P increase (except for humid temperate, H-t P+, where the increase is amplified) in part
because P is already higher during the historical period, so the P rate increases slower. Finally, classes subpolar
(SP) and humid cold (H-c) show an amplification of both ET and P increase, and for the case of class H-c, modu-
lation is detected by the statistical test. During summer, ET and P do not present a single type of modulation in
areas where P increases (P+), but humid cold (H-c) areas present an enhancement of ET increase (+98%), and the
subpolar (SP) class has an amplification of P increase (+20%), both detected by the statistical test.

For R, classes with P decrease due to climate change (P−) and have a reduction in negative trends (from −8 to −42%), with the exception of arid (A P−) and semiarid warm (SA-w P−) classes, which have an inversion from negative to positive R trends (but absolute values are small, even if modulation is detected for A P−). For summer means, R declines are attenuated for all classes with P−. For classes with P+, R increase is attenuated, except for humid temperate (H-t P+) class, following results for P. Summer values also closely follow results for P, except for semiarid temperate class (SA-c P+) which present an attenuation of an increase, humid cold (H-c) which reduces the R decline, and subpolar (SP), which attenuates R increase. Arid regions and semiarid/humid cold areas have detectable modulations according to the statistical test.

Tas shows a slowing of the yearly warming trend for almost every class, except for arid and semiarid temperate,
with P decline due to climate change (A P- and SA P-), and humid temperate with P increase (H-t P+). These
modulations agree with those of Tasmin, but for the Tasmax, the single class with acceleration of warming is
<figure><img src="image_16.png" /><figcaption>**Figure 8.** Trends of regional average grid cells according to climatic classes in 1980–2100 for REF and HSL, and hillslope flow modulation in % for (a, c, e, g) yearly values, (b, d, f, h) summer values. Summer is defined as JJA for northern grid cells and DJF for southern grid cells. Variables are (a, b) ET, (c, d) P, (e, f) R, and (g, h) Tas. Statistical significance of the hillslope flow modulation was tested with the Mann-Kendall test (*p* = 0.05). Significant hillslope flow modulations are written in bold with an *, and the color of the cell corresponds to the modulation type as shown at the bottom. The percentage of modulation is calculated as for Table 1a. Climate class names are the same as Figure 7.</figcaption></figure>
humid temperate with an increase in P (H-t P+, Figure S10 in Supporting Information S1). During summer, acceleration of warming becomes more common for Tas (for the Tasmax as well, Figure S10 in Supporting Information S1), but the sole detectable modulation by the statistical test is attenuation of warming, for yearly and summerly values.

These regional scale analyzes show that, where climate change decreases P rates, hillslope flow attenuates the related SM and ET declines, and the subsequent P and R declines. In regions where P increases, in contrast, hillslope flow amplifies the resulting increases of ET while it attenuates the ones of P and R. Eventually, hillslope flow is shown to weaken the projected intensification of hydrological extreme events. Our regional results also confirm the weak attenuation of warming found at the global scale.

**4. Discussion**

This study presents two coupled simulations with a single model; thus, our results may be affected by spurious effects due to internal model variability. In addition, our model represents the effect of groundwater on soil moisture and topography-driven subgrid variability in a very simplified way, with many shortcomings, including the lack of dynamics in both land use and the extension of lowland humid areas, and the lack of PFT differences in lowland and upland areas. All this forces us to consider these results as qualitative indicators of potential first order effects.

Our results can be summarized by Figure 9a on historical land average values for the main affected variables. In our simulations, hillslope flow increases SM, ET, and P at the global scale and reduces R, except in areas where the increase in P is greater than that in ET. Both offline and online studies show the same impact for SM, ET, and R if groundwater-SM interaction is activated, and online simulations present an increase in the P rate as well. Likewise, a study using LMDZOR with prescribed water table depths of one m (Wang et al., 2018) shows the same effects for SM, ET and P but with higher impact due to unrealistic water tables. The decrease in SHF and air temperature in our simulations (Figure 9a), due to enhanced ET, is also observed in online studies. For the air temperature, we detected a stronger effect on a larger area for Tasmax than for Tas when comparing HSL and REF. This effect has been observed in coupled simulations when irrigation is included (Thiery et al., 2017).

To assess whether our results carry some generality, we compare them with other assessments of groundwater impact on land surface fluxes, that is, offline simulations by Koirala et al. (2014), Decharme et al. (2019), and Miguez-Macho et al. (2007), and near-surface climate, that is, coupled land-atmosphere simulations by Lo and Famiglietti (2011), Leung et al. (2011), and Anyah et al. (2008).

Among the models compared here, only ORCHIDEE-HSL and ISBA-CTRIP (Decharme et al., 2019) include a representation of subgrid variability. Other modeling approaches, such as MAT-GW (Koirala et al., 2014) and CLM3.5 (Lo & Famiglietti, 2011), assume a capillary rise over a flat grid cell. Models that do not account for subgrid heterogeneity may show stronger effects than those with subgrid variability. For example, ET enhancement in MAT-GW is higher (9%) than its reference simulation compared to the enhancement in ORCHIDEE-HSL (4.2%) and in ISBA-CTRIP (1.1%). In any case, the effects on land surface fluxes and near-surface climate are consistent in direction across all the studies.

Other limitations in our results for land average values are related to modeling choices. First, we prescribe a static lowland fraction with many uncertainties (Tootchi et al., 2019), but this fraction is likely to be dynamic, according to the concept of variable contributing area (Beven & Kirby, 1979). Second, the deactivation of permafrost may induce a stronger signal in high latitude regions, as permafrost can hinder effective hillslope flow (Sergeant et al., 2021). Finally, we overlook losing streams; their inclusion could lead to lower ET rates in our simulations if the water table depth is below the soil column (Brunner et al., 2011; Rashid et al., 2019), especially in flat semiarid areas and/or where extensive groundwater pumping takes place (Jasechko et al., 2021), but they can also contribute to higher ET rates if the stream sustains a shallow water table Miguez-Macho et al. (2007).

For future climates, comparing our results is more difficult. To our knowledge, there are no studies that analyze the effect of hillslope flow and groundwater-SM interaction over long-term trends under climate change. Some studies assess the effect of climate change on groundwater storage, for instance, in Wu et al. (2020) for key mid-latitude aquifers with CESM-LE and CLM4.0. The coherence of our results with other studies for land average values is especially important because the former makes our analysis for future climate more robust.
a. Changes on historical averages

<figure><img src="image_17.png" /></figure>

b. Modulation of climate change trends

<figure><img src="image_18.png" /></figure>

**Figure 9.** Schematic of main impacts. Summary of the main impacts and feedbacks due to the inclusion of hillslope flow during the historical period (a). Red color refers to an increase, blue color to a decrease. Summary of the main hillslope flow modulations at the global and regional scales (b). Colors from the table of hillslope flow modulation are used for the arrows and percentage corresponds to modulation at global scale.

Figure 9b summarizes the modulations of long-term climate change trends by hillslope flow. At the regional scale, the modulation of ET trends depends strongly on the sign of P changes with climate change, which is overlooked in the offline climate change impact assessment of Condon et al. (2020) with an integrated surface-groundwater hydrologic model over the contiguous United States. Hillslope flow attenuates climate change trends of all other hydrological variables. Therefore, it has the potential to partially alleviate the intensification of droughts projected by current climate models (Cook et al., 2018; Douville et al., 2021) and already ongoing (Padrón et al., 2020). It can also attenuate the projected runoff increases in regions and seasons with increased P. These regional attenuations contrast with the amplifications that dominate the global scale response of hydrological variables (but R). It is noteworthy that these global scale amplifications, like the attenuation of global warming in our simulations, contradict the conclusions from the GLACE-CMIP5 experiment. In this framework, Berg et al. (2016) show that, on average over land, land-atmosphere feedbacks amplify anthropogenic warming and attenuate the subsequent
increases of P and ET, therefore amplifying aridity (P/PET) increases. These contrasted results at regional and global scales show the need to develop process-based understanding of climate change interactions with the water cycle.

For air temperature, we observe that hillslope flow slightly attenuates warming at both scales (Figure 9b) due to amplification of ET increase and attenuation of SHF increase. However, in some areas (Figure 8), we observe a small amplification of warming for the HSL simulation, which is difficult to explain. In arid regions, faster warming might be the result of decreasing rates of ET. Elsewhere, accelerating warming occurs even if there is an inversion of ET trend from negative to positive, for example, in the South Atlantic area of the USA (SA USA, Figure S8 in Supporting Information S1). This warming amplification could be linked to changes in net radiation at the surface, as has been reported in coupled simulations with irrigation (Thiery et al., 2017) for areas far away from irrigated regions.

Some modeling choices could limit our results for hillslope flow modulation on climate change trends. The first is the activation of permafrost, which could have an impact on the speed and strength of climate change trends in high latitude areas (Sergeant et al., 2021). The second is overlooking losing streams, which may become more common in the future due to widespread water table decline (de Graaf et al., 2019), and could prevent any impact on ET by the disconnection between groundwater and the lowland fraction. The third is land use change, which is not included in our simulation and has a complex relationship with groundwater, land surface fluxes, and climate (Zipper et al., 2019).

**5. Conclusions**

This study analyzes the influence of hillslope flow on mean climate and climate change projections. It relies on a new parameterization of subgrid-scale hillslope flow in the ORCHIDEE LSM and compares two transient coupled land-atmosphere simulations under historical and strong climate change scenarios. The impact on projections is split into mean values and long-term trends for several hydroclimatic variables (P, ET, R, SHF, and air temperature).

During the historical period, the redistribution of water along topography by hillslope flow increases SM, especially in the Northern Hemisphere, and it enhances water exchanges from the surface to the atmosphere (+4.2% in ET) and from the atmosphere to the surface (+1.2% in P). ET enhancement occurs mostly in transition areas and during dry seasons, while the P response is mainly due to local convection (as in Canada) and poleward movements of rain belts (as in the Sahel-Sudan band). As a result, R tends to decrease (-2.8%), except in areas where there is a strong P increase. The presence of the lowland fraction may also lead to changes in the R peak in some regions (such as Southern Hudson Bay). On the other hand, ET enhancement shifts the energy partitioning by decreasing SHF, which cools down the near-surface atmosphere. The effect is weak for Tas (-0.1 °C) but is stronger and wider for Tasmax (-0.2 °C). Finally, as a result of the SM increase and subsequent changes in soil thermal inertia, Tasmin can increase regionally in central Africa and some areas of the Amazonian basin.

For long-term trends over the 21st century, climate change is the leading driver, but hillslope flow can modulate climate change effects. Globally, hillslope flow enhances the increase in both P and ET and slightly slows global warming. At the regional scale, when climate change decreases P rates, hillslope flow attenuates the decline in both ET and P. This effect of hillslope flow occurs in all regions, regardless of regional aridity and mean temperature. Seasonally, hillslope flow modulation could persist, as in the Western Europe and Mediterranean region, WEM. Furthermore, because it attenuates most regional P declines, hillslope flow also attenuates the related R decrease. In contrast, when climate change increases P rates, hillslope flow enhances the increase in ET and attenuates the increase in both P and R. However, modulation at the seasonal scale may be different, depending on regional changes. Regional warming is slower when hillslope flow is included, but the signal is weak for yearly and seasonal mean values and is stronger for the Tasmax than for the Tas.

Eventually, the areas most prone to modulation of climate change trends by hillslope flow are humid areas, either warm like tropical South America and West Africa, temperate like Western Europe, or cold like Canada and northern Europe. Semiarid and cold areas (such as the northern Great Plains in North America) and humid and warm/temperate regions (such as tropical South America and Western Europe) are prone to slower warming than elsewhere when hillslope flow is included.
Our results highlight the role of hillslope flow in sustaining land surface-atmosphere water fluxes. This attenuating role (except for ET increases where P increases due to climate change) suggests that the inclusion of hillslope processes would weaken the projected intensification of droughts and aridity. The effect of hillslope flow on regional projections also suggests an attenuation in the evolution of extreme temperature events such as heat waves (Mu et al., 2021), and in hydrological extreme events. To our knowledge, only a few models include the effects of groundwater on soil moisture in CMIP6: CESM2 and its LSM CLM5 (Danabasoglu et al., 2020; Lawrence et al., 2019), NorESM2, which uses CLM5 (Seland et al., 2020), E3SM-1.1, which includes CLM4.5 (Golaz et al., 2019), and CNRM-CM6-1 with its LSM ISBA-CTRIP (Decharme et al., 2019; Voldoire et al., 2019). Our results and their implications call for including or improving the representation of hillslope flow and the related groundwater storage in LSMs such as ORCHIDEE and for incorporating these processes into state-of-the-art climate change projections.

## Data Availability Statement

The ORCHIDEE version used for the HSL simulation run for this paper is Hillslope-r6515, and is available at: https://forge.ipsl.jussieu.fr/orchidee/wiki/GroupActivities/CodeAvalabilityPublication/ORCHIDEE-Hillslope-r6515. MSWEP (Beck et al., 2017) is available at http://www.gloh2o.org/mswep/ under request, GPCP (Adler et al., 2018) is available at https://www.ncei.noaa.gov/data/global-precipitation-climatology-project-gpcp-monthly/access/, GPCC (U. Schneider, Jost, et al., 2017) is available at https://opendata.dwd.de/climate_environment/GPCC/html/download gate.html, Fluxcom (Jung et al., 2019) is available at https://www.bgc-jena.mpg.de/geodb/projects/Home.php after registration, GLEAM (Martens et al., 2017) is available at https://www.gleam.eu/ after registration. For analysis, we used standard packages from R, https://www.R-project.org/. Simulation data used in the analysis are available in Arboleda-Obando et al. (2022).

## Acknowledgments

This research is part of the PhD project of Pedro F. Arboleda, funded by the EUR IPSL CGS (ANR Investissements d'avenir, Grant Number: ANR-11-IDEX-0004-17-EURE-0006). The work also received support from projects IGEM (ANR Grant No. ANR-14-CE01-0018-01) and BLUEGEM (ANR-21-SOIL-0001). The simulations were done using the IDRIS computational facilities (Institut du Développement et des Ressources en Informatique Scientifique, CNRS, France), under the allocation 2018-[A0050107406]. We also thanks L. Band and an anonymous reviewer for their constructive comments which helped to improve the manuscript.

## References

Adler, R., Sapiano, M., Huffman, G., Wang, J.-j., Gu, G., Bolvin, D., et al. (2018). The global precipitation climatology project (GPCP) monthly analysis (new version 2.3) and a review of 2017 global precipitation. *Atmosphere*, 9(4), 138. https://doi.org/10.3390/atmos9040138

Al-Yaari, A., Ducharne, A., Cheruy, F., Crow, W. T., & Wigneron, J. P. (2019). Satellite-based soil moisture provides missing link between summertime precipitation and surface temperature biases in CMIP5 simulations over conterminous United States. *Scientific Reports*, 9(1), 1–12. https://doi.org/10.1038/s41598-018-38309-5

Anyah, R. O., Weaver, C. P., Miguez-Macho, G., Fan, Y., & Robock, A. (2008). Incorporating water table dynamics in climate modeling: 3. Simulated groundwater influence on coupled land-atmosphere variability. *Journal of Geophysical Research*, 113, D07103. https://doi.org/10.1029/2007JD009087

Arboleda-Obando, P., Ducharne, A., Jost, A., Colin, J., Nous, C., & Ghattas, J. (2022). Influence of hillslope flow on hydroclimatic evolution under climate change—Dataset. *Zenodo*. https://doi.org/10.5281/zenodo.6521865

Band, L. E., Patterson, P., Nemani, R., & Running, S. W. (1993). Forest ecosystem processes at the watershed scale: Incorporating hillslope hydrology. *Agricultural and Forest Meteorology*, 63(1–2), 93–126. https://doi.org/10.1016/0168-1923(93)0024-2

Barlage, M., Chen, F., Rasmussen, R., Zhang, Z., & Miguez-Macho, G. (2021). The importance of scale-dependent groundwater processes in land-atmosphere interactions over the Central United States. *Geophysical Research Letters*, 48, e2020GL092171. https://doi.org/10.1029/2020GL092171

Beck, H. E., Van Dijk, A. I., Levizzani, V., Schellekens, J., Miralles, D. G., Martens, B., & De Roo, A. (2017). Mswep: 3-hourly 0.25° global gridded precipitation (1979–2015) by merging gauge, satellite, and reanalysis data. *Hydrology and Earth System Sciences*, 21(1), 589–615. https://doi.org/10.5194/hess-21-589-2017

Beck, H. E., Zimmermann, N. E., McVicar, T. R., Vergopolan, N., Berg, A., & Wood, E. F. (2018). Present and future Köppen-Geiger climate classification maps at 1-km resolution. *Scientific Data*, 5, 1–12. https://doi.org/10.1038/sdata.2018.214

Berg, A., Findell, K., Lintner, B., Giannini, A., Seneviratne, S. I., Van Den Hurk, B., et al. (2016). Land-atmosphere feedbacks amplify aridity increase over land under global warming. *Nature Climate Change*, 6(9), 869–874. https://doi.org/10.1038/nclimate3029

Beven, K. J., & Kirby, M. J. (1979). A physically based, variable contributing area model of basin hydrology. *Hydrological Sciences Bulletin*, 24(1), 43–69. https://doi.org/10.1080/02626667909491834

Beven, K. J., Kirkby, M. J., Freer, J. E., & Lamb, R. (2021). A history of TOPMODEL. *Hydrology and Earth System Sciences*, 25(2), 527–549. https://doi.org/10.5194/hess-25-527-2021

Boucher, O., Servonnat, J., Albright, A. L., Aumont, O., Balkanski, Y., Bastrikov, V., et al. (2020). Presentation and evaluation of the IPSL-CM6A-LR climate model. *Journal of Advances in Modeling Earth Systems*, 12(7), 1–52. https://doi.org/10.1029/2019MS002010

Brunner, P., Cook, P. G., & Simmons, C. T. (2011). Disconnected surface water and groundwater: From theory to practice. *Ground Water*, 49(4), 460–467. https://doi.org/10.1111/j.1745-6584.2010.00752.x

Brutsaert, W. (2005). *Hydrology* (Vol. 43, No. 7). Cambridge University Press. https://doi.org/10.1017/CBO9780511808470

Burn, D. H., & Hag Elnur, M. A. (2002). Detection of hydrologic trends and variability. *Journal of Hydrology*, 255(1–4), 107–122. https://doi.org/10.1016/S0022-1694(01)00514-5

Campoy, A., Ducharne, A., Cheruy, F., Hourdin, F., Polcher, J., & Dupont, J. C. (2013). Response of land surface fluxes and precipitation to different soil bottom hydrological conditions in a general circulation model. *Journal of Geophysical Research: Atmospheres*, 118, 10725–10810. https://doi.org/10.1002/jgrd.50627

Carsel, R. F., & Parrish, R. S. (1988). Developing joint probability distributions of soil water retention characteristics. *Water Resources Research*, 24(5), 755–769. https://doi.org/10.1029/WR024i005p00755
Cherlet, M., Hutchinson, C., Reynolds, J., Hill, J., Sommer, S., & von Maltitz, G. (2018). *World Atlas of desertification* (3rd ed.). Publication Office of the European Union. https://doi.org/10.2760/06292

Cheruy, F., Ducharne, A., Hourdin, F., Musat, I., Vignon, E., Gastineau, G., et al. (2020). Improved near surface continental climate in IPSL-CM6A-LR by combined evolutions of atmospheric and land surface physics. *Journal of Advances in Modeling Earth Systems, 12*, e2019MS002005. https://doi.org/10.1029/2019MS002005

Cheruy, F., Dufresne, J. L., Ait Mesbah, S., Grandpeix, J. Y., & Wang, F. (2017). Role of soil thermal inertia in surface temperature and soil moisture-temperature feedback. *Journal of Advances in Modeling Earth Systems, 9*, 2906–2919. https://doi.org/10.1002/2017MS001036

Clark, M. P., Fan, Y., Lawrence, D. M., Adam, J. C., Bolster, D., Gochis, D. J., et al. (2015). Improving the representation of hydrological processes in earth system models. *Water Resources Research, 51*, 5929–5956. https://doi.org/10.1002/2015WR017096

Condon, L. E., Atchley, A. L., & Maxwell, R. M. (2020). Evapotranspiration depletes groundwater under warming over the contiguous United States. *Nature Communications, 11*(1), 873. https://doi.org/10.1038/s41467-020-14688-0

Cook, B. I., Mankin, J. S., & Anchukaitis, K. J. (2018). Climate change and drought: From past to future. *Current Climate Change Reports, 4*(2), 164–179. https://doi.org/10.1007/s40641-018-0093-2

Cuthbert, M. O., Gleeson, T., Moosdorf, N., Befus, K. M., Schneider, A., Hartmann, J., & Lehner, B. (2019). Global patterns and dynamics of climate–groundwater interactions. *Nature Climate Change, 9*(2), 137–141. https://doi.org/10.1038/s41558-018-0386-4

D'Orgeval, T., Polcher, J., & Rosnay, P. D. (2008). Sensitivity of the West African hydrological cycle in ORCHIDEE to infiltration processes. *Hydrology and Earth System Sciences, 12*, 1387–1401.

Danabasoglu, G., Lamarque, J. F., Bacmeister, J., Bailey, D. A., DuVivier, A. K., Edwards, J., et al. (2020). The Community earth system model version 2 (CESM2). *Journal of Advances in Modeling Earth Systems, 12*, e2019MS001916. https://doi.org/10.1029/2019MS001916

de Graaf, I. E. M., Gleeson, T., (Rens, L. P. H.), van Beek, L. P. H., Sutanudjaja, E. H., & Bierkens, M. F. P. (2019). Environmental flow limits to global groundwater pumping. *Nature, 574*(7776), 90–94. https://doi.org/10.1038/s41586-019-1594-4

de Rosnay, P., Polcher, J., Bruen, M., & Laval, K. (2002). Impact of a physically based soil water flow and soil-plant interaction representation for modeling large-scale land surface processes. *Journal of Geophysical Research, 107*(D11), 4118. https://doi.org/10.1029/2001JD000634

Decharme, B., Delire, C., Minvielle, M., Colin, J., Vergnes, J., Alias, A., et al. (2019). Recent changes in the ISBA-CTRIP land surface system for use in the CNRM-CM6 climate model and in global off-line hydrological applications. *Journal of Advances in Modeling Earth Systems, 11*, 1207–1252. https://doi.org/10.1029/2018MS001545

Douville, H., Raghavan, K., Renwick, J., Allan, R. P., Arias, P. A., Barlow, M., et al. (2021). Water cycle changes. In V. Masson-Delmotte, P. Zhai, A. Pirani, S. L. Connors, C. Pean, T. K. M. S. Berger, et al. (Eds.), *Climate change 2021: The physical science basis. contribution of working group I to the sixth assessment report of the intergovernmental panel on climate change*. Cambridge University Press.

Ducharne, A., Koster, R. D., Suarez, M. J., Stieglitz, M., & Kumar, P. (2000). A catchment-based approach to modeling land surface processes in a general circulation model: 2. Parameter estimation and model demonstration. *Journal of Geophysical Research, 105*(D20), 24823–24838. https://doi.org/10.1029/2000JD900328

Fan, Y., Clark, M., Lawrence, D. M., Swenson, S., Band, L. E., Brantley, S. L., et al. (2019). Hillslope hydrology in global change research and earth system modeling. *Water Resources Research, 55*, 1737–1772. https://doi.org/10.1029/2018WR023903

Fan, Y., Li, H., & Miguez-Macho, G. (2013). Global patterns of groundwater table depth. *Science, 339*(6122), 940–943. https://doi.org/10.1126/science.1229881

Fan, Y., & Miguez-Macho, G. (2011). A simple hydrologic framework for simulating wetlands in climate and earth system models. *Climate Dynamics, 37*(1–2), 253–278. https://doi.org/10.1007/s00382-010-0829-8

Fan, Y., Miguez-Macho, G., Weaver, C. P., Walko, R., & Robock, A. (2007). Incorporating water table dynamics in climate modeling: 1. Water table observations and equilibrium water table simulations. *Journal of Geophysical Research, 112*, D10125. https://doi.org/10.1029/2006JD008111

Feddersen, J. J. (2005). A revised Thorntwaite-type global climate classification. *Physical Geography, 26*(6), 442–466. https://doi.org/10.2747/0272-3646.26.6.442

Felfelani, F., Lawrence, D. M., & Pokhrel, Y. (2021). Representing intercell lateral groundwater flow and aquifer pumping in the community land model. *Water Resources Research, 57*, e2020WR027531. https://doi.org/10.1029/2020WR027531

Furusho-Percot, C., Goergen, K., Hartick, C., Kulkarni, K., Keune, J., & Kollet, S. (2019). Pan-European groundwater to atmosphere terrestrial systems climatology from a physically consistent simulation. *Scientific Data, 6*(1), 320. https://doi.org/10.1038/s41597-019-0328-7

Gaetani, M., Janicot, S., Vrac, M., Famien, A. M., & Sultan, B. (2020). Robust assessment of the time of emergence of precipitation change in West Africa. *Scientific Reports, 10*(1), 7670. https://doi.org/10.1038/s41598-020-63782-2

Gascoin, S., Ducharne, A., Ribstein, P., Carli, M., & Habets, F. (2009). Adaptation of a catchment-based land surface model to the hydrogeological setting of the Somme River basin (France). *Journal of Hydrology, 368*(1–4), 105–116. https://doi.org/10.1016/j.jhydrol.2009.01.039

Gedney, N., & Cox, P. M. (2003). The sensitivity of global climate model simulations to the representation of soil moisture heterogeneity. *Journal of Hydrometeorology, 4*(6), 1265–1275. https://doi.org/10.1175/1525-7541(2003)004<1265:TSOGCM>2.0.CO;2

Gleeson, T., Befus, K. M., Jasechko, S., Luijendijk, E., & Cardenas, M. B. (2016). The global volume and distribution of modern groundwater. *Nature Geoscience, 9*(2), 161–167. https://doi.org/10.1038/ngeo2590

Gleeson, T., Marklund, L., Smith, L., & Manning, A. H. (2011). Classifying the water table at regional to continental scales. *Geophysical Research Letters, 38*, L05401. https://doi.org/10.1029/2010GL046427

Gleeson, T., Wagener, T., Döll, P., Zipper, S., West, C., Wada, Y., et al. (2021). GMD perspective: The quest to improve the evaluation of groundwater representation in continental to global scale models. *Geoscientific Model Development Discussions, 14*, 7545–7571. https://doi.org/10.5194/gmd-2021-97

Golaz, J.-c., Caldwell, P. M., Van Roekel, L. P., Petersen, M. R., Tang, Q., Wolfe, J. D., et al. (2019). The DOE E3SM coupled model version 1: Overview and evaluation at standard resolution. *Journal of Advances in Modeling Earth Systems, 11*, 2089–2129. https://doi.org/10.1029/2018MS001603

Greve, P., Roderick, M. L., Ukkola, A. M., & Wada, Y. (2019). The aridity Index under global warming. *Environmental Research Letters, 14*(12), 124006. https://doi.org/10.1088/1748-9326/ab5046

Guimberteau, M., Drapeau, G., Ronchail, J., Sultan, B., Polcher, J., Martinez, J.-M., et al. (2012). Discharge simulation in the sub-basins of the Amazon using ORCHIDEE forced by new datasets. *Hydrology and Earth System Sciences, 16*(3), 911–935. https://doi.org/10.5194/hess-16-911-2012

Horton, R. M., Mankin, J. S., Lesk, C., Coffel, E., & Raymond, C. (2016). A review of recent advances in research on extreme heat events. *Current Climate Change Reports, 2*(4), 242–259. https://doi.org/10.1007/s40641-016-0042-x

Hourdin, F., Rio, C., Grandpeix, J. Y., Madeleine, J. B., Cheruy, F., Rochetin, N., et al. (2020). LMDZ6A: The atmospheric component of the IPSL climate model with improved and better tuned physics. *Journal of Advances in Modeling Earth Systems, 12*, e2019MS001892. https://doi.org/10.1029/2019MS001892
Jasechko, S., Seybold, H., Perrone, D., Fan, Y., & Kirchner, J. W. (2021). Widespread potential loss of streamflow into underlying aquifers across the USA. *Nature*, **591**(7850), 391–395. https://doi.org/10.1038/s41586-021-03311-x

Jung, M., Koirala, S., Weber, U., Ichii, K., Gans, F., Camps-Valls, G., et al. (2019). The FLUXCOM ensemble of global land-atmosphere energy fluxes. *Scientific Data*, **6**(1), 74. https://doi.org/10.1038/s41597-019-0076-8

Koirala, S., Yeh, P. J., Hirabayashi, Y., Kanae, S., & Oki, T. (2014). Global-scale land surface hydrologic modeling with the representation of water table dynamics. *Journal of Geophysical Research: Atmospheres*, **119**, 75–89. https://doi.org/10.1002/2013JD020398

Koster, R. D., Dirmeyer, P. A., Guo, Z., Bonan, G., Chan, E., Cox, P., et al. (2004). Regions of strong coupling between soil moisture and precipitation. *Science*, **305**(5687), 1138–1140. https://doi.org/10.1126/science.1100217

Koster, R. D., Suarez, M. J., Ducharne, A., Stieglitz, M., & Kumar, P. (2000). A catchment-based approach to modeling land surface processes in a general circulation model: 1. Model structure. *Journal of Geophysical Research*, **105**(D20), 24809–24822. https://doi.org/10.1029/2000JD900327

Koster, R. D., Sud, Y. C., Guo, Z., Dirmeyer, P. A., Bonan, G., Oleson, K. W., et al. (2006). Glace: The global land-atmosphere coupling experiment. Part I: Overview. *Journal of Hydrometeorology*, **7**(4), 590–610. https://doi.org/10.1175/JHM510.1

Krinner, G., Viovy, N., de Noblet-Ducoudré, N., Ogée, J., Polcher, J., Friedlingstein, P., et al. (2005). A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system. *Global Biogeochemical Cycles*, **19**, GB1015. https://doi.org/10.1029/2003GB002199

Lawrence, D. M., Fisher, R. A., Koven, C. D., Oleson, K. W., Swenson, S. C., Bonan, G., et al. (2019). The Community land model version 5: Description of new features, benchmarking, and impact of forcing uncertainty. *Journal of Advances in Modeling Earth Systems*, **11**, 4245–4287. https://doi.org/10.1029/2018MS001583

Leung, L. R., Huang, M., Qian, Y., & Liang, X. (2011). Climate-soil-vegetation control on groundwater table dynamics and its feedbacks in a climate model. *Climate Dynamics*, **36**(1–2), 57–81. https://doi.org/10.1007/s00382-010-0746-x

Lo, M.-H., & Famiglietti, J. S. (2011). Precipitation response to land subsurface hydrologic processes in atmospheric general circulation model simulations. *Journal of Geophysical Research*, **116**, D05107. https://doi.org/10.1029/2010JD015134

Lorenz, R., Argüeso, D., Donat, M. G., Pitman, A. J., Hurk, B. V. D., Berg, A., et al. (2016). Influence of land-atmosphere feedbacks on temperature and precipitation extremes in the GLACE-CMIP5 ensemble. *Journal of Geophysical Research: Atmospheres*, **121**, 607–623. https://doi.org/10.1002/2015JD024053

Lurton, T., Balkanski, Y., Bastrikov, V., Bekki, S., Bopp, L., Braconnot, P., et al. (2020). Implementation of the CMIP6 forcing data in the IPSL-CM6A-LR model. *Journal of Advances in Modeling Earth Systems*, **12**, e2019MS001940. https://doi.org/10.1029/2019MS001940

Markovich, K. H., Maxwell, R. M., & Fogg, G. E. (2016). Hydrogeological response to climate change in alpine hillslopes. *Hydrological Processes*, **30**(18), 3126–3138. https://doi.org/10.1002/hyp.10851

Martens, B., Miralles, D. G., Lievens, H., van der Schalie, R., de Jeu, R. A. M., Fernández-Prieto, D., et al. (2017). GLEAM v3: Satellite-based land evaporation and root-zone soil moisture. *Geoscientific Model Development*, **10**(5), 1903–1925. https://doi.org/10.5194/gmd-10-1903-2017

Martinez, J. A., Dominguez, F., & Miguez-Macho, G. (2016a). Effects of a groundwater scheme on the simulation of soil moisture and evaporation over southern South America. *Journal of Hydrometeorology*, **17**(11), 2941–2957. https://doi.org/10.1175/JHM-D-16-0051.1

Martinez, J. A., Dominguez, F., & Miguez-Macho, G. (2016b). Impacts of a groundwater scheme on hydroclimatological conditions over southern South America. *Journal of Hydrometeorology*, **17**(11), 2959–2978. https://doi.org/10.1175/JHM-D-16-0052.1

Martínez-De La Torre, A., & Miguez-Macho, G. (2019). Groundwater influence on soil moisture memory and land-atmosphere fluxes in the Iberian Peninsula. *Hydrology and Earth System Sciences*, **23**(12), 4909–4932. https://doi.org/10.5194/hess-23-4909-2019

Maxwell, R. M., & Condon, L. E. (2016). Connections between groundwater flow and transpiration partitioning. *Science*, **353**(6297), 377–380. https://doi.org/10.1126/science.aaf7891

Miguez-Macho, G., Fan, Y., Weaver, C. P., Walko, R., & Robock, A. (2007). Incorporating water table dynamics in climate modeling: 2. Formulation, validation, and soil moisture simulation. *Journal of Geophysical Research*, **112**, 2006JD008112. https://doi.org/10.1029/2006JD008112

Milly, P. C., & Dunne, K. A. (2016). Potential evapotranspiration and continental drying. *Nature Climate Change*, **6**(10), 946–949. https://doi.org/10.1038/nclimate3046

Mizuochi, H., Ducharne, A., Cheruy, F., Ghattas, J., Al-Yaari, A., Wigneron, J.-P., et al. (2021). Multivariable evaluation of land surface processes in forced and coupled modes reveals new error sources to the simulated water cycle in the IPSL (Institute Pierre Simon Laplace) climate model. *Hydrology and Earth System Sciences*, **25**(4), 2199–2221. https://doi.org/10.5194/hess-25-2199-2021

Mu, M., De Kauwe, M. G., Ukkola, A. M., Pitman, A. J., Guo, W., Hobeichi, S., & Briggs, P. R. (2021). Exploring how groundwater buffers the influence of heatwaves on vegetation function during multi-year droughts. *Earth System Dynamics*, **12**(3), 919–938. https://doi.org/10.5194/esd-12-919-2021

Mualem, Y. (1976). A new model for predicting the hydraulic conductivity of unsaturated porous media. *Water Resources Research*, **12**(3), 513–522. https://doi.org/10.1029/WR012i003p00513

Ngo-Duc, T., Laval, K., Ramillien, G., Polcher, J., & Cazenave, A. (2007). Validation of the land water storage simulated by organising carbon and hydrology in dynamic ecosystems (ORCHIDEE) with gravity recovery and climate experiment (GRACE) data. *Water Resources Research*, **43**, W04427. https://doi.org/10.1029/2006WR004941

O'Neill, B. C., Tebaldi, C., van Vuuren, D. P., Eyring, V., Friedlingstein, P., Hurtt, G., et al. (2016). The scenario model Intercomparison project (ScenarioMIP) for CMIP6. *Geoscientific Model Development*, **9**(9), 3461–3482. https://doi.org/10.5194/gmd-9-3461-2016

Oki, T., Nishimura, T., & Dirmeyer, P. (1999). Assessment of annual runoff from land surface models using total runoff integrating pathways (TRIP). *Journal of the Meteorological Society of Japan, Series II*, **77**(1B), 235–255. https://doi.org/10.2151/jmsj1965.77.1B_235

Padrón, R. S., Gudmundsson, L., Decharme, B., Ducharne, A., Lawrence, D. M., Mao, J., et al. (2020). Observed changes in dry-season water availability attributed to human-induced climate change. *Nature Geoscience*, **13**(7), 477–481. https://doi.org/10.1038/s41561-020-0594-1

Perkins, S. E. (2015). A review on the scientific understanding of heatwaves—Their measurement, driving mechanisms, and changes at the global scale. *Atmospheric Research*, **164–165**, 242–267. https://doi.org/10.1016/j.atmosres.2015.05.014

Perkins, S. E., & Alexander, L. V. (2013). On the measurement of heat waves. *Journal of Climate*, **26**(13), 4500–4517. https://doi.org/10.1175/JCLI-D-12-00383.1

Quesada, B., Vautard, R., Yiou, P., Hirschi, M., & Seneviratne, S. I. (2012). Asymmetric European summer heat predictability from wet and dry southern winters and springs. *Nature Climate Change*, **2**(10), 736–741. https://doi.org/10.1038/nclimate1536

Rashid, M., Chien, R.-Y., Ducharne, A., Kim, H., Yeh, P. J., Peugeot, C., et al. (2019). Evaluation of groundwater simulations in Benin from the ALMIP2 project. *Journal of Hydrometeorology*, **20**(2), 339–354. https://doi.org/10.1175/JHM-D-18-0025.1

Schneider, A., Jost, A., Coulon, C., Silvestre, M., Théry, S., & Ducharne, A. (2017). Global-scale river network extraction based on high-resolution topography and constrained by lithology, climate, slope, and observed drainage density. *Geophysical Research Letters*, **44**, 2773–2781. https://doi.org/10.1002/2016GL071844
Schneider, U., Finger, P., Meyer-Christoffer, A., Rustemeier, E., Ziese, M., & Becker, A. (2017). Evaluating the hydrological cycle over land using the newly-corrected precipitation climatology from the global precipitation climatology Centre (GPCC). *Atmosphere*, 8(12), 52. https://doi.org/10.3390/atmos8030052

Seland, Ø., Bentsen, M., Olivié, D., Toniazzo, T., Gjermundsen, A., Graff, L. S., et al. (2020). Overview of the Norwegian Earth System Model (NorESM2) and key climate response of CMIP6 DECK, historical, and scenario simulations. *Geoscientific Model Development*, 13(12), 6165–6200. https://doi.org/10.5194/gmd-13-6165-2020

Sen, P. K. (1968). Estimates of the regression coefficient based on Kendall's tau. *Journal of the American Statistical Association*, 63(324), 1379–1389. https://doi.org/10.1080/01621459.1968.10480934

Seneviratne, S. I., Corti, T., Davin, E. L., Hirschi, M., Jaeger, E. B., Lehner, I., et al. (2010). Investigating soil moisture-climate interactions in a changing climate: A review. *Earth-Science Reviews*, 99(3–4), 125–161. https://doi.org/10.1016/j.earscirev.2010.02.004

Seneviratne, S. I., Wilhelm, M., Stanelle, T., Van Den Hurk, B., Hagemann, S., Berg, A., et al. (2013). Impact of soil moisture-climate feedbacks on CMIP5 projections: First results from the GLACE-CMIP5 experiment. *Geophysical Research Letters*, 40, 5212–5217. https://doi.org/10.1002/grl.50956

Sergeant, F., Therrien, R., Oudin, L., Jost, A., & Anctil, F. (2021). Evolution of Arctic rivers recession flow: Global assessment and data-based attribution analysis. *Journal of Hydrology*, 601, 126577. https://doi.org/10.1016/j.jhydrol.2021.126577

Smerdon, B. D. (2017). A synopsis of climate change effects on groundwater recharge. *Journal of Hydrology*, 555, 125–128. https://doi.org/10.1016/j.jhydrol.2017.09.047

Swenson, S. C., Clark, M., Fan, Y., Lawrence, D. M., & Perket, J. (2019). Representing intrahillslope lateral subsurface flow in the community land model. *Journal of Advances in Modeling Earth Systems*, 11, 4044–4065. https://doi.org/10.1029/2019MS001833

Tafasca, S., Ducharne, A., & Valentin, C. (2020). Weak sensitivity of the terrestrial water budget to global soil texture maps in the ORCHIDEE land surface model. *Hydrology and Earth System Sciences*, 24(7), 3753–3774. https://doi.org/10.5194/hess-24-3753-2020

Tebaldi, C., Debeire, K., Eyring, V., Fischer, E., Fyfe, J., Friedlingstein, P., et al. (2021). Climate model projections from the scenario model Intercomparison project (ScenarioMIP) of CMIP6. *Earth System Dynamics*, 12(1), 253–293. https://doi.org/10.5194/esd-12-253-2021

Thiery, W., Davin, E. L., Lawrence, D. M., Hirsch, A. L., Hauser, M., & Seneviratne, S. I. (2017). Present-day irrigation mitigates heat extremes. *Journal of Geophysical Research: Atmospheres*, 122, 1403–1422. https://doi.org/10.1002/2016JD025740

Tootchi, A. (2019). *Development of a global wetland map and application to describe hillslope hydrology in the ORCHIDEE land surface model* (Unpublished doctoral dissertation). Sorbonne Université.

Tootchi, A., Jost, A., & Ducharne, A. (2019). Multi-source global wetland maps combining surface water imagery and groundwater constraints. *Earth System Science Data*, 11(1), 189–220. https://doi.org/10.5194/essd-11-189-2019

van Genuchten, M. T. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. *Soil Science Society of America Journal*, 44(5), 892–898. https://doi.org/10.2136/sssaj1980.03615995004400050002x

Voldoire, A., Saint-Martin, D., Sénési, S., Decharme, B., Alias, A., Chevallier, M., et al. (2019). Evaluation of CMIP6 DECK experiments with CNRM-CM6-1. *Journal of Advances in Modeling Earth Systems*, 11(7), 2177–2213. https://doi.org/10.1029/2019MS001683

Vörösmarty, C. J., Fekete, B. M., Meybeck, M., & Lammers, R. B. (2000). Global system of rivers: Its role in organizing continental land mass and defining land-to-ocean linkages. *Global Biogeochemical Cycles*, 14(2), 599–621. https://doi.org/10.1029/1999GB900092

Walko, R. L., Band, L. E., Baron, J., Kittel, T. G. F., Lammers, R., Lee, T. J., et al. (2000). Coupled atmosphere-biophysics-hydrology models for Environmental modeling. *Journal of Applied Meteorology*, 39(6), 931–944. https://doi.org/10.1175/1520-0450(2000)039(0931:CABHMF)2.0.CO;2

Wang, F., Ducharne, A., Cheruy, F., Lo, M. H., & Grandpeix, J. Y. (2018). Impact of a shallow groundwater table on the global water cycle in the IPSL land-atmosphere coupled model. *Climate Dynamics*, 50(9–10), 3505–3522. https://doi.org/10.1007/s00382-017-3820-9

Wu, W.-Y., Lo, M.-H., Wada, Y., Famiglietti, J. S., Reager, J. T., Yeh, P. J., et al. (2020). Divergent effects of climate change on future groundwater availability in key mid-latitude aquifers. *Nature Communications*, 11(1), 3710. https://doi.org/10.1038/s41467-020-17581-y

Zampieri, M., D'Andrea, F., Vautard, R., Ciais, P., de Noblet-Ducoudré, N., & Yiou, P. (2009). Hot European summers and the role of soil moisture in the propagation of Mediterranean drought. *Journal of Climate*, 22(18), 4747–4758. https://doi.org/10.1175/2009JCLI2568.1

Zipper, S. C., Keune, J., & Kollet, S. J. (2019). Land use change impacts on European heat and drought: Remote land-atmosphere feedbacks mitigated locally by shallow groundwater. *Environmental Research Letters*, 14(4), 44012. https://doi.org/10.1088/1748-9326/ab0db3

Zobler, L. (1986). *A world soil hydrology file for global climate modeling* (Tech. Memo. 87802). National Aeronautics and Space Administration.