Interactive comment on “ Observed groundwater temperature response to recent climate change ” by K .

Climate change is known to have a considerable influence on many components of the hydrological cycle. Yet, the implications for groundwater temperature, as an important driver for groundwater quality, thermal use and storage, are not yet comprehensively understood. Furthermore, few studies have examined the implications of climate-change-induced groundwater temperature rise for groundwater-dependent ecosystems. Here, we examine the coupling of atmospheric and groundwater warming by employing stochastic and deterministic models. Firstly, several decades of temperature time series are statistically analyzed with regard to climate regime shifts (CRSs) in the long-term mean. The observed increases in shallow groundwater temperatures can be associated with preceding positive shifts in regional surface air temperatures, which are in turn linked to global air temperature changes. The temperature data are also analyzed with an analytical solution to the conduction–advection heat transfer equation to investigate how subsurface heat transfer processes control the propagation of the surface temperature signals into the subsurface. In three of the four monitoring wells, the predicted groundwater temperature increases driven by the regime shifts at the surface boundary condition generally concur with the observed groundwater temperature trends. Due to complex interactions at the ground surface and the heat capacity of the unsaturated zone, the thermal signals from distinct changes in air temperature are damped and delayed in the subsurface, causing a more gradual increase in groundwater temperatures. These signals can have a significant impact on large-scale groundwater temperatures in shallow and economically important aquifers. These findings demonstrate that shallow groundwater temperatures have responded rapidly to recent climate change and thus provide insight into the vulnerability of aquifers and groundwater-dependent ecosystems to future climate change.


Introduction
Atmospheric climate change is expected to have a significant influence on subsurface hydrological and thermal processes (e.g., Bates et al., 2008;Green et al., 2011;Gunawardhana and Kazama, 2012).While the consequences for groundwater recharge and water availability have been scrutinized by many studies (e.g., Maxwell and Kollet, 2008;Ferguson and Maxwell, 2010;Stoll et al., 2011;Taylor et al., 2013;Kurylyk and MacQuarrie, 2013), there is still a lack of studies regarding the observation and analysis of the effects of recent climate change on shallow groundwater temperatures (Kløve et al., 2013).Groundwater temperature (GWT) is known to be an important driver for water quality (e.g., Green et al., 2011;Sharma et al., 2012;Hähnlein et al., 2013), and therefore it is a crucial parameter for groundwater resource quality management (Figura et al., 2011).
Furthermore, increasing groundwater temperatures can have a significant influence on groundwater and river ecology (e.g., Kløve et al., 2013).Numerous studies on the impact of recent or projected climate change on the thermal regimes of surface water bodies and the associated impact for cold-water fish habitats have already been conducted (e.g., Kaushal et al., 2010;van Vliet et al., 2011van Vliet et al., , 2013;;Wenger et al., 2011;Isaak et al., 2012;Wu et al., 2012;Jones et al., 2014)

