Shallow Groundwater Thermal Sensitivity to Climate Change and Land Cover Disturbances: Derivation of Analytical Expressions and Implications for Stream Temperature Modeling

Climate change is expected to increase stream temperatures and the projected warming may alter the spatial extent of habitat for cold-water fish and other aquatic taxa. Recent studies have proposed that stream thermal sensitivities , derived from short-term air temperature variations, can be employed to infer future stream warming due to long-term climate change. However, this approach does not consider the potential for streambed heat fluxes to increase due to gradual warming of the shallow subsurface. The temperature of shallow groundwater is particularly important for the thermal regimes of groundwater-dominated streams and rivers. Also, recent studies have investigated how land surface perturbations , such as wildfires or timber harvesting, can influence stream temperatures by changing stream surface heat fluxes, but these studies have typically not considered how these surface disturbances can also alter shallow groundwater temperatures and streambed heat fluxes. In this study, several analytical solutions to the one-dimensional unsteady advection–diffusion equation for sub-surface heat transport are employed to estimate the timing and magnitude of groundwater temperature changes due to seasonal and long-term variability in land surface temperatures. Groundwater thermal sensitivity formulae are proposed that accommodate different surface warming scenarios. The thermal sensitivity formulae suggest that shallow groundwater will warm in response to climate change and other surface perturbations, but the timing and magnitude of the subsurface warming depends on the rate of surface warming, subsurface thermal properties, bulk aquifer depth, and groundwater velocity. The results also emphasize the difference between the thermal sensitivity of shallow groundwater to short-term (e.g., seasonal) and long-term (e.g., multi-decadal) land surface-temperature variability , and thus demonstrate the limitations of using short-term air and water temperature records to project future stream warming. Suggestions are provided for implementing these formulae in stream temperature models to accommodate groundwater warming.


