Uncertainty assessments and hydrological implications of climate change in two adjacent agricultural catchments of a rapidly urbanizing watershed

Science of the Total Environment 473–474 (2014) 326–337

Contents lists available at ScienceDirect

Science of the Total Environment journal homepage: www.elsevier.com/locate/scitotenv

Uncertainty assessments and hydrological implications of climate change in two adjacent agricultural catchments of a rapidly urbanizing watershed S.K. Oni a,b,⁎, M.N. Futter b,c, L.A. Molot d, P.J. Dillon c, J. Crossman c a

Environmental and Life Sciences Graduate Program, Trent University, 1600 West Bank Drive, Peterborough, ON K9J 7B8, Canada Department of Aquatic Sciences and Assessment, Swedish University of Agricultural Sciences, Uppsala, Sweden c Environmental and Resource Studies, Trent University, 1600 West Bank Drive, Peterborough, ON K9J 7B8, Canada d Faculty of Environmental Studies, York University, 4700 Keele Street, Toronto, ON M3J 1P3, Canada b

H I G H L I G H T S • • • • •

Predictive uncertainty in runoff is often larger in managed than pristine catchments. CGCM3 projected warmer (4 °C) and longer growing season (26%) in the future. Behavioral, not optimum parameter sets should be considered in impacted watersheds. Model results show adjacent catchments differ in their snow and groundwater dynamics. Human activities exacerbate the differences in integrated hydrological responses.

a r t i c l e

i n f o

Article history: Received 20 October 2013 Received in revised form 6 December 2013 Accepted 7 December 2013 Available online 27 December 2013 Keywords: Statistical downscaling Phenology Climate change Lake Simcoe Hydrological change Uncertainty analysis

a b s t r a c t Lake Simcoe is the most important inland lake in Southern Ontario. The watershed is predominantly agricultural and under increasing pressure from urbanization, leading to changing runoff patterns in rivers draining to the lake. Uncertainties in rainfall–runoff modeling in tributary catchments of the Lake Simcoe Watershed (LSW) can be an order of magnitude larger than pristine watersheds, hampering water quality predictions and export calculations. Here we conduct a robust assessment to constrain the uncertainty in hydrological simulations and projections in the LSW using two representative adjacent agricultural catchments. Downscaled CGCM 3 projections using A1B and A2 emission scenarios projected increases of 4 °C in air temperature and a 26% longer growing season. The fraction of precipitation falling as snow will decrease. Spring runoff is an important event in LSW but individual HBV best calibrated parameter sets under-predicted peak flows by up to 32%. Using an ensemble of behavioral parameter sets achieved credible representations of present day hydrology and constrained uncertainties in future projections. Parameter uncertainty analysis showed that the catchments differ in terms of their snow accumulation/melt and groundwater dynamics. Human activities exacerbate the differences in hydrological response. Model parameterization in one catchment could not generate credible hydrological simulations in the other. We cautioned against extrapolating results from monitored to ungauged catchments in managed watersheds like the LSW. © 2013 Elsevier B.V. All rights reserved.

1. Introduction Long term monitoring and observations provide increasing evidence of climate change signals on ecosystem functioning (Campbell et al., 2005; Euskirchen et al., 2007; Menzel and Fabian, 1999; Tetzlaff et al., 2013; Oni et al., 2013a; Vincent and Mekis, 2006; Zhang et al., 2000). These trends are likely to continue as global climate models (GCM) project a ⁎ Corresponding author at: Department of Aquatic Sciences and Assessment, Swedish University of Agricultural Sciences, Uppsala, Sweden. E-mail address: [email protected] (S.K. Oni). 0048-9697/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.scitotenv.2013.12.032

warmer future with intensification of hydroclimatic cycles (Huntington, 2006). While there are uncertainties surrounding quantification of the magnitude of these changes, the evidence for climate-related changes is robust enough to devise mitigation and adaptation plans for the future. Climate drives the hydrological cycle; it is a first order control which drives other biogeochemical processes (Laudon et al., 2013; Oni et al., 2013a). Climate change has led to greater warming of the earth and could impact hydrological processes at both spatial and temporal scales. Possible impacts include changes in runoff duration and distributions, as well as magnitude and frequency of extreme events (Arnell, 2004; Barnett et al., 2005; Cunderlik and Simonovic, 2005). The impact of

S.K. Oni et al. / Science of the Total Environment 473–474 (2014) 326–337

327

