Global groundwater warming due to climate change

Aquifers contain the largest store of unfrozen freshwater, making groundwater critical for life on Earth. Surprisingly little is known about how groundwater responds to surface warming across spatial and temporal scales. Focusing on diffusive heat transport, we simulate current and projected groundwater temperatures at the global scale. We show that groundwater at the depth of the water table (excluding permafrost regions) is conservatively projected to warm on average by 2.1 °C between 2000 and 2100 under a medium emissions pathway. However, regional shallow groundwater warming patterns vary substantially due to spatial variability in climate change and water table depth. The lowest rates are projected in mountain regions such as the Andes or the Rocky Mountains. We illustrate that increasing groundwater temperatures influences stream thermal regimes, groundwater-dependent ecosystems, aquatic biogeochemical processes, groundwater quality and the geothermal potential. Results indicate that by 2100 following a medium emissions pathway, between 77 million and 188 million people are projected to live in areas where groundwater exceeds the highest threshold for drinking water temperatures set by any country. Earth’s climatic system warms holistically in response to the radiative imbalance from increased concentrations of greenhouse gases 1 . While the ocean absorbs most of this additional heat 2 , the terrestrial subsurface and groundwater also function as a heat sink. With a stable climate, seasonal temperature variation penetrates to a depth of 10–20 m, below which temperatures generally increase with depth in accordance with the geothermal gradient 3 . However, present-day borehole temperature–depth profiles frequently show an inversion (that is, temperature decreasing with depth) for up to 100 m due to recent, decadal surface warming 4 . Deviations from steady-state sub-surface temperatures in deep boreholes (for example, >300 m) have been used to

Aquifers contain the largest store of unfrozen freshwater, making groundwater critical for life on Earth.Surprisingly little is known about how groundwater responds to surface warming across spatial and temporal scales.Focusing on diffusive heat transport, we simulate current and projected groundwater temperatures at the global scale.We show that groundwater at the depth of the water table (excluding permafrost regions) is conservatively projected to warm on average by 2.1 °C between 2000 and 2100 under a medium emissions pathway.However, regional shallow groundwater warming patterns vary substantially due to spatial variability in climate change and water table depth.The lowest rates are projected in mountain regions such as the Andes or the Rocky Mountains.We illustrate that increasing groundwater temperatures influences stream thermal regimes, groundwater-dependent ecosystems, aquatic biogeochemical processes, groundwater quality and the geothermal potential.Results indicate that by 2100 following a medium emissions pathway, between 77 million and 188 million people are projected to live in areas where groundwater exceeds the highest threshold for drinking water temperatures set by any country.
Earth's climatic system warms holistically in response to the radiative imbalance from increased concentrations of greenhouse gases 1 .While the ocean absorbs most of this additional heat 2 , the terrestrial subsurface and groundwater also function as a heat sink.With a stable climate, seasonal temperature variation penetrates to a depth of 10-20 m, below which temperatures generally increase with depth in accordance with the geothermal gradient 3 .However, present-day borehole temperature-depth profiles frequently show an inversion (that is, temperature decreasing with depth) for up to 100 m due to recent, decadal surface warming 4 .Deviations from steady-state subsurface temperatures in deep boreholes (for example, >300 m) have been used to evaluate terrestrial heat storage and to estimate past, pre-observational surface temperature changes at a global scale 5 .Previous multi-continental synthesis studies on subsurface warming provide critical information on climate dynamics, but impacts on groundwater resources and associated implications are commonly ignored.
With the advent of the Gravity Recovery and Climate Experiment (GRACE) satellites, global datasets and global hydrological models, there is an emerging body of global-scale groundwater research [6][7][8][9] .However, global-scale groundwater studies so far have focused on resource quantity (for example, levels, recharge rates and gravity signals), whereas global-scale research into groundwater quality, including temperature, is rare.Furthermore, prominent syntheses of the relationship between anthropogenic climate change and https://doi.org/10.1038/s41561-024-01453-xset by drinking water standards (Fig. 1d) and (3) where discharge of warmed groundwater will have the most pronounced impact on river temperatures and aquatic ecosystems (Fig. 1e).Our model is global, and its resolution limits detailed capture of small-scale processes, producing conservative results based on tested hydraulic and thermal assumptions, including realistic advection from basin-scale recharge.More localized processes may lead to higher groundwater temperatures in areas with increased downward flow (for example, river-based recharge) or elevated surface temperatures (for example, urban heat islands) (Supplementary Note 1 provides details).