Introduction
The water temperature of streams and rivers is an important determinant of aquatic ecosystem health due to its influence on physicochemical conditions and because many freshwater fish species can only tolerate a certain temperature range (Caissie, 2006;Elliott and Elliott, 2010;Hannah and Garner, 2015;Webb et al., 2008).Also, river thermal diversity enhances ecosystem complexity by providing thermally suitable habitat in reaches that would otherwise be uninhabitable for certain species (Cunjak et al., 2013;Ebersole et al., 2003;Kurylyk et al., 2015;Sutton et al., 2007).The thermal regimes of streams and rivers are controlled by energy fluxes across the water surface and the streambed (Fig. 1 well as the internal structure of the stream or river network (Guenther et al., 2015;Hannah et al., 2004;Herbert et al., 2011;Leach and Moore, 2011;Poole and Berman, 2001).The total streambed heat flux is composed of conductive and advective heat fluxes, which both depend on subsurface temperatures (Caissie et al., 2014;Moore et al., 2005;St-Hilaire et al., 2000).Large rivers tend to be dominated by surface heat fluxes, but streambed advective heat fluxes induced by groundwater-surface water interactions can influence the thermal regimes of certain headwater streams or smaller rivers (Caissie, 2006).The significance of streambed advective heat fluxes generally varies spatially and temporally within a channel and depends on, among other things, the groundwater discharge rate and the degree of shading (e.g., Brown and Hannah, 2008;Leach and Moore, 2011;Story et al., 2003).Due to the thermal inertia of the subsurface soil-water matrix, groundwater-dominated streams and rivers typically exhibit attenuated thermal responses to diel and seasonal variations in air temperature compared to surface runoff-dominated streams and rivers (Caissie et al., 2014;Constantz, 1998;Garner et al., 2014;O'Driscoll and DeWalle, 2006;Tague et al., 2007).Kelleher et al. (2012) defined the thermal sensitivity of a stream as the slope of the linear regression between air and water temperatures.These regressions are typically performed on temperature data collected for a period of at least 1 year and averaged on a daily, weekly, or monthly basis.The stream thermal sensitivity is thus a measure of the short-term (e.g., seasonal) change in water temperature in response to a short-term change in air temperature (Kelleher et al., 2012;Mayer, 2012).
Because groundwater temperature exhibits less seasonal variability than surface water temperature, it is not surprising that extrapolated stream thermal sensitivities obtained from short-term temperature data will typically indicate that the temperature of groundwater-dominated streams will be relatively insensitive to climate change.As noted by Johnson (2003), care should be taken when using air temperature correlations to explain stream temperature dynamics, as air temperature is not the dominant controlling factor in stream temperature dynamics.Rather, the high correlation between stream and air temperature arises because both variables are influenced by incoming solar radiation, the primary driver of stream temperatures (Allan and Castillo, 2007).The approach of using short-term stream thermal sensitivities to estimate multi-decadal stream warming essentially employs future air temperature as a surrogate for future stream surface heat fluxes (Gu et al., 2015;Johnson et al., 2014;Mohseni and Stefan, 1999), but it ignores changes to streambed heat fluxes due to groundwater warming.Thus, the short-term relationship between air and water temperatures is not necessarily representative of the concomitant warming of the lower atmosphere and surface water bodies on inter-annual or multi-decadal timescales (Arismendi et al., 2014;Bal et al., 2014;Luce et al., 2014).
Furthermore, many studies have investigated the response of stream thermal regimes to land surface perturbations, such as wildfires and deforestation, for the first few years following the disturbance.However, very few studies have considered how these perturbations could increase the temperature of groundwater discharge to these streams and thereby produce enhanced or sustained stream warming.In general, the common approach of ignoring future increases in groundwater temperature, and streambed heat fluxes in stream temperature models may underestimate future stream warming and associated environmental impacts (e.g., habitat loss for coldwater fish; Snyder et al., 2015).
There is increasing evidence that the thermal regimes of shallow aquifers are sensitive to climate change, permanent deforestation, and wildfires.Observed shallow groundwater temperature warming has already been related to recent trends in air temperature (an indicator of climate change) in Taiwan (Chen et al., 2011), Switzerland (Figura et al., 2011(Figura et al., , 2014)), and Germany (Menberg et al., 2014).Empirical and process-based models of energy transport in shallow aquifers have been used to suggest that future climate change will continue to warm shallow groundwater bodies (e.g., Gunawardhana and Kazama, 2011;Kurylyk et al., 2013Kurylyk et al., , 2014a;;Taylor and Stefan, 2009) as reviewed in detail by Kurylyk et al. (2014b).Previous studies have also noted groundwater warming in response to deforestation due to the removal of the forest canopy (e.g., Alexander, 2006;Guenther et al., 2014;Henriksen and Kirkhusmo, 2000;Steeves, 2004;Taniguchi et al., 1998).Others have observed subsurface warming following wildfires.Burn (1998) found that the mean annual surface temperature at a burned site in southern Yukon, Canada, was 0.6 • C warmer than the surrounding surface thermal regime, and this surface thermal perturbation rapidly increased shallow subsurface temperatures.
In all cases (i.e., climate change, deforestation, and wildfires), the surface disturbance warms shallow aquifers by increasing the downward heat flux from the warming land surface.For example, climate change can influence surface thermal regimes and subsurface heat fluxes by altering convective energy fluxes from the lower atmosphere and causing increased net radiation at the ground surface (Jungqvist et al., 2014;Kurylyk et al., 2013;Mellander et al., 2007).The influence of wildfires or forest harvesting on surface thermal regimes can be complex.The removal of the forest canopy can decrease transpiration and thus increase the energy available to warm the land surface (Rouse, 1976).Lewis and Wang (1998) demonstrated that the majority of surface and subsurface warming caused by wildfires at sites in British Columbia and Yukon, Canada, could be attributed to decreased transpiration.Decreased surface albedo and consequent increased net radiation at the land surface can also arise due to wildfires (Yoshikawa et al., 2003).The increase in surface temperature as a result of a land cover disturbance will depend on the original vegetative state, climate, ground ice conditions, and potential for vegetative regrowth (Liljedahl et al., 2007).In the case of a wildfire or in post-harvest tree planting, the vegetation may eventually regenerate, and the surface-energy balance and temperature may return to the pre-fire conditions (Burn, 1998).Kurylyk et al. (2013Kurylyk et al. ( , 2014a) ) demonstrated that shallow groundwater warming may eventually exceed the magnitude of surface water warming and thus stream temperature models that do not consider this phenomenon may be overly conservative.The empirical method proposed by Kurylyk et al. (2013) for estimating the magnitude of groundwater warming requires measured land surface temperature and depth-dependent groundwater temperature for model calibration, but there is often a paucity of such temperature data available at the catchment scale.Also, the numerical modeling described by Kurylyk et al. (2014a) is time intensive and requires considerable data for model parameterization.These previous approaches for quantifying groundwater warming are site specific, and thus the results are not generally transferable to existing models that are used to investigate stream thermal regimes.
The intent of this study is to provide alternative, parsimonious approaches for investigating factors that influence the timing and magnitude of groundwater temperature changes in response to climate change or land cover disturbances.The specific objectives of this paper are twofold: 1. derive easy-to-use formulae to estimate the thermal sensitivity of groundwater to different surface-temperature changes (e.g., seasonal cycle or multi-decadal increases); 2. demonstrate how these formulae can be utilized to estimate how the groundwater thermal sensitivity in idealized environments is influenced by the depth, groundwater recharge rate, and subsurface thermal properties.
The illustrative examples (objective 2) will also be used to demonstrate the difference in the subsurface thermal response to short-term (seasonal) and long-term (multidecadal) surface-temperature trends.Consequently, the results will be employed to highlight the limitations of employing empirical stream temperature models with constant coefficients obtained from short-term temperature records to project future stream warming.The results will also demonstrate how stream temperature models can be improved to accommodate groundwater warming using these simple approaches.

Methods
There are several approaches for estimating future groundwater temperature warming in response to changes in land cover or climate.It is well known that mean annual ground surface temperature and shallow groundwater temperature are approximately equal to mean annual air temperature plus some thermal offset (e.g., 1-4 • C) due to the insulating effect of snow (Zhang, 2005).Meisner et al. (1988) employed this knowledge to estimate future groundwater temperatures by adding a thermal offset to projections of future mean annual air temperature.The approach employed by Meisner et al. (1988) utilized mean annual surface temperature as a proxy for groundwater temperature and thus implicitly assumed that the aquifer and ground surface are always in thermal equilibrium.The equilibrium assumption was also invoked in the empirical function described by Kurylyk et al. (2013).Such an approach does not consider the lag that occurs between an increase in surface temperature and its subsequent realization at some depth within the subsurface (Lesperance et al., 2010) and thus is only valid for very shallow groundwater (e.g., < 5 m) or for long timescales.Analytical solutions to subsurface heat transfer differential equations can also be applied to estimate the influence of future climate change on groundwater temperature (Gunawardhana and Kazama, 2011;Kurylyk and MacQuarrie, 2014;Menberg et al., 2014), although these approaches have most often been applied for deeper aquifers.Finally, numerical models of groundwater flow and coupled heat transport can be employed to investigate the thermal evolution of aquifers due to warming surface temperatures (e.g., Gunawardhana and Kazama, 2012;Kurylyk et al., 2014a).These numerical models are more flexible and can accommodate multidimensional groundwater flow and heat transport and inhomogeneities in subsurface thermal properties, but they require extensive subsurface field data for model parameterization.
Herein, we employ analytical solutions to a onedimensional, unsteady heat transport equation to estimate subsurface-temperature evolution due to climate change, permanent land cover changes, and wildfires.These solutions are physically based and account for the lag in the thermal response of groundwater to surface-temperature changes.Also, unlike the solution employed by Taylor and Stefan (2009), these solutions accommodate the subsurface thermal effects of vertically moving groundwater.The solutions provide an indication of expected groundwater warming due to climate or land cover changes, and the results can be incorporated into stream temperature models in the absence of site-specific hydrogeological modeling.These simple analytical solutions are particularly useful for performing parsimonious analyses when there is a paucity of subsurface data (e.g., hydraulic conductivity distribution) for parameterizing groundwater flow and energy transport models.Also, analytical solutions limit the degrees of freedom for a particular analysis and thus facilitate a comprehensive evaluation of possible interactions between model inputs and resultant solutions.As we demonstrate, the forms of these solutions can also be utilized to derive mathematical expressions for groundwater thermal sensitivity to surface-temperature perturbations.The analytical solutions discussed in this paper invoke assump-tions, and the limitations arising from these assumptions will be discussed later.

Advection-diffusion heat transport equation
Shallow subsurface heat transfer occurs primarily due to heat conduction and heat advection (Domenico and Schwartz, 1990), although the latent heat released or absorbed during pore water freeze-thaw can also be important in cold regions (Kurylyk et al., 2014b).The one-dimensional, transient conduction-advection equation for subsurface heat transport is (Stallman, 1963) where λ is the bulk thermal conductivity of the soil-water matrix (W m −1 • C −1 ), T is the temperature at any point in space or time ( • C), z is the depth below the surface (m; down is positive and the land surface occurs at z = 0), q is the vertical Darcy flux (m s −1 ; down is positive), c w ρ w is the volumetric heat capacity of pure water (4.18 × 10 6 J m −3 • C −1 ; Bonan, 2008), t is time (s), and cρ is the bulk volumetric heat capacity of the soil-water matrix (J m −3 • C −1 ).The first term on the left of Eq. ( 1) represents the divergence of the conductive flux, the second term on the left represents the divergence of the advective flux, and the term on the right represents the rate of change of thermal storage.Subsurface heat transport phenomena and the physical meaning of the terms in Eq. ( 1) are reviewed in more detail by Rau et al. (2014) and Kurylyk et al. (2014b).Equation ( 1) is often rewritten in the form (Carslaw and Jaeger, 1959) where D is the bulk thermal diffusivity (thermal conductivity divided by heat capacity) of the soil-water matrix (m 2 s −1 ), and U is the velocity of a thermal plume due only to heat advection (m s −1 ).Even in the absence of conduction, the thermal plume will not migrate at the same rate as the Darcy velocity due to differences in the heat capacities of water and the medium (Markle and Schincariol, 2007;Luce et al., 2013).An expression for U can be obtained via a comparison of Eqs. ( 1) and ( 2) Often an effective thermal diffusivity term, which accounts for the combined thermal homogenizing effects of heat diffusion and heat dispersion, is utilized in place of the bulk thermal diffusivity term D in Eq. ( 2).However, it is still common to ignore the subsurface thermal effects of dispersion, which are often minimal in comparison to heat conduction (Kurylyk et al., 2014b;Rau et al., 2014).Equation (2) represents vertical subsurface heat transport processes and accounts for the thermal effects of heat conduction induced by a thermal gradient and heat advection induced by groundwater flow.Analytical solutions to this equation can be developed and applied to consider inter-relationships between groundwater flow, surface-temperature changes, and subsurface thermal regimes.We consider four analytical solutions to Eq. (2) (Table 1) that vary based on the nature of the surface boundary condition.These are discussed in subsequent sections.

Analytical solution 1: harmonic surface-temperature changes
The diel or seasonal land surface-temperature cycle can be approximated with a harmonic function.Suzuki (1960) derived an analytical solution to Eq. ( 2) subject to a sinusoidal surface-temperature boundary condition: where A is the amplitude of the harmonic surfacetemperature cycle ( • C), T m is the mean surface temperature ( • C), p is the period of the surface-temperature cycle (s), ϕ is a phase shift to align the timing of the surface-temperature signal with the sinusoid (rad), d is a thermal damping term (m −1 ), and L is a lag term (m −1 ).Equation ( 5) thus states that the harmonic temperature signal at the surface retains its period within the subsurface but is exponentially damped and linearly lagged with depth.Stallman (1965) demonstrated that the exact expressions for d and L are Equations ( 5) to ( 7) are generally collectively referred to as Stallman's equation.No initial conditions are presented for the solution of Stallman (1965) as it assumes that the boundary condition has been repeating the harmonic cycle indefi-nitely.This solution also depends on a lower boundary condition (T = T m ) at infinite depth.Various forms of this solution have been applied/inverted to infer rates of groundwater flow due to subsurface temperature-time series arising from daily or seasonal harmonic variations in surface temperature (e.g., Anderson, 2005;Hatch et al., 2006;Rau et al., 2014).
Here, we employ the solution of Stallman (1965) in a forward manner to demonstrate why seasonal changes in air and surface temperature are not manifested in subsurface thermal regimes below certain depths, and thus why groundwaterdominated streams and rivers exhibit low thermal sensitivity to seasonal weather variability.In particular, we consider the ratio of the amplitude of the seasonal groundwater temperature cycle at any arbitrary depth to the amplitude of the surface-temperature boundary condition.This dimensionless parameter, herein referred to as the exponential damping factor , can be obtained from Eqs. ( 4) and ( 5) 2.3 Analytical solution 2: step change(s) in surface temperature due to land cover disturbances Taniguchi et al. (1999a) demonstrated how an analytical solution presented by Carslaw and Jaeger (1959) could be modified to calculate the groundwater warming arising from a sudden and permanent increase in surface temperature.This increase in surface temperature could arise due to rapid and large-scale timber harvesting or other changes in land use.Menberg et al. (2014) proposed that superposition principles could be employed to modify the solution by Taniguchi et al. (1999a) by considering a series of shifts in the surfacetemperature boundary condition.Herein we employ the technique by Menberg et al. (2014) and consider up to two sequential shifts in the boundary condition.The first shift, which warms the surface temperature, occurs at t = 0, and after a period of time (t = t 1 ), the surface temperature returns to its value prior to the initial warming (T 0 ).Such a boundary condition could approximate the sudden temporary increase in mean annual surface temperature due to a wildfire and the subsequent return to pre-fire surface temperatures due to vegetation regrowth (Burn, 1998).Alternatively, this boundary condition could represent the effect of clearcutting followed by industrial tree planting.The subsequent surface cooling due to gradual vegetative regrowth could also be represented with a series of shorter less intense cooling phases, but for the illustrative examples in the present study we assume one warming shift followed by one cooling shift of equal magnitude: Boundary condition : where T 0 is the uniform initial temperature ( • C), T is the magnitude of the surface-temperature shift ( • C), erfc is the complementary error function, and t 1 is the duration of the period characterized by warmer surface temperatures (s).This solution and the remaining three solutions presented later also require a lower boundary condition at infinite depth (T = T 0 ).Equation ( 11) can be employed to consider the subsurface warming due to a permanent step change in surface temperature (i.e., no subsequent cooling due to vegetative regrowth) by setting t 1 to infinity.In this case, only the first line on the right-hand side of Eq. ( 11) is retained.Even when t 1 is set to infinity, Eq. ( 11) differs slightly from the solution presented by Taniguchi et al. (1999a) because uniform initial temperatures are assumed in the present study (Eq.9).These initial conditions ignore the influence of the geothermal gradient and imply that the recent climate has been relatively stable.We employ these simplifying assumptions given that we are primarily interested in shallower depths (e.g., < 25 m) where the influence of the geothermal gradient is not as important.Also, the boundary conditions for this solution and the following solutions do not include seasonal temperature variations.Thus, these solutions are valid for predicting the evolution of mean annual groundwater temperature.

Analytical solution 3: linear increase in surface
temperature due to climate change Carslaw and Jaeger (1959) also presented an analytical solution to Eq. ( 2) subject to linearly increasing surface temperature.This solution was later adapted by Taniguchi et al. (1999b) and applied to study groundwater temperature evolution due to climate change.Herein, the analytical solution is presented in a slightly simpler form as thermally uniform initial conditions are assumed (i.e., initial conditions are given by Eq. 9): where β is the rate of the increase in surface temperature ( • C s −1 ).Equation ( 13) has been applied in an inverse manner to consider the complex relationships between past surfacetemperature changes, groundwater flow, and measured subsurface temperature-depth profiles (e.g., Miyakoshi et al., 2003;Taniguchi et al., 1999b;Uchida and Hayashi, 2005).It has also been applied to forward model future groundwater temperature evolution due to projected climate change (Gunawardhana and Kazama, 2011).Herein, the surface boundary condition (Eq.12) is fitted to mean annual air temperature trends produced by climate models.Because it is surface temperature, rather than air temperature, that drives shallow subsurface thermal regimes, this approach tacitly assumes that mean annual surface and air temperature trends are coupled.Thus, air temperature is being used as a proxy for surface temperature in this approach.Snowpack evolution may invalidate this assumption (Mellander et al., 2007), and thus it is best employed where snowpack effects are minimal.Snowpack evolution would typically retard the rate of groundwater warming (Kurylyk et al., 2013).

Analytical solution 4: exponential increase in surface temperature due to climate change
It may be inappropriate to assume a linear surface temperature rise as in Eq. ( 13), because many climate scenarios suggest that the rate of climate warming will increase over time.
Figure 2 presents the globally averaged IPCC (2007) multimodel air temperature projections for two different emission scenarios.The global air temperature series projected for the conservative emission scenario B1 is much better represented by a linear function than the air temperature series for the aggressive A2 emission scenario, which exhibits significant concavity.
In such cases, the boundary condition would be better represented as an exponential function (Kurylyk and MacQuarrie, 2014).The solution presented here is simpler than the original form given that the initial conditions are assumed to be thermally uniform (initial conditions = Eq.9):  (data from, IPCC, 2007).Details concerning the exponential and linear fits to the IPCC projections are given in Sect.3.3.1.Modified from Kurylyk and MacQuarrie (2014).
where T 1 ( • C), b ( • C), and c (s −1 ) are parameters for the surface-temperature boundary condition which can be fitted to climate model projections.Note that T 1 + b must equal T 0 for the boundary and initial conditions to converge at t = 0, z = 0.The original initial condition function proposed by Kurylyk and MacQuarrie (2014) superimposed linear and exponential functions, and thus the more complex form of the solution can also be applied to forward model future climate change impacts on deeper subsurface-temperature profiles.These temperature profiles can deviate from the geothermal gradient due to groundwater flow or recent surfacetemperature changes (Ferguson and Woodbury, 2005;Reiter, 2005).The alternate forms of the boundary conditions presented in Eqs. ( 10), (12), and ( 14) are illustrated in Fig. 3.Each of the listed analytical solutions to the one-dimensional, transient advection-diffusion equation is provided in Table 1 with details to highlight their differences.

Effective aquifer depth
The analytical solutions discussed above can be utilized to estimate the influence of surface warming at any desired (c) Linearly increasing surface temperature (Eq.12) (d) Exponentially increasing surface temperature (Eq.14) Climate change (conservative)  10).The difference between these is the duration of the period of warm surface temperatures (t 1 = ∞ in a).(c-d) The boundary conditions for GST due to long-term climate change for conservative (linear; Eq. 12) and aggressive (exponential; Eq. 14) climate scenarios.
depth.However, groundwater discharge to streams is sourced from different depths within the aquifer depending on the recharge location and the subsurface flow paths (Fig. 4a).Because the water table slope in unconfined aquifers is typically subdued in comparison to the land surface slope (Domenico and Schwartz, 1990), soil water that recharges the aquifer further upslope typically has a longer residence time and reaches greater depths relative to the land surface than soil water recharging the aquifer close to the discharge point.Groundwater flow in aquifers is often conceptualized as occurring in different "flow channels" or "flow tubes" (Domenico and Schwartz, 1990), and groundwater discharge is a thermal and hydraulic mixture of different groundwater flow channels coming from different depths and converging at the discharge point (Hoehn and Cirpka, 2006 and Fig. 4).Thus, when employing one-dimensional solutions to investigate the thermal evolution of groundwater discharge to streams and rivers, an effective depth z eff (m) must be considered that represents the bulk aquifer depth (i.e., accounting for all discharging groundwater flow channels) as a single point within the subsurface (Fig. 4).As a first estimate, this depth may be taken as the average unsaturated zone thickness.Figure 4b shows the conceptual model employed in this study.Above the effective depth, heat transport and water flow is assumed to be predominantly vertical as is often the case within the unsaturated zone, in overlying aquitards, or even in the upper portion of the aquifer (e.g., Kurylyk et al., 2014b).Within the aquifer (located at the effective depth), groundwater discharges horizontally towards a stream, and horizontal conductive heat transport is assumed to be negligible due to the relatively low horizontal thermal gradients in this zone.Heat advection and associated thermal dispersion near the discharge point is assumed to dominate vertical heat transfer and thus create a thermally uniform zone.Thus, the aquifer is treated as a thin, horizontally well-mixed thermal reservoir discharging to a surface water body (Fig. 4b).This approach is somewhat analogous to how contaminant hydrogeology studies have considered aquifers to be well-mixed reservoirs with respect to solute concentrations (e.g., Gelhar and Wilson, 1974).Vertical heat transfer continues below the aquifer (Fig. 4b).Limitations of this approach are discussed later.

Groundwater thermal sensitivity to long-term surface-temperature perturbations
Groundwater thermal sensitivity is herein defined as the change in groundwater temperature at some depth and time divided by the driving change in surface (z = 0) temperature at the same time.For example, if the surface temperature increases by 2 • C and the groundwater temperature has only increased by 1.4 • C at that same time, then the groundwater thermal sensitivity is 0.7 (1.4 • C/2 • C).The temperature changes at the surface and in the aquifer are measured with respect to the initial temperatures at those locations.This definition for groundwater thermal sensitivity S ( • C • C −1 ) can be expressed in the following manner: The groundwater thermal sensitivity is the analog to the stream thermal sensitivity defined by Kelleher et al. (2012), although the temperature changes are measured on a longer timescale for groundwater (e.g., multi-decadal vs. seasonal).Equation ( 16) represents the thermal sensitivity at any arbitrary depth within the aquifer.The bulk (i.e., the entire portion of the aquifer discharging to the stream or river) groundwater thermal sensitivity in Eq. ( 16) can be found by replacing z with z eff .

Groundwater thermal sensitivity to a step increase in surface temperature (land cover disturbance)
The groundwater thermal sensitivity S s (subscript denotes nature of boundary condition) to a step increase in surface temperature occurring at t = 0 followed by subsequent surface cooling at t = t 1 can be found by inserting Eqs. ( 9), (10), and (11) into Eq.( 16) In Eq. ( 17), sensitivities for all times greater than t 1 are calculated with respect to the initial temperature perturbation T .Interestingly, the groundwater thermal sensitivity is not dependent on the magnitude of the step change in surface temperature T or the initial temperature T 0 , provided that the initial temperature is uniform.Equation ( 17) has the same form as the well-known solute transport analytical solution proposed by Ogata and Banks (1961) to calculate normalized solute concentrations.
As in the case of Eq. ( 11), Eq. ( 17) can be simplified to represent the influence of a permanent step increase (i.e., no subsequent cooling) in surface temperature by setting t 1 to infinity and only considering the first line on the right-hand side of the equation.

Groundwater thermal sensitivity to gradual increases in surface temperature (climate change)
Equation ( 16) can also be applied to obtain an expression for the groundwater thermal sensitivity S L ( • C • C −1 ) due to a linear increase in the surface-temperature boundary condition by inserting Eqs. ( 9), ( 12), and (13) into Eq.( 16) and simplifying: Thus, S L is independent of the initial temperature T 0 and the rate of surface warming β.
The groundwater thermal sensitivity S E ( • C • C −1 ) to an exponentially increasing surface temperature can be obtained by inserting Eqs. ( 9), ( 14), and (15) into Eq.( 16).The resultant solution can be further simplified by canceling terms and by remembering that T 0 is the sum of T 1 and b: A spreadsheet is included in the Supplement that facilitates the calculation of the results for each of the analytical so-  lutions and groundwater thermal sensitivity equations.The user may vary input parameters such as depth, thermal properties, groundwater velocity, time, initial temperature, and the surface-temperature boundary conditions.

Subsurface thermal properties
These analytical solutions assume that subsurface thermal properties are homogeneous, but in reality the bulk thermal properties of unconsolidated soils depend on many factors, including the mineral constituents, porosity, total moisture saturation, and the pore water phase (Farouki, 1981;Kurylyk et al., 2014b).Water has a much higher thermal conductivity than air; thus, the saturated zone typically is characterized by a higher bulk thermal conductivity than the unsaturated zone (Oke, 1978).Despite the existence of subsurface thermal property heterogeneities, natural variability in soil thermal properties is orders of magnitude less than the natural variability in hydraulic properties (Domenico and Schwartz, 1990), and thus homogeneous assumptions are better justified for subsurface heat transport than for subsurface water flow.Table 2 lists the bulk thermal properties for unfrozen sand, clay, and peat at three water saturations (volume of soil water/pore volume).These values are used to represent the typical ranges of thermal conductivities experienced in common unconsolidated soils.The bulk thermal diffusivities of these soils do not vary significantly at pore water saturations above 0.  2) were employed, and a recharge Darcy velocity of 0.2 m yr −1 was assumed.The boundary condition parameters T m , A, ϕ, and p were assigned values of 10 • C, 15 • C, −4.355 radians, and 31 536 000 s (1 year), respectively, to represent typical surface-temperature conditions for a forested site in New Brunswick, Canada (e.g., Kurylyk et al., 2013).

