Potential surface temperature and shallow groundwater temperature response to climate change : an example from a small forested catchment in east-central New Brunswick ( Canada )

Global climate models project significant changes to air temperature and precipitation regimes in many regions of the Northern Hemisphere. These meteorological changes will have associated impacts to surface and shallow subsurface thermal regimes, which are of interest to practitioners and researchers in many disciplines of the natural sciences. For example, groundwater temperature is critical for providing and sustaining suitable thermal habitat for coldwater salmonids. To investigate the surface and subsurface thermal effects of atmospheric climate change, seven downscaled climate scenarios (2046–2065) for a small forested catchment in New Brunswick, Canada were employed to drive the surface energy and moisture flux model, ForHyM2. Results from these seven simulations indicate that climate change-induced increases in air temperature and changes in snow cover could increase summer surface temperatures (range−0.30 to +3.49C, mean+1.49C), but decrease winter surface temperatures (range −1.12 to+0.08C, mean −0.53C) compared to the reference period simulation. Thus, changes to the timing and duration of snow cover will likely decouple changes in mean annual air temperature (mean+2.11C) and mean annual ground surface temperature (mean+1.06C). Projected surface temperature data were then used to drive an empirical surface to groundwater temperature transfer function developed from measured surface and groundwater temperature. Results from the empirical transfer function suggest that changes in groundwater temperature will exhibit seasonality at shallow depths (1.5 m), but be seasonally constant and approximately equivalent to the change in the mean annual surface temperature at deeper depths (8.75 m). The simulated increases in future groundwater temperature suggest that the thermal sensitivity of baseflow-dominated streams to decadal climate change may be greater than previous studies have indicated.


Drivers and importance of ground surface temperature
The impact of climate change on ground surface temperature (GST) is of interest to a diversity of scientific disciplines.For example, hydrologists are concerned with the influence of surface freezing and thawing on infiltration and runoff rates (Williams and Smith, 1989), agricultural scientists have shown that seed germination is affected by surface and near-surface temperature (Mondoni et al., 2012), and geotechnical engineers have linked soil strength properties to surface/subsurface temperature (Andersland and Ladanyi, 1994).Increased GST could also enhance decay rates and CO 2 release from soils and thereby act as a positive feedback mechanism to climate change (Eliasson et al., 2005).Potential effects of changes in winter GST include: altered nutrient concentrations in soil water, enhanced winter root mortality, and decreased runoff quality (Mellander et al., 2007).
Published by Copernicus Publications on behalf of the European Geosciences Union.

B. L. Kurylyk et al.: Surface and groundwater temperature response to climate change
Increases in mean annual air temperature will not necessarily result in equivalent changes in mean annual GST (Mann and Schmidt, 2003;Mellander et al., 2007;Zhang et al., 2005).Other physical processes must be considered to predict the associated increase in mean annual GST.For example, the duration of the snow-covered period is expected to decrease in northern latitudes due to increases in late fall and early spring air temperature (AT), which will reduce the insulating effect of the snowpack (Zhang, 2005) and alter the dynamics of atmosphere-surface heat exchange.A reduction in the snow-covered period would also lead to an increase in the amount of radiation absorbed by the ground (Bonan, 2008).Futhermore, the length of the growing season is expected to increase with an upward shift in the AT regime.Under a deciduous canopy, an increase in early-spring and late-fall foliation/defoliation could affect both thermal and hydrological processes by altering the amount of net radiation and evapotranspiration at the land surface (Bonan, 2008).The net effect of these positive and negative climate change feedback signals can be studied with a process-oriented model capable of simulating surficial thermal and hydrological processes.
Very few local-scale studies have investigated the impact of future climate change on GST regimes.Mellander et al. (2007) used two climate scenarios for 2091-2100, generated with the Hadley global climate model (GCM) and downscaled with a regional climate model (RCM), to predict changes in soil temperature in northern Sweden.They simulated a decrease in the period of persistent snowpack of 73-93 days, an increase in annual soil temperatures of 0.9-1.5 • C at 10 cm depth, and an advance in spring soil warming of 15-19 days.Salzmann et al. (2007) used the data from ten RCMgenerated and two incremental climate scenarios to drive the surface energy balance model TEBAL and predicted a potential range of increased GST (mean +3.5 • C) for a mountainous permafrost region in Switzerland.Each of their RCM simulations were driven with the HadAm3H GCM forced by the A2 or B2 emission scenarios.They suggested that their study should be expanded by using multiple GCMs.

