Long-term trends in groundwater recharge and discharge in a fractured bedrock aquifer – past and future conditions

Abstract

In the context of climate change, it is important to understand possible future projections and historical trends of groundwater recharge, flow, and discharge to surface reservoirs. Knowledge of a vast range of possible conditions is required to fully appreciate the variability of the hydrologic cycle and hence the long-term vulnerability of groundwater-dependent habitats. This research investigates historical trends for a groundwater–surface water interacting system that supports a fragile ecosystem in southern Quebec. A transient model was developed using MODFLOW to simulate site-wide groundwater flow for the study area. The model was used to simulate past hydrogeological conditions (1900–2010) using a new data set of available precipitation (rain and snowmelt) and temperature. This data set was used to simulate the overall groundwater budget and to determine groundwater discharge (river baseflow and spring flow) in the study area. This allows for the quantification of century-long trends in flow data, as well as the extreme maximum and minimum flows over 110 years. Recharge was variable, ranging from 41 to 197 mm/year over the study period. Lower recharge rates from 1950 to 1965 induced marked effects on spring flow. Although the trend is not statistically significant, there appears to be, for the second half of the study period (1966–2010), a tendency towards a reversal to an increase for recharge, hydraulic heads, spring flow and baseflows. A longer time series would be necessary to confirm this tendency. The simulated historical trends are compared with flow projections for future scenarios (2041–2070). The confirmation that the natural system has been subjected to a wide range of climatic conditions over the last century helps to inform about its resilience. This study highlights the utility of groundwater flow modeling using historical climate data sets to gain a better understanding of long-term trends for climate change-related hydrogeological and ecohydrological studies.