Seasonal surface-temperature influences on groundwater temperature
Stallman's equation (Eqs.5-7) can be utilized to investigate how idealized subsurface environments respond to seasonal surface-temperature changes.Figure 5 shows temperaturedepth profiles for each month and temperature-time series for different depths in a soil column driven by a harmonic boundary condition at the surface (Eq.4).The results were obtained from Eq. ( 5) for sandy soil (thermal properties; Table 2) and for a downwards Darcy velocity (i.e., recharge) of 0.2 m yr −1 .This recharge value was chosen as a representative basin groundwater recharge (Döll and Fiedler, 2008;Healy, 2010).Stallman's equation generally matches seasonal groundwater temperature data reasonably well in shallow subsurface environments, except in locations where snowpack can make the surface temperature non-sinusoidal and the subsurface thermal envelope (Fig. 5a) asymmetrical (Lapham, 1989).Regardless, Eq. ( 5) and Fig. 5 both demonstrate that the seasonal subsurface-temperature variability is exponentially attenuated with depth and is barely discernible beyond a certain depth (e.g., 10-14 m).
The exponential damping factor is the ratio of the amplitude of the seasonal temperature cycle at an arbitrary depth z to the amplitude of the seasonal surface-temperature cycle (Eq.8).It is thus a measure of how the subsurface thermal regime responds to seasonal temperature variations, and it can be considered the seasonal counterpart to the groundwa-ter thermal sensitivities derived from the analytical solutions experiencing long-term surface-temperature variability.Figure 6 illustrates that the exponential damping factor (or seasonal thermal sensitivity) for a given depth decreases for the discharge scenario (black series; Fig. 6) in comparison to the recharge scenario (dashed-blue series).In a discharge scenario, the upward advective flux is impeding the downward propagation of the surface-temperature signal, and thus the surface signal is more quickly attenuated.
Figure 6a-c also indicate that the soil thermal properties greatly influence the subsurface thermal response to seasonal temperature variability.In particular, due to the significantly lower thermal diffusivity of partially saturated peat (Table 2), the surface-temperature signal is more quickly damped in the peat soil (Fig. 6c) in comparison to the results obtained for sand (Fig. 6a) and clay (Fig. 6b).However, in each of the nine scenarios presented in Fig. 6, the parameter is less than 0.2 (amplitude reduced by at least 80 %) when the depth is greater than 5 m, which indicates that groundwater discharge does not have to be sourced from a very deep aquifer to decrease the stream thermal sensitivity to seasonal air temperature changes.