Groundwater temperatures
We use gridded data to calculate transient subsurface temperaturedepth profiles across the globe (Methods).Besides past and current temperatures, we present potential (modest mitigation) and worst-case (no mitigation) projections to 2100 based on the Shared Socioeconomic Pathway (SSP) 2-4.5 or SSP 5-8.5 climate scenarios of phase 6 the Coupled Model Intercomparison Project (CMIP6) (ref.37).Results can be accessed and visually explored using an interactive Google Earth Engine app available at https://susanneabenz.users.earthengine.app/view/subsurface-temperature-profiles. Figure 2a-c displays maps of mean GWT at the depth of the water table and at 5 and 30 m below ground surface for 2020.
Comparison with measured data demonstrates a good accuracy of the model given the global scale with a root mean square error of 1.4 °C and a coefficient of determination of 0.75 (Fig. 2d).An in-depth discussion on model reliability, uncertainty and limitations is given in Supplementary Note 2.
The median GWT at the water table in 2020 was 21.0 °C (5.6 °C, 29.3 °C; 10th, 90th percentile; Fig. 2a).In comparison, using the same ECMWF re-analysis (ERA-5) data product, air temperatures in 2020 were lower at 17.6 °C (1.4 °C, 27.0 °C).This thermal offset is attributable to various processes and conditions including snow pack insulation in colder climates 38 and increased temperatures with depth following the geothermal gradient.
Simulated temperature-depth profiles are displayed at six example locations in Fig. 2e, including their seasonal envelope.Supplementary Note 3 provides a discussion of seasonality.Whereas all locations show an inversion of the temperature-depth profile, the depth at which this thermal gradient 'inflection point' (ref.4) is reached varies greatly based on the rate and duration of recent climate change.At the example location in Mexico, temperatures begin to increase with depth (as expected based on the local geothermal gradient) from approximately 10 m downwards, whereas at the example location in Brazil, the inflection point reaches a depth of 45 m (Fig. 2c).Globally, it has reached 15 (<1, 40) m (Extended Data Fig. 1a).Heat advection from vertical groundwater flow may also influence the depth of the inflection point 4 , but only heat diffusion is considered in our model as this is the dominant heat-transport mechanism at the modelled spatial scale (Methods).
To better assess the impact of recent climate change on groundwater temperatures at the water table depth, we compare annual mean GWTs from 2000 and 2020.Over this 20-year period, GWTs increased on average by 0.3 (0.0, 0.8) °C (Fig. 3a).We do not find any distinct large-scale patterns.However, some of the highest temperature increases occur in parts of Russia (for example, > + 1. 5 ∘ C north of Novosibirsk), while parts of Canada experienced cooling (for example, < −0. 5 °C in Saskatoon) between the two years.Both regions have shallow water tables, with GWTs tightly coupled to seasonal surface temperature variations and short-term intra-annual changes, rather than just the long-term surface temperature signals.As such, one hot summer can drastically alter the modelled GWT difference between 2000 and 2020.The influence of weather conditions for a given year is also notable in the depth profiles for six selected locations (Fig. 3d).Noticeable variations occur in the upper 5 m of mean temperature groundwater (for example, refs.10,11) concentrate on quantity leaving quality aspects unexplored 12 .Water temperature, sometimes known as the 'master environmental variable' (ref.13), is an understudied groundwater quality parameter in the context of climate change.
Whereas global studies of river and lake warming have been conducted 14,15 , there are no global assessments of climate change impacts on groundwater temperatures (GWTs).This is despite the high importance of groundwater, which represents the largest global reservoir of unfrozen freshwater 16 , providing at least part of the water supply for half the world 17 and close to half of the global irrigation demand 18 .It also sustains terrestrial and aquatic ecosystems 19 , particularly in the face of climate change 10 .Given the role of temperature as an overarching water quality variable and observational evidence of groundwater warming in different countries in response to recent climate change 4,20,21 , the potential impact of climate warming on groundwater temperatures at a global scale remains a critical knowledge gap.
Groundwater temperature influences a suite of biogeochemical processes that alter groundwater quality 22 .For example, an increase in temperatures reduces gas solubility and raises metabolism of organisms, with an increased rate of oxygen consumption and a shift in redox conditions 23 .Because many aquifers already possess low oxygen concentrations, a small change in temperature could trigger a shift from an oxic to a hypoxic or even an anoxic regime 24,25 .This switch can in turn facilitate the mobilization of redox-sensitive constituents such as arsenic, manganese and phosphorus 26,27 .Increases in soluble phosphorus in groundwater discharging to surface water can trigger harmful algal blooms 28 , and elevated arsenic and manganese contents in potable water supplies pose direct risks to human health 29 .Groundwater warming will also cause a shift in groundwater community composition with a challenge to biodiversity and the risk of an impaired cycling of carbon and nutrients 24,25 .Shallow soil and groundwater warming may also cause temperatures in water distribution networks to cross critical thresholds, with potential health implications such as the growth of pathogens such as Legionella spp. 30.
Diffusive discharge of thermally stable groundwater to surface water bodies modulates their temporal thermal regimes 30 .Also, focused groundwater inflows can create cold-water plumes that provide thermal refuge for stressed aquatic species 31 , including many prize cold-water fish.Accordingly, groundwater warming will increase ambient water temperatures in surface water bodies and the temperatures of groundwater-sourced thermal refuges.Spring ecosystems will also be affected.For example, crenobionts (true spring water species) have a very narrow temperature optimum and tolerance; hence, warming groundwater near the mouths of springs will lead to changes in their reproduction cycles, food web interactions and finally a loss of sensitive species 32 .
Groundwater warming can also have positive effects as the accumulated thermal energy can be recycled through shallow, low-carbon geothermal energy systems 33 .Whereas studies typically focus on recycling the waste heat from anthropogenic sources, particularly from subsurface urban heat islands 34 , the subsurface heat accumulating due to climate change also has the potential to sustainably satisfy local heating demands 35 .However, increased warming will make cooling systems less efficient 36 .
Here we develop and apply a global-scale heat-transport model (thermal diffusion) to quantify groundwater temperatures in space and time and their response to recent and projected climate change (Fig. 1a,b).Our objective is to reveal the potential magnitude and long-term implications of ongoing shallow groundwater warming and to identify 'hotspots' of concern.The model utilizes standard climate projections to drive global groundwater warming down to 100 m below ground surface but with a focus on temperatures at the depth of the water table.We discuss (1) where aquifer warming will influence the viability of shallow geothermal heat recycling in the shallow subsurface (Fig. 1c), (2) given how it impacts microbial activity and groundwater chemistry, where groundwater temperature may cross key thresholds https://doi.org/10.1038/s41561-024-01453-xrange profiles with temperature changes of 1.1 °C at the location in Australia, compared with 0.5 °C at the location in Nigeria.These effects of intra-annual and short-term interannual variations in weather are attenuated at greater depths (for example, 30 m).Long-term (climate change) effects penetrate deeper, although groundwater warming may be less pronounced with depth due to the time lag between surface and subsurface temperature signals (Fig. 3c).