, but
Published by Copernicus Publications on behalf of the European Geosciences Union.
the thermal sensitivity of shallow aquifers to climate change is a relatively unstudied phenomenon (e.g., Brielmann et al., 2009Brielmann et al., , 2011;;Taylor and Stefan, 2009;Kurylyk et al., 2013Kurylyk et al., , 2014a)).The thermal response of GWT to climate change is of particular interest to river temperature analysts, as the thermal regimes of base-flow-dominated streams or rivers and hydraulically connected aquifers are inextricable linked (Hayashi and Rosenberry, 2002;Tague et al., 2007;Risley et al., 2010).Furthermore, groundwater-sourced cold-water plumes within river mainstreams are known to provide thermal refuge for threatened cold-water fish (e.g., Ebersole et al., 2001;Breau et al., 2007), and questions have arisen regarding the sustainability of these groundwater-dependent ecosystems (GDEs) in a warming climate (Deitchman and Loheide, 2012).The current lack of knowledge regarding the thermal vulnerability of GDEs to the climate-changeinduced warming of shallow GWT has been highlighted as a research gap in several recent studies (e.g., Bertrand et al., 2012;Mayer, 2012;Kanno et al., 2014).
Thermal signals arising from changes in ground surface temperatures (GSTs) propagate downward into the subsurface, causing GWT to deviate from the undisturbed geothermal gradient.Heat transport theory has been applied for inverse modeling of temperature-depth profiles to infer paleoclimates based on measured deviations from the geothermal gradient (e.g., Mareschal and Beltrami, 1992;Pollack et al., 1998;Beltrami et al., 2006;Bodri and Cermak, 2007) and for forward modeling the impact of projected climate change on measured temperature-depth profiles (e.g., Gunawardhana and Kazama, 2011;Kurylyk and MacQuarrie, 2014).Such studies are often based on the assumption that long-term trends in GST will track long-term trends in surface air temperature (SAT), although this has been a matter of considerable debate (e.g., Mann and Schmidt, 2003;Chapman et al., 2004;Schmidt and Mann, 2004).For example, decreases in the duration of thickness of the insulating winter snowpack due to rising SAT can paradoxically lead to decreased winter GST (Smerdon et al., 2004;Zhang et al., 2005;Mellander et al., 2007;Mann et al., 2009;Kurylyk et al., 2013), and thus cause a decoupling of mean annual SAT and GST trends.
Heat advection due to groundwater flow may also perturb subsurface temperature-depth profiles, and it can be difficult to determine if deviations from a geothermal gradient have arisen from past climate change or from groundwater flow (Reiter, 2005;Ferguson and Woodbury, 2005;Ferguson et al., 2006).Thus, several analytical solutions have been proposed that account for subsurface thermal perturbations arising from a combination of climate change and vertical groundwater flow (e.g., Taniguchi et al., 1999a, b;Kurylyk and MacQuarrie, 2014).The solutions vary depending on the nature of the surface boundary conditions employed (e.g., linear, exponential, or step trends in temperature), which can be used to match measured or predicted GST trends for a region.These solutions do not account for horizontal groundwater flow, which can also perturb subsurface thermal regimes in certain environments (Ferguson and Bense, 2011;Saar, 2011).Numerical solution techniques can also be applied to account for inhomogeneous subsurface thermal properties, complex surface temperature evolution, and groundwater flow (e.g., Kooi, 2008).Figura et al. (2011) show that temperature variations in Swiss aquifers that are recharged by river water through bank infiltration can be related to changes in climate oscillations systems by applying a statistical regime shift analysis.Characterizing changes in time series of various climatic, physical and biological parameters with the concept of abrupt regime shifts has been the focus of numerous studies in the last 2 decades (e.g., Hare and Mantua, 2000;Overland et al., 2008).In this context, a regime is often defined as a period with quasi-stable behavior or with a quantifiable quasiequilibrium state (deYoung et al., 2004), and accordingly a rapid transition between states with differing average characteristics over multi-annual to multi-decadal periods is referred to as a regime shift (Bakun, 2004).
In this study, we demonstrate the direct influence of atmospheric temperature development on shallow GWT at two sites in Germany by analyzing time series of SAT and GWT with regard to abrupt changes in the long-term annual mean.Compared to previous studies, which used borehole temperature profiles for the analysis of temperature coupling between the atmosphere and the subsurface, the measured annual GWT time series in this study allow for an evaluation of this coupling on a shorter timescale with a higher temporal resolution.Furthermore, we compare different spatially averaged temperature time series from individual weather stations to global mean air temperature change bringing our observations in the context of global climate change.The magnitudes of the regime shifts and the time lags between the shifts in the chosen time series are evaluated under consideration of the different thermal processes in the subsurface and the site-specific hydrogeological settings.A standard analytical solution to the conduction-advection subsurface heat transfer equation is applied to investigate the physical thermal processes underlying the observed correlation between SAT regime shifts and GWT rise.