Impacts of land cover disturbances on groundwater temperatures
Beyond the depth of seasonal temperature fluctuations (Fig. 5), groundwater temperature will still be influenced by long-term surface-temperature perturbations.For instance, Recharge (q = 0.2 m yr -1 ) No flow (q = 0 m yr -1 ) Discharge (q = 2 m yr -1 ) Figure 6.Exponential damping factor (seasonal temperature sensitivity) (Eq.8) vs. depth for (a) sandy soil, (b) clay soil, and (c) peat soil.The thermal properties were taken from Table 2 assuming a volumetric water saturation of 50 %.Results are presented for Darcy velocities of 0.2 m yr −1 (recharge; downwards flow), 0 (conduction-dominated thermal regime), and −2 m yr −1 (discharge; upwards flow) and a period of 1 year.A higher discharge value was used in comparison to the recharge value given that discharge is typically concentrated over a smaller area than recharge.11) driven with the step boundary condition (Eq.10), with T = 2 • C and t 1 = infinity (solid lines) or 25 years (dashed lines).The subsurface thermal properties were taken from the 50 % saturated sand and peat values in Table 2, and the recharge rate was 20 cm yr −1 .The results shown in (b) were calculated with Eq. ( 17) using the same parameters as (a).
due to a sudden and permanent (t 1 = ∞, Eq. 10) mean annual surface temperature increase of 2 • C.This is approximately the long-term mean annual surface temperature increase observed by Lewis (1998) in response to deforestation.This is at the lower end of the range (1.6 to 5.1 • C) in the mean annual surface temperature increases noted by Taniguchi et al. (1998) following forest removal in Western Australia.The groundwater warming, rather than the temperature, is obtained by setting the initial temperature to zero (T 0 ; Eqs. 10 and 11).
Results are presented for sandy soil and peat soil as these two soils respectively exhibit the highest and lowest thermal diffusivities given in Table 2. Due to the nature of the surface thermal boundary condition, these groundwater warming se-