climate change would not be uniform throughout the world. For example, Canada has experienced a greater warming of approximately 0.3 °C (Vincent and Mekis, 2006; Zhang et al., 2000) above the global average of 0.7 °C reported by the International Panel on Climate Change (IPCC) by the end of the 20th century (IPCC, 2007). One possible explanation that has been offered is the high latitude location of Canada (e.g. Tetzlaff et al., 2013; Laudon et al., 2013). While efforts to better understand the potential impact of climate change on local watershed hydrology are ongoing, significant gaps in knowledge still exist as responses vary from region to region. This is due to several factors including differences in local runoff generation processes and landscape forms as well as varying degrees of landscape disturbances by human activities within watersheds (Franczyk and Chang, 2009; LSRCA, 2008; O'Connor et al., 2011). While scenarioneutral probabilistic ensembles are becoming more common in climate impact studies (e.g. Harris et al., 2010), the use of scenario-driven GCMs are still the most widely used approach to assess the plausible impacts of future climate on watershed hydrology and ecosystem functioning (Arnell, 2004; Oni et al., 2012a,b; Wilby and Dawson, 2013). However, there are a number of uncertainties associated with the GCM projections and downscaling to local conditions under IPCC emission scenarios (Wilby et al., 2002). Our preparation for future hydrological changes appears to be insufficient due to the larger uncertainty in our present hydrological model structures and forecasting capabilities (Allen and Ingram, 2002; Juston et al., 2013). This is more pronounced in snow-dominated catchments (Barnett et al., 2005; Lawrence and Haddeland, 2011; Oni et al., 2012b) and can be amplified by changing land uses. For example, urbanizing watersheds are characterized by rapid land use changes (Franczyk and Chang, 2009; Furusho et al., 2013) and associated landscape disturbances can shift the rainfall–runoff relationships away from natural processes upon which the conceptualization of most hydrological models are based (e.g. Seibert and McDonnell, 2010). As a result, the uncertainty in rainfall–runoff modeling in managed watersheds can be an order of magnitude larger than pristine watersheds or those propagated by GCMs. A continuous assessment of watershed hydrological responses is therefore necessary to constrain all sources of uncertainty in hydrological simulations before any credible future projections can be made in managed, rapidly urbanizing, and snowdominated watersheds such as Lake Simcoe. The Lake Simcoe watershed (LSW) is the most important inland lake watershed in Southern Ontario for socioeconomic reasons (agriculture, recreation and tourism, drinking water supply and angling activities etc.). However, the LSW has been under pressure from urban developments and intensive agriculture in the recent decades. This has led to changes in land use and changing rainfall–runoff relationships across the watershed (Oni et al., 2013b). Here we conducted a robust assessment to constrain the uncertainty in hydrological simulations and climate projections in two adjacent river catchments of LSW (Beaver and Pefferlaw Rivers) where impacts of urbanization are limited but where there have been large impacts from agriculture. It is often assumed in hydrological modeling that adjacent catchments with similar characteristics could have similar runoff generation processes and thus parameterization in one should generate credible runoff conditions in the other (Patil and Stieglitz, 2013; Oni et al., 2012a). We tested this hypothesis and the suitability of this approach of using calibrated parameter sets from one catchment to simulate hydrological processes in adjacent river catchments.

the southern end of the catchment (Fig. 1). There have been increasing pressures from urbanizations in some parts of the LSW (LSRCA, 2008; Oni, 2011), leading to a changing runoff ratio (Oni et al., 2013b). The Beaver and Pefferlaw River catchments were used in this study. They are adjacent catchments with characteristic branched headwater structures where streams of lower order drained into higher order streams and rivers that flow northward into the lake (Fig. 1). Both catchments have limited urban development and are largely impacted by mixed agricultural activities that include both intensive and non-intensive agriculture (Jin et al., 2013). The Beaver River drains ~282 km2 while the Pefferlaw River drains ~ 332 km2. The Pefferlaw River in this study is the combination of Pefferlaw Brook and Uxbridge Brook (Fig. 1). Both catchments have similar weather and land use patterns (Oni et al., 2013b) and share some similar physiographic features (LSRCA, 2012a,b). The headwaters of the Pefferlaw River (and some part of Beaver River) originate in the Oak Ridges Moraines (ORM) located at the southern part of the LSW (Fig. 1). The ORM was formed at the end of the last glaciation and the high infiltration capacity and total lack of surficial drainage are common to moraine features (Johnson, 1997). The moraine consists of an aquifer complex (confined and unconfined aquifers) that makes it an important recharge zone that contributes baseflow to the headwaters and provides a reliable supply of drinking water (LSRCA, 2012a,b). The peak of the moraine (Uxbridge wedge) is located in the Pefferlaw River catchment, making a steep elevation gradient of 397–220 m above sea level (masl) from ORM to the Lake Simcoe shoreline. In contrast, the Beaver River watershed has a low elevational gradient (LSRCA, 2012b). There are scattered rural communities in the central to downstream regions of the Beaver River watershed. Rural/urban developments in the Pefferlaw River watershed are concentrated around the headwater areas in close proximity to ORM. Both catchments have similar bedrock geology (soft limestone and shale bedrock) that are covered by thick fertile soils (Jin et al., 2013). However, soils in the Beaver River catchment are Brunisols (formed under imperfectly drained conditions) with some clay loam, sandy loam and organic muck in the headwaters (Johnson, 1997). Soils in the Pefferlaw River catchment are Gray Brown Luvisols that thin out toward the north, thereby limiting groundwater recharge in the central part of the catchment (LSRCA, 2012a). Downstream reaches of the Beaver River drain an extensive wetland complex (19% of the total catchment area) that is concentrated in the center of the catchment and forms a large, well defined recharge zone. In contrast, wetlands in Pefferlaw cover 16–18% of the total catchment area but the wetland system is more patchy and scattered (Jin et al., 2013). This makes the influence of wetland hydrology larger in the Beaver than Pefferlaw River. In addition to the wetlands, a large area of the Beaver River channel is underlain by organic sediment (Baulch et al., 2013). Both catchments have flow and weather monitoring stations within their catchment boundaries. Flow and weather monitoring stations at Udora in the Pefferlaw River catchment and Woodville in the Beaver River catchment were used in the hydrologic modeling and impact analysis presented herein (Fig. 1). Both the Udora (44° 15′ 45″ N, 79° 09′41″ W) and Woodville (44°24′ N, 78°58 W) stations are monitored by Environment Canada (EC). Available data included runoff (m3 s− 1), precipitation (mm d−1) as well as minimum, average and maximum daily air temperature (°C).

2. Material and method

The Statistical Downscaling Model (SDSM) was used to downscale Canadian Global Climate Model 3 (CGCM3) predictors to LSW conditions. Downscaling is important in order to bridge the void between global predictors and local predictands at a watershed scale (Oni, 2011). Using forty years of baseline predictand variables (precipitation and air temperature) from local weather stations (1961–2000) and the equivalent predictor variables from National Center for Environmental

2.1. Study site description Lake Simcoe (722 km2) is a dimictic lake located in south-central Ontario (44° 25′ N, 79° 20′ W). The lake drains a large terrestrial system (~ 2899 km2) via its headwater tributaries, which originate mainly at

2.2. Climate downscaling

328

S.K. Oni et al. / Science of the Total Environment 473–474 (2014) 326–337

Fig. 1. Map of entire Lake Simcoe watershed and its headwater hydrologic system. Beaver and Pefferlaw River catchments are highlighted with Oak Ridges Moraines at the southern end of the watershed. Pefferlaw River is gauged below the confluence of the Pefferlaw and Uxbridge Brook.

Prediction (NCEP) reanalysis data, statistical relationships were established in SDSM following the procedure described in Wilby et al. (2002). The SDSM combines a regression based approach with a stochastic weather generator to project plausible future conditions. In the first stage, quality control checks were performed on the predictand, NCEP and CGCM3 predictor variables. This was followed by screening of the predictand variables for statistical relationships with the NCEP reanalysis data. The NCEP (and the future CGCM3 predictor) variables used in the downscaling temperature series include specific humidity at 500 hpa and 850 hpa as well as mean temperature at 2 m. The predictor variables used in SDSM to downscale precipitation include surface meriodional velocity, surface vorticity, geopotential height at 500 hpa and specific humidty. Using multiple stepwise regressions, the established statistical relationship between predictand and NCEP reanalysis data was applied to future CGCM3 predictor variables. Both correlation matrices and scatter plots were used as measures of calibration performance. Stochasticity was added to the downscale variables using the stochastic weather generator in SDSM. In downscaling precipitation, conditional process was used due to intermediate forcing between regional and local conditions but unconditional process was used for temperature downscaling. More information about SDSM downscaling is available in Wilby et al. (2002).

Future climate projections were made under IPCC A1B and A2 emission scenarios represented in CGCM3. Due to inconsistencies and unavailability of long term historical weather data at any stations located within the LSW, predictand series from the nearby ESSA Ontario Hydro station was used in the downscaling and supplemented by data from the Barrie Water Pollution Control Centre station. Both the NCEP reanalysis data and/or the future CGCM3 predictor variables were made available by the Data Access Integration (DAI) portal of EC. More information on statistical downscaling in LSW is presented in earlier studies (Crossman et al., 2013; Oni et al., 2012a). The downscaled climate series were used as input variables to drive a rainfall–runoff model. Future hydroclimatic changes were based on thirty year split periods of 2041–2070 and 2071–2100 relative to the baseline period (1961–2000). We evaluated the change in growing season under the present day (using the baseline data) and future conditions (using the projected temperature data). Growing degree days (GDD) were calculated using Eq. (1). GDD ¼ ½ðTmax þ Tmin Þ–Tbase

ð1Þ

where Tbase was a 5 °C threshold temperature, and Tmax and Tmin were the corresponding daily maximum and minimum air temperature (°C). Growing season length (days) was calculated as the number of

S.K. Oni et al. / Science of the Total Environment 473–474 (2014) 326–337

days with GDD greater than 5 °C. Also, possible change in rainfall and snowfall were estimated. Using data from November to April inclusive, precipitation was partitioned into either rainfall (T N 0 °C) or snowfall (T b 0 °C) to evaluate potential future changes in precipitation patterns.

2.3. Hydrologic modeling 2.3.1. Model description The HBV (Hydrologiska Byråns Vattenbalansavdelning) is a physically-based semi distributed rainfall–runoff model used for runoff simulations and forecasting (Lindström et al., 1997). It was developed originally at the Swedish Meteorological and Hydrological Institute (SMHI). HBV has been widely used in hydrological studies throughout the world and has therefore led to a plethora of versions to address different hydrological issues. These include flood forecasting and spillway design (Bergström et al., 1992), nutrient and carbon budgets (Crossman et al., 2013; Oni et al., 2011; Futter et al., 2007), climate change impacts (Lawrence and Haddeland, 2011; Oni et al., 2012a,b), and land use change (Seibert and McDonnell, 2010) among others. The HBV light version was in used in this study (Seibert and McDonnell, 2010). The principal conceptual framework of the HBV model is based on natural processes in snow dominated catchments (Fig. 2), making the model a useful tool for understanding and detecting when watershed hydrology deviates from natural rainfall–runoff processes (Zégre et al., 2010). The model also works well in detecting hydrological changes in tropical regions (e.g. Gebrehiwot et al., 2013). Since the model conceptualization was based on a snow-dominated environment, the model structure was based on routines of snow accumulation and melt, soil moisture accounting and runoff response modules (Lindström et al., 1997; Oni et al., 2012b).

329

2.3.2. Model variables and parameterizations The model uses a parameter threshold temperature, TT (0 °C) to set the boundary between rainfall and snowfall in the precipitation series (Fig. 2; model parameters are described in Table 1). Snow accumulation and melt were simulated in HBV using the degree-day method (Lindström et al., 1997) represented by the parameter CFMAX (mm day−1 °C−1). The parameter SFCF accounts for either correction of evaporation from snowpack that was not properly represented in the model or correction factor for snowfall measurements (Seibert and McDonnell, 2010). The soil moisture accounting routine controls the model water balance and can derive its inputs from either or both precipitation (rainfall or snowfall) and snowmelt from the snow routine (Fig. 2). The parameter CWH denotes a threshold at which the snowmelt or rain water is retained in the snowpack and possible refreezing of the waters in the snowpack with the aid of refreezing coefficient CFR. Field capacity, FC (mm) sets the maximum soil moisture storage capacity. The limit of evapotranspiration, LP (mm), sets the percent fraction of FC that determines the value at which evapotranspiration becomes potential evapotranspiration (Lindström et al., 1997). The parameter BETA, β (−) regulates the soil moisture routine's contribution to discharge from precipitation (either as snowmelt or rainfall) as the soil moisture condition approaches field capacity (Seibert and McDonnell, 2010). The parameter PERC (mm day−1) determines the maximum percolation rate of the infiltrated waters to the lower soil box (Fig. 2). The runoff response function K1 (day−1) and K2 (day−1) determine model runoff generation as well as partitioning of runoff as overland/intermediate and baseflow (Fig. 2). Soil moisture SM (mm) is an output variable that represents the net soil water balance in the upper soil box, excluding water loss to evapotranspiration and percolation. Hydrological effective rainfall, HER (mm) (Fig. 2) is also an output variable that can be described as the available

Fig. 2. Conceptual representation of rainfall–runoff processes in the HBV hydrological model, showing interdependent of climate and hydrological processes. Threshold Temperature denotes a temperature where incoming precipitation is partitioned to either rainfall or snowfall and can be influenced by climate variability and change. Land use changes impact water flowpaths and transport between terrestrial–aquatic continuum.

330

S.K. Oni et al. / Science of the Total Environment 473–474 (2014) 326–337

Table 1 Model parameter descriptions, ranges used in the Monte Carlo runs and their best performing parameter values in Beaver and Pefferlaw Rivers. KS denotes Kolmogorov–Smirnov Statistics and significance estimates for behavioral model parameters (best 200) in Beaver and Pefferlaw Rivers at p b 0.05 (italics) and p b 0.01 (bold). Routine

Snow

Soil

Response

Notation

TT CFMAX SFCF CFR CWH FC β Cet PERC α K1 K2 MAXBAS

Parameter description

Threshold temperature Degree day factor Snowfall correction factor Refreezing coefficient Snowmelt coefficient Field Capacity Soil curve parameter Evapotranspiration correction factor Percolation to lower soil box Alpha Quick/subsurface recession rate Baseflow recession rate Triangular weighting functions

Range

BestParSet

Min

Max

Beaver

Pefferlaw

−0.3 2.56 0.59 0.38 0.20 151 3.90 0.5 0.68 0.001 0.04 0.008 2

−0.5 4.99 0.99 0.62 0.36 250 6.84 0.8 0.97 0.004 0.18 0.028 3

−0.5 4.00 0.63 0.47 0.23 203 4.08 0.75 0.96 0.002 0.12 0.015 2

−0.5 3.28 0.87 0.42 0.21 244 4.79 0.69 0.94 0.002 0.16 0.009 3

soil water that can initiate runoff as it penetrates and recharges upper and lower soil boxes. Soil moisture deficit, SMD (mm), is an estimate of soil dryness derived from the difference of FC and antecedent SM (Fig. 2). These are important hydrological variables necessary for coupling of the HBV to other landscape scale biogeochemical models (e.g. Futter et al., 2007; Jin et al., 2013; Oni et al., 2011, 2012a). More details about model parameter descriptions and calibration values are summarized in Table 1.

KS

Units

0.080 0.745 0.905 0.100 0.100 0.167 0.355 0.090 0.110 0.090 0.595 0.420 0.380

°C mm day−1°C−1

mm °C−1 mm day−1 mm day−1 day−1 day

the CDF of the most sensitive behavioral parameter sets between the two catchments were evaluated using Kolmogorov–Smirnov (KS) statistics. The future behavioral runs were made by coupling the downscaled climate series to HBV to generate future projections under A1B and A2 emission scenarios. 3. Results 3.1. Climate and growing season

2.3.3. Model calibration To calibrate the HBV model, data from 1992 to 2000 hydrological years were used in both catchments for comparison as data were available during this period in both catchments. The hydrological year in this study represents September 1 in a year to August 31 of the following year. In the first stage of model calibration, HBV was manually calibrated to present day conditions (1992–2000) in the two adjacent catchments. Model calibrations were first performed by changing parameters in each routine to produce satisfactory hydrographs that match the observed runoff. HBV requires daily time series of air temperature T (°C) and precipitation P (mm) to generate corresponding runoff that can be compared to field observations during the calibration period. Model goodness-of-fit was assessed from a combination of visual examination of the simulated and the observed runoffs and Nash–Sutcliffe statistics (Nash and Sutcliffe, 1970). An uncertainty analysis in the form of a generalized sensitivity analysis (Spear and Hornberger, 1980) was conducted on the model parameters to determine the parameter sets that mostly represent the integrated catchment behavior. This was achieved by subjecting the model parameter space to a Monte Carlo sampling procedure. This was done in two iterations (fifty thousand times) in each of the catchments. During each run, parameters were sampled randomly from the prior uniform distributions within the confinement of the parameter ranges (Table 1). Uniform distributions have been reported to be appropriate in catchments where there is no adequate information on the internal catchment behavior (Zégre et al., 2010). Two hundred top performing parameter sets were selected based on Nash Sutcliffe statistics (NS) as the behavioral runs. The NS statistic (Nash and Sutcliffe, 1970) measures the qualitative performance of the model and has been used widely in hydrological modeling as an index to measure model goodness-of-fit in reproducing runoff conditions. The NS values range from 1 to minus infinity with a value of 1 indicating a perfect fit and a value of 0 or less signifying that the mean of the observed runoff is a better predictor than the model simulation (Krause et al., 2005). The behavioral parameter sets were subjected to further analysis by checking the deviation of their cumulative frequency distributions (CDF) from rectangular distribution in each catchment. Differences in

The analysis of long term historical weather over the forty year baseline period showed that there were significant monotonic increases in mean annual air temperature [Mann–Kendall (MK) = 3.15; P b 0.01] and precipitation [MK = 2.90; P b 0.01]. The long term mean annual temperature in the region was ~6.5 °C with about a 0.7 °C increase in temperature and the long-term precipitation increased by 6.3% over the baseline period. Both A1B and A2 scenarios projected that the mean air temperature would continue to increase throughout the 21st century. During July, minimal differences were demonstrated between the baseline and future temperatures. During winter, however, large increases in temperature were projected for the future A1B scenario (Fig. 3a). A similar pattern was observed under the A2 scenario with little increase in temperature in summer and largest warming in the winter months (Fig. 3b). In summary, both IPCC scenarios projected more extreme increases in winter temperatures by up to 4 °C, a shorter winter and a slightly extended spring season. There was no significant change in annual precipitation under either scenario due to high monthly variability. Summer reductions were observed during the 2050s and the 2070s, but during winter changes were much more variable (Fig. 3c, d). No major change was observed under the A2 scenario from January to July, with a lower amount projected in autumn (Fig. 3d). Analysis of long term historical temperature series showed that average growing season length in LSW was about 201 days. After subjecting the projected temperature series to similar analysis, both A1B and A2 emission scenarios projected significant increases in the growing season length (Fig. 4). For example, A1B projected up to 20% increase in growing season length while A2 projected 26% by the end of the 21st century (Fig. 4a). These changes resulted from increases in the growing degree day temperature in the spring and autumn/early winter months (Fig. 4b). However, little difference was observed in winter growing degree days (January to March) as winter is expected to be a dormant period. When the projected precipitation was partitioned into snowfall and rainfall (using a TT value of 0 °C), both scenarios projected a shift in precipitation toward increasing rainfall but a decline in snowfall relative to the baseline condition (Fig. 5). The interannual variations are larger in rainfall than snowfall projections.

S.K. Oni et al. / Science of the Total Environment 473–474 (2014) 326–337

331

Fig. 3. Seasonal change in the future temperature (a, b) and precipitation (c, d) under CGCM3 A1B and A2 emission scenarios.

3.2. Hydrological modeling 3.2.1. Model calibration The overall runoff dynamics in the catchments were controlled mostly by the sensitivity of snow and runoff routine parameters (Table 1). The cumulative distribution function (CDF) of parameter TT which partitions precipitation into either snowfall or rainfall was not particularly sensitive in either catchment. The CDF followed similar behavioral patterns which were not significantly different from one another (Table 1). However, CFMAX and SFCF were highly sensitive as their CDFs were significantly different (Fig. 6). The CFMAX tends towards maximum range in Beaver River (Fig. 6a) while the SFCF was skewed towards the maximum range in Pefferlaw River (Fig. 6b). While FC showed some weak sensitivity in both catchments (Table 1), β is the only parameter among the soil routine parameters that showed high sensitivity (Fig. 6c) and was significantly different between the two catchments (Table 1). The CDF of parameter β skewed toward the

Fig. 4. a) Long term patterns of growing season length from 1961 to 2100 and b) seasonal evaluation of growing degree days between baseline and CGCM3 future climate scenarios.