Data and site description
For the analysis of shallow GWT, we use time series from four observation wells in porous and unconfined aquifers in Germany (Table 1, Fig. 1a and b).Two of the wells are installed in the surrounding area of Cologne outside the small villages of Dansweiler and Sinthern in agricultural areas.The other two wells are located in a rather densely vegetated forest, called Hardtwald, close to the city of Karlsruhe and are therefore named Hardtwald 1 and 2. The proximate surroundings of all four wells were undisturbed over the last decades, so that variations in GWT due to land use changes are unlikely.The distances from the observation wells to the nearest streams are several kilometers (Table 1); thus the influence of river water on the groundwater temperature in the wells can be excluded.The two study areas close to Karlsruhe and Cologne are located approx.240 km apart from each other and belong to different aquifer systems.Yet, the basic geological and hydrogeological settings of the two aquifers are rather similar (Tables 1 and 2).Table 2 lists some basic hydrogeological properties of the studied aquifers and the observation wells.The depth of water table differs considerably between the two well fields, and is approximately 17 m for the Cologne aquifer and 7 m near Karlsruhe.Variations in the depth of water table during the observation period are within ±1 m for the Dansweiler and Sinthern wells and more pronounced in the Hardtwald wells at ±3 m and are likely caused by a pumping station nearby.However, no statistically significant trend was observed over the last decades in the water level of the observation wells.Both aquifers are recharged by infiltration of meteoric water through the unsaturated zone with estimated recharge rates of 221 ± 4 5mm yr −1 for the Cologne aquifer and 228 ± 45 mm yr −1 for the aquifer near Karlsruhe (Table 2).A schematic cross section of the two aquifers near Cologne (left) and Karlsruhe (right) in Fig. 1c and d show the average depth of the water table below surface level and the depth of the underlying aquitard.Details on the wells' constructions are also depicted with the overall depth and the locations of the filter screens (black areas) that indicate the depth where the pumped water is captured.Furthermore, Fig. 1c and d shows the distance between the well pairs as well as the distances to the weather stations from which the SAT time series were obtained.
GWT in all observations wells was measured one to six times per year for a period of at least 32 years  during frequent water quality assessments by the local groundwater authorities.The measurement protocol, which is standardized by the environmental state agencies to ensure data quality and comparability, has undergone no significant changes in the last decades.During the specified procedure, water is pumped from the wells until the water temperature and other on-site parameters are constant.The temperature measurements are thereby conducted with a probe directly at the outlet, to minimize influences by ambient air temperatures.An examination of the time series for seasonal effects revealed that they contain certain minor seasonal effects with annual variations up to ±0.2 K, which indicates an impact of ambient air temperature on the GWT during the sampling.The natural temperature variations due to seasonal GST variations in depths of over 20 m (Table 2) are expected to be less than 0.1 K as can be demonstrated by Stallman's (1965) equation.In most years, at least two measurements per year were available, so that the arithmetic mean was adopted as an annual mean value to minimize such effects.It should also be noted that the measurement accuracy is in the range of ±0.1 K. Also changes in the measurement procedure, such as variations in the pumping rate or in the placement of the   pump within the well, as well as changes in the measurement equipment, can influence the measured GWT and were considered for the evaluation and interpretation of the data.
Annual SAT data are available from weather stations operated by the German Weather Service (DWD) outside the cities of Cologne and Karlsruhe in agricultural surroundings (Fig. 1a and b).Though located several kilometers from the observation wells, the SAT from these stations is expected to yield a good approximation for the development of SAT at the well sites.Furthermore, for the evaluation of abrupt shifts in the time series of SAT and GWT, the absolute temperature is only of minor importance, while the main focus is on the timing of the shifts and the temperature differences.For the comparison with air temperatures on a larger scale, we use time series of mean air temperature anomalies based on the reference period 1951-1980 from the NASA Goddard Institute for Space Studies (GISS) (e.g., Hansen et al., 2010).Of the spatially averaged temperature data sets available, we evaluate the annual global mean from land-surface air and sea-surface water temperature anomalies and the annual zonal mean for the Northern Hemisphere between 90 and 24 • N based on land-surface air temperature anomalies.

Regime shift analysis
There are several possibilities to statistically evaluate temperature changes in time series with rather simple functional forms.Seidel and Lanzante (2004) compared different ap-proaches (e.g., linear and flat steps models) and revealed that often time series of atmospheric temperatures can be represented more appropriately by models using breakpoints than by models assuming monotonic functions.Hence, we here apply a sequential t test analysis for regime shifts (STARS) to detect possible abrupt climate regime shifts (CRSs) in the temperature time series (Rodionov, 2004;Rodionov and Overland, 2005).The STARS method has been successfully used by recent studies to identify abrupt changes in the longterm mean of environmental time series (Marty, 2008;North et al., 2013) and GWT time series (Figura et al., 2011).STARS is a parametric test that can detect multiple regime shifts and needs no a priori assumption for the timing of possible shifts.Identification of a shift is based on the calculation of the regime shift index (RSI), which represents the cumulative sum of the normalized deviations from the mean value of a regime and thus reflects the confidence of a regime shift (Rodionov, 2004).For the regime shift analysis, several test parameters need to be adjusted to account for specific characteristics, such as the length of the tested time series.The target significance level in our analysis is set to 0.15, which corresponds to the p level of false positives.The actual p value of an identified shift between subsequent regimes is calculated separately with Student's t test.The cut-off length of the test corresponds to a low-pass filter, so that regimes with a shorter length are disregarded in the analysis (Rodionov and Overland, 2005).Here, we set the cut-off length to 10 years as atmospheric oscillations often occur at decadal intervals (Overland et al., 2008).Furthermore, the Huber weight parameter (set to 1 in our study) included in the STARS procedure improves the treatment of outliers by weighting them proportionally to their deviation from the mean value (Overland et al., 2008).As pointed out by Seidel and Lanzante (2004), atmospheric data tend to be highly temporally auto-correlated, so that especially in short time series, spurious regime shifts may be detected due to serial correlation (Rudnick and Davis, 2003).Therefore, we apply a pre-whitening procedure that removes the red noise component from the temperature time series prior to testing for a regime shift (Rodionov, 2006).To investigate the potential stationarity within detected regimes, the non-parametric Mann-Kendall test for the absence of trend is also applied to the temperature data (von Storch, 1995).