B. L. Kurylyk et al.: Shallow groundwater thermal sensitivity to climate change and land cover disturbances
ries exhibit a convex upward curvature.The results for the two depths (5 and 20 m) indicate that the lag between the surface and subsurface warming increases with increasing depth.For the sandy soil, the temperature at a depth of 20 m increases by 1.77 • C after 100 years, whereas at 5 m depth, this magnitude of warming was realized after only 14 years.Thus, for initially uniform conditions, deeper aquifers will generally remain colder longer than shallow aquifers, as it takes longer for the warming signal to be advected or conducted downwards.Furthermore, Fig. 7a also indicates that soils with a higher thermal diffusivity (i.e., sand) will initially transport the surficial warming signal through the subsurface more rapidly than soils with lower thermal diffusivity (i.e., peat).However, because the subsurface is slowly equilibrating with the new constant surface temperature, the solid series representing the results for the different depths and soils begin to converge as time increases.
In the case of vegetation regrowth, the surface-temperature warming due to the land cover disturbance would be temporary.As an illustrative example, Fig. 7a (dashed lines) shows the groundwater warming produced by Eq. ( 11) at two depths and for two soils due to a sudden 2 • C increase in surface temperature that persists for only 25 years (t 1 ; Eq. 10).If desired, the equation could be further enhanced to accommodate a gradual cooling phase, rather than the instant cooling employed in the present study, using the more general formula described by Menberg et al. (2014).In Fig. 7a, the dashed and solid lines overlap prior to the cooling phase occurring at 25 years.The dashed temperature curves after 25 years represent the thermal recovery period.The groundwater warming curve for a depth of 5 m and the more diffusive soil (sand) is sharp, whereas the groundwater warming curve for a depth of 20 m and the less diffusive soil (peat) is more diffused and lagged.For example, the maximum groundwater warming (0.88 • C) for the peat soil at a depth of 20 m occurs at 33 years, which is 8 years after the surface warming has ceased.Thus, thermal impacts to coldwater streams caused by deforestation may persist several years after vegetation regrowth has occurred, particularly if groundwater discharge to the stream is sourced from a deeper aquifer.However, these effects would likely not be significant as the warming signal would be strongly damped at such depths.
Figure 7b shows the aquifer thermal sensitivities in response to a sudden permanent (solid lines) or temporary (dashed lines) step increase in surface temperature, which correspond to the same warming scenarios as shown in Fig. 7a.As indicated in Eq. ( 17), these thermal sensitivity curves are similar to the groundwater warming curves (Eq.11 and Fig. 7a), but are scaled by a factor of T −1 .Hence, the thermal sensitivity curves due to a step increase in surface temperature are normalized with respect to the boundary temperature increase and are thus independent of the T value.The results presented in Fig. 7 clearly demonstrate that shallow groundwater will initially warm rapidly in response to permanent deforestation and then the rate of temperature increase will decrease with time.This arises due to the initially high thermal gradient and heat conduction arising from the abrupt surface step change in temperature.The resultant impacts of groundwater warming on streambed conductive and advective heat fluxes should be considered in models that simulate stream temperature warming due to deforestation -at least for streams where groundwater discharge has been shown to influence stream temperature.Of particular note, small headwater streams, which are often groundwater dominated, can warm more rapidly than larger streams in response to deforestation because, for natural vegetative conditions, smaller streams typically experience more shading than larger rivers (e.g., Caissie, 2006).
The results shown in Fig. 7 are presented for a recharge scenario (q = 0.2 m yr −1 ).This approach is conservative because recharge environments will typically warm more rapidly in response to rising surface temperatures than discharge environments, as conduction and advection are acting in parallel in the former case.The analytical solutions provided in this study for simulating subsurface warming due to long-term surface-temperature trends (Eqs.11, 13, and 15) are better suited for recharge environments than discharge environments, as groundwater discharge can bring up warm groundwater from deeper within the aquifer in accordance with the geothermal gradient.This phenomenon is not accounted for in the uniform initial conditions (Eq. 9).These solutions can be modified to allow for linearly increasing temperature with depth to account for the geothermal gradient (Kurylyk and MacQuarrie, 2014;Taniguchi et al., 1999a, b), but this adds complexity to the resultant sensitivity formulae.Also as previously noted, this study is primarily concerned with shallow aquifers where heat fluxes due to surface-temperature changes can dominate the influence of the geothermal gradient.

