Climate Change Yields Groundwater Warming in Bavaria, Germany

Thermodynamic coupling between atmosphere and ground yields increasing aquifer temperatures as a consequence of global warming. While this is expected to manifest as a gradual warming in groundwater temperature time series, such continuous long-term recordings are scarce. As an alternative, the present work examines the use of repeated temperature-depth profiles of 32 wells in southern Germany, that were logged during campaigns in the early 1990s and in 2019. It is revealed that the temperatures have increased in nearly all cases. We find a moderate to good depth-dependent correlation to trends in air temperature, which however is strongly influenced by local hydrogeological and climate conditions. While during the last three decades, air temperatures have increased by a rate of 0.35 K (10a)−1 on average, the temperature increase in the subsurface is decreasing with depth, with median values of 0.28 K (10a)−1 in 20 m and only of 0.09 K (10a)−1 in 60 m depth. Still, the slow and damped warming of the groundwater bodies are remarkable, especially considering naturally very minor temperature changes in pristine groundwater bodies and predictions of atmospheric temperatures. This entails implications for temperature-dependent ecological and hydro-chemical processes, and also for the heat stored in the shallow ground. Moreover, it is demonstrated that the annual heat gain in the groundwater bodies below 15 m due to climate change is in the range of one third of the state’s heat demand, which underlines the geothermal potential associated with the change in natural heat fluxes at the ground surface.


INTRODUCTION
Global warming is one of the most pressing challenges in the 21st century. Atmospheric climate and temperature variations are excessively studied, and a gigantic amount of climate data are continuously collected to delineate past warming trends and to predict future impacts of greenhouse gas accumulation. There is much less attention to the subsurface thermal regime, which is thermodynamically coupled with the atmosphere and thus also influenced by climate change. In fact, borehole temperature depth-profiles that deviate from the normal geothermal gradient served as early witnesses of regional atmospheric temperature increase (Wang and Lewis, 1992;Pollack and Chapman, 1993). This has in particular been studied in boreholes where vertical conduction is the dominant heat transport process and groundwater flow can be neglected. In contrast, temperature profiles recorded in wells are favored to infer the role and intensity of vertical groundwater flow (Bredehoeft and Papaopulos, 1965;Sorey, 1971;Taniguchi, 1993;Taniguchi et al., 1999;Li et al., 2019). The use of natural temperature variations in aquifers has also been recognized as a precious tracer to understand groundwater flow systems, in particular when interacting with surface waters (Constantz, 2008;Rau et al., 2010;Molina-Giraldo et al., 2011;Saar, 2011;Coluccio and Morgan, 2019;Kaandorp et al., 2019).
Temperature profiles measured in wells reflect spatially and temporally varying heat inputs in aquifers from the surface and thus can be used to examine thermal coupling at the ground surface (Gunawardhana and Kazama, 2011;Kurylyk et al., 2013;Burns et al., 2016;. Especially when the shallow subsurface is dominated by horizontal flow, changes in atmospheric temperatures and land use represent thermal signals that are conduced to the aquifer and become visible in well-logs. These changes are pronounced in cities, where accelerated heat flux from urban warming, sealed ground, and buried infrastructures yields large scale subsurface urban heat islands (Ferguson and Woodbury, 2004;Menberg et al., 2013;Zhu et al., 2015;Benz et al., 2016;Epting et al., 2017;Bayer et al., 2019;Hemmerle et al., 2019), and urban, industrial and waste sites are revealed to cause the most prominent local heat anomalies in Central European aquifers (Tissen et al., 2019). In less-disturbed rural areas, groundwater temperatures are reported to slowly increase as well, which is obviously the response of the shallow ground to recent climate change (Maxwell and Kollet, 2008;Bloomfield et al., 2013;Kurylyk et al., 2014;Menberg et al., 2014;Colombani et al., 2016). In comparison to atmospheric temperature recordings, however, continuously monitored longterm time series of groundwater temperatures are scarce, often strongly superimposed by local effects and thus our current picture of subsurface warming due to climate change is not very clear. As alternative means, repeated transient temperature logs in wells can be used Benz et al., 2018a). Conductive heat transport through the ground and aquifer is slow, however, the time period between two measurements has to be in the order of years and ideally decades.
In this study, the main objective is to reveal subsurface temperature variations. We report the findings of rare well temperature profiles measured again after more than 25 years. Focus is set on a large spatial coverage with several wells in order to reveal regional groundwater warming and to minimize the fingerprint of local effects. The wells are all located in the country of Bavaria in southern Germany, where climate change effects on groundwater are broadly discussed. As in many other countries, aquifers represent the major source of freshwater and thus any factors that impair groundwater reservoirs are of prime hydrological interest. Therefore, related work and recent reviews of groundwater climate change focus mainly on change of groundwater recharge and availability (Earman and Dettinger, 2011;Green et al., 2011;Stoll et al., 2011;Alam et al., 2019;Bloomfield et al., 2019;Zhang et al., 2020). Potential hydrochemical effects from climate-driven shifts in groundwater levels on contaminant mobilization are stressed by Jarsjö et al. (2020). The consequences of ongoing shallow groundwater warming are only rarely discussed but manifold (Riedel, 2019). Among these are fundamental changes to microbial activity and ecosystems in groundwater that rely on long-term stable thermal conditions (Kløve et al., 2014;Griebler et al., 2019), as well as potential threats on spring and cave ecosystems (Jyväsjärvi et al., 2015;Mammola et al., 2019). Leonhardt et al. (2017) emphasize for the springs in the Bavarian National Parks that climate change effects are most critical, and just recently therefore a program of long-term monitoring of biotic and abiotic data including temperature at 15 of these springs has been initiated by the authorities.
In the following, we describe the temperature profile measurements in Bavarian wells, which serve as the basis of this work. The presented data preparation is tailored to the purpose of this study, which is on quantifying the long-term thermal evolution of the aquifers on the country scale. This is examined by including available air temperature records and comparing the changes in the atmosphere with those found in the subsurface. Furthermore, in order to elucidate the heat accumulation associated with climate induced warming of the subsurface, also the change of the regional geothermal potential is assessed.