Analytical solutions
One common governing equation for studying transient subsurface heat transport is the one-dimensional conduction equation for homogeneous media, which equates the divergence of the conductive flux with the rate of the change of thermal energy in the medium (Carslaw and Jaeger, 1959;Domenico and Schwartz, 1990): where κ is the bulk thermal diffusivity of the subsurface (m 2 s −1 ), T is temperature ( • C), z is depth (m), and t is time (s).The governing heat transport equation becomes slightly more complex when advective heat transport (or "forced convection") due to groundwater flow is considered: where U (m s −1 ) is the thermal plume velocity under pure advection and a function of the Darcy velocity q (downwards or recharge is positive, m s −1 ), the bulk volumetric heat capacity of the soil-water matrix C (J m −3 • C −1 ), and the volumetric heat capacity of water C w (J m −3 • C −1 ): The governing conduction-advection Eq. ( 2) employs several limiting assumptions, including spatiotemporally constant groundwater velocity over the entire domain (including depths below the well screen), one-dimensional heat transport, homogenous thermal properties, constant pore-water phase, and isothermal conditions between the soil grains and pore water.Here we employ a distinct analytical solution to Eq. ( 2) to simulate the influence of a climate regime shift on GWT.We assume thermally uniform initial conditions and boundary conditions that are subject to a series of n step increases in GST: where T 0 is the initial uniform temperature ( • C) prior to the beginning of the regime shift, GST i is the step increase in GST for regime shift i ( • C), H is the Heaviside step function, and t i is the time (s) of the beginning of regime shift i.In this formulation, GST i refers to a step change in GST in comparison to the GST conditions immediately preceding that change (not necessarily in comparison to initial GST, T o ).We ignore short-term (e.g., annual) variations in SAT and GST and rather drive the subsurface heat transport models with temperatures averaged for a given climate regime and then instantaneously increased at the beginning of the next climate regime.The thermally uniform initial conditions are a reasonable assumption given that we begin by considering mean annual GWT at or near the water table following a relatively stable climate regime (i.e., prior to 1988, Fig. 2).Moreover, for the wells observed, the vadose zones and nearsurface aquifers are too shallow to realize the influence of any geothermal gradient.The isothermal condition assumption previously noted extends to the surface boundary, which implies that the groundwater recharge entering the semi-infinite domain at the ground surface has a temperature equal to the mean annual surface temperature for that climate regime.The transient conduction-advection heat transport model (TCA model) employed in this study is an analytical solution to the transient conduction-advection Eq. ( 2) subject to the initial and boundary conditions given in Eqs. ( 4) and (5).This solution was originally developed by Carslaw and Jaeger (1959) and subsequently employed by Taniguchi et al. (1999b) to study subsurface temperature evolution due to land cover changes in regions of significant groundwater flow.Because we assume initially thermally uniform conditions in the unsaturated zone and shallow groundwater, the resultant solution is simpler than in the original derivations.Unlike the original derivation, it is also presented here with superposition principles applied to allow for a series of regime shifts rather than one event.This superposition approach is valid given the linearity of the governing partial differential equation and the boundary and initial conditions (Farlow, 1982): where T (z, t) is the spatiotemporally varying subsurface temperature (GWT, • C), κ is the bulk thermal diffusivity of the subsurface (m 2 s −1 ), and erfc is the complementary error function.The Heaviside function indicates that the subsurface thermal influence of each regime shift i in the boundary condition is not realized until the time t exceeds t i .Comparisons between the model results and measured GWT indicate whether these simple analytical solutions are applicable for modeling the influence of observed and projected climate regime shifts in the wells considered in this study.
It should be noted that it is the GST rather than the SAT that drives subsurface thermal regimes and thus forms the boundary condition in Eq. ( 5).However, complete GST time series were not available for the locations considered in this study.Thus, in the present study, the magnitude and timing of the regime shifts in GST are obtained from the local SAT data as follows.In all cases, the timing of the GST regime shifts is assumed to correspond to the timing of the local SAT regime shifts for that location obtained from the statistical analysis.This approach is reasonable given the efficient heat transfer that occurs between the lower atmosphere and the ground surface (e.g., Bonan, 2008).The magnitude of the GST regime shift was set to be equal to the magnitude of the SAT.Measured SAT and GST data (not shown) indicate that this approach is valid as the measured magnitude of the climate regime shift in 1988 was 1.1 • C in both the SAT and GST data near Cologne.No GST data were available for the sites near Karlsruhe (Hardtwald sites, Fig. 1 and Table 2), but it is reasonable to assume that the magnitude of the GST changes track the magnitude of the SAT changes like in the case of the sites near Cologne.
Table 3 presents the assumed subsurface thermal properties for each well for both the saturated and unsaturated zones.A potential range in these values was estimated from literature values taking into account variations in lithology obtained from drilling logs as well as the variability of water content in the unsaturated zone ranging from dry to saturated conditions (VDI, 2010;Menberg et al., 2013).Because we consider temperature rise at different depths below the water table within the aquifer, the effective thermal diffusivities utilized in the analytical solution for each of the four locations were obtained from a weighted arithmetic average (weighted by zone thickness) of the saturated and unsaturated zone thermal diffusivities.For example, the unsaturated zone thickness was taken as the depth to the water table, and the saturated zone thickness was taken as the distance from the water table to a point along the well screen.Different points in the well screen were considered, as described in the results, because the vertical well capture zone flow dynamics may be complex depending on the nature of the pumping and heterogeneities in near-well hydraulic properties.This is particularly important for the Hardtwald wells, which have longer well screens than in the case of the Dansweiler or Sinthern wells (Fig. 1).
Regional recharge rates were extracted from Table 2 with a potential range to reflect the variability of recharge in this region over the last decades (Erftverband, 1995;W. Deinlein, personal communication, 2013).Similar thermal properties and recharge values are assumed for Hardtwald 1 and Hardtwald 2 based on their similar land cover and subsurface properties and the geographical proximity (about 200 m) between the wells.