Drivers and importance of shallow groundwater temperature
The subsurface thermal regime is driven by water and energy fluxes across the ground surface and the geothermal flux from the Earth's interior.Seasonal and decadal variations in GST can be propagated downwards via conduction and advection and thereby perturb the temperature of shallow (i.e.< 10 m) groundwater (Taylor and Stefan, 2009).Thus, atmospheric climate change may result in changes to seasonal and mean annual GST (as previously discussed), which could translate to shifts in the timing and magnitude of the seasonal groundwater temperature cycle.Because groundwater temperature is less variable than surface water temperature, groundwater discharge provides a thermal buffer for riverine systems during the summer and winter months (Caissie, 2006;Hayashi and Rosenberry, 2002;Webb et al., 2008).In eastern Canada, summer river temperatures are approaching the critical threshold for salmonids (Breau et al., 2007(Breau et al., , 2011;;Swansburg et al., 2004).Discrete cold-water plumes formed by groundwatersurface water interactions have been shown to provide critical thermal relief for stressed cold-water fishes during hightemperature events (Breau et al., 2011;Cunjak et al., 2005;Ebersole et al., 2003;Torgersen et al., 2012).As the climate warms, cold-water fishes may become increasingly dependent on these groundwater-sourced thermal refugia.
Recently, there have been a number of publications investigating the thermal sensitivity of rivers to environmental conditions (e.g.Kelleher et al., 2012;Mayer, 2012;O'Driscoll and DeWalle, 2006;Tague et al., 2007).Many of these studies have demonstrated that groundwater-dominated streams are less sensitive to AT variability on a seasonal basis; however, the response of groundwater temperature (and consequently the temperature of baseflow-dominated streams) to decadal climate change has not been well-studied.For example, Chu et al. (2008) stated "climate change induced differences in precipitation and temperature that may influence the magnitude and timing of groundwater discharge should be addressed in future analyses [of stream temperature and fish habitat]".Mayer (2012) acknowledged the dearth of information regarding the response of groundwater temperature to climate change and suggested that it posed a challenge for modelling future river temperatures.
Several researchers have attempted to address the impact of climate change on groundwater temperature.Taylor and Stefan (2009) employed an analytical solution to a simple conduction equation with a sinusoidal boundary condition to infer that the change in the mean annual groundwater temperature would be equivalent to the change in mean annual GST for Minnesota, USA.They suggested that the current seasonal amplitude of the groundwater temperature cycle would be relatively unchanged in a warmer climate.Gunawardhana et al. (2011), Gunawardhana andKazama (2011), andKurylyk andMacQuarrie (2013a) employed analytical solutions to a one-dimensional conductionadvection equation with an increasing GST to investigate the subsurface response to climate change in the Sendai Plain, Japan.This type of boundary condition does not allow for an investigation of seasonal trends in groundwater temperature.Gunawardhana and Kazama (2012) also used downscaled data from 15 GCMs for the Sendai Plain to drive a numerical model (VS2DH) of groundwater flow and heat transport to investigate subsurface thermal perturbations due to climate change.These simulations were performed on a coarse temporal resolution (1 yr) and also did not address seasonal effects.Others (e.g.Bense et al., 2009;Ge et al., 2011) have simulated the impact of rising surface temperatures on permafrost degradation, but the focus of these contributions was primarily the hydraulic evolution of suprapermafrost aquifers, rather than the thermal evolution.

Research objectives
The objective of this contribution is to provide answers to the following scientific questions: 1. will changes in mean annual and seasonal GST closely mimic climate change-induced changes in mean annual and seasonal AT?
2. can monthly GST data be utilised to predict monthly groundwater temperature in shallow aquifers by adopting an empirical approach?
3. how will shallow groundwater temperature respond to climate change on a monthly basis and what are the implications for salmonid thermal refugia?These questions will be answered in reference to a small forested catchment with available field data (AT, GST, and groundwater temperature), and in which cold groundwater discharge has been observed to provide thermal refugia for salmonids.Thus, the answers to these scientific questions will primarily be valid for our particular study site.This work differs from previous future GST studies (question 1) by employing multiple GCMs and by investigating climate change impacts on surface processes in a warmer climate, albeit still with seasonal snow cover.The answer to the second equation will provide a simplistic alternative to implementing analytical solutions, which are replete with restrictions, particularly for the GST boundary condition.Finally, the answer to the third question should provide surface water hydrologists and ecologists with valuable information concerning the thermal sensitivity of shallow aquifers and baseflow-dominated rivers and the resilience or sensitivity of cold-water fish habitat to a warming climate.

Study site
The geographic location selected for this study is the Catamaran Brook catchment (∼ 53 km 2 ) in east-central NB, Canada (46 • 52 N latitude, 66 • 06 W longitude). Catamaran Brook is a third-order tributary of the Little Southwest Miramichi River (Fig. 1).The Catamaran Brook catchment has a mixed coniferous (65 %) and deciduous (35 %) Acadian forest cover (Cunjak et al., 1990).Portions of the catchment were clear-cut in 1996 (Alexander, 2006).The annual average precipitation in the region is 1140 mm, with approximately 33 % falling as snow (Cunjak et al., 1993).The region experiences a humid continental climate with arid, cold winters and warm, humid summers (Cunjak et al., 1990).
As indicated in Fig. 2, the Catamaran Brook catchment is covered with a thin blanket of coarse sandy till, which is underlain by a dense clay till and Silurian and Devonian bedrock (Alexander, 2006).The water table is shallow and lies within the surficial sand and gravel deposit (Alexander, 2006).
Discharge in Catamaran Brook remains primarily groundwater-sourced throughout the year (Caissie et al., 1996).Thus, the summer water temperature of Catamaran Brook is generally considerably lower than that of the main stem of the Little Southwest Miramichi River, which has a high width to depth ratio and responds rapidly to solar radiation (Caissie et al., 2007).Juvenile Atlantic salmon (Salmo salar) have been documented seeking thermal refuge in the cold-water plume at the mouth of Catamaran Brook during high temperature events (Breau et al., 2007).During a heat wave in 2010, very dense fish aggregations (> 10 000 fish) were observed in nearby thermal refugia within the Little Southwest Miramichi River (T.Linnansaari, UNB, personal communication, 2012).

Reference period and projected climate data
Future climate scenarios are generated by GCMs driven by emission scenarios that invoke assumptions about future population growth, technology, and economic development (Nakicenovic and Swart, 2000).GCM simulations are performed on coarse computational grids (e.g.250 km by 250 km), thus their results should be downscaled to translate the data from the coarse grid to local climate conditions.Downscaling approaches have been thoroughly reviewed in past contributions (e.g.Wilby and Wigley, 1997;Xu, 1999).A simple downscaling approach is the daily translation (DT) method, which is in the family of "statistical" or "quantilequantile mapping" downscaling techniques (Teutschbein and Seibert, 2012).Statistical downscaling is predicated on the assumption that local climate conditions can be determined from large-scale climate variables using linear or non-linear transfer functions (Jeong et al., 2012a).In the DT method, a GCM is initially run for a reference period containing local observations.Scaling factors for precipitation and AT are then determined from the reference period simulations and local observations using empirical cumulative distribution functions.GCM simulations for a future time period (emission scenario) are then downscaled by applying the daily scaling factors.
Many more complex statistical downscaling methods have been developed; one of these is the hybrid multivariate linear regression (HMLR) method (Jeong et al., 2012a,b).
In this method, the local climate variables are obtained from the GCM simulations using multiple regression functions determined from reference period simulations.Because regression-based methods often have difficulty producing the observed variability in local climate predictands, a stochastic generator is used to increase the variance in the datasets.These two methods (i.e.DT and HMLR) were employed in this study.
Output from GCMs can also be dynamically downscaled with RCMs.RCMs produce results on a finer computational grid than GCMs, using GCM output as boundary values.However, RCMs tend to introduce additional biases (Wilby et al., 2000).For this reason, results from RCMs are often bias-corrected.This can be accomplished through techniques similar to those used for downscaling.In the present study, both of the dynamically downscaled climate series (CRCM 4.2.3 aev-A2 and CRCM 4.2.3 agx-A2, Table 1) were bias-corrected using the DT method addressed earlier.
This present study is part of a multi-disciplinary initiative investigating salmonid thermal refugia and their sensitivity to climate change.The climate scenarios utilized in this study (Table 1) were provided by collaborating climatologists.The HMLR-downscaled data were contributed by the Université du Québec Institut National de la Recherche Scientifique (D. Jeong, INRS, personal communication, 2011), while the other climate data series were produced from the third Coupled Model Inter-comparison Project database of GCM output (CMIP3, Meehl et al., 2007) and dynamically downscaled using the Canadian Regional Climate Model (CRCM4.2.3, de Elia et al., 2008;Huard, 2011) or statistically downscaled with the DT method (Huard, 2011).In total, seven projected climate scenarios (Table 1) were produced for the period of 2046-2065 using six GCMs and one RCM, two downscaling methods, and three emission scenarios.This period was selected for the present study due to the number of available downscaled climate scenarios.These particular climate scenarios were selected because they span the range of plausible future climatic conditions for eastcentral New Brunswick.These climate data provide the basis for predicting future GST and groundwater temperature.The term "Run ID" (Table 1) refers to a particular simulation performed within the RCM.In this case, the primary difference between the two RCM runs (Agx and Aev) is the GCM driver (CGCM3 and Echam5).
Figure 3 provides the projected changes in mean annual precipitation and mean annual AT for each of the seven scenarios listed in Table 1.The observed/reference climate data  were taken from the Environment Canada (EC) database of adjusted and homogenised Canadian climate data (EC, 2011).All of the scenarios predict a rise in mean annual AT (with range of +0.4-3.9 • C for the 2046-2065 period compared to the 1961-2000 period), but the projections for annual precipitation vary significantly in magnitude and direction (−12 to +49 %).It should be noted that the intent of this contribution is not to provide an in-depth analysis of the climate data but rather to drive the surface and subsurface models with multiple plausible climate scenarios produced with a variety of established methods.Thus, the intent is to examine the sensitivity of GST and groundwater temperature to external forcing rather than to make an assessment of the accuracy of predictions regarding the future state of surface and subsurface thermal regimes.

Physically based model of surface temperature
The physically based surface energy balance model ForHyM2 (Forest Hydrology Model v.2, Arp and Yin, 1992; Yin and Arp, 1993) was selected to simulate daily GST from the reference and projected climate data.This physically based model is capable of simulating the complex relationship between future climate change and the lower atmosphere-surface energy exchange.ForHyM was originally developed for simulating water fluxes through shallow forest soils; this model successfully reproduced data from a deciduous forest in Ontario, Canada, and a coniferous forest in Quebec, Canada (Arp and Yin, 1992).Yin and Arp (1993) later developed ForSTeM, to simulate soil temperature.ForSTeM was created to be coupled to ForHyM, and these two models with subsequent revisions are collectively referred to as ForHyM2.ForHyM2 has been applied at many other sites (e.g.Balland et al., 2006;Bhatti et al., 2000;Houle et al., 2002;Meng et al., 1995;Oja et al., 1995) to simulate hydrologic fluxes and/or soil temperatures; at all of these sites, model simulations closely agreed with field observations.ForHyM2 has also performed well at simulating snowpack depths under various forest canopies (Balland et al., 2006;Houle et al., 2002).ForHyM2 requires daily AT, precipitation, and relatively few site characteristics for input conditions.Thus, simulations can be performed without extensive field work or site parameterisation.GST simulations can be started without initial conditions, provided that the first time step is not during a period with snow on the ground.ForHyM2 simulations were conducted in a simplified manner (one-dimensional) with the entire catchment represented as a point; thus surface heterogeneities were ignored.The observed climate data and downscaled projected climate data were utilised to drive the ForHyM2 simulations on a daily time step.The GST was taken as the temperature of the forest floor surface.A more detailed description of the ForHyM2 model mechanics and input parameters is included in the supplementary material.

Comparison of uncalibrated ForHyM2 simulations with measured GST
To test the accuracy of performing uncalibrated GST simulations in ForHyM2, simulated GST results were generated for a reference period containing GST measurements.Alexander (2006) installed temperature loggers (VEMCO Minilog) in the Catamaran Brook catchment to record hourly AT, GST, and groundwater temperature (see temperature loggers, Fig. 2).AT observations were recorded in shelters at a height of 1.5 m above the ground surface.Groundwater temperature data were recorded in the monitoring well indicated in Fig. 1 at four depths: 1.5, 2.75, 5.25, and 8.75 m (Fig. 2).GST data were recorded near the monitoring well.The temperature data utilised in the present study were collected during 1 October 2001-30 September 2003.Additionally, measured precipitation data from the EC historical weather database for the Miramichi weather station (EC, 2010) were used as input for the ForHyM2 simulations.In general, the ForHyM2-simulated GST time series were in agreement with the GST observations (Fig. 4).The associated coefficient of determination (R 2 ) was 0.98, thus ForHyM2 was judged to have performed favourably for these site conditions and meteorological data considering that no model calibration was undertaken.

Empirically based estimation of groundwater temperature
Several previous studies investigating the impact of climate change on groundwater temperature have employed an analytical solution to the following governing equation for one-dimensional heat transport in regions of significant groundwater flow (Domenico and Schwartz, 1990): where λ is the thermal conductivity of the subsurface (M L t −3 T −1 ), T is the spatiotemporally varying groundwater temperature, q is the vertical Darcy flux (L t −1 ), C and C w are the volumetric heat capacities of the medium and water, respectively (M L −1 t −2 T −1 ), z is depth (L), and t is time.Analytical solutions to Eq. ( 1) have been derived for GST that is linearly or exponentially increasing on a decadal basis (Taniguchi et al., 1999;Kurylyk and MacQuarrie, 2013a) or periodically varying on a seasonal basis (Goto et al., 2005;Stallman, 1965).However, these boundary conditions are inappropriate for the present study.A linearly or exponentially increasing GST boundary condition ignores seasonal variations in temperature.Furthermore, a periodically varying (i.e.sinusoidal) GST boundary condition is a poor approximation of the annual GST cycle in seasonally snow-covered catchments due to the insulating effect of the winter snowpack (Lapham, 1989;Zhang, 2005).There are many other limitations associated with Eq. ( 1), including spatiotemporally constant groundwater flux and thermally homogeneous subsurface properties.Taylor and Stefan (2009) utilised a solution to a simplified form of Eq. ( 1) that ignored advective heat transport due to groundwater flow.However, advective heat transport could be significant in the Catamaran Brook basin, particularly in the shallow aquifer.The subsurface heat transport modelling capabilities in ForHyM2 (Supplement) were also not utilised for the present study because the model equations ignore advective heat transport due to groundwater flow.
In light of the limitations associated with the ForHyM2 model and previously published analytical solutions, an empirical GST-to-groundwater temperature transfer function was adopted for simulating the monthly groundwater temperature response to rising GST: where GWT i is the groundwater temperature for month i, GST (i−L) is the GST for month (i − L), and D, L, and B are empirical parameters that are temporally constant, but depthdependent.Although this equation is empirically based, the equation parameters relate to physical processes.Analytical solutions to the transient conduction equation have demonstrated that the seasonal groundwater temperature cycle is damped and lagged with respect to the seasonal GST cycle.These damping and lagging effects increase with depth (Bonan, 2008, p. 134;Taylor and Stefan, 2009).The unitless D parameter produces the damping effect of the subsurface thermal diffusivity, while the L parameter (units of months) produces the lagging effect between a GST change and its subsurface realisation.The B term ( • C) accounts for shallow heat transfer phenomena that may not be included in the other two parameters, which are conduction based.These phenomena include: latent heat effects due to freezing and thawing, vadose zone heat transfer processes in the vapour phase, and groundwater advection.
There are admittedly limitations to adopting an empirical approach for relating GST and groundwater temperature.However, linear and non-linear regression-based ATstream temperature transfer functions have been developed for examining the response of surface water temperature to future warming climates (Kelleher et al., 2012;Mayer, 2012;Mohseni et al., 2003, and references therein).Furthermore, GST is the actual driver for shallow groundwater temperature, whereas AT is merely used as a surrogate for radiation, the primary driver of river temperature (Allen and Castillo, 2007).Thus, the application of this GSTgroundwater temperature transfer function should be at least as insightful as its surface water counterparts.
The measured groundwater temperature and GST collected by Alexander (2006) were used to estimate the depthdependent values of D, L, and B by minimising the root mean square error (RMSE) between the measured and simulated monthly groundwater temperatures for each depth.It should be noted that two of the data loggers were installed in the sand and gravel aquifer, while the other two were installed in the clay aquitard (Fig. 2).Future projections of groundwater temperature at each depth were obtained by driving the GST-groundwater temperature transfer function with the future projections of GST obtained from ForHyM2.The flow of data from the GCMs to the GST-groundwater temperature transfer function is illustrated in Fig. 5.

Projected climate change impacts on GST and snow cover
The most pronounced increase (2.64 • C) in mean annual GST, compared to the reference period simulated GST (1961GST ( -2000)), was for the MIROC 3.2 HIRES-A1B climate data, while the only projected decrease in the mean annual GST (−0.15 • C) was obtained for the CSIRO Mk3.0-B1 climate data (Fig. 6).It should be noted that changes in mean annual AT for these particular scenarios were +3.96 and +0.39 • C, respectively (Fig. 3).Thus, it appears that the projected increases in mean annual AT are damped at the ground surface.
In addition to the projected variations in mean annual GST, the ForHyM2 results also suggest that changes will occur to  seasonal GST.As indicated in Figs. 6 and 7, the trends in seasonal GST do not necessarily reflect the trends in mean annual GST.For example, the projected changes in mean annual GST are generally positive, whereas the trends in seasonal GST are typically positive for the spring, summer, and fall, but negative for the winter.Additionally, Table 2 indicates that the changes in seasonal GST will not necessarily closely follow the changes in seasonal AT, although this appears to be the case in the summer.Figure 8 indicates that the snow-covered period is expected to decrease under the various climate scenarios.The predicted reduction in the average snow-covered period could potentially range from 13 (CSIRO Mk3.0-B1) to 49 days (CRCM 4.2.3 aev-A2).Additionally, Fig. 8 suggests that the mean annual snowpack depth will also likely decrease.Estimated reductions in mean annual snowpack depth range from 35 (CSIRO Mk3.0-B1) to 77 % (CRCM 4.2.3 aev-A2).The simulated snowpack depths for the other five climate scenarios are within the range of the CSIRO Mk3.0-B1 and CRCM 4.2.3 aev-A2 results.

GST to groundwater temperature transfer function
The best fits for the L, D, and B parameters in the GSTgroundwater temperature transfer function for each of the four depths are presented in Table 3. Figure 9 gives plots for the measured and simulated monthly groundwater temperature at each of the four depths.The resultant correlation coefficient (R value) and RMSE between the measured and simulated data are also indicated.
According to the analytical solution employed by Taylor and Stefan (2009), the lag term (L) and the natural logarithm of the damping term (D) should be linearly related to the depth.Figure 10 indicates that the relationship between the depth and the L term or the natural logarithm of the D term can be reasonably approximated by a line with a zero intercept.

Projected climate change impacts on groundwater temperature
For the sake of clarity, only the groundwater temperature results from the reference period simulation  and the maximum (MIROC 3. CSIRO-Mk3.0-B1climate data is relatively unchanged from the reference period.
Figure 12 shows the change in average monthly groundwater temperature projected for each of the climate scenarios.As expected, the variability in the changes of monthly groundwater temperature decreases with depth for each of the climate scenarios.

Relationship between climate-induced changes in AT and GST
Figures 3 and 6 demonstrate that the magnitude of the increases in projected mean annual AT for the seven climate scenarios for this study site (+0.4 to +3.9 • C) are larger than the simulated changes in mean annual GST (−0.15 to +2.64 • C).Thus, the findings of the present study concur with Mann and Schmidt (2003), who proposed the debated notion that decadal mean annual GST changes will not necessarily track mean annual AT changes (see Chapman et al., 2004;Schmidt and Mann, 2004).Decadal snow-cover evolution can decouple mean annual GST and mean annual AT trends by altering the winter thermal resistance between the lower atmosphere and the ground surface.Figure 4 indicates that the variability in simulated and measured GST is considerably less in winter than in the other seasons due to the insulating effect of the snowpack.This effect is simulated in ForHyM2 by the inclusion of a thermal resistance layer during snow-covered periods (Supplement).Figure 8 suggests that the snowpack insulating effect will likely be reduced in the coming decades due to a reduction in the snow-covered period and the average snowpack depth.Thus, even when snow cover exists in a warmer climate, the equivalent thermal resistance will be limited by snowpack thinning.This will increase the winter heat transfer between the lower atmosphere and the ground surface and result in colder winter GST as indicated in Fig. 7a.This decrease in simulated winter GST limits increases in mean annual GST and may be an important negative climate change feedback mechanism for the subsurface thermal regime.Table 2 indicates that changes in average summer GST (−0.30 to +3.49 • C) will closely replicate changes in average summer AT (−0.30 to +3.49 • C) for this location.Table 2 also suggests that increases in fall GST will be slightly less than increases in fall AT (average fall ratio, GST/ AT = 93 %) and that increases in spring AT will be damped at the ground surface (average spring ratio, GST/ AT = 72 %).The spring damping effect is caused by a lingering snowpack.Thus, changes in seasonal GST will likely follow changes in seasonal AT during the warmer period, but not during the colder period.

The empirical GST to groundwater temperature function
The results presented in Fig. 9 indicate that an empirical relationship of the form of Eq. ( 2) is flexible enough to reproduce groundwater temperature at coarse spatial (∼ 2 m) and temporal (monthly) resolutions.As expected, the RMSE values were higher at shallow depths where the seasonal variability in groundwater temperature is greater.However, the relatively constant correlation coefficients (R values) suggest that the function is consistent in addressing the variability in measured groundwater temperature to depths of approximately 9 m. Figure 9 suggests that the function has a problem with simulating groundwater temperature in April.We expect that this error arises due to advective heat transport associated with the spring snowmelt and major recharge event during this period and/or with latent heat absorbed during the thawing of the upper few cm of soil.The coefficients of determination (R 2 ) in Fig. 10 illustrate that the lag term (L) and the damping factor (D) generally varied with depth as expected.The slope of the natural logarithm of the damping factor vs. depth relationship can be used to infer a soil thermal diffusivity of 7.1 × 10 −7 m 2 s −1 (see Eq. 4 in Taylor and Stefan, 2009) which is in agreement with the typical saturated soil thermal diffusivity (7.0 × 10 −7 m 2 s −1 ) reported by Bonan (2008).Thus, the behaviour of the D and L terms concurred with our expectations based on classic heat conduction theory.As discussed, B is an encompassing parameter to account for any heat transfer process other than saturated zone conduction.Near-surface phenomena (e.g.freezing, evaporative losses to the atmosphere, thermal conductivity heterogeneities due to variable saturation) and groundwater flow can significantly perturb a temperature-depth profile.These phenomena exhibit seasonality, which would contribute to the errors in the estimated groundwater temperature for particular months.Table 3 indicates that the values of B decreased with increasing depth, which we expected given our proposition that B accounts for near-surface phenomena and advective heat transport.Clearly, the effects of near-surface phenomenon would decrease with depth.Advective heat transport would also decrease with depth given the presence of the clay aquitard (Fig. 2).In general, the behaviour of the equation parameters (i.e.D, L, and B) is suggestive that they produce the physical effects that we postulated.

Shallow groundwater temperature respond to climate change
Figures 11 and 12 indicate that the MIROC 3.2 HIRES-A1B and CSIRO Mk3.0-B1 climate data resulted in the highest and lowest groundwater temperatures, respectively.These results were consistent for each month and depth.This is not surprising given that these two climate scenarios also resulted in the highest and lowest mean annual AT (Fig. 3) and mean annual GST (Fig. 6).Thus in general, increased mean annual AT will result in increased mean annual GST and mean annual groundwater temperature; however, the response of monthly groundwater temperature to rising AT (and consequently rising GST) may be complex.For example, Fig. 12 indicates that at 1.5 m depth the maximum changes in the groundwater temperature for the MIROC 3.2 HIRES-A1B scenario (+3.43 • C) will occur in May, and the minimum changes in groundwater temperature (+1.26 • C) will occur in March.Thus, climate change-induced increases in shallow groundwater temperature in the Catamaran Brook catchment likely exhibit seasonality.Figure 12 demonstrates that this seasonality will decrease with depth.At a depth of 8.75 m, the increase in monthly groundwater temperature is approximately constant and equivalent to the rise in mean annual GST.This study has also demonstrated that, due to increased summer GST and decreased winter GST (Fig. 7), the range in the annual cycle of groundwater temperature will increase.For example, the difference between the maximum (September) and minimum (February) average monthly groundwater temperature at a depth of 1.5 m was 13.8 • C for the reference period simulation, and the differences projected for the seven climate scenarios ranged from 14.0 to 17.2 • C.These findings contrast with those of Taylor and Stefan (2009) who proposed that the amplitude of shallow groundwater temperature would not be significantly affected by a warming climate.
One motivation for this contribution was to examine the effect of rising groundwater temperature on groundwatersourced salmonid thermal refugia.Stream temperature researchers have defined the thermal sensitivity of a stream as the slope of the AT-stream temperature relationship (e.g.Mayer, 2012).Following this approach, we define the subsurface thermal sensitivity as the ratio of the climate change-induced increase in summer (June-August) groundwater temperature to the increase in summer AT.We are primarily concerned with summer temperatures as groundwater discharge is most critical for the preservation of salmonid populations during this time of the year (Breau et al., 2007;Cunjak et al., 2005).Table 4 presents the simulated subsurface thermal sensitivities for each climate scenario at each of the four depths.The sensitivity values are negatively correlated with depth because simulated summer GST increases closely followed summer AT increases (Table 2), but simulated changes in mean annual GST were significantly damped compared to changes in mean annual AT.With increased depth, the subsurface thermal effects of seasonal GST are reduced, and the mean annual GST becomes the primary driver of groundwater temperature.The mean summer subsurface thermal sensitivity values for Catamaran Brook (0.66-0.88,Table 4) are somewhat higher than stream temperature sensitivities obtained by previous researchers examining AT-surface water temperature relationships for shorter periods (Kelleher et al., 2012;Mayer, 2012;O'Driscoll and DeWalle, 2006).For example, Kelleher et al. (2012) found that the average thermal sensitivity of smaller streams (0.52) was less than that of larger (stream order > 3) streams (0.83) for stream systems in Pennsylvania, USA.Mayer (2012) used weekly AT and stream temperature data and found that the summer thermal sensitivities varied from 0.25 to 0.58 for six streams in Oregon, USA.Kelleher et al. (2012, and references therein) proposed that thermal sensitivities of stream water temperature to AT variations derived from high-frequency temperature data (e.g.seasonal, monthly, weekly, or daily) could be employed to estimate stream temperature sensitivity to climate change.This approach may not be valid in groundwater-dominated streams, because it is fundamentally based on the premise that groundwater temperature will respond in the same manner to both high-frequency and low-frequency AT variations.Groundwater temperature may be insensitive to daily or even seasonal AT variations, but it is sensitive to decadal AT variations, as we have demonstrated by the subsurface thermal sensitivity values in Table 4. Thus extrapolated seasonally derived thermal sensitivities obtained for baseflowdominated streams may underestimate their true thermal sensitivity to decadal climate change.
The changes in shallow (1.5 or 2.75 m) summer groundwater temperature projected for the MIROC 3.2 HIRES-A1B data are approximately 3 • C (Fig. 12).Changes to the temperature of groundwater discharge on this order would negatively impact riverine ecosystems by increasing the temperature of groundwater-sourced thermal refugia and increasing local surface water temperature in groundwaterdominated streams and rivers by increasing the heat flux at the river bed (see Caissie et al., 2007;Hebert et al., 2011).Although groundwater temperature will respond to decadal climate change, groundwater-dominated tributaries will generally continue to remain colder than the main stems and thereby induce riverine thermal heterogeneity.Catchments for tributaries that provide thermal refugia should be protected from deforestation (Alexander, 2006;Bourque and Pomeroy, 2001) and aggregate extraction (Markle and Schincariol, 2007), which have been shown to increase shallow groundwater and river temperatures.

Limitations of the approach
Our modelling approach for estimating future groundwater temperature from GST assumes that the empirical parameters (B, D, and L) will be unaffected by a warming climate.It seems reasonable to assume that this will be the case for D and L based on the physical processes they represent.The subsurface damping and lagging effects are primarily controlled by the period of the seasonal GST cycle, the soil diffusivity, and depth (Bonan, 2008), and none of those controls will likely be significantly affected by atmospheric or surficial climate change in this catchment.Climate change-induced variations in the timing or magnitude of precipitation may impact soil moisture and consequently thermal diffusivity; however, these effects would be more noticeable in a catchment with a deeper water table.It is more likely that the B parameter would be affected by a warming climate.For example, increased precipitation and recharge rates in this catchment (Kurylyk and MacQuarrie, 2013b) could potentially lead to accelerated advective heat transport due to increased groundwater flow.Increases in extreme precipitation events could potentially alter the B parameter due to changes in hydrological processes such as intensified subsurface lateral flow (horizontal advection) and heterogeneous surface ponding.Furthermore, at the latitude considered for this study, latent heat effects arising from pore-water phase changes may actually increase in a warmer climate, as winter GST are projected to decrease.As previously discussed, the importance of B decreases with depth, thus the results for deeper depths (e.g.8.75 m) may be more reliable.Also, this empirical function assumes that the subsurface is in dynamic thermal equilibrium with the surface; however, it could potentially take years for the rise in mean annual GST to be fully realised in deeper (> 10 m) aquifer systems.These complex subsurface climate change feedbacks could be simulated with physically based groundwater flow and heat transport models that can accommodate freezing and thawing effects.We anticipate that this research will motivate others to perform additional studies in other aquiferriver systems.

Conclusions
To investigate the relationships between climate change and GST in the Catamaran Brook catchment, simulations were conducted using the physically based surface flux balance model ForHyM2.In general, the ForHyM2 results indicate that the changes in summer GST will mimic changes in summer AT for this catchment.However, the model results also indicate that, due to complex snow-cover evolution effects, winter GST will likely decrease (−1.12 to +0.08 • C).This will result in mean annual GST changes (−0.15 to +2.64 • C) that are significantly damped with respect to mean annual AT changes (+0.39 to +3.96 • C).
Measured groundwater temperature and GST were utilised to develop an empirical relationship between the surface and subsurface thermal regimes for this catchment.This empirical function was then driven with ForHyM2-produced GST projections to estimate the impact of atmospheric climate change on groundwater temperature.Results from these simulations indicated that shallow (1.5 m) groundwater temperature in this catchment is very sensitive to increases in summer AT.High-frequency GST signals, such as seasonal variations in GST, are damped at greater depths; however, lowfrequency GST signals, such as decadal climate change, are retained at these depths.Thus, even groundwater temperature at deeper depths will respond to increasing AT and GST.The results presented in this contribution are site-specific; however, our approach could be employed to investigate surface and subsurface thermal processes in other catchments with seasonal snow cover.
This study has also demonstrated the limitations inherent in predicting future climate change impacts using a single projected climate series based on one emission scenario, simulated with one GCM, and downscaled using only one approach.Climate modelling involves many assumptions and, as such, an array of climate scenarios should be considered.
We have also demonstrated that baseflow-dominated streams may exhibit more sensitivity to climate change than previous contributions have indicated.Salmonids are threatened by rising river temperatures in eastern North America, and the ability of groundwater to buffer rising river temperatures may be overestimated.Future physically based groundwater flow and heat transport modelling will be conducted to further investigate the potential increase in the temperature of groundwater discharge and the subsequent impact to cold-water fish habitat.

Fig. 2 .
Fig. 2. Cross-sectional view of the portion of the Catamaran Brook catchment containing the monitoring well and temperature loggers (adapted from Alexander, 2006).

Fig. 3 .
Fig. 3. Projected changes in average annual precipitation vs. projected changes in mean annual AT for east-central New Brunswick.The climate data from the observed period (1961-2000) were chosen for a baseline for which to compare the future period (2046-2065) period simulations.The 1961-2000 observed data series and each projected climate series (2046-2065) have their own distinct set of statistics (e.g. annual mean and standard deviation).

SimulatedFig. 4 .
Fig. 4. Measured and ForHyM2-simulated GST and ForHyM2simulated snowpack for the Catamaran Brook catchment for the October 2001-September 2003 period.Measured and simulated GST data are relatively constant during the snow-covered period.

Fig. 5 .
Fig. 5.Translation of GCM-simulated data into daily GST-simulations.Daily GST were then averaged on a monthly basis and entered as input to the empirical GST-groundwater temperature (GWT) transfer function.

Fig. 6 .
Fig. 6.ForHyM2-simulated mean annual GST for each climate scenario (2046-2065).Simulated mean annual GST based on the historical climate data (1961-2000) is indicated by the data point to the far left.Error bars indicate one standard deviation in mean annual GST.

Fig. 7 .Fig. 8 .
Fig. 7. ForHyM2-simulated average seasonal GST for each climate scenario (2046-2065) and the reference period.The y axes show various ranges in GST but are drawn to the same scale.Winter = December-February, spring = March-May, summer = June-August, and fall = September-November.

Fig. 9 .Fig. 10 .
Fig. 9. Comparison of observed and simulated monthly groundwater temperature at different depths: (A) 1.5, (B) 2.75, (C) 5.25, and (D) 8.75 m.The y axes are not all drawn to the same scale.

Fig. 11 .Fig. 12 .
Fig. 11.Results from the empirical GST-groundwater temperature transfer function driven by the ForHyM2-produced GST temperature for the reference period and two projected climate scenarios.Of the seven climate scenarios, the CSIRO Mk3.0-B1 and the MIROC 3.2 HIRES-A1B climate scenarios consistently produced the lowest and highest monthly groundwater temperature data, respectively.

Table 1 .
Details for each of the seven climate scenarios employed in the current study.

Table 2 .
Changes in seasonal * AT and ForHyM2-simulated GST for each climate scenario.

Table 3 .
Depth-dependent parameter values for the GSTgroundwater temperature transfer function.

Table 4 .
Simulated summer subsurface thermal sensitivities for each depth, averaged for the seven climate scenarios.