minimum range in Pefferlaw River and towards the maximum in Beaver River (Fig. 6c). Both the MAXBAS and recession (K1 and K2) parameters in the runoff response routine module were also sensitive (Fig. 6). MAXBAS showed strong sensitivity in the Pefferlaw River by tending toward maximum range but almost non-sensitive in Beaver River (Fig. 6d). Both of the recession parameters which partition the generated runoff into overland/subsurface and baseflow were sensitive in both catchments with their behavioral parameter set significantly different from each other (Table 1). The parameter K1 was skewed toward minimum range in Beaver River (Fig. 6e) but K2 was skewed toward minimum range in Pefferlaw River (Fig. 6f). The significant difference between the values of the most sensitive parameters demonstrates that despite their close proximity, the hydrological behavior of the two study catchments is markedly different. 3.2.2. Runoff simulations The HBV applications presented here adequately simulated the inter-annual hydrological conditions of the two adjacent catchments (Fig. 7). Model performances in simulating runoff conditions in both catchments were similar with a NS value of 0.53 in Beaver River (Fig. 7a) and 0.55 in Pefferlaw River catchment (Fig. 7b). The models captured baseflow well except during 1995 when the observed

Fig. 5. Long term shift in mean monthly precipitation (November–April) to either rainfall and snowfall under future CGCM3 A1B and A2 emission scenarios relative to the historical baseline (1961–2000). Precipitation was partitioned to either rainfall when TT N 0 °C or snowfall when TT b 0 °C respectively. The error bars represent standard error.