Regime shifts in air and groundwater temperatures
At least two climate regime shifts (CRSs) could be detected in the later decades of all analyzed time series (Fig. 2).The time series of global mean temperature change and zonal mean temperature change in 90-24 • N show significant (STARS, p < 0.005) positive shifts in 1977in , 1987in , 1997in and 1977in , 1988in and 1998 (Table 4 (Table 4).The observation of shifts in air temperature change in these years is in good agreement with the observation of decadal shifts in atmospheric oscillation indices in the late 1970s, late 1980s and late 1990s (Overland et al., 2008)    local SAT data from Cologne and Karlsruhe.However, this is not surprising as previous studies observed that the CRS in the late 1970s was most prominent in the North Pacific region (Hare and Mantua, 2000;Overland et al., 2008), and less accentuated in Europe.The same applies to the CRS in the late 1990s (Overland et al., 2008;Swanson and Tsonis, 2009), which is reflected by the differing RSI values in Fig. 2. While the high RSI for the CRS in 1997 in the global mean temperature change indicates a significant shift, the RSIs for the late 1990s CRS in the German SAT time series are much lower than the RSIs in the late 1980s.Figura et al. ( 2011) correlated the abrupt increase in SAT in Switzerland with a change in the Arctic Oscillation (AO) that has a strong influence on air temperatures in Europe.However, no such change in the AO Index was found in the late 1990s, suggesting that the CRS in the German SAT is also coupled to the general air temperature increase in the Northern Hemisphere.Two regime shifts were detected in the GWT time series for the four wells near Cologne and Karlsruhe.These shifts correspond to the CRS in the atmosphere with a certain time lag (Fig. 2, Table 4).The regime shifts in GWT time series are all statistically significant (p < 0.01), except for the second regime shift in the late 1990s in Dansweiler.Two prominent outliers in the third regime of the time series influence the statistical significance for this shift, while the RSI value is calculated under consideration of the outliers according to the Huber weight parameter.Furthermore, the RSI values in Fig. 2 for the second shifts in Dansweiler and Sinthern are not the final values, as the 10-year cut-off length of the STARS test in the last regime has not yet been reached.In general, the time series of GWT show a more gradual increase than the SAT time series.In particular, the GWT in the Sinthern well appears to exhibit a linear trend rather than a step increase, which is discussed later.The GWT time series exhibit considerable inter-annual variability, which appears to be more significant in the Hardtwald wells than in Dansweiler and Sinthern.Potential reasons for these rather large fluctuations in annual GWT are related to the uncertainties associated with the measurements as mentioned in the method section.Other possible factors that influence the inter-annual variability could be the pumping station close to the wells in the Hardtwald, where groundwater is extracted at irregular intervals, and impacts by undetected land use changes in the close surroundings.

Statistical analysis of time lags and magnitude of temperature change
The time lags between the regime shifts in SAT and GWT are listed in Table 4.The regime shifts in global mean temperature change and the zonal mean in 90-24 • N occur simultaneously, except for the regime shift in the late 1980s that has a time lag of 1 year.However, as annual mean values are used for the analysis, the accuracy of the shift detection is limited to ±2 year, so that the shifts occur within the uncertainty range.The same applies to the first regime shifts in the local SAT time series in Cologne and Karlsruhe.Changes in local SAT are also expected to be temporally and spatially highly heterogeneous due to the variability of local climate and the complexity of atmospheric circulation systems (Hansen et al., 2010).The observed CRSs in shallow GWT lag behind the abrupt increase in local SAT by 1-4 years (Table 4).In Karlsruhe the time lag is generally small with 1 year for all shift events, while the time lags in Cologne vary between 2 and 4 years.This difference in the time lags reflects the specific hydrogeological site conditions with the unsaturated zone in Cologne (17 m) being significantly thicker than in Karlsruhe (7 m, Table 2).The thermal properties in the unsaturated zone differ significantly from those in the saturated zone (Table 3).Thus the propagation of the thermal signal in Cologne is retarded due to the lower thermal diffusivity than in Karlsruhe.
The magnitudes of the temperature increase between two subsequent regimes in the zonal mean SAT change are considerably higher than in the global mean SAT change (Fig. 3), because the global temperature data set contains ocean temperature measurements, and ocean temperatures are known to respond more slowly to climatic forcing due to the ocean's large thermal inertia (Hansen et al., 2010).The abovementioned temporal and spatial heterogeneity of the CRS also accounts for the higher increase in SAT in the German time series, which is above the average of the zonal mean in 90-24 • N. The significant abrupt increase in the long-term mean of SAT with the late 1980s CRS of close to 1 • C was likewise observed in Swiss SAT by Figura et al. (2011).
The magnitudes of the increases in the long-term means of GWT are lower and damped by up to 70 % compared to the shift magnitude in SAT (Fig. 3).This damping arises from the fact that, due to the thermal inertia of the subsurface, the GWT has not yet fully equilibrated with the GST at the time when the regime shift is observed in the GWT.The magnitudes of the regime shifts in Fig. 3 also reveal that the damping in the time series from the Hardtwald wells is more pronounced than the damping in Dansweiler and Sinthern.This likely occurs due to depth of groundwater extracted for temperature measurements.For example, the depth to the midpoint of the well screen is higher for the Hardtwald wells than it is for the Dansweiler and Sinthern wells (Table 2).This will be investigated in more detail with the TCA model.

Stationarity within the regimes
In order to investigate the stationarity within the identified regimes the Mann-Kendall test for the absence of trend was performed for the individual regimes.The resulting p values are listed in Table 5, in which high p values close to 1 indicate stationary conditions.No significant trends could be found within the individual regimes of the examined SAT time series, suggesting that the temperature increase in the last decades can be attributed completely to the detected CRS.
In the GWT time series, the p values of the Mann-Kendall test are generally lower (median of 0.20, Table 5) than the p values of the SAT time series (median of 0.53), indicating that the SAT time series are more stationary than GWT time series.This more gradual increase in GWT reflects the effects of subsurface heat transport dynamics, which convert the sharp surface temperature signal to a more diffuse subsurface temperature signal.A significant trend (p < 0.05) with a slope of 0.13 • C was detected in the third regime (2001)(2002)(2003)(2004)(2005)(2006) in the Sinthern well.However, it has to be noted that this regime is quite short, and thus the trend analysis may be biased by the last two rather high temperature values in 2005 and 2006.In the regimes before 1991, the p values of the time series in Dansweiler and Sinthern are 0.05 and 0.06, respectively, and thus close to the critical p value of 0.05 indicating a more gradual increase rather than abrupt changes.For the wells near Karlsruhe no significant trends were found in GWT within the regimes, which indicates that the temperature increase in the time series can be linked to the regime shifts.   2 and 3.The GWT data at the lower and higher ends of the temperature envelope are obtained with the ranges in well screen depth, thermal diffusivity, heat capacity, and recharge rates (Tables 2 and 3).
To compare the performance of the regime shift analysis to an approach with linear temperature increase, the RMSE values for the statistical step function model and a linear model were calculated for each time series (not shown).This analysis revealed that the RMSE of the step function fit for all GWT and SAT time series is slightly lower than the RMSE of the linear fit, indicating that the step function model performs slightly better.Thus, it can be stated that, with the exception of the potentially biased last regime in Sinthern, all regimes in the time series of GWT are statistically stationary, which corroborates the feasibility of the application of regime shift analyses on GWT time series in addition to the low p values of STARS (Table 4).

Analytical model
Predicted GWTs were obtained from the analytical solution in Eq. ( 6) (TCA model) with the thermal properties and recharge rates given in Table 3 and the magnitude and timing of the regime shifts given in Table 4. Due to the availability of GWT data in each well, model runs were started in 1970.Figure 4 shows the measured GWT, assigned GST boundary condition, and predicted GWT for each of the four wells.The range of predicted GWT (shaded area, Fig. 4) is derived from the range of thermal properties, well screen depths, and recharge values utilized as input parameters to the model (Table 3).In particular, the upper boundaries of the temperature envelopes in Fig. 4 were obtained with Eq. ( 6) using depths to the tops of each well screen (Table 2), maximum thermal diffusivities (Table 3), and minimum heat capacities (Table 3; see Eq. 3).The lower boundaries of the temperature envelopes were obtained using depths to the bottom of the well screens, minimum thermal diffusivities, and maximum heat capacities.Finally, the best estimates (red series, Fig. 4) for the predicted GWT data for each well were obtained using depths equal to the midpoints of the well screens and mean thermal diffusivity and heat capacity.In all cases, the thermal properties were taken as the weighted arithmetic average of the unsaturated and saturated zone thermal properties as described in Sect.3.2.
Note that the GST data simulated for the Hardtwald wells are characterized by a wider range in the predicted temperature envelopes.This range is primarily due to the longer well screens in the case of the Hardtwald wells than for the Sinthern and Dansweiler wells (Table 2).Hereafter, when we refer to the TCA model results we ignore the range in the modeling results and only allude to the specific results obtained using the midpoint of the well screens and mean thermal properties given in Table 3 (i.e., red series, Fig. 4).
The TCA model-predicted trends in GWT generally concur with the long-term trends exhibited in the measured data for Dansweiler, Hardtwald 1, and Hardtwald 2; however, the TCA model underpredicts the rise in the Sinthern GWT data.These differences suggest that, although they were assumed to be equal, the magnitude of the GST regime shifts in Sinthern may have been greater than that in Dansweiler, or thatdue to subsurface heterogeneity -the pumped water may be predominantly sourced from above the Sinthern well screen.Modeling results (not shown) indicate that, for the mean thermal properties and recharge values, the z value used in Eq. ( 6) would have to be approximately 9 m for the predicted and observed GWT trends to generally concur.Furthermore, the recharge rates in this well may have been greater than the obtained regional recharge rates for this area.Higher recharge would lead to higher heat advection, which would reduce the lag between a GST signal and its realization in the subsurface (see range in predicted Sinthern GWT, Fig. 4).Similarly, higher thermal diffusivity would generally lead to higher GWT in Sinthern, as the Sinthern GWT is still adjusting to the GST regime shifts in the data shown in Fig. 4. Finally, the last few years of measured GWT data are not available for the Sinthern well.GWT data in the nearby Dansweiler well decreased during this period; thus the visual fit between the measured and predicted Sinthern GWT would likely improve if these data were available.
Our approach does not reproduce inter-annual variability in GWT due to the nature of the GST boundary condition, which is constant for a given climate regime (Fig. 4).Interannual variability in GWT could theoretically be reproduced by considering a series of "GST regimes" that only last 1 year; however, the objective of the present study was to examine the subsurface thermal influence of climate regime shifts -not inter-annual SAT or GST variability.Finally, it is interesting to note that the abrupt regime shifts applied in the simplified boundary condition manifest themselves as gradual changes in the predicted GWT evolution in the deeper wells due to the influence of the heat capacity and thermal inertia of the subsurface.These findings demonstrate that observed gradual increases in shallow GWT are not necessarily suggestive of gradual trends in GST.The effects of the abrupt GST regime shifts are discernible in the upper edge of the temperature envelopes in Hardtwald 1 and 2 (i.e., the GWT signal is diffused, but the impact of the piecewise boundary condition is still discernible).This is due to the fact that these particular results were obtained for depths to the top of the well screens of only 10 m (Table 2).
With the exception of the anomalous Sinthern data, the general agreement between the predicted and observed trends in GWT data (Fig. 4) indicates that the TCA model can yield first-order approximations of the thermal sensitivity of these shallow aquifers to past or future climate regime shifts by conforming the boundary condition to climate model projections.The boundary condition form employed in this study could be matched to future climate projections by considering a series of short GST regimes, or alternatively a boundary condition that considers a gradual rise in GST could be employed (e.g., Kurylyk and MacQuarrie, 2014).The form of the analytical solution indicates that if a new long-term stable climate is achieved, the GWT will eventually rise an equivalent magnitude to the changes in GST, which are often in turn assumed to follow changes in SAT.In the absence of snowpack evolution or land cover changes, any perceived damping in GWT changes in comparison to SAT changes based on statistical analyses likely results from the lagged subsurface thermal response to the boundary condition.
There are limitations associated with employing analytical solutions to simple one-dimensional heat transport equations.Several assumptions associated with the conductionadvection equation have been previously noted.For example, the governing equation, and hence the analytical solution, assumes that the water phase is constant.This assumption is justified in the present study considering that no permafrost thaw (which retards soil warming; Kurylyk et al., 2014b) is occurring.Also, the solution assumes homogeneous thermal properties; however, we considered heat transport in both the saturated and unsaturated zones.The thermal diffusivities of the unsaturated zone for the wells considered in the present study were up to 30 % lower than the saturated zone thermal diffusivities (Table 3).We considered both zones by employing a weighted arithmetic average (based on zone depths) for the effective thermal diffusivity.Also, recharging water may be at a temperature different than the mean annual surface temperature, particularly if the recharge mechanism is snowmelt.However, snowmelt-induced recharge is minimal at the observation areas in this study.In general, due to these limitations, the results presented in Fig. 4 should be considered first-order approximations of the sensitivity of these shallow aquifer thermal regimes to climate regime shifts.

Implications for future river temperatures and groundwater-dependent ecosystems
Although the wells analyzed in this study were not located near streams, the timing and magnitude of the measured GWT rise can provide insight into the potential warming of alluvial aquifers feeding ecologically important rivers.Gaining rivers and streams can be strongly influenced by the thermal regimes of surrounding aquifers (e.g., Tague et al., 2007;Kelleher et al., 2012), and this is often particularly true during the dry, warm season when base flow can provide the majority of the river or stream discharge.Thus, deterministic models of future water temperature of base-flow-dominated streams and rivers should explicitly account for the future thermal regimes of aquifers.Various studies have demonstrated that the thermal regimes of rivers respond to a warming climate, and these studies have generally tacitly ignored GWT rise due to climate change.The results of this study however contradict this assumption by indicating that shallow GWT will respond to SAT warming and that the lag time between SAT warming and the associated increase in shallow GWT can be rather short (< 5 years).Similar results were obtained by Kurylyk et al. (2014a), who employed a numerical model of groundwater flow and energy transport driven by downscaled climate scenarios to demonstrate a potential damping and short lagging of future groundwater discharge temperature rise in response to air temperature changes.
Given the expected warming of rivers across the globe (van Vliet et al., 2011(van Vliet et al., , 2013)), researchers have rightfully proposed that cold-water fish will begin to increasingly rely on the occurrence and distribution of suitable cold-water refugia (e.g., Brewer, 2013).Our results suggest that GDEs and groundwater-sourced cold-water refugia will also warm in response to climate change.The magnitude and timing of the GWT warming however will depend on several factors, including the timing and magnitude of the SAT warming, changes in precipitation (and thus recharge and advection), the depth of the overlying soil and and the presence or absence of seasonal snowpack.

Conclusions
By applying a sequential t test analysis for regime shifts (STARS) to time series of air and groundwater temperatures, we empirically demonstrated that groundwater temperatures in shallow aquifers show temperature changes that correspond to positive shifts in local SAT in Germany, which in turn can be traced back to increasing global SAT.This observed direct coupling of atmospheric and groundwater temperature development through the unsaturated zone implies that climate warming affects not only aquifers recharged by riverbank infiltration (Figura et al., 2011) but also a large number of shallow aquifers on a wide spatial scale.The regime shifts in GWT occur with a certain time lag to the CRS depending mainly on the thermal properties and thickness of the overlying soil.The magnitude of these regime shifts in GWT compared to the shifts in SAT is damped by the thermal propagation of the temperature signal into the subsurface, leading to a more gradual increase in GWT.This damping perceived in the statistical analyses is predominantly an artifact of the lagged subsurface thermal response.However, despite the extenuation of the temperature signal in the subsurface and the mixing of shallow groundwater during pumping, significant temperature shifts were found in the extracted groundwater.
Process-oriented modeling was also performed with an analytical solution to the conduction-advection equation.In three of the four observation wells, the simulated decadal GWT trends generally concurred with the measured decadal GWT trends.However, inter-annual variability could not be reproduced due to the simplistic nature of the boundary condition, which equals the long-term mean surface temperature.This agreement indicates that the solution to the conductionadvection equation can also be applied to obtain first-order estimates of the influence of future climate change on subsurface thermal regimes.
Our results indicate that increasing SATs are prone to have a substantial and swift impact, not only on soil temperatures but also on large-scale, shallow groundwater temperatures in productive and economically important aquifers.Furthermore, this study has demonstrated that long-term series of pumped groundwater temperature can be analyzed using stochastic approaches to examine the relationship between local and global climate change and local groundwater temperature evolution.

Figure 1 .
Figure 1.(a, b): locations of the four observation wells and two weather stations used in the present study.(c, d): conceptual sketch of the well settings in the aquifers close to Cologne (left) and Karlsruhe (right).The black zones in the wells indicate the location of the filter screens.Please note the different scales in the subsurface.

Figure 2 .
Figure 2. Time series of temperature data with long-term means (dashed lines) and observed regime shifts with calculated RSI values (grey bars).

Figure 3 .
Figure3.Left: magnitude of regime shifts in time series of atmospheric and groundwater temperatures.Right: relative damping of the regime shift magnitude in groundwater temperatures compared to regional atmospheric regime shift, calculated as 100 minus the ratio of T in GWT to T in SAT in percent.

Figure 4 .
Figure 4. Measured GWT, predicted GWT, and assigned GST boundary conditions for the TCA model (Eq.6) for each well versus the year.Red lines indicate GWT results obtained using the mean well screen depth, thermal properties and recharge rates presented in Tables2 and 3.The GWT data at the lower and higher ends of the temperature envelope are obtained with the ranges in well screen depth, thermal diffusivity, heat capacity, and recharge rates (Tables2 and 3).

Table 2 .
Hydrogeological data of the four observation wells.

Table 3 .
Range of thermal properties and recharge values utilized in the analytical solution.

Table 4 .
Time lags and final p values of the observed regime shifts in air and groundwater temperature.The specific years indicate the first year of the new regime.Time lags are defined as the period between the occurrence of a regime shift in local SAT and the corresponding successive shift in GWT.

Table 5 .
Results (p values) of the Mann-Kendall test for the absence of a trend for all regimes in SAT and GWT time series.