Dans un contexte de changements climatiques, il est important d’analyser les projections futures possibles et les tendances historiques des écoulements d’eau souterraine (recharge, flux, apports aux eaux de surface). La prise en compte d'une vaste gamme de conditions possibles est nécessaire afin d’apprécier pleinement la variabilité du cycle hydrologique et, conséquemment, la vulnérabilité à long terme des habitats qui dépendent des eaux souterraines. Cette recherche examine les tendances historiques d’un système aquifère-eaux de surface qui soutient un écosystème fragile dans le sud du Québec. Un modèle transitoire a été développé en utilisant MODFLOW pour simuler l'écoulement des eaux souterraines dans la zone d'étude. Le modèle a été utilisé pour simuler les conditions hydrogéologiques passées (1900–2010) en utilisant des données de précipitation disponible au ruissellement, à l’infiltration et à l'évapotranspiration (pluie plus eau de fonte de la neige) et de température. Ces données ont permis de simuler le bilan hydrique souterrain afin de déterminer les apports d’eau souterraine (débit de base des cours d'eau et débit aux résurgences) dans la zone d'étude. Ce travail mène à la quantification des tendances séculaires et des extrêmes d’écoulement souterrain sur une période de 110 ans. La recharge est variable, fluctuant de 41 à 197 mm/an au cours de la période d'étude. La baisse de recharge de 1950 à 1965 induit des effets marqués sur les débits aux résurgences. Bien que non significatives, les tendances pour la deuxième portion de la période d'étude (1966–2010) semblent indiquer une inversion vers une légère augmentation de la recharge, des charges hydrauliques, des débits aux résurgences et des débits de base des cours d'eau. Une série temporelle plus longue serait nécessaire pour confirmer cette tendance. Les tendances historiques simulées sont comparées aux projections de scénarios futurs (2041–2070). La confirmation que le système à l’étude a été soumis à un large éventail de conditions climatiques au cours du dernier siècle contribue à informer sur sa résilience. Cette étude met en évidence l'utilité de la modélisation des écoulements souterrains à partir de données climatiques historiques afin de mieux comprendre les tendances à long terme dans le cadre d’études hydrogéologiques et écohydrologiques prenant en compte les changements climatiques.

Introduction

Groundwater–surface water interactions can be strongly influenced by climate-induced changes in the hydrologic cycle. Local or regional studies that quantify groundwater flow and the connections with surface waters are increasingly important to inform land management decisions in order to maintain source water availability and to sustain healthy ecosystems. For example, biodiversity is enriched by groundwater discharge features such as springs and seeps which support species locally and also feed into downstream lakes and rivers, providing abundant, rich habitats (Wood et al. 2005; Meyer et al. 2007; Winter 2007).

In the last decade, many investigations about the prediction of groundwater-surface interactions in future climate conditions have been conducted in Canada. These studies typically include investigation of trends in hydraulic head and discharge (e.g. baseflow, spring flow, peatland–aquifer interactions) stemming from different future recharge scenarios (e.g. Scibek and Allen 2006; Jyrkama and Sykes 2007; Scibek et al. 2007; Toews and Allen 2009; Sulis et al. 2011; Bourgault et al. 2014; Kurylyk et al. 2014; Levison et al. 2014a; Levison et al. 2014b; Rivard et al. 2014) or historical data sets (e.g. Chen et al. 2004; Rivard et al. 2009).

Recharge is the atmospheric forcing that connects changes in climate to groundwater reservoirs influenced by variations in precipitation (volume, intensity and timing) and temperature. The range of predicted and simulated recharge changes across Canada is quite large and there is much uncertainty associated with the estimated recharge volumes (Kurylyk and MacQuarrie 2013). More field observations are required to fill the existing knowledge gap on the relationship between groundwater recharge and long-term climate change. Estimates of future recharge in several provinces are as follows (not necessarily all for the same future period, using various approaches): +11 to +25% in the Grand Forks aquifer, British Columbia, Canada (Scibek et al. 2007); +53% in the Grand River watershed, Ontario (Jyrkama and Sykes 2007); –59 to +15% in the Chateauguay River watershed, southern Quebec (Croteau et al. 2010); –4 to +15% at Covey Hill, Quebec (Levison et al. 2014b); –6 to +58% in the Otter Brook watershed, New Brunswick (Kurylyk and MacQuarrie 2013); and +14 to +45% in the Annapolis Valley, Nova Scotia (Rivard et al. 2014). Reductions in recharge are usually considered the “worst-case scenario” from an ecohydrological standpoint (e.g. Bourgault et al. 2014), since less recharge can induce declines in groundwater levels which can trigger drying of streams, springs and seeps.

Holman et al. (2012) note that the uncertainty of future conditions needs to be addressed by running a broad range of scenarios incorporating both climate and future socio-economic and land-use changes. In order to help understand and interpret what will happen in the context of climate change, it is useful to investigate historical trends of groundwater recharge, flow and discharge, in addition to future projections. Historical trends provide insight into the hydrologic cycle over a variety of observed temperature and precipitation patterns and extremes. This is particularly important to better understand the long-term sustainability of water resources, as well as the long-term vulnerability and resilience of groundwater-dependent habitats and species (Johnson et al. 2004; Sommer and Froend 2011).

Long-term data sets that can be used to support groundwater modeling investigations of climate change include direct observations of temperature and precipitation and/or groundwater levels, as well as reconstructed paleoclimatic data. Rivard et al. (2009) investigated long-term (up to 50 years) baseflow and groundwater level data, finding no overall trends for national scale recharge decline in Canada. However, regional recharge reduction trends were observed (e.g. in Atlantic Canada). Chen et al. (2004) related climate variability of the Winnipeg region to groundwater levels recorded as early as 1961, finding that increasing temperatures predicted for the future could reduce net recharge and impact groundwater levels. The interannual variability of recharge dynamics has also been related to climate indices (e.g. Tremblay et al. 2011). Brutsaert (2008) determined long-term groundwater storage trends in two large watersheds in Illinois, based on observed stream flow records (since the 1940s) and groundwater level records (since the 1960s). With a method developed using tree ring indices, groundwater levels, precipitation and stream flow data, Hanson et al. (2004, 2006) showed how hydrologic responses induced by climate variations can be better understood over the long term using reconstructed time series.

Past climatic extremes, and general trends encountered using long-term historical data, can provide insight and additional lines of evidence into how groundwater–surface water systems will react to climate-driven perturbations into the future. This type of data and historical investigation of trends have not been widely investigated and thus should be more frequently included in groundwater investigations for a more complete understanding of hydrological systems (Iwatsuki et al. 2010; Ouysse et al. 2010; Holman et al. 2012).

The aim of this research is to investigate historical trends of recharge and discharge in a fractured bedrock aquifer over a 110-year period, from 1900 to 2010. This is done using a new data set of available precipitation (rain and snowmelt) and temperature (Poirier et al. 2014) in a groundwater flow model representing a headwater ecosystem located in the Covey Hill region of southern Quebec, Canada. The simulated historical trends of heads, river baseflows and spring flow rates are compared with flow projections simulated for climate change scenarios for the 2041–2070 period. A previous investigation at this research site, in contrast, focused on peatland–aquifer interactions (and other aspects of the hydrological system) for “current” conditions (2007–2010), and employed a sensitivity analysis to investigate a range of potential future conditions (i.e. cases specifically for +10%, –15% and –30% of calibrated recharge; Levison et al. (2014b). Spring flows, which are now considered in the current study, were not previously addressed in Levison et al. (2014b).

Study area

The research area is the 200-km² Covey Hill Natural Laboratory (Larocque et al. 2006), located on the Quebec–New York State border 65 km southwest of Montreal (Figure 1). This field site has been used for several previous hydrogeological and ecological investigations (e.g. Fournier 2008; Pellerin et al. 2009; Gagné 2010; Lavoie et al. 2013; Levison et al. 2014a; Levison et al. 2014b). The Natural Laboratory is within the Chateauguay River watershed where numerous hydrogeological characterization studies have recently been conducted (Nastev et al. 2008; Croteau et al. 2010; Lavigne et al. 2010; Tremblay et al. 2010).

Figure 1. (a) Location of the study area (Covey Hill Natural Laboratory). The Hemmingford (HMG) weather station is shown on the inset map; (b) regional cross-section. Note that the delineated “watershed limits” correspond to the gauging station watersheds. The Système d'information hydrogéologique (SIH) wells are not represented in this figure. Modified with permission from Levinson et al. (2014b).

Covey Hill, the northernmost extension of the Adirondack Mountains, rises to an elevation of about 345 m above sea level (asl) and is comprised of Cambrian sandstone of the Potsdam group (Covey Hill Formation), deformed and fractured during the Appalachian orogeny (Globensky 1986). Horizontal to sub-horizontal bedding planes dip 1 to 5° (Clark 1966). Groundwater resides in the sandstone aquifer that is generally unconfined and flows primarily in sub-horizontal and sub-vertical fractures (Nastev et al. 2008). Bedrock is exposed near the hilltop, and other areas are covered by thin, permeable and sandy Saint-Jacques till (Lasalle 1981). Glaciolacustrine sediments occur below 220 m asl (Parent and Occhietti 1988) and sandy beach deposits are between 80 and 100 m asl near the foot of the hill (Tremblay et al. 2010). At the base of Covey Hill, littoral sediments, composed of highly permeable sands and gravels from the erosion of the rock substrate by the Champlain Sea, are abundant on the northern side (for additional details see Levison et al. 2014a; Levison et al. 2014b).

One of the oldest peatlands in Quebec is near the summit of Covey Hill at about 300 m asl (Figure 1). It is a 51-ha sphagnum-dominated bog. Organic matter started accumulating approximately 13,000 years BP (Lavoie et al. 2013). At various locations within the Natural Laboratory, there is continuous monitoring of groundwater levels and river flows. Covey Hill is an important recharge zone for the sandstone aquifer (Croteau et al. 2010). Groundwater flows in a radial direction from the hilltop generally towards the north. The water table is shallow, located between 2 and 15 m below the ground surface. Head data are available from nine instrumented private monitoring wells, in three Geological Survey of Canada wells, and in two bedrock piezometers (4.5 and 15 m depth) located near the peatland. The provincial water well database, the Système d'information hydrogéologique (SIH), provides 371 heads distributed throughout the Natural Laboratory (Ministère du Développement durable, de l’Environment et de la Lutte contre les changements climatiques – Québec [MDDELCC] 2010). Annual variations in groundwater levels increase with elevation on Covey Hill (Levison et al. 2014b).

Within the Natural Laboratory, the Allen River, Outardes River and Schulman stream (see Figure 1) are gauged, and their 2007–2010 baseflows ranged from 0.07–0.40, 0.02–0.40 and 0.001–0.024 m³/s, respectively (Levison et al. 2014b). The endangered Allegheny Mountain Dusky Salamanders (Desmognathus ochrophaeus) reside in groundwater discharge zones (small, low-flow springs) on the northeastern face of the hill, and permanent streams are home to the Northern Spring Salamanders (Gyrinophilus porphyriticus), a threatened species (Larocque et al. 2013; Levison et al. 2014a). The groundwater springs discharge from the sandstone aquifer. Five springs on the northeastern face of Covey Hill at elevations ranging from 140 to 170 m asl have been instrumented since 2007 with temperature loggers (Levison et al. 2014a). Springs that flow continuously are found at lower elevations, while intermittent springs (flowing, for example, from autumn to early summer and drying out during the hot summer months) are located farther up the hill (Levison et al. 2014a).

Methods and approach

Available meteorological data and climate change scenarios

Precipitation and temperature data starting in 1961 are available from the Hemmingford weather station (see inset map, Figure 1) which is located about 11 km from the study area (Environment Canada 2010). There were no other weather stations in the vicinity of the Covey Hill Natural Laboratory prior to this period. Snow typically falls between November and March. The long-term average annual precipitation from 1961 to 2010 is 872 mm (Environment Canada 2010). The average annual temperature is 6.4°C and the average monthly temperature ranges from –9.6°C in January to 20.6°C in July with a winter minimum of –30.5°C and a summer maximum of 30.3°C. The Oudin et al. (2005) equation is used to calculate potential evapotranspiration (PET). Oudin et al. (2005) provides PET estimates based on mean daily air temperature and on extraterrestrial radiation, estimated following Morton (1983). Using data from the onset of the Hemmingford weather station, the average annual PET for the 1961–2010 period is 637 mm.

Gridded meteorological data for the 1900–2010 period were generated by the Centre d’expertise hydrique du Québec (Poirier et al. 2014), an agency of the Quebec Ministry of Environment (Ministère du Développement durable, de l’Environnement et de la Lutte contre les changements climatiques – MDDELCC). The spatial interpolation method relies on measured values at climatological stations. These stations, mainly owned by the MDDELCC, are operated by the Ministry’s branch of atmospheric monitoring activities (Direction du suivi de l’état de l’environnement [DSÉE] – Service de l’information du milieu atmosphérique [SIMAT]) through a climate-monitoring program. The interpolation is performed at a daily time step, using simple isotropic kriging with monthly mean variograms at a 0.1° (~10-km) resolution over southern Quebec (Poirier et al. 2014). The available data cover the 1900–2010 period. This study makes use of kriged data of mean temperature and vertical inflows (VIs). The VIs are defined as the sum of precipitation and snowmelt available daily for runoff, percolation or evaporation. In the absence of snow on the ground and when the precipitation is liquid, VIs and precipitation data are the same. Snow accumulation and melting is estimated based on the snowmelt module in the HYDROTEL model, a distributed hydrological model (Fortin et al. 2001; Turcotte et al. 2007). This data set has been used in at least one additional case study in the Lanoraie region of southern Quebec (e.g. Bourgault et al. 2014). It is important to underline the fact that the first weather station was installed around 1920 and the first station measuring snow on the ground was installed around 1965. Before these dates, VI values are less reliable since they are estimated using data from stations situated farther away.

From 1900 to 2010, the average annual temperature varied from 4.1 to 8.3°C. Minimum and maximum annual VIs vary from 606 mm (1964) to 1302 mm (1918). The 1900–2010 VI data were combined with the daily temperature values to calculate monthly net precipitation (P_net = VI – PET) for the entire period. In order to verify the validity of the VI and temperature values, a comparative analysis was conducted using this data set and the meteorological data recorded at the Hemmingford weather station from the onset of the records in 1961 until 2010. Annual and monthly averages of daily temperature and precipitation as well as annual P_net values were used for this comparison.

For the investigation of future climate change, 10 projections for the 2050 horizon (2041–2070) were selected from the 25 dynamically downscaled simulations available for the Covey Hill area (see Levison et al. 2014a). Levison et al. (2014b), in contrast, conducted a sensitivity analysis using the calibrated groundwater model to investigate the impacts on the flow system from changes in recharge ranging from –30% to +10%. The simulations are from the Canadian Regional Climate Model (CRCM) (Music and Caya 2007) and the North American Regional Climate Change Assessment Program (NARCCAP), and were supplied by the Ouranos Consortium on Regional Climatology and Adaptation to Climate Change (Table 1). The 10 simulations account for 85% of the future climate variability projected for the study site as established by a cluster analysis carried out on the range of available regional climate model (RCM) scenarios. PET for the future scenarios were calculated using the simulated temperature data (as for the 1900–2010 period) and the Oudin et al. (2005) formula. P_net values were calculated using future time series of daily precipitation and PET data. For the groundwater modelling, described in further detail below, the winter recharge (derived from P_net values) is transferred to the spring.

Table 1. RCM runs used for this study (see Mearns et al. 2012 for model acronym details).

RCM	GCM	Member	Domain	Emissions scenario
CRCM4.2.3	CGCM3	5	AMNO	A2
CRCM4.2.3	CGCM3	2	AMNO	A2
CRCM4.2.3	ECHAM5	1	AMNO	A2
CRCM4.2.3	ECHAM5	2	AMNO	A2
CRCM4.2.3	Arpège UnifS2	–	AMNO	A1B
CRCM4.2.0	CGCM3	4	AMNO	A2
HRM3	HADCM3	–	QC	A2
CRCM	CCSM	–	N. Amer.	A2
ECP2	GFDL	–	N. Amer.	A2
RCM3	CGCM3	–	N. Amer.	A2

Notes:

RCM, Regional Climate Model; GCM, General Circulation Model; CRCM, Canadian Regional Climate Model version; CGCM3, Coupled Global Climate Model version 3; ECHAM5, Max Planck Institute for Meteorology Global Climate Model version 5; AMNO, North America; Arpège UnifS2, France National Centre for Meteorological Research Global Circulation Model; HRM3, Hadley Center Regional Model version 3; HADCM3, Hadley Center Coupled Model version 3; QC, Quebec; CCSM, Community System Climate Model; N. Amer., North America; ECP2, Experimental Climate Prediction Center Regional Spectral Model; GFDL, Geophysical Fluid Dynamics Laboratory; RCM3, Regional Climate Model version 3.

Groundwater flow model development

A groundwater flow model was constructed using MODFLOW (Harbaugh 2005), assuming that the sandstone aquifer behaves as an equivalent porous medium, to simulate groundwater flow in the fractured bedrock and interactions between the aquifer, peatland and streams (see Levison et al. 2014b). In the current study, the groundwater model has been used to: (1) simulate hydrogeological conditions from 1900 to 2010 using the historical VI and temperature data, and (2) quantify water dynamics in response to climate change using 10 RCM scenarios for 2041–2070.

The model domain is divided into eight recharge zones and four hydraulic conductivity zones (Figure 2). The model is 96 m thick and is discretized in 16 layers. The layer thickness increases from 0.25 m at the surface to 30 m at the base. There are 9698 cells of 135 × 135 m, which are refined to 67.5 m over and around the peatland. The model is 173 km² and covers the area of the Natural Laboratory. It is limited to the north by the Noire River and to the northwest by the Outardes River. The Noire River is represented as a mixed boundary condition (River package in MODFLOW) and is located in the top two layers. Along the northern boundary of the model, this allows some of the groundwater to flow out of the aquifer as river baseflow. Below this, in the deeper layers, the boundary condition is set as a specified head which is considered to be groundwater flow to the regional aquifer. The eastern no-flow boundary is a flow line based on the piezometric map, located approximately 9 km parallel to and east of the Allen River. No-flow boundaries are also used for the southern and southwestern limit, set on the watershed limits of the Allen and Outardes Rivers. The base of the model is a no-flow boundary (see Levison et al. 2014b for a three-dimensional illustration of the model domain). Calibration targets were the available head measurements (bedrock observation wells, SIH database, private monitoring wells) and the baseflows estimated for the three gauging stations (Allen and Outardes Rivers and Schulman stream). Details on the calibration method and results are presented in Levison et al. (2014b).

Figure 2. The groundwater model domain with recharge and hydraulic conductivity zones (from Levison et al. 2014b). Recharge zone 1 covers the area of the peatland.

The Allen River and the Outardes River (both gauged) are represented using MODFLOW's River package in the top two layers. The Schulman stream (gauged) and other unnamed small permanent streams and tributaries are represented using MODFLOW's Drain package in the top two layers. The simulated springs are each represented using a single Drain node set at the elevation of the model surface. As such, water only flows from the springs when the aquifer head reaches the elevation of the spring. Springs (Drain nodes) were created at the locations of the five instrumented springs, as well as at four other locations at similar elevations where aquifer heads reach the ground surface in the model. These nine springs were simulated to determine flow properties, although it can be expected that there are additional unidentified springs within the large study area.

For the 1900–2010 period, transient-state simulations were performed using a seasonal stress period (91 days) to simulate varying flow conditions, including groundwater discharge to rivers and to springs. The spin-up period length for the 110-year analysis was 20 years, using the average P_net values for the first 20 years of the century. Monthly stress periods (30.4 days) were used to simulate groundwater discharge to rivers and to springs from climate model projections for the period 2041–2070. The results were compared to the 1971–2000 reference period. A 20-year spin-up period was also used for these simulations. A longer stress period for the 1900–2010 simulations was chosen to reduce model calculation time.

The groundwater recharge was first calibrated for the steady-state model and was therefore inferred from the observed hydrological data (hydraulic heads and baseflows). For the 1900–2010 period, multipliers were used to modify the value of calibrated recharge for each transient stress period, using the ratio of the P_net for the stress period of interest to the average P_net for the calibrated period, as described in Levison et al. (2014b). This reduces potential evapotranspiration to actual evapotranspiration and takes into consideration runoff. This method assumes that seasonal (or monthly) recharge is distributed analogously to the P_net ratios. It is assumed that a lower net precipitation will lead to a reduction in both runoff and infiltration. This might not hold true if rainfall intensity increases, for example. More detailed approaches exist to simulate long-term recharge values, such as using a physically-based recharge model (Rivard et al. 2014), or a daily weather generator coupled with a recharge model where actual evapotranspiration is simulated based on vegetation characteristics and soil water content (cf. Herrera-Pantoja and Hiscock 2008; Herrera-Pantoja et al. 2012). However, such approaches are more labour intensive and beyond the scope of this work.

For the 1900–2010 simulations (data from VIs and temperature), the multipliers were calculated for each 91-day season using the ratio of the P_net for the season of interest to the average P_net for the calibrated steady-state period, following Levison et al. (2014b). For the 2041–2070 (data from climate models) simulations, these multipliers were calculated for each 30.4-day month. In the model, the winter recharge (or recharge for below-zero temperature months) is set to zero and transferred to the spring (or first month with an average temperature above 0°C). For durations where the P_net is negative (i.e. for some summers or summer months), the recharge in the model is also set to zero.

Results and discussion

Comparison of vertical inflows data with measured data (1961–2010)

As noted previously, a comparative analysis was conducted using this data set and the meteorological data recorded at the Hemmingford weather station from the onset of the records in 1961 until 2010. There is a good fit between the Hemmingford (HMG) weather station data and the gridded temperature data (Centre d'expertise hydrique du Québec – CEHQ), even though the gridded data slightly overestimate the station data (Figure 3a). Annual average temperature is comparable for both data sets (HMG = 6.4°C; CEHQ = 6.7°C). Interannual averages of yearly minimum and maximum temperatures from 1961 to 2010 are also similar, with no statistical difference between data sets according to a Wilcoxon–Mann–Whitney test (minimum: HMG = –23.8°C; CEHQ = –23.4°C; maximum: HMG = 30.3°C; CEHQ = 30.8°C). The coldest month of the year is January (HMG = –9.6°C; CEHQ = –9.3°C) and the warmest is July (HMG = 20.6°C; CEHQ = 20.9°C). With both the HMG and CEHQ temperature data, the largest variability occurs during the winter months (December, January and February). The greatest discrepancy between average monthly data is 0.4°C for the months of June, July and September. However, there is no statistically significant difference between the gridded and station data for any month of the year (Wilcoxon–Mann–Whitney test; α = 0.05).

Figure 3. Comparison of Poirier et al. (2014) data with measured data from the Hemmingford (HMG) weather station from 1961 to 2010: (a) monthly temperature, and (b) monthly precipitation. Error bars represent minimum and maximum monthly values over the 1961–2010 period. CEHQ = Centre d'expertise hydrique du Québec.

On a yearly basis, the gridded VI data compare well to the weather station precipitation data. Interannual averages of yearly total precipitation/VI over the 1961–2010 period are similar (HMG = 926 mm; CEHQ = 933 mm) and there is no statistical difference between data sets (Wilcoxon–Mann–Whitney test; α = 0.05). However, on a monthly basis, the gridded VI data show marked differences from the measured precipitation data at the Hemmingford weather station (Figure 3b). This is due to the fact that the VI and weather station data do not represent the same phenomena during the winter months. In the winter months, there will only be VIs if there is water available for recharge or runoff (i.e. if precipitation falls as rain due to above-freezing temperatures and/or if there is melting of the snowpack). VIs do not account for the solid precipitation stored in the snowpack until temperatures rise above freezing, during the spring snowmelt. Winter precipitation is in this way transferred to spring months. March and April VIs should therefore roughly equal the sum of December, January and February precipitation, minus the portion of the precipitation that would have fallen as rain and the portion of the snowpack that would have melted on days with above-freezing temperatures.

It should thus be expected that VIs will be generally lower than station precipitation for the months of December, January and February, and considerably higher for the months of March and April. As a corollary, both data sets should be comparable from May through to November. The greatest HMG precipitation is measured in July, with an average monthly precipitation of 99 mm. In contrast, the highest VIs are in April, when snowmelt is most important, with available water amounts ranging from 24 to 399 mm. Statistically significant differences between data sets are found from December to April, while no significant differences are found for the other months of the year (Wilcoxon–Mann–Whitney test; α = 0.05).

Apart from the differences expected for the winter months (where vertical inflows are lower than precipitation because of storage in snow cover) and spring months (where the vertical inflows are greater than precipitation due to the melting of the snow cover), the correspondence between the gridded data and the weather station data was judged satisfactory, thus justifying the use of the vertical inflows data for the P_net calculation used to model groundwater flow over the 1900–2010 period. It is noteworthy that although this may be true for the 1961–2010 period, it may not hold true for the first half of the century. During the 1900–1960 period, the spatial interpolation on which the VIs are based was done using a smaller number of weather stations. The VI data thus have a lower reliability for the first half of the century due to low weather station density (Poirier et al. 2014). The quality of the VI data increases markedly after 1950, with the increasing number of weather stations.

1900–2010 simulated results

The VI and temperature data sets supplied by the CEHQ for the 1900–2010 period were used to calculate P_net values, i.e. the water available for surface runoff and groundwater recharge (Figure 4). The VI time series does not show any trend during the 1900–2010 period series (Mann–Kendall test; p = 0.783). Although significant trends in precipitation data were observed for certain parts of Canada, notably in the northernmost regions, similar trends were not detected in southern Quebec (Zhang et al. 2000). The temperature time series shows a statistically significant increasing trend (Mann–Kendall test; p < 0.001). This is consistent with trends observed across Canada over the twentieth century (Zhang et al. 2000; Nalley et al. 2013). A statistically significant increasing trend (Mann–Kendall test; p < 0.001) is also detectable in the PET data set over the 1900–2010 period. Because the PET interannual fluctuations are low in contrast to the fluctuations of the VI, the P_net time series evolves in a similar manner to that of the VI time series. P_net varies greatly from year to year, ranging from a minimum of 208 mm (1964) to a maximum of 754 mm (1918). The trend in temperatures produces a P_net time series with a slightly decreasing statistically significant trend (Mann–Kendall test; p = 0.011). When looking at the 10-year moving average, short-period trends become apparent. P_net generally decreases from 1900 to 1950 and a sharper decline is observed from 1950 to 1965. An increase occurs between 1965 and 1975, followed by a decrease from 1975 to 1985, and there is a slight increase after 1985. Although they could be due to the inherent temporal variability and lower reliability of the kriged VI and temperature data, especially from 1900 to 1965, these “intra-century” trends are similar to those presented in other studies for southern Canada (Zhang et al. 2000; Hamilton et al. 2001).

Figure 4. Annual vertical inflow (VI), evapotranspiration (PET), net precipitation (P_net), P_net 10-year moving average and temperature (T) from 1900 to 2010.

Figure 5 shows annual recharge (used in the MODFLOW model), average simulated hydraulic heads for the entire model domain, average simulated spring flow rates, and average simulated groundwater discharge to rivers (baseflows) from 1900 to 2010. The results for these four variables show similar trends. Until 1950, the results indicate strong temporal variability. This may be explained by the relatively low reliability ​​of the VI data set before 1950 (Poirier et al. 2014).

Figure 5. Annual recharge from 1900 to 2010 compared to (a) average simulated hydraulic heads for the study area, (b) average simulated spring flow and (c) average simulated baseflows for the Allen River, the Outardes River and the Schulman stream.

Average annual recharge varies from 41 to 197 mm (Figure 5). Recharge values fluctuate significantly from year to year, but the range includes the average annual recharge of 113 m calibrated by Levison et al. (2014b) for the recent past. Recharge generally decreases from 1900 to 1950, with a sharper decline from 1950 to 1965. Recharge then increases over the 1965–1975 decade, before decreasing again over the next decade (1975 to 1985) to finally stabilise after 1985. The downward trend over the century is significantly different from zero according to the Mann–Kendall trend test (α = 0.05). These results are the same as those observed above for the P_net time series simply because recharge is derived from P_net values (see the Methods and approach section). When looking more closely at the intra-century fluctuations it can be seen that the 1900–1965 downward trend is significant, but there is no trend statistically different from zero for the second portion of the century (1966–2010; Table 2).

Table 2. Long-term trends for calculated recharge and for simulated heads, spring flows and river/stream baseflows.

	Trends 1900–2010
	1900–1965	1966–2010
Recharge (m)	↓^*
	↓^*	null
Heads (m)	↓^*
	↓^*	↑
Spring flow (m^3/s)	↓^*
	↓^*	↑
Baseflow – Allen River (m^3/s)	↓^*
	↓^*	↑
Baseflow – Outardes River (m^3/s)	↓^*
	↓^*	↑
Baseflow – Schulman stream (m^3/s)	↓^*
	↓^*	↑

^* Statistically significant results.

As expected, the trends for average simulated heads are similar to those for recharge. The average head over the 1900–2010 period is 176.3 m, ranging from a minimum of 172.8 m to a maximum of 179.0 m. Average heads decrease significantly over the 1900–2010 period and also for the 1900–1965 period. However, average heads increase, although not statistically significantly, from 1966 to 2010 (Table 2). In Levison et al. (2014b), the average head over the model domain for the spring season in 2010, for example, was 178.1 m. Simulated recharge changes of –30%, –15% and +10% induced average head changes of –4.2 m, –1.9 m and +1.1 m, respectively (Levison et al. 2014b).

The simulated average outflow to springs varies from 0.00013 to 0.00195 m³/s, which is greater than the 0.00009 m³/s value estimated in the field in May 2010 [it is assumed that the flow rates during snowmelt are higher than this; see Figure 5b, and Levison et al. (2014a)]. In the absence of reliable flow records, these results are considered reasonable. For the VI data set, 2010 ranks in the 75th percentile for wettest years. However, only considering the VIs for March and April (snowmelt periods), 2010 ranks in the 25th percentile. Thus, 2010 can be considered a “wetter” year overall, but it did not have a very snowy winter nor important snowmelt compared to other years. The low recharge period from 1950 to 1965 induces a marked effect on spring flows, which decrease significantly. The results of the trend tests are identical to those of the average heads, i.e. there is a significant downward trend over the entire century and over the 1900–1965 period, and a non-significant upward trend from 1966 to 2010 (Table 2). The individual springs were dry for a minimum of 3.8% to a maximum of 89.4% of the simulated seasons. Sustained spring flows are ecologically important since they provide habitats for many species, including the Allegheny Mountain Dusky Salamanders.

The average annual baseflows simulated for the Allen River range from 0.031 to 0.130 m³/s. For the Outardes River they range from 0.027 to 0.071 m³/s, while for the Schulman stream they range from 0.005 to 0.013 m³/s (Figure 5c). The order of magnitude and the range of baseflows are considered to be correctly simulated. This is assessed considering that the baseflow separation of field data reported in Levison et al. (2014b) captures peak baseflows (e.g. 0.40 m³/s) while the simulated data (1900–2010) are based on seasonal stress periods (that is, an average value of recharge is applied for the entire season, which has the effect of dampening the simulated peak baseflow). Field-derived flow rates for the recent past (2007–2010) are 0.05 to 0.40 m³/s (Allen River), 0.02 to 0.40 m³/s (Outardes River) and 0.001 to 0.024 m³/s (Schulman stream) determined by baseflow separation (cf. Levison et al. 2014b). Similarly to the spring flows, the low recharge period from 1950 to 1965 induces a marked effect on baseflows, which decrease significantly. The results of the trend tests are comparable to those for average heads and spring flows (Table 2). It is important to note that none of the three rivers dry out during the 110 years simulated. These results are particularly important as the Covey Hill rivers are primary habitats for the Northern Spring Salamander and many other species.

Comparison with climate change scenarios

The transient-state model was used to simulate hydraulic heads as well as spring flows and river baseflows under 10 climate scenarios for the baseline (or reference) period (1971–2000) and the 2050 horizon (2041–2070 period). The results were compared to determine a statistically significant difference between these two periods. A significant difference is interpreted as a possible impact of climate change on a variable. Results are summarized in Table 3. The climate scenario simulation results were also compared with the simulated results obtained using observed precipitation and temperature from the HMG weather station for 1971–2000, as well as with the simulation results obtained with the CEHQ gridded data VI and temperature data from 1900 to 2010.

Table 3. Groundwater flow simulations with climate change scenarios. Future projections for calculated recharge and for simulated heads, spring flows and river/stream baseflows.

	Climate scenario
	REF	FUT
Average annual recharge (mm)	128	136^*
Average heads for the entire study domain (m)	176.5	177.3^*
Average spring flow (m^3/s)	0.0011	0.0013^*
Average baseflow – Allen River (m^3/s)	0.094	0.103^*
Average baseflow – Outardes River (m^3/s)	0.074	0.083^*
Average baseflow – Schulman Stream (m^3/s)	0.014	0.016^*

Notes:

REF, reference period; FUT, future period.

^* Statistically significant results for the ensemble mean.

The average annual recharge estimated from the HMG data for the reference period is 124 mm (0.124 m/year; Figure 6a). The average recharge obtained for the 1900–2010 period is slightly lower (117 mm). The simulated recharge variability over the past century is greater than that of the observed period (1971–2000; period when observed climate data is available from the HMG weather station), ranging from a minimum of 41 mm to a maximum of 198 mm, with values from the observed period varying from 69 mm/year to 170 mm/year. For the reference period, climate simulations satisfactorily reproduce observed period data (average of 128 mm), but their dispersion is slightly greater than that of the observed period values. The spread of the reference climate simulations is also greater than that from 1900 to 2010. A majority of the future climate simulations (seven out of 10) predict an increase in recharge when compared to the reference simulations, but these differences are not statistically significant. However, the ensemble mean indicates a statistically significant increase for the 2050 horizon (Wilcoxon–Mann–Whitney test; α = 0.05).

Figure 6. Comparison with climate change conditions for: (a) recharge, (b) average hydraulic heads, (c) total spring outflows, (d) baseflows to the Allen River, (e) baseflows to the Outardes River and (f) baseflows to the Schulman stream. The significant differences (Wilcoxon–Mann–Whitney test; α = 0.05) between the reference and future periods are marked with *. OBSERVED is the simulation using the observed climate data, ENS MEAN is the ensemble mean, REF is the reference period and FUT is the future period. The 1st, 25th 50th 75th and 99th percentiles are plotted.

The average hydraulic head for the 1971–2000 period using the HMG data is 176.6 m (Figure 6b). Although the difference is not significant, this value is slightly higher than the average value of 176.3 m over the last century. As might be expected, the variability over the century is greater than that simulated using the observed data from 1971–2000. The variability of hydraulic heads over the past century is also greater than that of reference climate simulations, although the latter are able to reproduce average heads with an ensemble mean of 176.6 m. Future heads from all climate simulations are fairly variable (ranging from an average value of 176.1 m for the CRCM_CCSM model to an average value of 178.5 m for the CRCM4.2.3_CGCM3#5 model). The ensemble average projects an increase of 0.78 m. The difference between the reference and future period is statistically significant (Wilcoxon–Mann–Whitney test; α = 0.05) for nine out of 10 models.

During the 1971–2000 period, results obtained with the HMG data indicate an average spring flow rate of 1.08 × 10⁻³ m³/s, ranging from a minimum of 0.67 × 10⁻³ m³/s to a maximum of 1.51 × 10 ⁻³ m³/s (Figure 6c). The 0.96 × 10⁻³ m³/s average value of spring flows over the 1900–2010 period is slightly lower than the spring flows simulated with the HMG data, but the interannual variability for this period is greater, ranging from 0.13 × 10⁻³ to 1.95 × 10⁻³ m³/s. Climate simulations for the reference climate scenarios satisfactorily reproduce the mean and interannual variability observed between 1971–2000, with an ensemble average of 1.09 × 10⁻³ m³/s for this period. All but two of the climate models show increased spring flow for the future. The ensemble average projects an increase of 0.23 × 10⁻³ m³/s, or 21% in spring flow. This could be beneficial to the spring-dwelling Mountain Dusky Salamander, as well as other species that depend on these springs. The differences between the reference and future periods are statistically significant (Wilcoxon–Mann–Whitney; α = 0.05) for all models except CRCM4.2.3_ECHAM#2 and HRM3_HADCM3.

Baseflows for the three gauged rivers show similar results to the spring flows (Figure 6d–f). Results for the 1971–2000 period obtained with the HMG data indicate average values of 0.092, 0.073 and 0.014 m³/s for the Allen River, the Outardes River and the Schulman stream, respectively. The average baseflows over the 1900–2010 period are slightly lower than these values (Allen: 0.087 m³/s; Outardes: 0.071 m³/s; Schulman: 0.013 m³/s). In all three cases, the interannual variability over the past century is greater than the variability between 1971 and 2000. Climate simulations adequately reproduce both the central tendency and the annual variability during the reference period, with ensemble averages of 0.094, 0.074 and 0.014 m³/s for the Allen River, the Outardes River and the Schulman stream, respectively. Future climate simulations project increased baseflows for the three streams (Allen: 0.008 m³/s or 8%; Outardes: 0.008 m³/s or 11%; Schulman: 0.002 m³/s or 14%). This increase in baseflows could be potentially beneficial to stream salamanders’ habitat. The differences between the reference and future periods are statistically significant (Wilcoxon–Mann–Whitney; α = 0.05), with the exception of models CRCM4.2.3_ECHAM#2, CRCM_CCSM and HRM3_HADCM3.

In general, reference climate scenarios adequately reproduce the variability of observed climate data for the 1971–2000 period with respect to average heads, spring flows and baseflows for the two rivers and one stream studied. The ability of the climate models to reproduce observed conditions could increase the level of confidence in the ability of future projections to represent actual future conditions for these variables. Because it spans a longer time frame, the spread of values for the 1900–2010 data set is greater than that for the reference period (either with HMG data or with climate scenarios), as may be expected. Moreover, most 1971–2000 simulation results are found in the upper portion of the 1900–2010 data range while the lowest – or driest – values simulated over the century are underrepresented in the reference simulations. However, 1900–2010 data before 1950 are considered less reliable due to low weather station density used for data interpolation.

The confirmation that the Covey Hill natural system has thus been subjected to a wide range of climatic conditions over the last century may represent good news in terms of the system’s resilience, or its capacity to adapt to stress and change. Indeed, the system has shown little change despite the more variable (and sometimes drier) conditions than what has been observed in the last 30 to 40 years. The fact that most future simulations are in the upper portion of the range (or even beyond the range) of values ​​simulated over the past century should not a priori represent a concern for the Covey Hill natural system, from a hydrological standpoint. However, changes in socio-economic factors and/or land-use (not investigated in this research) could also contribute to altering the system.

Conclusions

In order to gain a better understanding of hydrological responses to future climate forcing, this study compared a range of hydrological reactions to past century-long and future climate scenarios at the Covey Hill Natural Laboratory in southern Quebec, where a groundwater-surface water interacting system supports a fragile ecosystem. The results provide a broader view about the hydrological conditions that an endangered species (salamanders) has survived (in the past) and to which it might be exposed (in the future). Specifically, this study determined historical trends, using a 1900 to 2010 vertical inflow data set and temperature data as input to a transient-state groundwater flow model. The simulated historical trends were compared with flow projections for the 2041–2070 period.

The modeling of the groundwater system over the last century highlighted the responses to the intrinsic climate variability during this period. The P_net analysis showed a decreasing trend of potential recharge from 1900 to 2010, which impacted the simulated groundwater levels. However, this trend tends to reverse itself in recent decades. This suggests that the increase in precipitation projected for 2050 would in effect offset the increase in temperature (and evapotranspiration) and therefore may not cause significant decreases in available water for recharge, and consequently groundwater levels on Covey Hill. The hydrosystem appears to be highly resilient and climate might not negatively impact groundwater-sustained habitats on Covey Hill in the near future. These a priori positive results are important, not because they limit the need to protect the land, but because they reduce the uncertainty on possible future conditions. The past simulated conditions vary due to uncertainty in the VI data, but this variability is not as great as the variability for the future simulated conditions. Therefore, the impact of the VI uncertainty (especially before 1950) is considered to not be a major factor.

Many research studies have now been conducted that investigate climate change-induced trends in groundwater levels and groundwater–surface water interactions, through numerical modeling of future scenarios compared with recent past conditions. Findings from this study confirm the importance of maintaining the long-term operation of ecosystem monitoring sites, and of performing paleohydrological investigations, to understand the range of climatic conditions to which a natural system is adapted and to better constrain future predictions and expected extremes based on already encountered conditions.

Acknowledgements

This project was funded by the Ouranos consortium on regional climatology and adaptation to climate change, as part of the “Fonds vert” for the implementation of the Quebec Government Action Plan 2006–2012 on climate change and its measure 26 (Grant #554007-107). The authors would like to thank the Nature Conservancy of Canada and the Covey Hill landowners for providing access to properties on Covey Hill.