Groundwater temperature data is recorded in Bavaria by logging temperatures in observation wells using a water level meter with a temperature sensor. Data was gathered during the period from 1992 to 1994 and in 2019. During the 1992–4 period, wells were sampled every 3 months on an irregular basis within one year, so that up to four measurements were obtained for each well. These measurements were performed by the Landesamt für Umwelt Bayern (LfU Bayern), which is the environmental state office of the federal state of Bavaria in Germany. In 2019, the same set of wells were measured in April and July/August field campaigns using a water level meter type 120 – LTC from HT Hydrotechnik GmbH. For the 2019 dataset, precision and accuracy of the temperature sensor is 0.1 K. For the 1992–4 sampling campaign, the precision is also 0.1 K. The accuracy for this dataset is assumed to be in the same magnitude, assuming the use of standard devices even though no detailed information on the applied loggers is preserved. The vertical resolution of sampling was highest close to the surface where temperature variations are commonly strongest. Temperatures were recorded in 1 m steps above 20 m, in 2 m steps between 20 and 40 m and in 10 m steps below 40 m for the 1992–4 campaign. In 2019, temperature was sensed at 0.5 m intervals above 10 m, at 1 m intervals between 10 and 40 m, at 2 m intervals between 40 and 60 m and at 5 m intervals below 60 m. For further analysis, all measured temperatures were then interpolated linearly at 1 m resolution.

Initially the data set measured in 1992–4 covered 347 temperature-depth profiles (209 of the wells are shallower than 10 m). From these, 95 wells had a complete seasonal record of four seasonal measurements. Temperature data from these 95 wells have been made available by the LfU Bayern for this study. Out of these, 46 wells were measured again in 2019. The other wells were either shallower than 15 m, have been dismissed from the observation well network (renaturalised), could not be found in the field or could not be accessed (no access permission for wells on private property).

For this study we applied a threshold for seasonal temperature variation of 0.5 K below 20 m. 12 wells that exceeded this threshold were cut for further analysis. These wells showed inconsistencies or clearly revealed the effect of disturbing processes either in the 1992–4 or in the 2019 measurements. In one of the wells abnormal changes of more than 4 K between the 1990s and 2019 period together with an inconsistent signal in 2019 were related to strong super positioning of local anthropogenic disturbances and thus these measurements were also excluded. At one well, the groundwater table was below 100 m depth, and in one well the probe got stuck in both 2019 measurements after less than 2 m below groundwater level. These two wells have also been excluded for further analysis.

The 32 remaining observation wells are located in different hydrogeological settings and distributed all over Bavaria. A map showing the location of the observation wells and air temperature stations is given in Figure 1A. Temperature-depth profiles are available as Supplementary Material to this article and are additionally archived in World Data Center PANGAEA.Temperature-depth profiles are also available at: https://doi.pangaea.de/10.1594/PANGAEA.923529.

To derive atmospheric temperature trends, we used arithmetic annual mean values from 23 air temperature stations in Bavaria. These data are provided by the Climate Data Center (CDC) run by the Deutscher Wetterdienst (DWD). According to recommendation of the World Meteorological Organization (World Meteorological Organization, 2017) air temperature data was split into two 30-year periods. Herein after referred to as reference period 1960–1989, and study period 1990–2019. Air temperature stations were selected for having consistent, gap-free records from 1960 onwards. Detailed descriptions and metadata for each weather station can be found on the accessed FTP server by the DWD.

Temperature variations in the subsurface also change the total thermal energy stored in the subsurface. From the perspective of shallow geothermal energy utilization, this energy determines the theoretical geothermal potential (Bayer et al., 2019). This represents the heat in place and can be calculated based on the caloric equation of state for a water saturated homogeneous solid:E  = [n⋅cw+ (1−n)⋅cs]⋅V⋅ Twhere n is the porosity of the solid, c_w and c_s in kJ (m³K)⁻¹ are the volumetric heat capacities of water and solid, V (m³) is the reservoir volume and T (K) is the temperature of the reservoir. This equation can be condensed by distinction of a component related to material properties, m= n⋅cw+ (1−n)⋅cs, and a part describing the thermal field of the reservoir, τ =V⋅T. If we assume a laterally homogeneous reservoir where temperature only varies with depth (z), τ becomes  τ=A⋅∫T(z) dz =A⋅θ, where A is the area and θ is the integral of the temperature over depth. Thus, the heat in place is E=m⋅A⋅θ. Accordingly, the change of the thermal energy stored in the subsurface with constant material properties can be described as ΔE(τ)=m⋅A⋅Δθ.