Impacts of climate change on groundwater temperatures
Equations ( 13) and ( 15) can be employed to investigate the sensitivity of groundwater temperatures to long-term gradual surface-temperature changes such as those experienced during climate change.The IPCC (2007) multi-model results (Fig. 2) are globally averaged results, and these data will be used to form the surface boundary conditions for the illustrative examples.

Exponential and linear boundary conditions
The IPCC air temperature anomalies for this century produced by the conservative emission scenario B1 were fitted to a linear surface-temperature function (Fig. 2).The best fit between the linear function and the projected B1 air temperature warming was obtained with a slope β of 5.41 × 10 −10 • C s −1 (1.7 • C per century, see Eq. 12).Also,  15) where T 1 , b, and c are −1.59• C, 1.59 • C, and 3.68 × 10 −10 s −1 , respectively (to match the IPCC A2 projections).The subsurface thermal properties were for 50 % saturated soil (Table 2), and the recharge rate was 20 cm yr −1 .The aquifer thermal sensitivities shown in (c) and (d) were calculated with Eqs. ( 18) and ( 19), respectively.
the exponential function was employed to represent the IPCC multi-model results obtained using the more aggressive, nonlinear A2 emission scenario (Fig. 2).The optimal exponential fit was obtained with fitting parameters b and c of 1.59 • C and 3.67 × 10 −10 s −1 , respectively (Eq.14).The root mean square error (RMSE) values for the exponential and linear fits are presented in Fig. 2. The fitting parameter T 1 (T 0 − b) can be adjusted to obtain the desired initial temperature, and herein we consider the subsurface warming (rather than the temperature per se) by setting initial temperatures to 0 • C (i.e., T 1 = −b).

Groundwater warming due to climate change
Equation ( 13) was employed to illustrate how an idealized, shallow aquifer would respond to a slow linear surface temperature rise (Fig. 3c). Figure 8a shows the groundwater warming results at different depths and for different soils calculated with Eq. ( 13) by applying a 0.017 • C yr −1 linear surface warming as the boundary condition (B1; Fig. 2).The starting date is the year 2000.Similar to the results presented above for land cover disturbances, the surface warming is more rapidly propagated to shallower depths (i.e., 5 m vs. 20 m) and for more thermally diffusive soils (sand vs. peat).After 100 years, the 1.7 • C surface warming produced a 1.6 • C increase in groundwater temperature for the sandy soil at a depth of 5 m (solid red series), but only a 0.94 • C increase for the peat soil at a depth of 20 m (dashed-black series; Fig. 8a).
Figure 8b shows the groundwater warming results produced with the analytical solution with an exponentially increasing surface temperature (Eq.15).The soil thermal properties and recharge rates are identical for the results shown in Fig. 8a and b, and thus the only difference between the two figure panels is the surface-temperature boundary condition.The results shown in Fig. 8a and b for a given soil type and depth (i.e., same color and line type) begin to significantly diverge after approximately 30 years because the IPCC A2 multi-model projections exhibit more extreme warming than the B1 projections after 2030 (Fig. 2).

Groundwater thermal sensitivity due to climate change
Figure 8c and d show the groundwater thermal sensitivity (Eqs.18 and 19) results due to the linear surface warming and the exponential surface warming shown in Fig. 8a and b, respectively.Although the surface warming scenario shown in Fig. 8b is much more pronounced than that shown in Fig. 8a, it is interesting to note that the groundwater thermal sensitivity results for these warming scenarios are very similar (Fig. 8c and d), since the thermal sensitivity is essentially the thermal effect divided by the driving cause.Due to the lag between the surface warming and the subsurface thermal response, the subsurface thermal regime will never reach equilibrium with the surface thermal regime when the boundary condition represents continuous surface temperature increases.Hence, the groundwater thermal sensitivities will never attain unity unless a stable surface-temperature regime is eventually established.However, Fig. 8c and d indicate that the groundwater thermal sensitivity increases with time as the magnitudes of both the surface and subsurface temperature warming increase, and thus the relative impact of the lag decreases.For example, after 100 years, the thermal sensitivity of the sandy soil at a depth of 5 m is about 0.90 for both the B1 linear warming scenario (Fig. 8c) and the A2 exponential warming scenario (Fig. 8d).Thus, shallow groundwater at this depth and for these conditions would be expected to warm by approximately 90 % of the surface temperature increase within 100 years.

Implications for groundwater-dominated streams and rivers
The consideration of groundwater temperature in stream temperature modeling is especially relevant in small streams where surface heat fluxes do not dominate the total energy budget.In fact, small streams are generally very dependent on groundwater inputs and temperatures, and their low thermal capacity (shallow depth and volume) makes them very vulnerable to any surface or subsurface-energy flux modifications (e.g., Matheswaran et al., 2014).This has been shown in many timber harvesting studies, where the smallest streams have experienced the greatest increase in stream temperature following forest removal (e.g., Brown and Krygier, 1970).Thus, quantifying future changes in shallow groundwater flow and temperatures is essential for a better understanding of the future thermal regimes of groundwater-dominated rivers and associated impacts to aquatic organisms (Kanno et al., 2014).The results presented in Fig. 8 demonstrate the limitations inherent in inferring future stream warming from stream thermal sensitivities obtained from seasonal stream and air temperature data.For instance, the seasonal groundwater thermal sensitivity ( ) values presented in Fig. 6 indicate that groundwater temperature beyond 10 m depth generally exhibits minimal sensitivity to seasonal variations in weather.Thus, stream thermal sensitivities obtained from seasonal air and stream temperature data are typically low for groundwater-dominated streams (Kelleher et al., 2012).However, as Fig. 8c and d illustrate, groundwater warming at depths greater than 10 m may still be significant in response to long-term surface-temperature changes, such as would be experienced under climate change.Due to the inter-relationships between the thermal regimes of stream and aquifers and the differences between the thermal sensitivities of shallow aquifers to short-term (Fig. 6) and longterm (e.g., Figs.7b and 8) surface-temperature changes, it is not generally valid to extrapolate thermal sensitivities for groundwater-dominated streams obtained from sub-annual data to project long-term stream warming.
These results also demonstrate the potential limitations of using relatively short (e.g., < 25 years) records of interannual air and water temperature data to obtain estimations of future stream warming (e.g., Luce et al., 2014).Results for the present study (Fig. 8c and d) indicate that even at a timescale of 25 years, the thermal sensitivities of relatively shallow (e.g., 10 m) groundwater reservoirs may be low compared to thermal sensitivities that could be attained after 100 years of surface warming.These results suggest that what is interpreted as a damped stream thermal sensitivity to inter-annual air temperature variability in the case of groundwater-dominated streams may actually be a delayed thermal sensitivity due to the lag in the warming of groundwater and associated streambed heat fluxes.We acknowledge, however, that employing thermal sensitivities derived from inter-annual temperature data to project future stream warming is preferable to considering thermal sensitivities from seasonal temperature data (Luce et al., 2014).The appropriateness of these methods depends on the depth to the aquifer, the degree of groundwater contribution to the stream/river, the subsurface thermal properties, and the timescale of interest.