332

S.K. Oni et al. / Science of the Total Environment 473–474 (2014) 326–337

Fig. 6. Cumulative frequency distribution functions (CDF) of behavioral model parameters (top 200) in Beaver and Pefferlaw River catchments.

Fig. 7. HBV simulations of present day runoff conditions in a) Beaver River and b) Pefferlaw River from Sept. 1, 1992 to Aug.31, 2008 hydrological years, using the best performing parameter set from the ensemble of Monte Carlo model runs.

baseflow was abnormally high, and underrepresented by the model. This high baseflow phenomenon was unique to the Beaver River and was not observed in the Pefferlaw River or in other years in the Beaver River. On the seasonal scale, spring runoff dynamics differed between the two catchments. The spring peak flow occurred during April in Beaver River but was extended over March and April in the Pefferlaw River (Fig. 8). The results also showed that representing the integrated hydrologic behavior of semi-natural, snow-dominated catchments with a single best performing parameter set under-predicted spring runoff events (Fig. 8). Despite the acceptable simulations of long term runoff dynamics in the catchments (Fig. 7), best parameter sets missed the spring peak by 15% in Beaver (Fig. 8a) and by up to 32% in Pefferlaw River (Fig. 8b). The mean from the ensemble of behavioral parameter sets performed better in representing the peak spring and late autumn events in the Beaver River (Fig. 8a) but slightly overestimated baseflow conditions in summer months. In the Pefferlaw River, however, the 90th percentile of behavioral runs outperformed both the best parameter set and the behavioral mean in representing the runoff conditions in Jan–May, including the snowmelt-driven spring runoff peaks (Fig. 8b). Representative behavioral parameter sets (Fig. 8) were coupled with CGCM3 climate series (Fig. 3) to project plausible hydrological