Accumulated energy
The overall increase in GWT can be quantified as accumulated energy (Methods).By 2020, a net energy amount of 14 × 10 21 J has already been absorbed by the terrestrial subsurface (Fig. 4a; 119 (45, 202) MJ m −2 ) since the beginning of the industrial revolution.In comparison, 436 × 10 21 J or about 25 times more has been absorbed by the oceans over a similar time period 39 .A review of Earth's energy imbalance identifies a total heat gain of 358 × 10 21 J for the time period 1971-2018 only, attributing about 6% of that to land areas including permafrost regions (21 × 10 21 J, that is, a similar magnitude as our estimate) 40 .In a similar range is the 23.8 × 10 21 J that was stored in the continental landmass since 1960 following a recent study; 90% is from heat storage 41 .

Implications for drinking water quality
Whereas groundwater warming offers benefits for geothermal heating systems, the accumulated heat also threatens water quality.In many developing countries or in poor and rural areas within developed countries, groundwater may be consumed directly without treatment or storage.It may also indirectly impact temperatures of drinking water within pipes 42 .In these regions in particular, the changes in water chemistry or microbiology that are associated with groundwater warming have to be carefully considered.
According to the World Health Organization, only 18 of 125 countries have temperature guidelines for drinking water 43 .These temperature guidelines, which are often aesthetic guidelines, range from 15 °C to 34 °C, with a median of 25 °C.Figure 4b shows where annual maximum groundwater temperatures at the geothermal gradient inflection point, that is, the most conservative depth as it is the coldest point in the temperature-depth profile, are above these thresholds in 2020.At this time, more than 29 million people live in areas where our modelled maximum

Implications for groundwater-dependent ecosystems
The ecosystems most dependent on groundwater are those in the aquifers themselves.A temperature increase may threaten groundwater biodiversity and ecosystem services 44,45 .Also, the increased metabolic rates of microbes caused by warming will accelerate the cycling of organic and inorganic matter, additionally fuelled by the increasing importance of dissolved organic carbon to the subsurface 46 .Combined with decreasing groundwater recharge as projected for many North African, southern European and Latin American countries 47 , this may transform oxic subsurface environments into anoxic 24 .
Groundwater warming also threatens many riverine groundwaterdependent ecosystems and the industries (for example, fisheries) that they support 48 .To capitalize on past continental-scale research related to groundwater, river temperature and ecosystems, we compare our modelled spatial patterns of groundwater warming in the conterminous United States to a recent distributed analysis of 1,729 stream sites 49 .The amplitude and phase of seasonal temperature signals in these surface water bodies were used to reveal the thermal influence and source depth of groundwater discharge to these streams, with about 40% classified as groundwater dominated.Our results show that GWT at the water table for the groundwater-dominated stream sites increased by 0.1 (0.0, 0.4) °C between 2000 and 2020 and 1.3 (0.3, 2.6) °C and 1.9 (0.4, 4.5) °C between 2000 and 2100 following SSP 2-4.5 and SSP 5-8.5 median projections, respectively (Fig. 4c,f and Extended Data Fig. 3g).Twenty-fifth percentile projections reveal 0.7 (−0.1, 1.5) °C and 1.0 (0.0, 2.9) °C and 75th percentile projections 2.0 (0.  and 7).
The warming groundwater will inevitably raise the ambient temperature of surface water systems thermally influenced by groundwater discharge.Furthermore, such groundwater warming will even more strongly impact the thermal regimes of groundwater-fed thermal refuges (for example, at the outlets of springs or groundwater-dominated tributaries flowing into rivers) and cause them to more regularly cross critical temperature thresholds for resident species seeking relief from thermal stress.Given the connection between aquifer thermal regimes and river sediment temperatures 50 , groundwater warming also threatens the thermal suitability of benthic ecosystems and spawning areas for fish 51 , posing a major risk to fisheries and other dependent industries.