Addressing groundwater warming in stream temperature models
The present study demonstrates the importance of surfacetemperature forcing on groundwater temperature, particularly for shallow aquifers.The potential influence of shallow groundwater warming on stream temperatures is not generally considered in existing empirical stream temperature models.The equations proposed in this study can be used to develop an approach to approximate the timing and magnitude of groundwater temperature warming in response to long-term surface-temperature changes.As described below, this information may be integrated within existing stream temperature models that consider streambed heat fluxes.The upper boundary condition for the equations presented in this study is the ground surface temperature.Thus, the projected trends in catchment land surface temperature due to future climate change or land cover disturbances must be obtained prior to utilizing these equations.In the case of climate change without related snowpack changes, mean annual surface-temperature trends are often assumed to follow mean annual air temperature trends (see Mann and Schmidt, 2003).This simplification facilitates the boundary condition generation because air temperature trends can be readily ob-   Taniguchi, 1993) tained from the output of climate models.However, in the case of land cover changes (e.g., urbanization) or snowpack evolution, mean annual air temperature trends may be decoupled from mean annual surface-temperature trends (Mann and Schmidt, 2003;Mellander et al., 2007).In this situation, a simple surface heat-flux balance model can be applied to calculate the surface-temperature changes due to changes in the climate and/or land cover.A detailed discussion on appropriate techniques for simulating these relationships can be found in Mellander et al. (2007), Kurylyk et al. (2013), andJungqvist et al. (2014).
Once the surface-temperature trends are obtained, they can then be fitted to the appropriate boundary condition function (Fig. 3).The associated analytical solution (Table 1) and groundwater thermal sensitivity formula can be utilized to perform simulations of future subsurface warming and/or groundwater thermal sensitivity due to the surfacetemperature change.It should be noted that these solutions only calculate increases in mean annual groundwater temperature and do not account for seasonality.It is generally reasonable to assume that the amplitude and timing of the seasonal groundwater cycle will not be greatly influenced by climate change (Taylor and Stefan, 2009), provided snowpack conditions or the seasonality of soil moisture will not change significantly (Kurylyk et al., 2013).
In addition to the surface-temperature boundary condition terms, the analytical solutions must be parameterized with subsurface thermal properties, vertical groundwater flow information, and effective aquifer depth.Subsurface thermal properties can be obtained from information regarding the soil type and typical water saturation of the sediment overlying the aquifer (Table 2).Vertical groundwater flow rates can be obtained from field measurements (e.g., using heat as a hydrologic tracer, Gordon et al., 2012;Lautz, 2010;Rau et al., 2014) or from regional or local groundwater recharge and discharge maps.Potential changes in groundwater recharge (Crosbie et al., 2011;Kurylyk and MacQuarrie, 2013;Hayashi and Farrow, 2014) and groundwater discharge (Kurylyk et al., 2014a;Levison et al., 2014) due to changes in climate or land cover could also be considered.The aquifer effective depth can be crudely estimated as the average unsaturated zone or aquitard thickness overlying the aquifer (e.g., Fig. 4).Such information may be available from well data, geophysical surveys, or regional maps of the groundwater table depth (Fan et al., 2013;Snyder, 2008).Further research is required to assess approaches for more accurately determining the effective aquifer depth.A reasonable range of the input variables to these equations should be considered to generate an envelope of predicted groundwater warming (see Fig. 4 of Menberg et al., 2014).Such a range could incorporate uncertainties arising from, for example, heterogeneities in soil thermal properties and interannual variability in groundwater recharge (Hayashi and Farrow, 2014).Table 3 lists alternative options for parameterizing the equations presented in this study.The parameter values used in the present study are representative of conditions often observed.
To determine the influence of warming groundwater on stream temperatures, the future groundwater thermal sensitivity can be applied to estimate the resultant changes to streambed heat fluxes.There are different approaches available for estimating streambed heat fluxes from subsurface temperatures depending on whether the total streambed energy flux or the apparent sensible flux is being considered (e.g., Caissie et al., 2014;Moore et al., 2005), but in either case, the streambed fluxes depend on subsurface temperature, particularly the temperature immediately below the stream.These changes in streambed heat fluxes can then be combined with simulated changes in stream surface heat fluxes, and the resultant change in stream temperature can be obtained in a deterministic stream temperature model.Such an approach to estimate long-term evolution of stream temperatures would be more realistic than considering a stream temperature model driven by air temperature only, as both surface and streambed heat fluxes can be important in stream temperature dynamics.

Limitations
The unsteady heat advection-diffusion equation utilized in this study (Eq.2) assumes one-dimensional groundwater flow and heat transport, spatiotemporally invariant groundwater flow, isothermal conditions between the soil grains and water at every point, and homogeneous thermal properties.Flashy groundwater flow regimes with very short residence times (e.g., aquifers with large fractures) may invalidate the assumption of thermal equilibrium between the subsurface environment and the mobile water.In such settings, recharge seasonality may exert strong control on the temperature of groundwater discharge (Luhmann et al., 2011).Horizontal groundwater flow can perturb subsurface thermal regimes, at least in regions with significant horizontal thermal gradients (Ferguson and Bense, 2011;Reiter, 2001), and there may be a vertical discontinuity in vertical water flow across aquifers due to horizontal discharge to surface water bodies.Aquifers that exhibit considerable lateral hydraulic heterogeneities may be characterized by flow regimes that are not well represented by the conceptual model (Fig. 4).
Herein, we propose that the average depth to the groundwater table may be a reasonable approximation for the effective depth (z eff ).This approach assumes that the groundwater temperature at the bottom of the vertical flow tubes is fully mixed and that there is no modification of the temperature signal as the groundwater flows horizontally towards the discharge location (Fig. 4).This assumption may be violated in very shallow aquifers with slow groundwater flow (i.e., low horizontal advection and dispersion) due to vertical conductive heat fluxes from the surface in the vicinity of the discharge location.
In shallow aquifers, groundwater velocity varies seasonally and is driven by the seasonality of precipitation, but subsurface hydraulic storage properties tend to damp the seasonality of groundwater flow in comparison to precipitation.Equation (2) also assumes that no soil thawing occurs as a result of the surface-temperature change, but latent heat absorbed during soil thaw can significantly retard subsurface warming (Kurylyk et al., 2014b).Ignoring soil thaw is reasonable, except in permafrost regions, because in ephemerally freezing regions the dynamic freeze-thaw process only influences the seasonality of groundwater temperature, and does not significantly influence the change in mean annual groundwater temperature in response to long-term climate change (Kurylyk et al., 2014a).
At very shallow depths (e.g., < 3 m), the subsurface thermal regime can be considered to be in equilibrium with the mean annual surface temperature.Because the lag between surface and subsurface warming is negligible in this case, the solutions presented in this study are not overly useful at very shallow depths.Also, at greater depths (e.g., 25 m), the influence of the geothermal gradient should be explicitly considered.In such cases, the equations proposed in this study can be modified to incorporate a geothermal gradient (Kurylyk and MacQuarrie, 2014;Taniguchi et al., 1999a, b).Despite these limitations, the analytical solutions presented here can be employed to obtain reasonable estimates of the evolution of mean annual groundwater temperature due to climate change and land cover disturbances for a broad range of aquifer depths.For example, Menberg et al. (2014) applied a modified form of Eq. ( 11) to calculate groundwater warming trends that generally concurred with measured  groundwater warming trends recorded at forested and agricultural sites in Germany.We anticipate that other studies may also benefit from these approaches.