We also tested a different approach where the arithmetic mean difference profile is linearly fitted to each observation profile. This approach is based on standard statistical measures, as related spatial temperature profile analysis is scarce, there is no earlier work known following the same approach. In this model, profiles are calibrated for the lowest root-mean squared error (RMSE) between the reference profile and the measured ΔT profiles at each station. This allows for the possibility to get both a spatial and a total depth coverage at each observation well. The resulting profiles were integrated over the depths between 15 and 100 m and afterward interpolated via inverse distance weighting on a 1 km × 1 km grid with a total of 70,610 grid cells for the state of Bavaria.

Recent trends in air temperature are calculated from annual mean values of 23 air temperature stations, uniformly distributed over Bavaria (Figure 1A). Figure 1B shows the temperature record since 1960 normalized by the corresponding mean temperature of the reference period from 1960 to 1989 for each station. Annual mean air temperature records for individual stations exhibit a high year-to-year variability of up to 3 K, but also a high consistency of the trendlines among each other. This indicates that relative air temperature changes behave uniformly within the studied area and are forced by a large-scale climatic regime. Calculating a linear regression over longer periods allows to deduce the decadal temperature change rate (ΔT₁₀), which is the slope of the linear regression per decade. For the study period of 1990–2019, the linear regression of these air temperature stations yields a slope of 0.35 ± 0.11 K (10a)⁻¹. Compared to the linear regression slope of the previous 30 years period of 0.14 ± 0.07 K (10a)⁻¹, the rate of temperature change for the study period is significantly elevated. The change in slope of the linear regression is consistent with the median temperature change of 1.06 K between 1990 and 2019 and the subsequent 30 years (1960–1989), which results in a decadal temperature change of also 0.35 K (10a)⁻¹.

To reveal subsurface temperature variations, we will first describe recent groundwater temperature variations in the wells individually. Recorded temperature profiles reflect a limited section of the subsurface. This section is limited by the groundwater table toward the top and the drilling depth of the observation well to the bottom. From the seasonal measurements 1992–4 period, we can confer seasonal temperature variation by depth. For depths shallower than 15 m, seasonal variations are higher than 0.1 K and thus measurable. Supplementary Figure S1 shows the mean values of the ranges (min-max differences) per depth at the selected 32 stations.

Vertical temperature changes below 15 m commonly follow the local geothermal gradient according to the basal heat flux from the interior of the Earth. They also conserve thermal signatures over a longer time period, e.g., from land-surface temperature changes (e.g., in urban areas (Banks et al., 2009) as well as from changes in regional groundwater temperatures – which closely reflects mean annual air temperatures in the recharge area (Burns et al., 2017)). Changes in repeated temperature profiles measured over timeframes of more than a year can therefore reveal long-term temperature variations, such as those related to global warming. The temperature difference, ΔT(z), between the 1992–4 and the 2019 period is calculated from mean temperature profiles for each time period (cf., Supplementary Figure S2). This yields ΔT-profiles for each individual station depicted in Figure 2A. To deal with alternating sampling dates in the 1992–4 period and for better comparability, the temperature change is normalized to the decadal temperature change rate (ΔT₁₀) with respect to the time span between the periods. The decadal temperature change rates per depth are depicted as boxplots in Figure 2B. At a depth of 20 m, ΔT₁₀ varies between −0.01 and 0.49 with a mean of 0.28 K (10a)⁻¹. Toward greater depths, the mean ΔT₁₀ decreases to 0.16, 0.09, 0.07, 0.05 K (10a)⁻¹ at 40, 60, 80 and 100 m depth, respectively.

The proposed rates are in line with recent single depth groundwater observations from the western adjacent federal state of Baden-Württemberg and observations in Austria in the south. For Baden-Württemberg, Riedel (2019) reported ΔT₁₀ rates of 0.1–0.4 K (10a)⁻¹ based on a state-wide groundwater quality dataset for the period from 2000 to 2015. Riedel (2019) also found a quasi-depth dependency with temperatures in spring waters averaging at a higher ΔT₁₀ rate of 0.3 K (10a)⁻¹ over a rate of 0.2 K (10a)⁻¹ in groundwater. In Austria, Benz et al. (2018b) observed a slightly higher temperature change rate of 0.7 ± 0.8 K over the 20-year period from 1994 to 2013 (ΔT₁₀ ∼ 0.37 ± 0.42 K (10a)⁻¹). In both studies, temperatures change rates are not linked to depths. Despite this it has to be noted that the average measurement depth in the Austrian wells is relatively shallow at 7 ± 4 m below ground surface.