Summary and model application
In summary, global climate change is leading to increased atmospheric and surface water temperatures, both of which have already been assessed across spatial scales ranging from local to global.Here we contribute to the global analyses of environmental temperature change and of groundwater resources through the presentation of projected groundwater temperature change to 2100 at a global scale.Our analyses are based on reasonable hydraulic and thermal assumptions providing conservative estimates and allow for both the hindcasting and forecasting of groundwater temperatures.Future groundwater temperature forecasts are based on both SSP 2-4.5 and 5-8.5 climate scenarios.We provide global temperature maps at the depth of the water table, 5 and 30 m below land surface, and these highlight that places with shallow water tables and/or high rates of atmospheric warming will experience the highest groundwater warming rates globally.Importantly, given the vertical dimension of the subsurface, groundwater warming is inherently a three-dimensional (3D) phenomenon with increased lagging of warming with depth, making aquifer warming dynamics distinct from the warming of shallow or well-mixed surface water bodies.
To facilitate more detailed future analyses, the temperature maps are included in a Google Earth Engine app at https://susanneabenz. users.earthengine.app/view/subsurface-temperature-profiles.The gridded GWT output could be integrated with global river temperature models 52 to more holistically understand future warming in aquifers and connected rivers.Whereas the warming of Earth's groundwater poses some opportunities for geothermal energy production, it increasingly

Diffusive heat transport
We hindcast monthly subsurface temperatures (and therefore also groundwater temperatures (GWTs) based on the assumption of local equilibrium) from the surface to a depth of 100 m for the years 2000 to 2020.We also force our model with future projections following SSP 2-4.5 and SSP 5-8.5 up to the year 2100.Subsurface temperatures in the shallow crust are generally controlled by one-dimensional (1D) (vertical) diffusive heat transport.Heat advection due to water flow plays a lesser and often inconsequential role in controlling subsurface temperatures [54][55][56] , particularly at larger spatial scales that average out focused groundwater flows in faults and fractures and groundwater exchange with surface water bodies.We adopt our 1D diffusion-dominated approach rather than a 3D numerical model of coupled groundwater flow and heat transfer as there are presently neither the parameterization data nor the computing power to enable such a coupled, 3D water and thermal transport model at a global scale.Also, whereas the influence of heat advection on steady-state or transient, subsurface temperature-depth profiles can be detected with precise temperature loggers and yields valuable insight into vertical groundwater fluxes when heat is used as a groundwater tracer 57 , the rate of shallow groundwater warming is often not thought to be strongly influenced by typical basin-average, vertical groundwater flux rates.Accordingly, heat advection has been ignored in some past local-scale groundwater warming studies (for example, ref. 58).However, to further investigate the thermal effects of multi-dimensional flow, we run a suite of scenarios and find that advection only exerts a minor influence on groundwater warming rates for typical groundwater flow conditions (Supplementary Note 1), enabling us to employ our approach.Appropriate initial conditions can be far more important for reliable simulation of temperature-depth profiles than the inclusion of heat advection 59 .To ensure our initial conditions are not influenced by any preceding climate change, we initiate our model in 1880 when the industrial revolution had not yet increased atmospheric greenhouse gasses and the climate was relatively stable.As default initial setting, we define a temperature-depth profile that increases linearly with depth z from the surface T S in accordance with the geothermal gradient a: T(z) = T S + az (ref.55).In permafrost regions, warming above critical thresholds requires latent heat to thaw ground in addition to the sensible heat to raise the temperature.As we do not include latent heat effects, model results are not presented for permafrost regions 60 .
We use the following analytical solution to the transient 1D heat diffusion equation for a semi-infinite homogeneous medium subject to a series of n step changes in surface temperature 55 : where j is a step change counter (counting by month), t is time, T S (t) is the time series of the ground surface temperature, D is the thermal diffusivity and erfc is the complementary error function.This equation is often used in an inverse manner to reconstruct pre-observational ground surface temperature history from observed, deep temperature-depth profiles, demonstrating its utility for investigating the response of subsurface thermal regimes to surface warming.We run our model in Google Earth Engine (GEE) 61 , and the results are presented in the form of a Google Earth Engine app openly accessible at https://susanneabenz.users.earthengine.app/view/subsurface-temperature-profiles.The application presents zoomable maps of annual mean, maximum and minimum GWT at different depths and seasonal variability (maximum minus minimum) for selected years and climate scenarios.All datasets were created at a native 5 km resolution at Earth's surface.However, Google Earth Engine automatically rescales images shown on the map based on the zoom level of the user.Charts that represent temperatures at a given location at a 5 km scale are created by clicking on the map and can be exported in CSV, SVQ or PNG file formats.For all analyses showing annual mean data at the water table depth, we first calculate monthly temperatures at the associated monthly groundwater level before averaging the results.