Summary and conclusions
Stream temperature models often ignore the potential for future groundwater warming.This simplifying assumption is employed because mean annual groundwater temperature is relatively constant (or thermally insensitive) on the intraannual or short inter-annual timescales that it is typically measured.We have suggested in this study that although seasonal surface-temperature changes are damped in the shallow subsurface, long-term changes in surface temperatures can be propagated to much greater depths.This phenomenon has been known for some time in the field of thermal geophysics (e.g., Lesperance et al., 2010), but it is generally overlooked in stream temperature modeling.Due to the difference in the subsurface thermal response to seasonal and multi-decadal surface-temperature changes, it may be inappropriate to infer multi-decadal warming of groundwater-dominated streams based on linear regressions of short-term air and water temperature data.
Previous studies have identified the potential importance of considering shallow groundwater temperature warming when projecting future stream temperature (Kurylyk et al., 2013(Kurylyk et al., , 2014a)).These studies have employed methods that either require extensive surface and subsurface-temperature Hydrol.Earth Syst.Sci., 19, 2469-2489, 2015 www.hydrol-earth-syst-sci.net/19/2469/2015/ data collection or detailed numerical modeling.In many cases, these methods may be prohibitive.Several analytical solutions and associated groundwater thermal sensitivity equations are presented in this study as alternative approaches for estimating a range for the potential timing and magnitude of future groundwater warming in response to climate change or land cover disturbances.These are most applicable to idealized environments, but the methods can be employed to obtain first-order approximations of future groundwater warming in natural environments (see Menberg et al., 2014).The subsurface warming scenarios can be considered within existing stream temperature models to investigate whether groundwater warming is an important consideration for the future thermal regime of a particular stream (Snyder et al., 2015).
The present study has highlighted the importance of shallow groundwater sensitivity to surface warming.Although groundwater warming has been inferred from subsurface temperature-depth profiles at many sites, few long-term data sets of directly measured groundwater temperature exist to corroborate the methods proposed herein (Menberg et al., 2014).The initiation of long-term shallow groundwater temperature monitoring sites would provide a better understanding of the processes linking atmospheric and subsurface warming (e.g., Caldwell et al., 2014).
The Supplement related to this article is available online at doi:10.5194/hess-19-2469-2015-supplement.

Figure 1 .
Figure 1.Heat fluxes at the water surface and streambed for the cross section of a gaining stream or river (modified from Caissie, 2006).

B
. L. Kurylyk et al.: Shallow groundwater thermal sensitivity to climate change and land cover disturbances

Figure 2 .
Figure2.IPCC multi-model globally averaged air temperature anomaly projections for the twenty-first century relative to the air temperature data for 1980-1999 for emission scenarios B1 and A2 (data from, IPCC, 2007).Details concerning the exponential and linear fits to the IPCC projections are given in Sect.3.3.1.Modified fromKurylyk and MacQuarrie (2014).

Figure 3 .
Figure 3. (a-b) The boundary conditions for ground surfacetemperature (GST) disturbances due to land cover changes.Both (a) and (b) represent the boundary condition given in Eq. (10).The difference between these is the duration of the period of warm surface temperatures (t 1 = ∞ in a).(c-d) The boundary conditions for GST due to long-term climate change for conservative (linear; Eq. 12) and aggressive (exponential; Eq. 14) climate scenarios.
Figure 4. (a) Groundwater flow and heat transport in a two-dimensional cross section of an aquifer-stream system.(b) Conceptual model of the physical processes shown in (a).Dashed arrows indicate heat transport, and solid arrows indicate water flow.

Figure 5 .
Figure 5. (a) Temperature-depth profiles for each month obtained from Stallman's equation (Eqs.5-7) for homogeneous soil subject to harmonic seasonal surface-temperature variation.(b) Temperature-time series generated with Stallman's equation for depths of 0, 1, 5, and 10 m.In (a) and (b), the thermal properties for sand at 50 % saturation (Table2) were employed, and a recharge Darcy velocity of 0.2 m yr −1 was assumed.The boundary condition parameters T m , A, ϕ, and p were assigned values of 10 • C, 15 • C, −4.355 radians, and 31 536 000 s (1 year), respectively, to represent typical surface-temperature conditions for a forested site in New Brunswick, Canada (e.g.,Kurylyk et al.,  2013).
Fig. 7a (solid lines)  shows the groundwater warming produced with Eq. (11) at different depths and for different soils Figure 7. (a) Groundwater temperature warming due to a permanent (solid lines) or temporary (dashed lines) step increase in surface temperature vs. the time since the surface warming began.(b) Groundwater thermal sensitivity vs. time for each of the eight scenarios presented in (a).The results shown in (a) were obtained with Eq. (11) driven with the step boundary condition (Eq.10), with T = 2 • C and t 1 = infinity (solid lines) or 25 years (dashed lines).The subsurface thermal properties were taken from the 50 % saturated sand and peat values in Table2, and the recharge rate was 20 cm yr −1 .The results shown in (b) were calculated with Eq. (17) using the same parameters as (a).

Figure 8 .
Figure 8. Groundwater temperature warming due to a linear trend (a) and an exponential trend (b) in surface temperature vs. the time since the surface warming began.(c) and (d) groundwater thermal sensitivity vs. time for each of the six scenarios presented in (a) and (b), respectively.The results shown in (a) were obtained with Eq. (13) with β = 5.41 × 10 −10 • C s −1 based on the IPCC B1 projections and setting T 0 = 0 • C. The results shown in (b) were obtained with Eq. (15) where T 1 , b, and c are −1.59• C, 1.59 • C, and 3.68 × 10 −10 s −1 , respectively (to match the IPCC A2 projections).The subsurface thermal properties were for 50 % saturated soil (Table2), and the recharge rate was 20 cm yr −1 .The aquifer thermal sensitivities shown in (c) and (d) were calculated with Eqs.(18) and (19), respectively.

Table 1 .
Details regarding the four analytical solutions employed in this study.

Table 2 .
Bulk thermal properties of some common soils and their dependence on saturation. *

Table 3 .
Parameters for equations considered in this study.