Global patterns and dynamics of climate-groundwater
interactions

M. O. Cuthbert, T. Gleeson, N. Moosdorf, K. M. Befus, A. Schneider, J.

Hartmann, B. Lehner

► To cite this version:

M. O. Cuthbert, T. Gleeson, N. Moosdorf, K. M. Befus, A. Schneider, et al.. Global patterns and dynamics of climate-groundwater interactions. Nature Climate Change, 2019, 9, pp.137-141. 10.1038/s41558-018-0386-4 . insu-03746422

HAL Id: insu-03746422

https://insu.hal.science/insu-03746422v1

Submitted on 29 Sep 2025

**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>
ORCA - Online Research @
Cardiff

This is an Open Access document downloaded from ORCA, Cardiff University's institutional
repository:https://orca.cardiff.ac.uk/id/eprint/117553/

This is the author's version of a work that was submitted to / accepted for publication.

Citation for final published version:

Cuthbert, M. O., Gleeson, T., Moosdorf, N., Befus, K. M., Schneider, A., Hartmann, J. and Lehner, B. 2019. Global patterns and dynamics of climate-groundwater interactions. Nature Climate Change 9, pp. 137-141. 10.1038/s41558-018-0386-4

Publishers page: https://doi.org/10.1038/s41558-018-0386-4

Please note:

Changes made as a result of publishing processes such as copy-editing, formatting and page numbers may
not be reflected in this version. For the definitive version of this publication, please refer to the published
source. You are advised to consult the publisher's version if you wish to cite this paper.

This version is being made available in accordance with publisher policies. See
http://orca.cf.ac.uk/policies.html for usage policies. Copyright and moral rights for publications made
available in ORCA are retained by the copyright holders.

<figure><img src="image_2.png" /></figure>
1 Global patterns and dynamics of climate-groundwater interactions

2

3 **Authors:** M.O. Cuthbert<sup>1,2,3*</sup>, T. Gleeson<sup>4</sup>, N. Moosdorf<sup>5</sup>, K.M. Befus<sup>6</sup>, A. Schneider<sup>7</sup>,

4 J. Hartmann<sup>8</sup> and B. Lehner<sup>9</sup>

5 <sup>1</sup>School of Earth and Ocean Sciences, Cardiff University, Main Building, Cardiff, CF10 3AT,

6 UK

7 <sup>2</sup>Water Research Institute, Cardiff University, The Sir Martin Evans Building, Museum Ave,

8 Cardiff CF10 3AX, UK

9 <sup>3</sup>Connected Waters Initiative Research Centre, UNSW Sydney, NSW 2052, Australia

10 <sup>4</sup>Department of Civil Engineering and School of Earth and Ocean Sciences, University of

11 Victoria, PO Box 1700, ECS 316, 3800 Finnerty Road, Victoria, British Columbia, Canada

12 V8W 2Y2

13 <sup>5</sup>Leibniz Centre for Tropical Marine Research (ZMT), Fahrenheitstrasse 6, 28359 Bremen,

14 Germany

15 <sup>6</sup>Department of Civil and Architectural Engineering, University of Wyoming, Laramie,

16 Wyoming, USA 82071

17 <sup>7</sup>Sorbonne Université, CNRS, EPHE, Milieux environnementaux, transferts et interactions

18 dans les hydrosystèmes et les sols, METIS, F-75005 Paris, France

19 <sup>8</sup>Institute for Geology, Center for Earth System Research and Sustainability (CEN),

20 Universität Hamburg, Hamburg, Germany

21 <sup>9</sup>Department of Geography, McGill University, 805 Sherbrooke Street West, Montreal H3A

22 0B9, Quebec, Canada

23

24 *Correspondence and request for data sets can be sent to Mark Cuthbert at

25 cuthbertm2@cardiff.ac.uk

26

27
28 **Summary Paragraph**

29 Groundwater is the largest available store of global freshwater¹, upon which more than two
30 billion people rely². It is therefore important to quantify the spatiotemporal interactions
31 between groundwater and climate. However, current understanding of the global scale
32 sensitivity of groundwater systems to climate change³,⁴ – as well as the resulting variation in
33 feedbacks from groundwater to the climate system⁵,⁶ – is limited. Here, using groundwater
34 model results in combination with hydrologic datasets, we examine the dynamic timescales
35 of groundwater system responses to climate change. We show that nearly half of global
36 groundwater fluxes could equilibrate with recharge variations due to climate change on
37 human (~100 year) timescales, and that areas where water tables are most sensitive to
38 changes in recharge are also those that have the longest groundwater response times. In
39 particular, groundwater fluxes in arid regions are shown to be less responsive to climate
40 variability than in humid regions. Adaptation strategies must therefore account for the
41 hydraulic memory of groundwater systems which can buffer climate change impacts on water
42 resources in many regions, but may also lead to a long, but initially hidden, legacy of
43 anthropogenic and climatic impacts on river flows and groundwater dependent ecosystems.
44
45

44

45
Groundwater flow systems exist in dynamic balance with the climate, connecting interacting zones of recharge (i.e. the replenishment of water in the subsurface) and discharge (the loss of groundwater from the subsurface), with multiple feedbacks. As climate varies, changes in the quantity and location of natural groundwater recharge lead to changes in groundwater storage, water table elevations and groundwater discharge¹. These changes in time and space play a central role in controlling the exchange of moisture and energy across the Earth's land surface⁵,⁶ and connect processes critical to, for example, hydro-ecology, as well as carbon and nutrient cycling⁷. Climate-groundwater interactions may also have played a key role in the evolution of our own and other species⁸ and continue to be critical in setting the availability of water for abstraction by humans in coupled food-water-energy systems¹. Recent global mapping of water table depths⁹ and the critical zone¹⁰ suggest where interactions of climate and groundwater may be most tightly coupled. However, they do not resolve where groundwater systems are most sensitive to changes in climate and vice versa, or the timescales over which such changes may occur.

Here, we derive and combine global scale analytical groundwater model results and other hydrologic data sets to provide the first global assessment of the sensitivity of groundwater systems to changes in recharge in both space and time (Figure 1), and discuss their utility as an emergent constraint in understanding and modelling groundwater interactions with climate and other Earth systems at the global scale.