Ground surface temperatures
We use two distinct ground surface temperature time series: (1) one for the analysis of current (2020) temperatures based primarily on the ERA-5 data 62 and (2) one for the analysis of projected changes based on CMIP6 data 37 .On the basis of available computational power and data, we are not able to utilize monthly temperatures for the entire time period between the years 1880 and 2100.Instead, we present monthly temperatures from 1981 onwards and annual mean temperatures for 1880.The threshold 1981 is selected as ERA-5 data were available in Google Earth Engine from this point on when developing the model.
As these data are input into the analytical step function model (equation ( 1)), we supplement them with mean temperatures of the early 1980s (that is, three-year mean 1981 to 1984) to reduce artefacts of the sudden onset of seasonal signals in our data.An example of the ground surface temperature time series is shown in Supplementary Fig. 11.
For the analysis of current GWT, we use monthly mean soil temperature at 0-7 cm depth for the years 1981 to 2022 based on the ERA-5-Land monthly average reanalysis product 62 to form the ground surface temperature boundary condition for equation (1).These data have a native resolution of 9 km at the surface and are available through the GEE data catalogue.We also used annual ground temperature anomalies of 1880 of the top layer following the Goddard Institute for Space Studies (GISS) atmospheric model E 63 .This dataset gives the temperature difference between 1880 and 1980 in a horizontal resolution of 4° × 5° (approximately 444 km × 555 km at the equator) and can be extracted from https://data.giss.nasa.gov/modelE/transient/Rc_ij.1.11.html.To obtain absolute temperatures of 1880, we subtract the anomalies from three-year mean temperatures (1981 to 1984) of the ERA-5 data.
Future projections of ground surface temperatures are based on monthly soil temperatures closest to the surface for scenarios SSP 2-4.5 and SSP 5-8.5 from the CMIP6 programme available from 2015 to 2100.Model selection and methodology follow previous work 64 , but were updated to CMIP6 based on availability.In total we use nine models: BCC Where available, we used data from the variant label r1i1p1f1; however, for GISS-E2-1-G and HadGEM3-GC31-LL, these were not available, and we had to use r1i1p1f2 or r1i1p1f3 instead.Furthermore NorESM2-MM was missing data for January 2015; thus, we replaced them with data from December 2014 from the historic scenario.Data were collected from the World Climate Research Programme at https://esgf-node.llnl.gov/search/cmip6/.In addition, monthly data of the historic scenario were prepared for January 1981 to December 2014 and the annual mean data for 1880.To account for the difference between the CMIP6 models and ERA-5 reanalysis, we adjust the CMIP6 outputs based on mean temperatures T from ERA-5 between 1981 and 2014 (that is, the overlap between ERA-5 and the CMIP6 historic scenario) for each of the CMIP6 models separately as follows: Temperatures are determined for each model before being presented as the median and the 25th and 75th percentiles.

Thermal diffusivity
For our analysis we use the ground thermal diffusivity D: where λ (W m −1 °C −1 ) is the bulk thermal conductivity and C V ( J m −3 °C −1 ) is the bulk volumetric heat capacity.Ground thermal conductivity and volumetric heat capacity for various water saturation values are derived following previous examples 35,65 .This method links λ and C V values for different soil and/or rock types following the VDI 4640 guidelines 66 to a global map of soil and/or rock type.This map is based on grain size information of the unconsolidated sediment map database (GUM) 67 .
Where there is no available sediment class, we link to soil type in GUM.When this is also not available, we rely on the global lithological map database (GLiM) 68 .All required datasets were uploaded to Google Earth Engine in their native resolution.For assigned values, refer to Supplementary Table 1.
We acknowledge that the distribution of subsurface thermal properties is heterogeneous.However, specific heat capacity and thermal conductivity for rocks are both well constrained to within less than half an order of magnitude 69,70 compared with the many orders of magnitude for hydraulic conductivity 71 .We also note that water saturation can change the individual thermal properties and have accordingly run our model for six example locations with three different diffusivity values: (1) a dry soil, (2) a moist soil (default) and (3) a water saturated soil (Supplementary Fig. 12).The influence of water saturation on thermal diffusivity can be complex as both the heat capacity and thermal conductivity increase with water content (equation ( 3)).Overall, for locations with unconsolidated material in the shallow subsurface, groundwater warming rates increase with water saturation.However, the effect is nonlinear and the overall impact of water saturation on the thermal diffusivity is negligible for relative saturation values between 0.5 and 1 (ref.72).A map of the diffusivity utilized here is given in Supplementary Fig. 13a.