Groundwater Temperature Data
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.

Air Temperature Data
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.

Theoretical Geothermal Potential
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 . This represents the heat in place and can be calculated based on the caloric equation of state for a water saturated homogeneous solid: where n is the porosity of the solid, c w and c s in kJ (m 3 K) −1 are the volumetric heat capacities of water and solid, V (m 3 ) 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 · c w + (1 − n) · c s , 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 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.

RESULTS AND DISCUSSION Recent Variations in Air Temperature
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 10 ), 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) −1 . Compared to the linear regression slope of the previous 30 years period of 0.14 ± 0.07 K (10a) −1 , 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 , which results in a decadal temperature change of also 0.35 K (10a) −1 .

Recent Variations in Groundwater Temperature
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 temperatureswhich 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 10 ) with respect to the time span between the periods. The decadal temperature change rates per depth are depicted as boxplots in Figure 2B. Riedel (2019) also found a quasi-depth dependency with temperatures in spring waters averaging at a higher ΔT 10 rate of 0.3 K (10a) −1 over a rate of 0.2 K (10a) −1 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 10 ∼ 0.37 ± 0.42 K (10a) −1 ). 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.

Correlation Between Air and Groundwater Temperature Change Rates
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 10 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 10 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.

Shallow Geothermal Potential of Recent Temperature Variations
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 3 K) −1 ). 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).

CONCLUSION
The scope of this study was to extract regional long-term trends of groundwater temperature, and for this purpose, well profiles remeasured after up to 27 years were compared to each other. In order to remove as much as possible momentary, local variability, we compared profiles that were obtained as averaged logs from repeated campaigns throughout each year. The total number of suitable wells was 32 which cover broadly the entire state of Bavaria, Germany, and these offered a substantial insight into the shallow thermal regime in the subsurface. Clearly, each profile was different, influenced by the local hydrogeological, climatic and potentially anthropogenic conditions. However, short term variability and seasonality was shown to be marginal for these wells beneath 15 m, and this is where long-term climate change left an observable thermal imprint. Nearly all wells revealed increasing groundwater temperatures. The increase and the discrepancy among the well profiles was largest at shallow depth, and a maximum warming of 1.5 K was observed. Only one well showed a slight cooling by 0.2 K, but on average, 0.7 K warmer groundwater is found at 15 m depth which declines with depth along the profile. Over the past ∼30 years, the warming climate transferred heat in the subsurface at a rate that is estimated roughly more than one quarter of the state's annual heating demand. Rising ground (water) temperatures increase shallow geothermal system efficiency, and accelerated ground heat flux represents a form of thermal recharge of shallow geothermal reservoirs. This ground heat gain thus can be considered as favourable for geothermal use. In contrast, groundwater ecosystems used to nearly static thermal conditions need to adjust and may be under increasing stress. Even if the changes are considered marginal so far, and impacts are considered acceptable, the groundwater temperature will further increase as delayed response to the past temperature changes, and are likely to rise further in response to future atmospheric warming.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.