We have characterized the mode of groundwater-climate interactions as being either principally bi-directional or uni-directional using an improved formulation of the water table ratio (*WTR*)¹¹,¹² mapped globally at high resolution (Figure 1a, Figure S1-2). The *WTR* is a measure of the relative fullness of the subsurface and thus the extent of the water table's interactions with topography. Values of *WTR*>1 indicate a topographic control on water table
conditions broadly correlating to shallow (<10 metres below ground level, m bgl) water table depths (*WTDs*) globally (see Methods and Figure S3). This is indicative of a prevalently bi-directional mode of groundwater-climate interaction (Figure 1c) where the climate system can both give to the groundwater system in the form of recharge, and receive moisture back via evapotranspiration if *WTDs* are shallow enough.

The land surface in such regions rejects a proportion of the potential recharge, and groundwater can have a limiting control on land-atmosphere energy exchanges⁵; a tight two-way coupling between groundwater and surface water is also common. In contrast, in ‘recharge controlled’ areas where *WTR*<1, water tables are more disconnected from the topography and, while groundwater may still receive recharge from the land-surface, the extent of two-way interaction between climate and groundwater is limited and the mode of interaction is predominantly uni-directional (Fig. 1c).

We find that regions where *WTR*>1 cover around 46% of the Earth's land area (see Methods, Figure 1a,b) and contribute to the large, but until recently underestimated, extent of groundwater-vegetation interactions globally¹⁰,¹³,¹⁴. Consistent with previous regional analyses and the form of the governing equation (see Methods), our results indicate that bi-directional interactions are more likely to occur in areas with high humidity, subdued topography and/or low permeability. In contrast, regions with *WTR*<1 are more common in drier climates or more mountainous topography¹¹.

In order to assess the large scale temporal sensitivity of climate-groundwater interactions we have used an analytical groundwater solution to quantify groundwater response times (*GRTs*) globally and at high resolution. *GRT* is a measure of the time it takes a groundwater system to re-equilibrate to a change in hydraulic boundary conditions¹⁵. For example, the *GRT* estimates the time to reach an equilibrium in baseflow to streams (or other boundaries) after a change in recharge rate, potentially from climate or land use change. Our results indicate that
groundwater often has a very long hydraulic memory with a global median $GRT$ of nearly 6000 yrs, or approximately 1200 yrs when hyper-arid regions, where recharge is <5 mm/y, are excluded (Figure 1d,e). Only 25% of Earth's land surface area has response times of less than 100 yrs (herein called 'human timescale'). However, this is equivalent to nearly 44% of global groundwater recharge flux, calculated by aggregating contemporary recharge over the land area where $GRT<100$ y, expressed as a proportion of the total global recharge. Around 21% by area have uni-directional climate-groundwater interactions and response times on human timescales, mostly associated with high permeability geology suggesting a strong lithological control (Figure 2a).

The remainder (4%) in areas with bi-directional climate-groundwater interactions are mostly located in the humid, lowland, tropical regions with unconsolidated sediments (e.g. Amazon and Congo Basins, Indonesia), low-lying coastal areas (e.g. Florida Everglades, Asian megadeltas) or in high latitude, low topography humid settings (e.g. northeastern Canada, parts of northern Europe).

A powerful advantage of using analytical groundwater equations such as the *WTR* is that they allow us to directly assess the spatial sensitivity of the mode of climate-groundwater interactions. By taking the derivative of *WTR* with respect to recharge (Figure S4) we have a measure of the sensitivity of the relative fullness of the subsurface to changes in recharge (see Methods). Our results indicate that the mode of climate-groundwater interaction is very insensitive to relative changes in recharge (Figure 2b, Figure S5), with only 5% of the Earth's land surface switching mode for a 50% relative change in recharge rate. This represents a large change in natural groundwater recharge in the context of projections for the coming century¹⁶. However, when absolute recharge rate changes are considered, more sensitivity is apparent and a pattern emerges (Figure S6-7) that indicates the strong inverse relationship between the spatial and temporal sensitivity of groundwater systems to changes in recharge
that we observe (Figure 3b). At small, local scales our calculations may have relatively large
uncertainties, stemming from the uncertainties in global data sets used for the analysis
particularly for hydraulic conductivity (see Methods). However, at the larger scales
considered here, Monte Carlo Experiments (MCE) indicate that, once the variance in each
parameter is combined, the global estimates have relatively small standard deviations
(Figures 1-2, S2).