Geothermal gradient
When advection is absent, the geothermal gradient a (°C m −1 ; equation ( 1)) is the rate of temperature change with depth due to the geothermal heat flow Q (W m −2 ) and thermal conductivity λ (W m −1 °C −1 ): with global values for λ derived as described earlier, and the mean heat flow Q available as a global 2° equal area grid (about 222 km at the equator) 73 .Due to their resolution, these data do not incorporate fractures and major faults, and we thus are not able to estimate groundwater temperatures at these locations properly.The grid was uploaded to GEE in its native resolution for analysis (Supplementary Fig. 13b).

Water table depth
Much of our analysis and interpretation focuses on the future projection of temperatures at the water table depth.We therefore use the results of a previously published global groundwater model 74,75 with a 30 sec grid (about 1 km at the equator) to obtain the mean water table depth for 2004 to 2014.These data are available as monthly averages that we uploaded to GEE in their native resolution.In temperate climates, the model underestimates the observed water table depth by 1.5 m, and we therefore set the minimum water table depth to 1.5 m as was done in a previous study 35 .Still, whereas the global-scale hydro(geo)logical model of Fan et al. 74,75 can reveal large-scale patterns, it is of limited use for small-scale analysis and must be used with caution.Hence we run additional information for best-and worst-case scenarios where we add or subtract 10 m to the depth of the water table (Supplementary Note 4).
To calculate mean annual GWTs at the water table, temperatures for each month were determined at the corresponding water table depth by setting z in equation ( 1) to this depth.Future changes of water table elevation are challenging to predict, and we therefore base our analysis on the assumption that future water table elevations are unchanging.If we assume that the water table will rise, then warming would be more extreme; should the water table lower, warming as projected here is overestimated.A more detailed discussion, modelling water table changes of ± 10 m, can be found in Supplementary Note 4.However, we note that a modelled temperature-depth profile (equation (1)) is not impacted by the choice of the water table depth, and thus the results at 10 and 30 m are independent of the water table model.

Model evaluation
To assess the performance of our GWT calculations, we use two datasets of measured GWT or borehole temperatures.First, we compare our data to (multi-)annual mean shallow GWTs introduced in Benz et al. 35 .These data comprise more than 8,000 individual locations, primarily in Europe, where GWTs were measured at least twice between 2000 and 2015 at less than 60 m depth.Measurements are filtered based on their seasonal radius, a measure describing if a well was observed uniformly over the seasons and mean temperatures are therefore free of seasonal bias 76 .Second, we compare our data to temperature-depth profiles from the Borehole Temperatures and Climate Reconstruction Database at https://geothermal.earth.lsa.umich.edu/core.html.For these data, an exact date and depth of measurement are known.We filter the database based on time of measurement and depth of the first measurement, using only data taken after the year 2000 and starting at less than 30 m depth, resulting in 72 borehole measurements.To evaluate the model, we compare it to the observed groundwater temperatures described above.We compare the shallow (multi-)annual mean temperatures to mean temperatures at 30 m depth (the middle between 0 m and 60 m, the maximum depth of the observations) between 2000 and 2015.
For the dataset of one-time borehole temperature-depth profiles, we compare the shallowest data points to temperatures from our model at the same depth (rounded to the nearest metre), month and year.

Example locations
We use six locations distributed over all latitudes as examples in many of our figures, with locations in Australia (longitude 149.12°, latitude −35.28°),Brazil (−47.92°,−15.77°),China (116.39°,39.90°), Mexico (−99.12°,19.46°), Norway (10.74°, 59.91°) and Nigeria (7.49°, 9.05°).For convenience, each point is at the location of the capital city.However, as our model is not able to adequately describe the impact of urban heat on measured groundwater temperatures, groundwater at these locations is expected to be warmer, potentially by several degrees.Our focus is on the rate of warming in response to climate change.

Depth of the geothermal gradient 'inflection point'
To find the depth d i down to which annual mean temperature-depth profiles T(z) are inverted (that is, decrease with depth as opposed to increase following the geothermal gradient 4 ), we find the maximum depth where T(d i ) > T(d i+1 ).Given our computational resources, we test this at a resolution of 1-m steps for the first 10 m, then in 5-m steps down to 50 m depth and lastly in 10-m steps down to the maximal depth of 100 m.