S.K. Oni et al. / Science of the Total Environment 473–474 (2014) 326–337

Fig. 8. Integrated catchment versus model runoff behaviors using ensemble model runs in (a) Beaver and (b) Pefferlaw River catchments.

333

an annual mean temperature of only 6.5 °C. Both the winter and spring season will be affected as the rates of snow accumulation and melting are functions of temperature (Barnett et al., 2005). This is contrast to rainfall–dominated watersheds that were only sensitive to, and driven by, rainfall and evapotranspiration (Gebrehiwot et al., 2013; Lawrence and Haddeland, 2011). Model projection of a shorter and warmer winter regime indicates a possible reduction in the areal extent of snow cover. A decrease in snowfall depth or areal snow cover could suggest a loss of snow-held water storage in the winter and could shift the timing and magnitude of spring melt runoff peaks (Barnett et al., 2005). This can be explained by the shift in projected precipitation patterns away from snowfall dominance that characterizes LSW to possible rainfall dominance in the future. As suggested by the IPCC, a declining snow cover from the 1960s is evidenced from their continental-scale compilation of analyses of remotely sensed satellite images in high latitude regions (McCarthy et al., 2001). A similar analysis of historical satellite images is required in the LSW to give us further insights on the extent of spatiotemporal retraction of historical snow cover, a proxy for the plausible future trajectory as climate changes. Future work in LSW could utilize an ensemble of climate models to constrain the uncertainty in climate predictions. The use of ensemble climate projections that combine different climate and hydrological models (e.g. Barnston et al., 2003; Block et al., 2009) has been shown to outperform projections based on single climate model (Oni et al., 2012b). However, some credence could be given to the downscaled series as a previous detailed ANOVA conducted by Crossman et al. (2013) has shown that there was no significant difference between the NCEP, present day CGCM3 predictors and observed predictands. Also, incorporating the temporal evolution of land use changes into the modeling process might reduce the uncertainty in the hydrological projections around LSW. 4.2. Change in growing season

responses of the two adjacent catchments under A1B and A2 emission scenarios (Fig. 9). Future projections under both emission scenarios demonstrated that the largest changes could be in the form of a reduction in spring runoff towards the end of the century (and most pronounced in the Beaver River). These could lead to significant loss of seasonality in both catchments (Fig. 9). The annual average of the long term runoff projections (2040–2099 hydrologic year) were compared for the two adjacent catchments (Fig. 10). Results indicated that spring runoff could decrease by up to 35% in the Beaver and up to 25% in the Pefferlaw River catchment by the end of the century (Fig. 10). 4. Discussion 4.1. Hydroclimatic change Changes in climatic conditions have large implications for snow dominated watersheds (Barnett et al., 2005; Cunderlik and Simonovic, 2005). The impacts will be larger on water resources in high latitude locations (Tetzlaff et al., 2013) including Canada where the climate has already become warmer and predominantly wetter across the country in the last 20th century (Vincent and Mekis, 2006; Zhang et al., 2000). The historical climate signal is evident in LSW as long term weather data over the historical baseline period showed that the watershed is evolving towards warmer conditions. The 0.7 °C increase in annual temperature observed over the 40-year baseline period is large and may contribute in part to the water quality impairment in tributary rivers in LSW observed in recent decades (e.g. Eimers et al., 2005; Jin et al., 2013; Whitehead et al., 2011; Winter et al., 2007). Both the A1B and A2 scenarios considered here suggest that the increase in annual temperatures could reach 3.5 °C by 2100. This falls in the range of the IPCC ensemble projections i.e., a 1.8 to 4 °C global increase in temperature by 2100 (IPCC, 2007), and such an increase will have significant implications for a snow-dominated watershed with

Growing season is a function of air temperature and results suggest that there could be more days with air temperature above 5 °C in the future. Growing season can serve as a second order control of hydrological processes in heavily agriculture impacted watersheds such as Lake Simcoe. The average growing season length in LSW was estimated as 201 days in our baseline period. Since agricultural activities are large around LSW and constitute up to 51% of total land use in Pefferlaw and 65% in Beaver River (Baulch et al., 2013; Jin et al., 2013), the projected longer growing season by up to 26% in the future would have large implications on the net primary production, evapotranspiration regimes, net soil moisture balances and the extent of winter snow cover or density in the watershed. The implication of this is the alteration in the freeze–thaw cycles. Spring could be extended into late winter as the last frost day could come earlier, or autumn could encroach into early winter as the first frost day of the next hydrological year comes later. This could lead to possible increases in concentration of dissolved organic matter (DOM) in the streams due to the strong influence of adjacent agricultural landscapes surrounding the river coupled with ambient conditions that could lengthen the growing season (Wohlfart et al., 2012). Oni et al. (2013b) demonstrated strong positive concentration–runoff relationships during the growing than dormant season. The projected future increases in growing degree days could, however, have a large impact on the crop yield in this agriculture dominated watershed. For example, Lobell and Asner (2003) reported a 17% decline in corn/soybeans yield for every degree rise in growing season temperature in the United States. A more recent observation by You et al. (2009) similarly noted a 3–10% decline in wheat yield for every degree rise in the growing season temperature in China. The widely reported decline in crop yield with rising temperature is explainable as many temperate crops require prolonged cold winter treatment to gain flowering competence needed for effective yield in the following spring and summer seasons (Amasino, 2004; Craufurd and Wheeler, 2009).

334

S.K. Oni et al. / Science of the Total Environment 473–474 (2014) 326–337

Fig. 9. Integrated runoff behavior and projections in Beaver and Pefferlaw River catchments in 2040–2070 (column 1) and 2071–2100 (column 2) windows following ensemble indices derived in Fig. 8.

This has been termed “vernalization” (Chouard, 1960), as many plants have winter memory effects (Kim et al., 2009; Sung and Amasino, 2005). The process of vernalization is not only needed by cereal crops but also by perennial plants including fruit trees, to break dormancy before flowering and ensure an efficient fruit yield. Change in growing season length could therefore lead to earlier flowering/maturity (Craufurd and Wheeler, 2009) or increased occurrence of cases where maximum crop yields comes earlier in the spring/summer. In addition to the altered vernalization processes (Sung and Amasino, 2004), there could also be increased pest infestation and virulence of diseases due to more favorable conditions for their proliferation as winter/spring seasons get warmer (Rosenzweig et al., 2001). These show the possible implications of increasing growing season length associated with future

warming of winter/spring in the LSW but there is a paucity of literature to directly quantify the effects in the watershed. 4.3. Uncertainty assessments Considering the level of both quantifiable and unquantifiable uncertainties associated with climate projections, impacts on watershed hydrological processes could be larger than presently anticipated if the climate of tomorrow is driven towards future projections. There are several limitations and uncertainties associated with the climate change modeling and hydrologic impact analyses presented here. These include, but are not limited to, the number of climate models used and the uncertainties inherent in the downscaling process in SDSM. The regression-based approach of SDSM could only account for a fraction

S.K. Oni et al. / Science of the Total Environment 473–474 (2014) 326–337