The global pattern of *GRT* (Figure 1d) indicates a propensity for longer hydraulic memory in
more arid areas. Despite the expected scatter due to geomorphological and lithological
heterogeneity, there is a power law relationship between median *GRT* and groundwater
recharge (*R*) such that *GRT* ∝ 1/*R*^y with y ~ 2 (Figure 3a). This discovery is not directly
expected from the form of the governing equations but is rather an emergent property of
groundwater system interactions with the Earth's land surface and climate system. The
principal control on the observed power law is the distribution of perennial streams
(Figure S8) to which the *GRT* is most sensitive, and which itself is strongly controlled by
climate (Figure S9-11). How to characterize, quantitatively, this climatic control on the
perennial stream distributions is a pertinent question for further hydro-geomorphological
research.

We should not therefore expect *GRTs* to be static nor consider them as 'time constants'
despite being mathematically equivalent to other diffusion processes. Rather, *GRTs* will
evolve in time as both climate and geology vary the geometry and hydraulic properties of
groundwater flow systems. This will occur over long but diverse timescales associated with
changing river geometries.

Despite its importance, most global climate, Earth system, land surface and global hydrology
models exclude groundwater or do not allow groundwater to flow between model grid cells¹⁸-
²⁰. While our results suggest that the spatial distribution of the mode of climate-groundwater
interactions may be rather static over century long timescales, we have shown that nearly a
half of the world's groundwater flux is responsive on 100 y timescales. Hence in order to
capture the important mass and energy transfers correctly, which may affect regional
precipitation and temperature dynamics⁵,⁶, lateral flow circulation of groundwater must be
incorporated into the next generation of global models rather than assuming within-grid-cell
hydrological closure of the water budget as is currently often assumed²¹⁻²³. Our *GRT*
calculations provide direct estimates of spin-up times to improve groundwater-enabled global
models, without having to use the currently employed methods of extrapolation²². Given the
long *GRTs* present over much of the Earth's land-surface, defining initial conditions with an
equilibrium water table calculated for present-day climate conveniently, but wrongly,
assumes stationarity in groundwater levels and fluxes. Since groundwater is known to be the
part of the hydrological system that takes longest to achieve equilibrium²⁴, new approaches
that incorporate the existence of long term transience should continue to be developed²⁵.

The global distribution of *GRTs* suggests that widespread, long-term transience in
groundwater systems persists in the present day due to climate variability since at least the
late Pleistocene in many semi-arid to arid regions (Figure 3a). This is consistent with
observations of larger than expected groundwater gradients, given the current low recharge,
that have been observed in present day arid zones²⁵. While groundwater residence time and
groundwater response time are fundamentally different concepts, we also note the
correspondence between high *GRT* and significant volumes of fossil-aged groundwater
storage in arid regions²,²⁶. The outcome of this result is that groundwater discharge to oases,
rivers or wetlands in otherwise dry landscapes will be particularly intransient in comparison
to climate change, in as much as climate controls the variations in groundwater recharge.
However, our results also indicate that groundwater response times tend to be greater in
regions where water tables are most sensitive to changes in recharge (Figure 3b). This
follows from the fact that both the groundwater response time and the derivative of the water table ratio share a strong dependence on the square of the distance between perennial streams (L, compare Equations 10 and 14).

Away from these more arid contexts, the responsiveness of groundwater systems has recently been demonstrated to be as important as climate controls for the development of hydrological drought²⁷. For example, low *GRT* systems tend to enhance the speed of propagation of meteorological drought through to hydrological drought whereas higher *GRT* systems attenuate climate signals to a greater extent but also show greater lags in recovery from drought. Thus, even within relatively small geographic areas, geological variations can lead to very different drought responses even under similar climate variability. By way of a specific example, increasing lags between meteorological and hydrological drought indicators have been observed between the two most significant aquifers in the UK²⁸ in a manner consistent with what would be expected from our estimates of *GRT* (i.e. Cretaceous Chalk limestone - *GRTs* of months to years, Permo-Triassic sandstone – *GRTs* of years to 100s years, Figure 1d).

Our analysis therefore provides a new framework for understanding global water availability changes under climate change. First, the discovery of a power law relating groundwater recharge and *GRT* suggests that important areas of groundwater discharge in naturally water scarce parts of the world are likely to be more resilient to climate fluctuations than humid areas. However, where groundwater response times are higher, water tables also tend to be most sensitive to changes in recharge in the long term. Hence, accounting appropriately for groundwater-climate interactions within analyses of global water scarcity in the context of climate change is thus of great importance when explicitly considering the contribution of groundwater storage changes²⁹. Second, the long memory of groundwater systems in drylands also means that abrupt (in geological terms) changes in recharge or widely
distributed groundwater abstraction will leave longer legacies. There may also be initially 'hidden' impacts on the future of environmental flows required to sustain streams and wetlands in these regions. It is critical therefore that climate change adaptation strategies which shift reliance to groundwater¹ in preference to surface water should also take account of lags in groundwater hydrology³⁰ and include appropriately long timescale planning horizons for water resource decision making. Third, robust assessments of the impact of climate change on hydrological drought require estimates of 'groundwater responsiveness'²⁷.

The timescale of such responses can be directly informed by our results and improve the decision making process with regard to adaptation strategies to changing drought frequencies under climate change.

206

## 207 **Figure Captions**

**Figure 1.** Global distributions of water table ratios (*WTR*) and groundwater response times (*GRT*) with their conceptual interpretation as metrics of climate-groundwater interactions. (a) Global map of log(*WTR*) with hyper-arid regions of recharge (*R*) < 5 mm/y shaded grey¹⁷. (b) Frequency distribution of global values of log(*WTR*). (c) Conceptual model for *WTR* as a metric for either bi-directional or uni-directional groundwater-climate interactions - *WTR* is dependent on *R*, terrain rise (*d*), distance between perennial streams (*L*) and the saturated thickness of the aquifer (*b*). (d) Global map of log(*GRT*). (e) Frequency distribution of global values of *GRT* - median 5727 yrs (standard deviation, σ = 376 yrs), or 1238 yrs (σ = 92 yrs) when hyper-arid regions are excluded. (f) Conceptual model of *GRT* as a metric of the temporal sensitivity of groundwater-climate interactions.

217

**Figure 2.** Global distributions of the temporal and spatial sensitivity of the mode of climate-groundwater interactions. (a) Temporal sensitivity: percentage of uni-directional and bi-directional groundwater systems, by area globally, that will re-equilibrate significantly to changes in recharge on the timescale of <100 y or >100 y. (b) Spatial sensitivity: percentage of the global area that would change mode from bi-directional to uni-directional climate-groundwater interactions, or vice versa, for a relative change of 50% in recharge, given an unlimited amount of time. Mapped values use the baseline parameter set (see Methods). The median percentage
224 coverage of Earth's landmass for each category from the Monte Carlo Experiments is labelled in the key with
225 standard deviations in percentage coverage shown in brackets. Grey areas represent contemporary recharge
226 <5 mm/y (ref¹⁷).

227

228 **Figure 3.** Global quantitative inter-relationships between climate and the temporal (GRT) and spatial (WTR) sensitivity of groundwater-climate interactions. (a) Globally, median GRT values scale approximately with the inverse of recharge (R) squared. Relationships between recharge and aridity index categories are shown on the top axes as derived in Figure S12. Box extents are at 25-75% percentiles, with Tukey whiskers and outliers. Histograms within each box represent median GRT values from each MCE realisation. (b) The sensitivity of climate-groundwater interactions in time (GRT) and space (dWTR/dR) are log-correlated. Each point uses median values for a geographic location from the MCE realisations. Inset plots are frequency distributions of the slope and r² derived from linear regressions carried out for each realisation indicating consistency in the relationship across the uncertainty range.

237

238 **References**

239 1 Taylor, R. G. *et al*. Ground water and climate change. *Nature Climate Change* **3**, 322-329 (2013).

241 2 Jasechko, S. *et al*. Global aquifers dominated by fossil groundwaters but wells vulnerable to modern contamination. *Nature Geoscience* **10**, 425-429 (2017).

243 3 Döll, P. Vulnerability to the impact of climate change on renewable groundwater resources: a global-scale assessment. *Environmental Research Letters* **4**, 035006 (2009).

246 4 Green, T. R. *et al*. Beneath the surface of global change: Impacts of climate change on groundwater. *Journal of Hydrology* **405**, 532-560 (2011).

248 5 Maxwell, R. M. & Kollet, S. J. Interdependence of groundwater dynamics and land-energy feedbacks under climate change. *Nature Geoscience* **1**, 665-669 (2008).
250 6 Maxwell, R. M. & Condon, L. E. Connections between groundwater flow and transpiration partitioning. *Science* **353**, 377-380 (2016).

252 7 Cole, J. J. et al. Plumbing the global carbon cycle: integrating inland waters into the terrestrial carbon budget. *Ecosystems* **10**, 172-185 (2007).

254 8 Cuthbert, M. et al. Modelling the role of groundwater hydro-refugia in East African hominin evolution and dispersal. *Nature Communications* **8** (2017).

256 9 Fan, Y., Li, H. & Miguez-Macho, G. Global patterns of groundwater table depth. *Science* **339**, 940-943 (2013).

258 10 Xu, X. & Liu, W. The global distribution of Earth's critical zone and its controlling factors. *Geophysical Research Letters* **44**, 3201-3208 (2017).

260 11 Gleeson, T., Marklund, L., Smith, L. & Manning, A. H. Classifying the water table at regional to continental scales. *Geophysical Research Letters* **38** (2011).

262 12 Haitjema, H. M. & Mitchell □ Bruker, S. Are water tables a subdued replica of the topography? *Groundwater* **43**, 781-786 (2005).

264 13 Koirala, S. et al. Global distribution of groundwater□vegetation spatial covariation. *Geophysical Research Letters* **44**, 4134-4142 (2017).

266 14 Schenk, H. J. & Jackson, R. B. Mapping the global distribution of deep roots in relation to climate and soil characteristics. *Geoderma* **126**, 129-140 (2005).

268 15 Carr, E. & Simpson, M. Accurate and efficient calculation of response times for groundwater flow. *Journal of Hydrology* **558**, 470-481 (2017).

270 16 Döll, P., Kaspar, F. & Lehner, B. A global hydrological model for deriving water availability indicators: model tuning and validation. *Journal of Hydrology* **270**, 105-134 (2003).

273 17 Döll, P. & Fiedler, K. Global-scale modeling of groundwater recharge. *Hydrology and Earth System Sciences* **12**, 863-885 (2008).
275 18 Döll, P., Douville, H., Güntner, A., Schmied, H. M. & Wada, Y. Modelling
276 freshwater resources at the global scale: Challenges and prospects. *Surveys in
277 Geophysics* **37**, 195-221 (2016).

278 19 Sood, A. & Smakhtin, V. Global hydrological models: a review. *Hydrological Sciences Journal* **60**, 549-565 (2015).

280 20 Wood, E. F. *et al*. Hyperresolution global land surface modeling: Meeting a grand challenge for monitoring Earth's terrestrial water. *Water Resources Research* **47** (2011).

283 21 Koirala, S., Yeh, P. J. F., Hirabayashi, Y., Kanae, S. & Oki, T. Global □ scale land surface hydrologic modeling with the representation of water table dynamics. *Journal of Geophysical Research: Atmospheres* **119**, 75-89 (2014).

286 22 Milly, P. C. *et al*. An enhanced model of land water and energy for global hydrologic and earth-system studies. *Journal of Hydrometeorology* **15**, 1739-1761 (2014).

288 23 Schaller, M. F. & Fan, Y. River basins as groundwater exporters and importers: Implications for water cycle and climate modeling. *Journal of Geophysical Research: Atmospheres* **114** (2009).

291 24 Ajami, H., McCabe, M. F., Evans, J. P. & Stisen, S. Assessing the impact of model spin□up on surface water□groundwater interactions using an integrated hydrologic model. *Water Resources Research* **50**, 2636-2656 (2014).

294 25 Schulz, S. *et al*. Improving large-scale groundwater models by considering fossil gradients. *Advances in Water Resources* **103**, 32-43 (2017).

296 26 Befus, K. M., Jasechko, S., Luijendijk, E., Gleeson, T. & Cardenas, M. B. The rapid yet uneven turnover of Earth's groundwater. *Geophysical Research Letters* **44**, 5511–5520 (2017).
299 27 Van Lanen, H. A. J., et al. Hydrological drought across the world: impact of climate and physical catchment structure. *Hydrology and Earth System Sciences* **17**, 1715-1732 (2013).

302 28 Bloomfield, J. P., & Marchant, B. P. Analysis of groundwater drought building on the standardised precipitation index approach. *Hydrology and Earth System Sciences* **17**, 4769-4787 (2013).

305 29 Damkjaer, S. & Taylor, R. The measurement of water scarcity: Defining a meaningful indicator. *Ambio* **46**, 513-531 (2017).

307 30 Alley, W. M., et al. Flow and storage in groundwater systems. *Science* **296**, 1985-1990 (2002).

310 Acknowledgements:

311 312 Funding is gratefully acknowledged by: MOC for an Independent Research Fellowship from the UK Natural Environment Research Council (NE/P017819/1); JH from the German Science Foundation DFG (Cluster of Excellence “CliSAP,” EXC177, Universität Hamburg) and Bundesministerium für Bildung und Forschung Project PALMOD (Ref 01LP1506C); NM from The German Federal Ministry of Education and Research (BMBF) (Grant #01LN1307A); AS from Agence Nationale de la Recherche (ANR grant ANR-14-CE01-00181-01) and the French national programme LEFE/INSU; BL from Natural Sciences and Engineering Research Council of Canada (NSERC) (Discovery Grant RGPIN/341992).

321 Author Contributions

322 The idea for the paper was conceived by MOC and TG. Analyses were by all authors. The manuscript was written by MOC with input from all authors.

324
325 **Data Availability statement**

326 Once the paper is accepted the main data outputs (i.e. *WTR* and *GRT* digital maps) will be
327 made freely available for download via an online data repository and a link included in the
328 final published version of the paper.

329

330
331 METHODS

332 Derivation of Equations

333 Governing groundwater flow equations

334 The governing equations were formulated by considering an ideal homogeneous, horizontal
335 unconfined aquifer bounded at one end ($x = L/2$) by a stream assumed to be a constant head
336 boundary and at the other ($x = 0$) by a no-flow boundary representing a flow divide
337 (Figure S13). The one-dimensional (Boussinesq) equation of groundwater flow for such an
338 aquifer receiving homogeneous recharge can be given as follows:

$$
339 \quad \frac{\partial}{\partial x} \left( K h \frac{\partial h}{\partial x} \right) = S \frac{\partial h}{\partial t} - R(t) \qquad (1)
$$

340 where K is hydraulic conductivity [LT⁻¹], S is storativity [-], *h(x,t)* is hydraulic head [L], *t* is
time [T], *x* is distance [L] and *R(t)* is groundwater recharge [LT⁻¹].

342 If changes in transmissivity due to fluctuations in groundwater heads are assumed to be
343 negligible, Equation (1) may be linearised as follows:

$$
344 \quad T \frac{\partial^2 h}{\partial x^2} = S \frac{\partial h}{\partial t} - R(t) \qquad (2)
$$

345 where *T* is transmissivity [L²T⁻¹], and *T* = *KH*, with *H* the average saturated thickness [L].

346 The lateral boundary conditions are as follows:

$$
347 \quad \frac{\partial h(0,t)}{\partial x} = 0, h\left(\frac{L}{2}, t\right) = b \qquad (3)
$$

348 The parameter *L* is thus a characteristic length equivalent to the distance between perennial
349 streams which act as fixed head groundwater discharge boundaries.

350 Water table ratio (WTR) derivations
For steady state flow, where $h(x,t)$ becomes $h(x)$, the solution to Equation 1 for the stated boundary conditions is:

$$h(x) = \left( b^2 + \frac{R}{K} \left( \frac{L^2}{4} - x^2 \right) \right)^{0.5} \qquad (4)$$

At the flow divide, $x = 0$, therefore:

$$h(0) = \sqrt{b^2 + \frac{RL^2}{4K}} \qquad (5)$$

For steady state flow, the solution to the linearised form, Equation 2, for the stated boundary conditions is:

$$h(x) = \frac{R}{2T} \left( \frac{L^2}{4} - x^2 \right) + b \qquad (6)$$

At the flow divide, $x = 0$, therefore:

$$h(0) = \frac{RL^2}{8T} + b \qquad (7)$$

The WTR is defined¹² as the ratio of the head at the flow divide above the fixed head boundary (i.e. $h_0 - b$) to the maximum terrain rise above the fixed head boundary, $d$ [L].

This yields a new, non-linearised, form of the WTR, from Equation 5 as follows:

$$WTR_{NL} = \frac{\sqrt{b^2 + \frac{RL^2}{4K} - b}}{d} \qquad (8)$$

For the linearised form, from Equation 7, and as originally given by ref ¹², the WTR is:

$$WTR_L = \frac{RL^2}{8Td} = \frac{RL^2}{8KHd} \qquad (9)$$

Equations 8 and 9 become equivalent for combinations of small $L$ or $R$, or large $K$.

All maps and analysis presented in this paper use the non-linear form of the WTR (Equation 8) with the exception of Figure S1 where the two versions are compared, and
calculated using the *L* parameters derived using a minimum river discharge threshold of 0.1 m³/s. A comparison of global maps and frequency distributions for the linear and non-linear forms are shown in Figure S1-2. The frequency distribution comparison (Figure S2) shows that the new non-linear formulation has a narrower and more symmetric distribution with a median closer to zero than the linearised form. This is indicative of its better physical representation such that the extent of higher *WTRs* is limited by the feedback between higher water table elevation and concomitant increases in transmissivity inherent in the non-linear Boussinesq equation (Equation 1).

The *WTR* is a measure of the relative fullness of the subsurface and thus the extent of the water table's interactions with topography. We have therefore used the *WTR* to characterize the dominant mode of groundwater-climate interactions as being either principally bi-directional or uni-directional based on whether they are 'topographically controlled' (*WTR*>1) or 'recharge controlled' (*WTR*<1), respectively. This is a reasonable approximation since a global comparison with water table depths (*WTDs*) (Figure S3) indicates that *WTR*>1 broadly correlates to shallow (<10 metres below ground level) water table conditions. This condition is indicative of a prevalently bi-directional mode of groundwater-climate interaction where the climate system can both give to the groundwater system in the form of recharge, and receive moisture back where local variations in *WTDs* enable evapotranspiration to occur from groundwater directly. In contrast, areas with *WTR*<1 show increasingly large *WTDs* well beyond plant rooting depths leading to predominantly uni-directional climate-groundwater interactions where the groundwater system receives recharge from the climate system but there is more limited potential for feedback in the other direction.

The sensitivity of the *WTR* to changing recharge is given by differentiating Equation 8 with respect to *R*:
$$394 \qquad \frac{dWTR_{NL}}{dR} = \frac{L^2}{8Kd} \left(b^2 + \frac{RL^2}{4K}\right)^{-0.5} \qquad (10)$$

This equation represents the sensitivity of the maximum head to recharge relative to the topography which can be understood as the sensitivity of the ‘fullness’ of the subsurface to changes in recharge.

Following from Equations 8, we calculate the recharge required for the WTR to equal 1 for every grid cell as:

$$400 \qquad R_{WTR=1} = \frac{4K}{L^2} (d^2 + 2db) \qquad (11)$$

The difference between $R$ and the values given in Equation 11 then gives an expression for the change in recharge ($\Delta R$) needed to effect a change in the WTR across the transition between topography control (bi-directional climate-groundwater interactions) and recharge control (unidirectional climate-groundwater interactions) modes. In absolute terms this is:

$$405 \qquad \Delta R_{abs} = R - R_{WTR=1} \qquad (12)$$

and in relative terms it becomes:

$$407 \qquad \Delta R_{rel} = \frac{R - R_{WTR=1}}{R} \qquad (13)$$

### Groundwater response time (GRT) definition

The groundwater response time is, in general terms, a measure of the time it takes a groundwater system to respond significantly (as defined below) to a change in boundary conditions¹⁵,³¹⁻³⁵ and is defined here as follows:

$$412 \qquad GRT = \frac{L^2 S}{\beta T} \qquad (14)$$

where $\beta$ is a dimensionless constant, $T$ is transmissivity [$L^2 T^{-1}$], $S$ is storativity [-] and $L$ is the distance between perennial streams [L]. To illustrate why this equation defines a time of
response consider a groundwater mound such as that shown in Figure S13. Let the initial shape of the mound (of maximum height A), due to some steady recharge, be given by:

$$h(x, 0) = A \cdot \cos\left(\frac{\pi x}{L}\right) \qquad (15)$$

If recharge suddenly ceases (i.e., a step change) then it can be shown, in the manner of ref³³, that the solution to the linearised Equation 2 without recharge (i.e. R(t)=0) is:

$$h(x, t) = h(x, 0) \cdot \exp\left(-\frac{t}{GRT}\right) \qquad (16)$$

for β is equal to π².

Hence, for this case, the GRT controls the timescale for the groundwater levels to decay exponentially to reach 63% re-equilibrium after a change in boundary (recharge) conditions (i.e., an “e-folding” timescale). This value for β was chosen in order to be consistent with mathematically equivalent uses of ‘time constants’ (often denoted as τ), in other branches of science.

As outlined by ref³⁴, comparing the timescale of a particular forcing to the GRT can be a useful measure of the degree of transience a groundwater system will manifest in terms of variations in lateral groundwater flow. However, there is an important difference to note in the case of a step change in conditions, as used to define GRT in Equation 14, in comparison with a periodic variation in the forcing recharge (of period P). For the step change case outlined above, both heads and fluxes decay exponentially after the change in recharge. However, in the periodic case, where GRT >> P, variations in recharge lead to very stable groundwater fluxes (including at the downstream lateral boundary) but large temporal changes in groundwater head across much of the aquifer³⁵. Thus, it is important to distinguish between the control of GRT on the degree of transience in either heads or fluxes, depending on the nature of the boundary conditions.
438 Spatial input data and manipulation

439 *Global mapping of the distance between perennial streams (L)*

440 The distance between perennial streams (L) was calculated using a globally consistent river
441 network provided by the HydroSHEDS database³⁶ which was derived from the 90 m digital
442 elevation model of the Shuttle Radar Topography Mission (SRTM). For this study, we
443 extracted the global river network from the HydroSHEDS drainage direction grid at 500 m
444 pixel resolution by defining streams as all pixels that exceed a long-term average natural
445 discharge threshold of 0.1 cubic meters per second, resulting in a total global river length of
446 29.4 million kilometers. Smaller rivers with flows below this threshold were excluded as they
447 are impaired by increasing uncertainties in the underpinning data. However, the sensitivities
448 of the most important results of this paper to the chosen threshold are considered in our
449 uncertainty analysis below. Estimates of long-term (1971-2000) discharge averages have
450 been derived through a geospatial downscaling procedure³⁷ from the 0.5° resolution runoff
451 and discharge layers of the global WaterGAP model (version 2.2, 2014) a well-documented
452 and validated integrated water balance model¹⁶,³⁸. Only perennial rivers were included in the
453 assessment; intermittent and ephemeral rivers were identified through statistical discharge
454 analysis (lowest month of long-term climatology is 0) and extensive manual corrections
455 against paper maps, atlases and auxiliary data, including the digital map repository of
456 National Geographic³⁹. L was calculated for every pixel of the landscape (Figure S8) by
457 identifying the shortest combined Euclidean (straight-line) distance between two river
458 locations at opposing sides of the pixel. Neighbourhood low pass filters (5x5 kernel size)
459 were applied to remove outlier pixels and speckling. All calculations were performed in
460 ESRI© ArcGIS environment using custom-made scripts.

461 *Global mapping of the water table ratio (WTR), groundwater response times (GRT) and*
462 *other expressions*
Global WTR maps were created from the above equations using: the recharge rate (R in m/y), based on ref¹⁷, a minimum saturated thickness of the aquifer (b) set to 100 m (refs⁴⁰,⁴¹), the distance between two perennial streams (L, in m, as described above), intrinsic permeability values (m²) reported in ref⁴⁰ were converted to hydraulic conductivity (m/s) by assuming standard temperature and pressure (1 x 10⁷ multiplication factor) and then converted to units of m/y. The maximum terrain rise between rivers (d, in m) was based on the range of elevations in the 250m GMTED2010 data set⁴².

The GRT was mapped using the same L data and hydraulic conductivity values as for the WTR calculations. Transmissivity (T, m²/y) was calculated by multiplying the hydraulic conductivity with a fixed saturated thickness of 100 m (refs⁴⁰,⁴¹). It was assumed that storativity (S) for unconfined aquifers is dominated by the specific yield and that this can be approximated by mapped porosity values⁴⁵. Owing to the significant uncertainties in these assumptions for calculating T and S values the parameters were subjected to a Monte Carlo analysis as described below.

Each of the data sets was prepared to match a global equal-area projection with a grid size of 1 km x 1 km, and the calculations of the data sets were performed in ArcGIS. To avoid mathematical problems, for zero values of d and R, 1 and 0.00001 were added, respectively. For WTR estimates, regions where contemporary groundwater recharge was estimated as < 5 mm/y (ref¹⁷) were excluded from the analysis due to the increasingly large relative uncertainties in recharge below this range, and the resulting unrealistic sensitivity of the resulting WTR estimates. For deriving the frequency distributions and comparisons of parameters from the range of derived geo-spatial data sets, point values were taken from each raster of interest for 10,000 randomly distributed locations across the Earth's land-surface. Global distributions of the parameters d, K and S are given in Figure S10 and relationships
487 between *R* and *L*; *d* and *WTR*; and *R* and *WTR* are explored in Figure S9, and Figures S11, 488 respectively. All areal calculations ignore the Antarctic landmass.

489 Although we have made best use of coherent available global datasets at high (1 km) 490 resolution for the calculations, our results are intended for appropriate large scale 491 interpretation, not detailed local analysis.

492 **Justification of the model assumptions**

493 Our calculations are based on mapped surface lithology only and, as such, they represent a 494 first estimate of the response of unconfined groundwater across the global land surface. The 495 more complex responses of regional or local confined aquifers, which may be locally 496 important to discerning groundwater-climate interactions, are not considered. However, such 497 confined aquifers only cover around 6-20% of the Earth's surface<sup>43</sup>, are often located in more 498 arid parts of the world and are, by definition, inherently less connected to the land surface and 499 climate-related processes.

500 Using 1-D analytical solutions to the groundwater flow equations gives a powerful advantage 501 over the use of more complex models in enabling the sensitivity of the key parameters 502 controlling patterns and timescales of climate-groundwater interactions to be analysed 503 analytically. This, for example, allows us to sample the entire parameter space directly rather 504 than a restricted subset via a limited ensemble of more computationally expensive numerical 505 model runs. Equation 1 assumes the validity of the Dupuit-Forchheimer approximation 506 whereby the water table is assumed to be a true free surface governed by effective hydraulic 507 parameters and that water pressure in the direction normal to the flow is approximately 508 hydrostatic. This is a good approximation when the ratio of the lateral extent of the average 509 saturated depth is more than approximately 5 times its depth<sup>12</sup>, i.e. H/(L/2) < 0.2 (see 510 Figure S13). Calculating the maximum saturated depth *h*<sub>max</sub> as the smaller of *d*+*b* or *h*<sub>0</sub>, and
approximating the average saturated depth as $(h_{max} + b)/2$, we find that the criterion $H/(L/2) < 0.2$ is met in 96% of our global grid calculations. Locations which fail this test are all in mountainous regions where Equation 1 cannot account accurately for steep hillslope groundwater hydraulics and hence our results may be less reliable in such areas.

The GRT is a parameter which consistently appears in solutions to the groundwater flow equations and has been used for decades³² as a robust estimate for the timescale of re-equilibration of a groundwater system following a change in boundary conditions⁸,¹⁵,³⁰⁻³⁵,⁴⁴⁻⁴⁸.

Thus it is an appropriate metric for long term transience which is currently impossible to model in state of the art coupled groundwater-surface water models, which are limited to short run times even for regional scale analyses due to their massive computational demands.

More realistic aquifer geometries and initial water table configurations lead to behaviours which are more complex than the case of a simple exponential decay⁴⁶, and non-uniform flow fields (strong convergence or divergence) can also lead to variations in GRT (refs⁴⁴,⁴⁷,⁴⁸). We have therefore included these factors in an uncertainty analysis as outlined below.

While the models used here cannot represent the detailed process interactions in the way that a distributed fully coupled 3-D model would, they have a strong theoretical basis and show consistency with other large scale studies based on very different model assumptions and data sets. Justification for the approach of using *WTR* as a proxy for the mode of climate-groundwater interaction is given in at least four ways. First, at global scale, similarities of *WTR* with shallow *WTD* globally⁹ are strong (Figure S3), given the very different model assumptions and data sets employed in the two studies. Second, at a continental scale for the contiguous US a recent study compared the results of a physically based, 3-D, fully coupled surface water-groundwater model validated against water table depth data, against the *WTR* metric⁴¹. The results show scatter as expected due to variations in the derivation of the comparative characteristic length scales used in the comparison. However, general trends and
geographic patterns at a regional scale compare well between the WTD computed by the fully
coupled model and the calculated WTRs. Third, also at a continental scale for the contiguous
US, a systematic relationship has been shown between *WTR* and mean stream junction angles
which are indicative of a strong coupling between surface and subsurface⁴⁹. Lastly,
comparisons of *WTR* calculations against a more complex 3-D regional groundwater flow
model, has indicated that the *WTR* is a robust indicator of groundwater's connection to the
land surface as it is a strong predictor of the propensity for local versus regional flow
conditions⁵⁰. Our analyses thus allow us to make a robust first global scale estimate of the
sensitivity of climate-groundwater interactions, while enabling the range of uncertainty to be
fully and directly appreciated.

546 Uncertainties and Monte Carlo experiments

547 We ran 10,000 Monte Carlo experiments (MCE) at 10,000 randomly distributed locations
548 across the Earth's land-surface to investigate the range of uncertainty due to parameter
549 uncertainties as well as model structural simplifications.

550 Hydraulic conductivity (K) was allowed to vary log-normally within uncertainty ranges
551 defined in refs⁴⁰,⁵⁰, this parameter having by far the highest uncertainty of any others used in
552 our calculations. Groundwater recharge (R) values were taken from ref¹⁷ but allowed to vary
553 through a normal distribution with a standard deviation of 22% of this baseline, chosen
554 according to the difference with a contrasting global recharge distribution⁵²,⁵³ commonly used
555 in other global hydrological calculations. Storativity (S) was sampled from a normal
556 distribution with standard deviations of 25% of the mapped value⁵³. Although the absolute
557 error in the DEM used is only 1-2 m, we allowed the maximum terrain rise (d) to vary
558 normally with a standard deviation of 10% to allow for uncertainties due to gridding. The
559 minimum saturated thickness of the aquifer (b) was allowed to vary log-normally around
560 100 m with a standard deviation of 0.3 orders of magnitude. Sampled distributions were cut off at zero to stop meaningless negatives being included in the calculations.

562 Parameter uncertainty in the distance between perennial streams (L), calculated from the variation in L for an order of magnitude change in discharge threshold used to define the stream network (from 0.1 to 1 m³/s), gives a median uncertainty of a factor of 1.9. However, there is also additional uncertainty to L due to the choice of the one-dimensional groundwater flow solutions applied, which ignore non-uniform (i.e. convergent or divergent) flow fields which are common in real catchments. In order to account for the maximum likely range of possible uncertainty, we have compared the 1-D analytical solutions used here to cases of radial flow which represent an extreme 2-D non-uniform flow end-member for natural groundwater flow systems. By equating the distance between perennial streams (L) to be equal to the radius of the flow domain for the equivalent radial solutions, we can estimate the impact of this choice on both *WTR* and *GRT*. For *WTR*, by replacing Equation 6 with Eq 30.11 from ref⁵⁴, the average error is approximately a factor of 2. For the *GRT*, comparison of recession timescales for 1-D and radial flow cases (e.g. Appendix A of ref ⁴⁶) indicates a similar level of uncertainty due to non-uniform flow as for the *WTR*. We therefore added a log-normal variation in L with a standard deviation of 0.3 orders of magnitude to accommodate the likely range of combined parameter and structural uncertainty.

578

579 References (additional for Methods, numbering following on from main text references)

580 31 Domenico, P. A. & Schwartz, F. W. *Physical and chemical hydrogeology*. Vol. 506 (Wiley New York, 1998).

582 32 Downing, R., Oakes, D., Wilkinson, W. & Wright, C. Regional development of groundwater resources in combination with surface water. *Journal of Hydrology* 22, 155-177 (1974).
585 33 Erskine, A. & Papaioannou, A. The use of aquifer response rate in the assessment of groundwater resources. *Journal of Hydrology* **202**, 373-391 (1997).

587 34 Currell, M., Gleeson, T. & Dahlhaus, P. A new assessment framework for transience in hydrogeological systems. *Groundwater* **54**, 4-14 (2014).

589 35 Townley, L. R. The response of aquifers to periodic forcing. *Advances in Water Resources* **18**, 125-146 (1995).

591 36 Lehner, B., Verdin, K. & Jarvis, A. New global hydrography derived from spaceborne elevation data. *Eos, Transactions American Geophysical Union* **89**, 93-94 (2008).

593 37 Lehner, B. & Grill, G. Global river hydrography and network routing: baseline data and new approaches to study the world's large river systems. *Hydrological Processes* **27**, 2171-2186 (2013).

596 38 Alcamo, J. et al. Development and testing of the WaterGAP 2 global model of water use and availability. *Hydrological Sciences Journal* **48**, 317-337 (2003).

598 39 National Geographic, *Atlas of the World*. Ninth Edition edn, National Geographic (2010).

600 40 Gleeson, T., Moosdorf, N., Hartmann, J. & Beek, L. A glimpse beneath earth's surface: GLobal HYdrogeology MaPS (GLHYMPS) of permeability and porosity. *Geophysical Research Letters* **41**, 3891-3898 (2014).

603 41 Condon, L. E. & Maxwell, R. M. Evaluating the relationship between topography and groundwater using outputs from a continental-scale intergrated hydrology model. *Water Resources Research* **51**, 6602-6621 (2015).

606 42 Danielson, J. J. & Gesch, D. B. Global multi-resolution terrain elevation data 2010 (GMTED2010). Report No. 2331-1258, (US Geological Survey, 2011).
<table>
  <tr>
    <td>608</td>
    <td>43</td>
    <td>de Graaf, Inge E.M., et al. A global-scale two-layer transient groundwater model: Development and application to groundwater depletion. <i>Advances in water resources</i> 102, 53-67 (2017).</td>
  </tr>
  <tr>
    <td>609</td>
    <td>44</td>
    <td>Rousseau Gueutin, P. et al. Time to reach near steady state in large aquifers. <i>Water Resources Research</i> 49, 6893-6908 (2013).</td>
  </tr>
  <tr>
    <td>610</td>
    <td>45</td>
    <td>Rushton, K. R. &amp; Redshaw, S. C. <i>Seepage and groundwater flow: Numerical analysis by analog and digital methods</i>. (Chichester, 1979).</td>
  </tr>
  <tr>
    <td>611</td>
    <td>46</td>
    <td>Cuthbert, M. Straight thinking about groundwater recession. <i>Water Resources Research</i> 50, 2407-2424 (2014).</td>
  </tr>
  <tr>
    <td>612</td>
    <td>47</td>
    <td>Cuthbert, M. et al. Understanding and quantifying focused, indirect groundwater recharge from ephemeral streams using water table fluctuations. <i>Water Resources Research</i> 52, 827-840 (2016).</td>
  </tr>
  <tr>
    <td>613</td>
    <td>48</td>
    <td>Walker, G. R., Gilfedder, M., Dawes, W. R. &amp; Rassam, D. W. Predicting Aquifer Response Time for Application in Catchment Modeling. <i>Groundwater</i> 53, 475-484 (2015).</td>
  </tr>
  <tr>
    <td>614</td>
    <td>49</td>
    <td>Seybold, H., Rothman, D. H., &amp; Kirchner, J. W. Climate's watermark in the geometry of stream networks. <i>Geophysical Research Letters</i> 44, 2272-2280 (2017).</td>
  </tr>
  <tr>
    <td>615</td>
    <td>50</td>
    <td>Gleeson, T., &amp; Manning, A. H. Regional groundwater flow in mountainous terrain: Three-dimensional simulations of topographic and hydrogeologic controls. <i>Water Resources Research</i> 44, W10403 (2008).</td>
  </tr>
  <tr>
    <td>616</td>
    <td>51</td>
    <td>Gleeson, T., et al. Mapping permeability over the surface of the Earth. <i>Geophysical Research Letters</i> 38, L02401 (2011).</td>
  </tr>
  <tr>
    <td>617</td>
    <td>52</td>
    <td>De Graaf, I. E. M., et al. A high-resolution global-scale groundwater model. <i>Hydrology and Earth System Sciences</i> 19, 823-837 (2015).</td>
  </tr>
</table>
632 53 Wada, Y., et al. Global depletion of groundwater resources. *Geophysical Research*

633 *Letters* **37**, L20402 (2010).

634 54 Bruggeman, G.A., Analytical Solutions of Geohydrological Problems. *Developments*

635 *in Water Science* **46**, 3-959 (Elsevier, 1999)
a

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

log(Water Table Ratio)

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

BI-DIRECTIONAL
groundwater climate
interactions

c

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

b

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

UNI-DIRECTIONAL
groundwater climate
interactions

d

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

Groundwater Response Time (y)

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

e

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

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

<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /></figure>

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