Accumulated energy
To quantify shallow subsurface accumulated energy I ( J m −2 ), we compare mean annual temperature-depth profiles down to 100 m depth to the initial conditions T(z) = T S (t = 1,880) + az by solving the following integral in 1-m steps: https://doi.org/10.1038/s41561-024-01453-x This analysis utilizes annual mean subsurface temperatures T(z) for 2020 or 2100 for the current and projected analyses, respectively.The volumetric heat capacity C V (z) of the unsaturated zone (for z above the water table) and the saturated zone (for z below the water table) uses discrete values given in Supplementary Table 1.

Drinking water temperature thresholds
To assess the impact of groundwater warming on drinking water resources, we compare annual maximum groundwater temperatures to thresholds for drinking water temperatures summarized by the World Health Organization 43 .We do so for temperatures at the depth of the thermal gradient inflection point, the coldest point in the temperature profile and thus a best-case scenario, and for the depth of the water table to capture the 6% to 20% of wells that are no more than 5 m deeper than the water table 77 .To quantify populations at risk of exceeding the threshold, we compare the resulting maps with population counts.
For temperatures in 2022, we use the 2015 United Nations-adjusted population density from the Population of World Version 4.11 Model 78 .
For future scenarios, we rely on the global population projection grids for 2100 from the SSPs 79,80 .These data are available through the socioeconomic data and applications centre.

Impact on surface water bodies
Temperatures in surface water bodies are strongly influenced by atmospheric heat fluxes, but groundwater discharge and other processes can decouple temperatures in the atmosphere and water column.In the United States, 1,729 stream sites have been analysed by Hare et al. 49 to determine the dominance of groundwater discharge and to ascertain the relative depth (shallow or deep) of the associated aquifers.We use these sites to extract changes in mean annual groundwater temperature at the depth of the water table from our results to assess the impact of groundwater warming on these surface water bodies.
8-on average, there is 68 (13, 133) MJ m −2 of heat in the global subsurface saturated zone in 2020.By comparing the accumulated aquifer thermal energy in the United States (about 45 MJ m −2 ) with local residential heating demands (about 35,000 MJ per household in 2015 following the US Energy Information Administration 2015 Energy Consumption Survey), we find that, if recycled, the energy accumulated below an average home (250 m 2 for the floor area in new single-family houses following the 2015 'Characteristics of new housing' report, US Department of Commerce) in 2020 would fulfil about four months of heating demands.However, by 2100, global heat storage in the saturated zone is projected to increase to 233 (75, 363) MJ m −2 following SSP 2-4.5 and 352 (105, 536) MJ m −2 following SSP 5-8.5 median projections (Extended Data Figs.8 and 9 for 25th and 75th percentile projections).

Fig. 1 |
Fig. 1 | Processes and impacts related to groundwater temperature changes.a-e, Increases in surface air and ground surface temperatures (a) drive increases in groundwater temperatures (b) that, in turn, impact the geothermal potential for shallow geothermal energy systems (c), groundwater chemistry

bFig. 2 |
Fig. 2 | Current groundwater temperatures.a-c, Map of modelled mean annual temperatures at the depth of the water table (a), at 5 m below ground surface (b) and at 30 m below ground surface (c) in 2020.d, Comparison of modelled and observed groundwater temperatures.Blue markers are (multi-) annual mean temperatures observed between 2000 and 2015 at an unspecified depth against modelled temperatures of the same time period at 30 m depth.Grey markers are temperatures of a single point in time versus modelled temperatures of the same time and depth.A histogram of the errors (observed minus modelled temperatures) is shown in the upper left corner.e, Modelled temperature-depth profiles showing mean annual temperatures and the seasonal envelope for the locations displayed in a.Please note that we use bulk thermal properties, and the water table depth is thus not an input parameter into our model.

Fig. 3 |Fig. 4 |
Fig. 3 | Change in groundwater temperatures between 2000 and 2020 and between 2000 and 2100 following SSP 2-4.5.a-d, Recent (2000 to 2020) changes.e-h, Projected (2000-2100) changes.a,e, Map of the change in annual mean temperature at the depth of the water table.The line in the legend indicates 0 °C.b,c,f,g, Temperature change 5 m below the land surface (b,f) and 30 m below the land surface (c,g).d,h, Change in temperatures between 2000 and 2020 (d) and difference between 2000 and 2100 (h) as depth profiles for selected locations (symbols in a and e).Lines in h indicate median projections, whereas 25th to 75th percentiles (pct.) are presented as shading.