To infer the correlation between air and groundwater temperature variations, the air temperature at each observation well is calculated by inverse distance weighting of the five nearest air temperature stations. Figure 3 shows the temperature change per decade for each of the 32 observation wells for air temperature and the mean value of 10 m increments of the individual temperature profiles. For the respective depths vs. air temperatures, Pearson correlation coefficients (r) indicate a moderate but robust correlation with values between 0.4 and 0.6. P-values exhibit that statistical significance is increasing with depth despite a high-variability. In general, higher temperature change rates in air temperature are reflected by higher temperature change rates in the subsurface. However, local variations are high both in the penetration depth of the signal and in absolute values. Note that air temperature change rates have a low standard deviation of 0.8 K, which makes it hard to identify distinct variations between temperature stations, as the atmospheric temperature signal is more homogeneous than that in the subsurface. The diffuse signal present in the subsurface is hereby hard to infer with the minor variations in air temperatures (Figure 3A). The correlation metrics are highest for the ΔT₁₀ at 100 m and inferred air temperature. Assuming a process where the thermal field of the subsurface is mainly altered by a change in ground surface conditions (e.g., for shallow and unconfined aquifers), we would expect the contrary with decreasing metrics at greater depth induced by changing ground surface conditions. Obviously, correlating surface or air temperatures directly to subsurface temperatures presumes that changes in the temperature profile are a result of downward vertical heat transfer (advective or diffusive). If this is not the case, temperatures especially at greater depths could be altered by lateral heat transfer associated with the regional groundwater flow regime (Taniguchi et al., 1999; Zhu et al., 2015; Burns et al., 2016; Bense et al., 2020). In this case the temperature profile would carry an integrated signal dominated by annual mean temperatures of the recharge area. However, despite the potential influence of local hydrogeological conditions, it is also possible that the observed discrepancy is an artifact of the statistical analysis as the number of observations decreases with depth (from 24 to 16).

The substantial number of wells included in this study allows robust average estimates. Compared with the corresponding aboveground air temperatures, as expected, the decadal temperature increase is dampened in the subsurface. The median values for Bavaria are ΔT₁₀ = 0.35 K for the atmosphere, 0.28 K at 20 m and 0.09 K at 60 m depth. This means that at the deepest level, the changes, however, are close to measurement accuracy. Temperatures in air and subsurface show a moderate to good correlation, with higher scattering of groundwater than air temperatures. Major tendencies of air temperature are also reflected in the subsurface temperature record. This supports a good statistical basis of the measured data, even without detailed analysis of the site-specific heat transport processes.

The changes in subsurface temperatures also constitute an increase in thermal energy stored in the subsurface in response to recent shifts in the thermal regime of surface conditions. This additional heat can be accessed via shallow geothermal systems and thus increases the theoretical geothermal potential of the shallow subsurface. To accurately quantify this additional heat, both geological material properties and the thermal field have to be characterized at a high level. However, thermal properties such as heat capacity do not vary over a broad range in geological media (Stauffer et al., 2013). If we assume bulk material properties and superimpose an arithmetic mean ΔT profile over the depth interval between 15 and 100 m below surface, we can calculate a rough estimate of the magnitude of the heat flow into the subsurface. Table 1 lists the energy stored per year for a minimum, a maximum and median (50%) assumption of porosity (0.5–0.15) and volumetric heat capacity of the solid (1900–2,500 kJ (m³K)⁻¹). The energy stored per year varies between 155 (±125) and 212 (±171) PJ between the minimum and maximum scenarios (±1σ). For the median scenario the thermal energy stored by climate change in the subsurface per year equals 184 (±149) PJ. Compared to the primary energy demand (1944 PJ) and the heating demand (669.7 PJ) in 2017 (Ebert and Voigtländer, 2019), the energy stored per year in the subsurface for the median scenario equals 9.5% of the annual primary energy demand or 27.5% of the heating demand of Bavaria. This simplistic calculation, however, does not consider spatial variations and relies on the arithmetic mean differences of all profiles.

For an alternative approach where the arithmetic mean difference profile is linearly fitted to each observation profile, stored energies for the minimum and maximum scenario are 208 and 285 PJ with RMSEs of 40 and 54 PJ, accordingly. The median scenario yields a thermal energy input of 248 PJ with an RMSE of 47 PJ, which would equal 12.8% of the Bavarian primary energy demand or 37% of the heating demand. Supplementary Figures S3, S4 display the profile fitting and the spatial distribution of interpolated groundwater and air temperatures. Note that the spatial coverage of the observation wells with respect to the area of Bavaria is insufficient to produce fully reliable numbers, but it provides a first rough estimate of the order of magnitude of energy stored in the subsurface annually. Still, the presented spatial approach offers the opportunity to employ more precise regionalized material parameters for heterogeneous hydrogeological facies. In essence, we found the thermal energy stored by climate change in the subsurface per year to be in the magnitude of 10% of the total Bavarian primary energy demand, or one third of the heating demand. This continuous heat transfer represents an enormous replenishing resource that is fueled by climate change. The results from this local study are expected be applicable in similar magnitudes on a global scale as observed atmospheric temperature trends in Bavaria are in line with the global trend and the thermal connection between the surface and shallow subsurface has been proven to be robust also on a global scale (Benz et al., 2017).