Fig. 10. Changes in runoff shifts in Beaver and Pefferlaw River using criteria established in Fig. 9. Future runoff conditions were projected in Beaver River using behavioral mean of ensemble runoff projections and 90% percentile of behavioral ensemble in Pefferlaw River from 2040 to 2099 hydrologic year.

of total variation in the climate model being downscaled and this could be more pronounced on precipitation than temperature (Dibike and Coulibaly, 2005; Salathé, 2003; Wilby et al., 2002). However, the main focus of this study is on the uncertainty in hydrological modeling in managed watersheds than climate and this will be discussed further. Uncertainty in hydrological modeling is an important learning process about catchment behavior, model structure and input data towards a more robust projection of future conditions (Juston et al., 2013). While model structures would definitely contribute to uncertainty in our hydrological projections and could be constrained by comparing with other models in a future study, parameter uncertainty analysis and behavior are considered here. The best parameter sets derived from static calibration period might not work well in transitional catchments such as in LSW. This could be more pronounced in the future as the catchment drifts from snowfall to rainfall dominance and present a potential watershed management issue to the relevant stakeholders. 4.3.1. Integrated catchment behavior The parameter spaces in the HBV rainfall–runoff model can be multidimensional, leading to equifinality where several parameter spaces produce similar optimal model behaviors (Beven, 2006). This could lead to complementary parameter behavior. While the ranges of some sensitive parameters are narrower than those reported in other studies (Hamilton et al., 2000; Seibert and McDonnell, 2010; Seibert, 1997; Zégre et al., 2010), results are difficult to compare directly due to differences in local climate and fundamental physiographic characteristics of the catchments (e.g. elevation, geology, landscape features, dominant land covers etc.). All these contribute in part to the differences in overall catchment behavior. The similar behavioral patterns and non-sensitivity of the parameter TT indicates that the two adjacent catchments considered here were under the same weather influence under present day condition. However, it is not fully clear if the present day assumption of TT will hold for the rainfall–snowfall patterns in the future term. Other paired catchment studies have also observed similar TT behavioral patterns, although the parameter was highly sensitive in their catchments (Hamilton et al., 2000; Seibert, 1997). Despite similar weather influence, the model suggested that there could be differences in the rate of snowmelt accumulation and melt patterns between the two catchments. The CDF of CFMAX toward the minimum range, suggesting a greater tendency for winter snow accumulation in Pefferlaw River catchment than in the Beaver, perhaps because of possible higher melt rates in the Beaver River. This might explain why the uncertainty in capturing the spring peak is larger in Pefferlaw with an extended melt regime (March–April) than in the Beaver River which only peaked in April (Fig. 8). This could also explain why Pefferlaw tends towards higher SFCF values as areal distribution of snow cover could be less

335

uniform due to the steeper elevation gradient and interactions of orographic effects at the southern end of the catchment. The terrestrial area of Pefferlaw River is bound by the ORM at the south end of the catchment (Fig. 1). However, the behavioral patterns and sensitivity of these snow related parameters could change in the future as precipitation patterns shift from snow to more rainfall dominance (Lawrence and Haddeland, 2011). With the exception of the β parameter, the sensitivity of soil routine parameters (Table 1) was not as strong as snow and runoff response routine parameter sets. These observations were in line with earlier studies (except FC) elsewhere where response routine parameters were also shown to be more sensitive (Hundecha and Bárdossy, 2004; Seibert, 1999). The behavior of K1 suggested more quick/subsurface flows in the Pefferlaw compared to the Beaver River. This is explainable as Beaver River catchment has relatively low gradient that support less overland flows (LSRCA, 2012a,b). This could also explain the larger uncertainty in simulating spring floods in Pefferlaw than Beaver River. However, the behavioral pattern of K2 that controls the baseflow rate showed that both catchments differed significantly in terms of their groundwater dynamics. Pefferlaw River showed high sensitivity and tended towards the minimum range in contrast to virtual nonsensitivity of the parameter in the Beaver River catchment. This is explainable as baseflow conditions in Pefferlaw River is influenced by 1) the aquifer complex of ORM with continuous discharge of groundwater to the headwater stream networks and 2) a broad shallower aquifer towards the downstream reaches. However, the growing rural/urban communities around the headwater and increasing extraction of sand and gravel deposits in the ORM might increase groundwater sensitivity to disturbances in this region. The shallow aquifer towards the downstream is an important control on the river baseflow conditions and can be very sensitive to landscape disturbances due to its proximity to the surface (Johnson, 1997). However, agricultural and other human activities are concentrated in this downstream Algonquin shallow aquifer region. These can explain the sensitivity of baseflow recession and provide a hint on the overall future sensitivity of groundwater dynamics to land use changes in Pefferlaw River catchment. The non sensitivity of K2 (baseflow) in the Beaver River is explainable as the catchment has low gradient with a large wetland complex (about 19%) that aggregated around the center of the catchment (Baulch et al., 2013; Jin et al., 2013). This gives Beaver River catchment large storage capacity and slow flowing streams (LSRCA, 2012b). The catchment soils also have larger infiltration capacity to recharge the groundwater; a dominant driver of streamflow in the catchment (LSRCA, 2012b). The influence of wetland hydrology in the Beaver River could therefore be large and could explain the higher baseflow in part of the record (1995) that HBV could not account for (Fig. 7). The 1995 hydrological year was an unusually wet year in LSW (Oni et al., 2013b) coupled with the presence of large, well defined discharge zones around the wetland complex (Baulch et al., 2013). However, the uncertainty in simulating this baseflow condition was not pronounced in the Pefferlaw River. Most headwaters of Beaver headwaters are fed by several discharge springs with only a small part originating from the ORM (LSRCA, 2012b). The stronger headwater connectivity and continuous contributions of groundwater to the Beaver River could elevate the baseflow in wet years but could also buffer the sensitivity of the catchment to extreme weather especially in the dry years. This might explain why future emission scenarios projected higher baseflow conditions in the Beaver than in the Pefferlaw River. However, loss of welldefined seasonality appeared to be pronounced in the two adjacent river catchments (Fig. 9). 5. Conclusion While there are increasing efforts in predictive water quality modeling and nutrient budget calculations in the watershed, the uncertainty in hydrological representation of its tributary rivers are high due to

336

S.K. Oni et al. / Science of the Total Environment 473–474 (2014) 326–337