Extended Data Fig. 1 |Extended Data Fig. 3 |Extended Data Fig. 6 |Extended Data Fig. 7 |. 9 |
Depth to the inflection point.Shown is the depth down to which we can trace the impact of climate change in form of inverted temperature-depth profiles, that is temperature is decreasing with depth and not increasing with depth as expected based on the geothermal gradient.a and b, The depth to the geothermal inflection point in 2020 and 2100 following SSP 2-4.5.c, The depth to the geothermal inflection point in 2100 following SSP 5-8.5.https://doi.org/10.1038/s41561-024-01453-xWater table depth Change 2000 to 2100 (SSP 2-4.5, 75 th percentile projections) No data (permafrost) Extended Data Fig. 2 | Change in groundwater temperatures following SSP 2-4.5, 25th and 75th percentile projections.a-f, Map of the change in annual mean temperature between 2000 and 2100 following SSP 2-4.5 at the depth of the water table (under consideration of its seasonal variation).Temperatures in 2000 are based on the historic CMIP6 scenario.The line in the legend indicates 0 ∘ C. b and e, Annual mean groundwater temperature 5 m below the surface.c and f, Annual mean groundwater temperature 30 m below the surface.a-c, Annual mean groundwater temperature 25th percentile projected changes.d-f, Annual mean groundwater temperature 75th percentile projected changes.Change in groundwater temperatures between 2000 and 2100 and implications following SSP 5-8.5.a, Map of the change in annual mean temperature between 2000 and 2100 following SSP 5-8.5 (median projections) at the depth of the water table (under consideration of its seasonal variation).Temperatures in 2000 are based on the historic CMIP6 scenario.The line in the legend indicates 0 ∘ C. b, temperature change 5 m below the surface, and c, 30 m below the surface.d, Change in temperatures between 2000 and 2100 as depth profiles for selected locations.Lines indicate median projections whereas 25th to 75th percentile are presented as shading.e, Accumulated heat down to 100 m depth.The line in the legend indicates 0 MJ per m 2 .f, Map showing locations where maximum monthly GWTs at the thermal gradient inflection point (that is coldest depth) in 2100 are above guidelines for drinking water temperatures (DWTs).g, GWT changes between 2000 and 2100 at stream sites with a groundwater signature.Extended Data Fig. 4 | Change in groundwater temperatures following SSP5-8.5,25th and 75th percentile projections.a and d, Map of the change in annual mean temperature between 2000 and 2100 following SSP5-8.5 at the depth of the water table (under consideration of its seasonal variation).Temperatures in 2000 are based on the historic CMIP6 scenario.The line in the legend indicates 0 ∘ C. b and e, Annual mean groundwater temperature 5 m below the surface.c and f, Annual mean groundwater temperature 30 m below the surface.a to c, Annual mean groundwater temperature 25th percentile projected changes.d to f, Annual mean groundwater temperature 75th percentile projected changes.Implication of groundwater warming for SSP 2-4.5 25th and 75th percentile projections.a and d, Accumulated heat down to 100 m depth for SSP 2-4.5 25th and 75th percentile projections, respectively.The line in the legend indicates 0 MJ per m 2 .b and e, Locations where maximum monthly GWTs at the thermal gradient inflection point (that is coldest depth) in 2100 are above guidelines for drinking water temperatures (DWTs) for SSP 2-4.5 25th and 75th percentile projections, respectively.c and f, GWT changes between 2000 and 2100 at stream sites with a groundwater signature for SSP 2-4.5 25th and 75th percentile projections, respectively.Implication of groundwater warming for SSP 5-8.5 25th and 75th percentile projections.a and d, Accumulated heat down to 100 m depth for SSP 5-8.5 25th and 75th percentile projections, respectively.The line in the legend indicates 0 MJ per m 2 .b and e, Locations where maximum monthly GWTs at the thermal gradient point (that is coldest depth) in 2100 are above guidelines for drinking water temperatures (DWTs) for SSP 5-8.5 25th and 75th percentile projections, respectively.c and f, GWT changes between 2000 and 2100 at stream sites with a groundwater signature for SSP 5-8.5 25th and 75th percentile projections, respectively.https://doi.org/10.1038/s41561-024-01453-xAccumulated heat in the saturated zone (defined as below the water table down to 100 m depth) and maximum temperatures (based on monthly GWTs) at the depth of the geothermal inflection point showing exceedence of guideline thresholds for drinking water temperatures (DWTs) for 25th and 75th percentile SSP projections.a and b, Accumulated heat in the saturated zone for SSP 2-4.5 25th and 75th percentile projections, respectively.c and d, Locations where maximum temperatures exceed guideline thresholds for drinking water temperatures (DWTs) for SSP 2-4.5 25th and 75th percentile projections, respectively.e and f, Accumulated heat in the saturated zone for SSP 5-8.5 25th and 75th percentile projections, respectively.g and h, Locations where maximum temperatures exceed guideline thresholds for DWTs for SSP 5-8.5 25th and 75th percentile projections, respectively.
GWT exceeded 34 °C.If water is extracted at the depth of the water table, this increases to close to 31 million (Extended Data Fig.10). https://doi.org/10.1038/s41561-024-01453-x Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.© The Author(s) 2024 https://doi.org/10.1038/s41561-024-01453-x