impact of multiple stressors. Getting the hydrology right is therefore necessary for sustainable management of Lake Simcoe water resources under changing climate and land use conditions. Both A1B and A2 emission scenarios projected a warmer (up to 4 °C) and extended growing season (up to 26%) in the future. This will have large impacts on agricultural activities that dominate the LSW. Both scenarios also projected possibility of a shift in precipitation patterns toward rainfall dominance and a decline in annual runoff. This will have a large impact on biogeochemical processes and water balance in the watershed. Best parameter set missed the spring event by 15% in Beaver and up to 32% in Pefferlaw despite robust parameter optimization strategy. Utilizing behavioral parameters (as oppose to best parameter) is necessary for credible representation of present day hydrology of a managed watershed and to constrain the uncertainties in the future projections. Our previous modeling analysis in LSW had assumed that calibration functions in the Beaver River catchment for example, could be used to predict runoff conditions in the ungauged adjacent White Creek catchment (Jin et al., 2013; Oni et al., 2011, 2012a, 2013b). The assumption was based on the similarity of weather and land use patterns between Beaver River and White Creek catchments. Although Patil and Stieglitz (2013) suggested that streamflow predictions could be achieved in an ungauged basin by applying the sensitive parameter sets in one catchment to the neighboring catchments under the condition of spatial proximity and/or physical similarity. This does not appear to hold in the two adjacent agricultural catchments considered here as the sensitive parameter sets exhibited contrasting behavioral patterns. The contrasting patterns of the most sensitive parameters demonstrates that while the two adjacent catchments have similar land use patterns and shared some physiographic features, the differences in the integrated hydrological behavior of the two study catchments could be exacerbated by human activities within LSW. This indicates that parameter functions from one model can not necessarily be used to create models that are representative of the hydrological conditions in the other. However, this needs to be validated with other nested catchments of the LSW with varying degrees of landscape disturbances. Future works should utilize an ensemble modeling approach that will incorporate the temporal evolution of land use changes into the modeling processes. This will further constrain the uncertainty in hydrological or biogeochemical projections as climate and land use changes in the watershed. Acknowledgments We acknowledge the Data Access Integration (DAI, see http:// quebec.ccsn.ca/DAI/) team for providing the predictor data used in this study. The DAI data download gateway is made possible through collaboration among the Global Environmental and Climate Change Centre (GEC3), the Adaptation and Impacts Research Division (AIRD) of Environment Canada, and the Drought Research Initiative (DRI). We thank Jan Seibert for making the HBV light version freely available. We would also like to thank Environmental Canada (EC) as well as the Lake Simcoe Region Conservation Authority (LSRCA) for providing historical flow and weather data as well as other catchment data. MNF was supported by the MISTRA Future Forests program. We thank Hamdi Jarjanazi of MOE for making the map used in Fig. 1 and two anonymous reviewers for their comments that improved the manuscript. This research was supported by a Natural Science and Engineering Research Council (NSERC) strategic grant to PJ Dillon. References Allen MR, Ingram WJ. Constraints on future changes in climate and the hydrologic cycle. Nature 2002;419:224–32. Amasino R. Vernalization, competence, and the epigenetic memory of winter. Plant Cell 2004;16:2553–9. [Online]. Arnell NW. Climate change and global water resources: SRES emissions and socioeconomic scenarios. Glob Environ Chang 2004;14:31–52.

Barnett TP, Adam JC, Lettenmaier DP. Potential impacts of a warming climate on water availability in snow-dominated regions. Nature 2005;438:303–9. Barnston AG, Mason SJ, Goddard L, Dewitt DG, Zebiak SE. Multimodel ensembling in seasonal climate forecasting at IRI. Bull Am Meteorol Soc 2003;84:1783–96. Baulch HM, Futter MN, Jin L, Whitehead PG, Woods DT, Dillon PJ, et al. Phosphorus dynamics across intensively monitored subcatchments in the Beaver River. Inland Waters 2013;3:187–206. Bergström S, Harlin J, Lindström G. Spillway design floods in Sweden: I. New guidelines. Hydrol Sci J 1992;37(5):505–19. Beven K. A manifesto for the equifinality thesis. J Hydrol 2006;320:18–36. Block PJ, Souza Filho FA, Sun L, Kwon HH. A streamflow forecasting framework using multiple climate and hydrological Models1. J Am Water Res Assoc 2009;45:828–43. Campbell JL, Mitchell MJ, Groffman PM, Christenson LM, Hardy JP. Winter in northeastern North America: a critical period for ecological processes. Front Ecol Environ 2005;3: 314–22. Chouard P. Vernalization and its relations to dormancy. Annu Rev Plant Physiol 1960;11: 191–238. Craufurd P, Wheeler T. Climate change and the flowering time of annual crops. J Exp Bot 2009;60:2529–39. Crossman J, Futter M, Oni SK, Whitehead P, Jin L, Butterfield D, et al. Impacts of climate change on hydrology and water quality: future proofing management strategies in the Lake Simcoe watershed, Canada. J Great Lakes Res 2013;39:19–32. Cunderlik JM, Simonovic SP. Hydrological extremes in a southwestern Ontario river basin under future climate conditions. J Hydrol Sci 2005;50:654. Dibike YB, Coulibaly P. Hydrologic impact of climate change in the Saguenay watershed: comparison of downscaling methods and hydrologic models. J Hydrol 2005;307: 145–63. Eimers MC, Winter JG, Scheider WA, Watmough SA, Nicholls KH. Recent changes and patterns in the water chemistry of Lake Simcoe. J Great Lakes Res 2005;31:322–32. Euskirchen E, McGuire A, Chapin FS. Energy feedbacks of northern high‐latitude ecosystems to the climate system due to reduced snow cover during 20th century warming. Glob Chang Biol 2007;13:2425–38. Franczyk J, Chang H. The effects of climate change and urbanization on the runoff of the Rock Creek basin in the Portland metropolitan area, Oregon, USA. Hydrol Process 2009;23:805–15. Furusho C, Andrieu H, Chancibault K. Analysis of the hydrological behaviour of an urbanizing basin. Hydrol Process 2013. http://dx.doi.org/10.1002/hyp.9706. Futter M, Butterfield D, Cosby B, Dillon P, Wade A, Whitehead P. Modeling the mechanisms that control in-stream dissolved organic carbon dynamics in upland and forested catchments. Water Resour Res 2007;43:W02424. Gebrehiwot SG, Seibert J, Gärdenäs AI, Mellander P-E, Bishop K. Hydrological change detection using modeling: half a century of runoff from four rivers in the Blue Nile Basin. Water Resour Res 2013;49. http://dx.doi.org/10.1002/wrcr.20319. Hamilton A, Hutchinson D, Moore R. Estimating winter streamflow using conceptual streamflow model. J Cold Reg Eng 2000;14:158–75. Harris G, Collins M, Sexton D, Murphy J, Booth B. Probabilistic projections for 21st century European climate. Nat Hazards Earth Syst Sci 2010;10:2009–20. Hundecha Y, Bárdossy A. Modeling of the effect of land use changes on the runoff generation of a river basin through parameter regionalization of a watershed model. J Hydrol 2004;292:281–95. Huntington TG. Evidence for intensification of the global water cycle: review and synthesis. J Hydrol 2006;319:83–95. IPCC. The physical science basis. contribution of working group I to the fourth assessment report of the intergovernmental panel on climate change. In: Solomon S, Qin D, Manning M, Chen Z, Marquis M, Averyt KB, Tignor M, Miller HL, Solomon S, Qin D, Manning M, Chen Z, Marquis M, Averyt KB, Tignor M, Miller HL, editors. Climate change 2007: the physical science basis. Cambridge, United Kingdom and New York, NY, USA: Cambridge University Press; 2007. p. 996. Jin L, Whitehead PG, Baulch HM, Dillon PJ, Butterfield DA, Oni SK, et al. Modelling phosphorus in the Lake Simcoe catchments and lake system: scenario analysis to assess alternative management strategies. Inland Waters 2013;3:207–20. Johnson FM. Landscape ecology of the Lake Simcoe basin. Lake Reservoir Manage 1997;13(3):226–39. Juston JM, Kauffeldt A, Montano BQ, Seibert J, Beven KJ, Westerberg IK. Smiling in the rain: seven reasons to be positive about uncertainty in hydrological modelling. Hydrol Process 2013;27:1117–22. Kim DH, Doyle MR, Sung S, Amasino RM. Vernalization: winter and the timing of flowering in plants. Annu Rev Cell Dev 2009;25:277–99. Krause P, Boyle D, Bäse F. Comparison of different efficiency criteria for hydrological model assessment. Adv Geosci 2005;5:89–97. Laudon H, Tetzlaff D, Soulsby C, Carey S, Seibert J, Buttle J, et al. Change in winter climate will affect dissolved organic carbon and water fluxes in mid‐to high latitude catchments. Hydrol Process 2013;27:700–9. Lawrence D, Haddeland I. Uncertainty in hydrological modelling of climate change impacts in four Norwegian catchments. Hydrol Res 2011;42:457–71. Lindström G, Johansson B, Persson M, Gardelin M, Bergström S. Development and test of the distributed HBV-96 hydrological model. J Hydrol 1997;201:272–88. Lobell DB, Asner GP. Climate and management contributions to recent trends in US agricultural yields. Science 2003;299:1032. LSRCA. Lake Simcoe and Its Watershed. Toronto: Ministry of the Environment; 2008116. LSRCA. Pefferlaw River subwatershed plan. Lake Simcoe Region Conservation Authority. Newmarket; 2012a354. LSRCA. Beaver River subwatershed plan. Lake Simcoe Region Conservation Authority. Newmarket; 2012b336. McCarthy JJ, Canziani OF, Leary NA, Dokken DJ, White KS. Climate change 2001: impacts, adaptation, and vulnerability: contribution of Working Group II to the third assessment

S.K. Oni et al. / Science of the Total Environment 473–474 (2014) 326–337 report of the Intergovernmental Panel on Climate Change. Cambridge University Press; 2001. Menzel A, Fabian P. Growing season extended in Europe. Nature 1999;397:659. Nash JE, Sutcliffe J. River flow forecasting through conceptual models part I—a discussion of principles. J Hydrol 1970;10:282–90. O'Connor EM, McConnell C, Lembcke D, Winter JG. Estimation of total phosphorus loads for a large, flashy river of a highly developed watershed—seasonal and hysteresis effects. J Great Lakes Res 2011;37(Supplement 3):26–35. Oni SK. The effects of changing land use and climate on the hydrology and carbon budget of Lake Simcoe watershed, Ontario. Environmental and Life Sciences. PhD. Peterborough, Canada: Trent University; 2011310. Oni SK, Futter MN, Dillon PJ. Landscape-scale control of carbon budget of Lake Simcoe: a process-based modelling approach. J Great Lakes Res 2011;37:160–5. Oni S, Futter M, Molot L, Dillon P. Modelling the long term impact of climate change on the carbon budget of Lake Simcoe, Ontario using INCA-C. Sci Total Environ 2012a;414:387–403. Oni SK, Dillon PJ, Metcalfe RA, Futter MN. Dynamic modelling of the impact of climate change and power flow management options using STELLA: application to the Steephill Falls Reservoir, Ontario, Canada. Can Water Resour J 2012b;37:125–48. Oni SK, Futter M, Bishop K, Köhler S, Ottosson-Löfvenius M, Laudon H. Long term patterns in dissolved organic carbon, major elements and trace metals in boreal headwater catchments: trends, mechanisms and heterogeneity. Biogeosciences 2013a;10:2315–30. Oni SK, Futter MN, Molot LA, Dillon PJ. Adjacent catchments with similar patterns of land use and climate have markedly different DOC and runoff dynamics. Hydrol Process 2013b. http://dx.doi.org/10.1002/hyp.9681. Patil S, Stieglitz M. Modelling daily streamflow at ungauged catchments: what information is necessary? Hydrol Process 2013. http://dx.doi.org/10.1002/hyp.9660. Rosenzweig C, Iglesias A, Yang X, Epstein PR, Chivian E. Climate change and extreme weather events; implications for food production, plant diseases, and pests. Glob Chang Hum Health 2001;2:90–104. Salathé Jr EP. Comparison of various precipitation downscaling methods for the simulation of streamflow in a rainshadow river basin. Int J Climatol 2003;23:887–901. Seibert J. Estimation of parameter uncertainty in the HBV model. Nord Hydrol 1997;28:247–62. Seibert J. Regionalisation of parameters for a conceptual rainfall–runoff model. Agr Forest Meteorol 1999;98:279–93.

337

Seibert J, McDonnell JJ. Land-cover impacts on streamflow: a change-detection modelling approach that incorporates parameter uncertainty. Hydrol Sci J 2010;55:316–32. Spear R, Hornberger G. Eutrophication in Peel Inlet—II. Identification of critical uncertainties via generalized sensitivity analysis. Water Res 1980;14:43–9. Sung S, Amasino RM. Vernalization and epigenetics: how plants remember winter. Curr Opin Plant Biol 2004;7:4–10. Sung S, Amasino RM. Remembering winter: toward a molecular understanding of vernalization. Annu Rev Plant Biol 2005;56:491–508. Tetzlaff D, Soulsby C, Buttle J, Capell R, Carey SK, Laudon H, et al. Catchments on the cusp? Structural and functional change in northern ecohydrology. Hydrol Process 2013;27: 766–74. Vincent LA, Mekis É. Changes in daily and extreme temperature and precipitation indices for Canada over the twentieth century. Atmos Ocean 2006;44:177–93. Whitehead P, Jin L, Baulch H, Butterfield D, Oni S, Dillon P, et al. Modelling phosphorus dynamics in multi-branch river systems: a study of the Black River, Lake Simcoe, Ontario, Canada. Sci Total Environ 2011;412–413:315–23. Wilby RL, Dawson CW. The Statistical DownScaling Model: insights from one decade of application. Int J Climatol 2013;33:1707–19. Wilby RL, Dawson CW, Barrow EM. SDSM—a decision support tool for the assessment of regional climate change impacts. Environ Model Software 2002;17: 145–57. Winter JG, Eimers MC, Dillon PJ, Scott LD, Scheider WA, Willox CC. Phosphorus inputs to Lake Simcoe from 1990 to 2003: declines in tributary loads and observations on lake water quality. J Great Lakes Res 2007;33:381–96. Wohlfart T, Exbrayat J, Schelde K, Christen B, Dalgaard T, Frede H, et al. Spatial distribution of soils determines export of nitrogen and dissolved organic carbon from an intensively managed agricultural landscape. Biogeosciences 2012;9:4513–25. You L, Rosegrant MW, Wood S, Sun D. Impact of growing season temperature on wheat productivity in China. Agr Forest Meteorol 2009;149:1009–14. Zégre N, Skaugset AE, Som NA, McDonnell JJ, Ganio LM. In lieu of the paired catchment approach: Hydrologic model change detection at the catchment scale. Water Resour Res 2010;46:W11544. Zhang X, Vincent LA, Hogg W, Niitsoo A. Temperature and precipitation trends in Canada during the 20th century. Atmos Ocean 2000;38:395–429.