Modulation of linear and nonlinear hydroclimatic dynamics by mountain glaciers in Canada and Norway: Results from information-theoretic polynomial selection

Abstract

Historical streamflow and climate datasets were analyzed for low- and high-frequency hydroclimatic variability. Four glacial/non-glacial catchment pairs were considered, two from the southern Canadian Rocky Mountains and two from arctic coastal Norway. Analyses were performed using daily data, providing high seasonal resolution and facilitating the identification of possible nonlinear hydroclimatic processes. Spearman rank correlation, and an information theory-based polynomial selection method, were employed in parallel. The latter permits straightforward identification of highly nonlinear relationships, simultaneous consideration of multiple models and estimation of the probability of a given relationship, as distinct from conventional p-values. Highly nonlinear (parabolic) atmospheric teleconnections to the Arctic Oscillation and El Niño-Southern Oscillation were confirmed in Norway and Canada, respectively, and their corresponding hydrologic effects were detected; conversely, little evidence for deviation from linearity was found for long-term monotonic trends. Presence or absence of watershed glacial cover was found to fundamentally alter streamflow responses to climate variability and change. In particular, for ecologically highly relevant late-summer low flows, glaciers induced: (1) stronger negative long-term trends than observed for non-glacial basins, presumably reflecting net mass balance declines seen in nearby glaciers, and (2) parabolic teleconnections, largely absent in non-glacial basins, reflecting parabolic air temperature teleconnections and the presence or absence of glacial ice available for melting.

Des ensembles de données historiques sur les débits et sur le climat ont été analysés pour en dégager la variabilité hydroclimatique à basse et à haute fréquence. Quatre paires de bassins versants glaciaires/non glaciaires ont été étudiées, deux situées dans le sud des Rocheuses canadiennes et deux dans les régions arctiques des côtes de la Norvège. Des analyses ont été menées à l’aide de données quotidiennes. Elles ont fourni une résolution saisonnière élevée et ont facilité l’identification de processus hydroclimatiques non linéaires possibles. Le coefficient de corrélation de rang de Spearman et une méthode de sélection polynomiale fondée sur la théorie de l’information ont été employés en parallèle. Cette dernière permet la détermination directe de relations hautement non linéaires, la prise en considération simultanée de modèles multiples et l’estimation de la probabilité d’une relation donnée, le tout étant distinct des valeurs de p traditionnelles. Les téléconnexions atmosphériques (paraboliques) hautement non linéaires associées à l’oscillation arctique et à El Niño-oscillation australe ont été confirmées en Norvège et au Canada, respectivement, et leurs effets hydrologiques correspondants ont été détectés; inversement, peu de preuves de l’écart par rapport à la linéarité ont été recueillies pour les tendances monotones à long terme. Il a été constaté que la présence ou l’absence de couverture de dépôts glaciaires dans le bassin a fondamentalement altéré les réactions des cours d’eau à la variabilité et au changement climatiques. En particulier, pour les débits faibles de fin d’été hautement pertinents d’un point de vue écologique, les glaciers ont provoqué : (1) des tendances négatives à long terme plus fortes que celles observées pour les bassins non glaciaires, ce qui reflète probablement les diminutions nettes du bilan de masse observées dans les glaciers à proximité et (2) des téléconnexions paraboliques, largement absentes dans les bassins non glaciaires, ce qui reflète les téléconnexions paraboliques associées aux températures de l’air et la présence ou l’absence de glace glaciaire disponible pour la fonte.

View correction statement:

Erratum

Introduction

Continental-scale “water towers” (e.g. Qiu 2010) like the Rocky Mountains, Alps, Andes and Himalayas feed major and often transboundary river systems such as the Columbia, Yukon, Athabasca, Danube, Rhine, Indus, Ganges, Brahmaputra, Yangtze and Yellow (Huang He). Alpine glaciers and ice fields lie at the hydrologic heart of many of these mountain ranges, impressing a distinct signature upon water quantity, water quality and freshwater and coastal marine ecosystems (e.g. Ward 1994; Dorava and Milner 2000; Fleming 2005; Mark et al. 2005; Hood and Berner 2009; Neal et al. 2010; Jacobsen et al. 2012). Due to orographic precipitation and glacial melt, mountain headwaters often make inordinately large contributions to the flow of large rivers, and typically impart hydrological and biogeochemical signals that are clearly observable far downstream (e.g. Milner et al. 2009; Whitfield and Spence 2011; Jacobsen et al. 2012; Fleming and Sauchyn 2013). Growing global human populations and economies are expected to increase water scarcity by amplifying water demand, and potentially also by impacting water supply through increased net atmospheric greenhouse gas emissions which are believed to help drive global climatic changes (e.g. Vörösmarty et al. 2000; Arnell 2004; Murtaugh and Schlax 2009). As water systems operate closer to the margins, a given climatic fluctuation may be more likely to push the supply-demand balance from surplus to deficit, likely increasing the practical consequences of variability in climatic forcing. More broadly, a variety of mechanisms and timeframes of climatic variability and change exist, and these may all affect the discharge of northern alpine rivers as well as associated biogeochemical properties, hydroelectric generation capacity, nutrient fluxes, freshwater habitat quantity and quality, coastal ocean dynamics, and fish abundance, for example (e.g. Milner et al. 2009; Neal et al. 2010; Whitfield and Spence 2011; Jacobsen et al. 2012; Thomson et al. 2012; Fleming and Weber 2012). Glacier-influenced rivers in particular may be unusually sensitive to climatic shifts and therefore serve as important indicator systems (e.g. Milner et al. 2009).

The comparative hydroclimatic dynamics of glacial vs. non-glacial rivers is an emerging theme in mountain hydrology. Of course, it has long been understood that discharge in glacier melt-driven rivers is related to glacial variability and change (e.g. Pellicciotti et al. 2010; Marshall et al. 2011), and that glaciers can modify the severity of year-to-year streamflow variations (e.g. Henshaw 1933; Meier 1969; Fountain and Tangborn 1985; Moore 1992; Fleming and Clarke 2005). However, a growing body of empirical literature – jointly studying both glacial rivers, and non-glacial control catchments – is additionally demonstrating that the presence or absence of a glacier in a watershed’s headwaters can affect the basic nature of that river’s response to both coherent modes of climatic variation and longer-term climate changes (Neal et al. 2002; Fleming and Clarke 2003; Lafrenière and Sharp 2003; Stahl and Moore 2006; Fleming et al. 2006; Fleming et al. 2007; Hodgkins 2009; Brabets and Walvoord 2009; Casassa et al. 2009; Moore et al. 2009; Li et al. 2010; Baraer et al. 2012; Dahlke et al. 2012; Fleming and Weber 2012). These glacial modulating effects can be profound, even leading in some cases to selective teleconnectivity, or to opposite climate change responses from one catchment to the next within a given region. Although the underlying physical processes are complex, to first order, these differences between the hydroclimatic responses of nearby and otherwise closely comparable glacial and non-glacial rivers can be attributed to the presence or absence of an additional term in the water balance, i.e. a permanent ice mass potentially available for melt generation. Such issues are obviously crucial to understanding hydroclimatic heterogeneity, e.g. reliable regionalization. There is additionally growing interest in implications to salmon ecology: these keystone species often tend to return to fresh water to spawn during late summer and early autumn, when the distinguishing hydrologic characteristics of glacial rivers are most pronounced. Climatic variability and change may therefore influence corresponding freshwater habitat quantity and quality in different ways between nearby glacial and non-glacial rivers; for a review, see Moore et al. (2009). Clearly, implications for water supply availability and hydropower generation are also of interest. However, many of the foregoing results are very recent, and much remains to be learned about the comparative hydroclimatic dynamics of glacial and non-glacial rivers.

Here, novel statistical analysis techniques are applied in a case study of several paired catchments in Canada and Norway. The primary goal is to assess nonlinear impacts from both low- and high-frequency climatic variability, and whether these differ between glacial and non-glacial basins. Low-frequency climatic variability is taken to mean long-term monotonic trends, as potentially associated with anthropogenic climate changes; high-frequency variability refers here to organized patterns of year-to-year ocean-atmosphere circulation, with a focus on two dominant northern hemisphere modes (Wallace and Thompson 2002), El Niño-Southern Oscillation (ENSO) and the Arctic Oscillation (AO). Several considerations motivate this emphasis on nonlinear effects. Rigorous statistical analyses of historical datasets for long-term trends often use methods robust to nonlinearity, such as Spearman correlation and Mann-Kendall tests. However, linear trend assessments remain common, in part because they provide a more intuitive sense of change rates; further, recent work has questioned whether departures from nonlinearity are significant for monotonic trends observable in short, noisy historical records (e.g. Fleming and Weber 2012). A more complete understanding of the hydroclimatic dynamics of mountain rivers, and the refinement and streamlining of standard trend analysis techniques, could benefit from further examination of whether hydroclimatic trends are linear or nonlinear. More pressing, however, is the assessment of highly nonlinear teleconnections to ENSO and the AO. It is well recognized that teleconnections are generally nonlinear, and modern methods in statistical hydroclimatology typically accommodate that fact (e.g. Hoerling et al. 1997; Fleming et al. 2007). However, it is almost universally assumed that such relationships between large-scale climate indices and local-scale hydrometeorological responses are monotonic. In reality, neural network-based methods capable of accommodating full nonlinearity, applied to gridded continental- and hemispheric-scale meteorological data, have revealed that wintertime temperature and precipitation responses to both ENSO and the AO are significantly non-monotonic (Wu et al. 2005; Hsieh et al. 2006). Specifically, these responses appear to be parabolic, with extreme positive- and negative-state climate events yielding similar temperature and precipitation responses. Some local- to regional-scale analyses of climate (Pozo-Vázquez et al. 2005; B. Taylor, Environment Canada, pers. com. 2012) and lake ice cover (Bai et al. 2012) have generated results consistent with this wider finding. While the exact physical mechanism for non-monotonic teleconnections appears to remain unclear, their specific form might be related to quadratic terms in the governing equations for atmospheric circulation (Hsieh et al. 2006). Additionally, Hsieh et al. (2006) suggested these nonlinear teleconnections propagate perturbations farther, into regions where classical linear teleconnections are insignificant. Apparently, however, no attempt has yet been made to assess the streamflow implications of such highly nonlinear teleconnections. Detection of parabolic streamflow teleconnections could radically alter the standard conceptual model of, and the standard suite of methodological tools used to investigate, how water resources respond to climatic variability. Implications include improving the ability to infer climate variability impacts on fish populations, or to perform early-season water supply forecasting, which often relies on climate indices.

The approach taken was to organize datasets and analyses into data groups, each containing a glacial river, a comparable non-glacial river and representative precipitation and temperature stations. Each data group constitutes an approximate paired-catchment experiment, facilitating assessment of how glacial and non-glacial rivers differentially respond to similar variations in regional climatic forcing. Four such data groups (for a total of eight rivers) were used. In particular, two data groups were taken from each of two regions: the southern Canadian Rocky Mountains, and arctic coastal Norway. These areas are closely comparable in some respects (northern alpine basins with significant maritime influences and mixed glacial and nival streamflow regimes), yet highly distinct in others (western vs. eastern hemispheres and different latitudes, elevations and degrees of continentality). These distinctions permit meaningful evaluation of region-to-region heterogeneity in hydroclimatic dynamics, which is important given that some differential glacial/non-glacial phenomena are regionally variable (e.g. Stahl and Moore 2006). A polynomial selection technique, based on the Akaike Information Criterion (AIC), was applied to relatively high-resolution (daily) data. This approach provides approximate estimates of the probability that a hydroclimatic relationship is present, and allows assessment of strongly nonlinear ENSO and AO impacts on local-scale climate and hydrology.

Data

Study area and data sources

Two widely separated mountain regions were considered as case studies: the Canadian Rocky Mountains and adjacent ranges in southeast British Columbia and southwest Alberta, and the coastal ranges of arctic Norway in the vicinity of Nordland (Figure 1). Canadian temperature (T) and precipitation (P) time series were taken from the second-generation Adjusted and Homogenized Historical Climate Data (AHCCD). These are historical Meteorological Service of Canada (MSC) climate station data that have been examined and homogenized for climate studies (Mekis and Vincent 2011; Vincent et al. 2012). Canadian streamflow (Q) data were selected from the Reference Hydrometric Basin Network (RHBN), a subset of Water Survey of Canada (WSC) hydrometric stations identified as suitable for hydroclimatic studies (Brimley et al. 1999; Harvey et al. 1999). Daily observations of 2-m temperature and precipitation records for the selected locations in Norway were taken from the Norwegian Meteorological Institute’s climate station database. Norwegian streamflow data were obtained from the national hydrological database system HYDRA II, maintained by the Norwegian Water Resources and Energy Directorate (NVE). In general, catchments with <~1% glacial cover by area were considered non-glacial, and catchments with >~10% glacial cover by area were considered glacial; potentially ambiguous intermediate levels of glaciation were avoided to the extent possible.

Figure 1. Location map. Elevation contours provided only for approximate visual reference; interval is ~600 m. Flat areas at east side of North American and west side of Scandinavian inset maps are the northern Great Plains and Norwegian Sea, respectively. Note differing scales for the two inset maps.

Two data groups were studied for each region. In Canada (Figure 1; Table 1), data group C1 included streamflow data from the glacially influenced Illecillewaet River, and from a non-glacially influenced location on the Kootenay River. Temperature and precipitation data were taken from the Glacier National Park climate station. Climate is mixed maritime and continental, with mild summers and cold winters. On average, precipitation peaks in winter and is associated primarily with frontal systems, but summertime convective activity is additionally important. Streamflow shows a strong seasonal peak in summer, corresponding to snow and (for glacial rivers) ice melt, with very low flows during winter. The C1 basins are located on the western flanks of the Rocky Mountains and form significant tributaries to the Columbia River. Data group C2 included streamflow data from the glacially influenced Mistaya River, and from non-glacial Cataract Creek. T and P data were taken from the Banff climate station. These rivers are located on the eastern flank of the Rocky Mountains and are headwater tributaries to the Saskatchewan River. Being located on the leeward side of the western North American Cordillera, weather and climate for these watersheds are generally similar to nearby data group C1, but slightly more continental.

Table 1. Hydrometric and climate station summary. Elevations are given in meters above sea level (masl). MAT and TAP denote the averages, over the period of record used in the analyses, for mean annual air temperature and total annual precipitation, respectively.

	Name	ID	Latitude	Longitude	Area (km^2)	% glacier	Runoff (mm)	Elevation (masl)	MAT (°C)	TAP (mm)
C1	Illecillewaet	08ND013	51.01N	118.08W	1170	9.8	1424	n/a	n/a	n/a
C1	Kootenay	08NF001	50.89N	116.04W	420	0.6	365	n/a	n/a	n/a
C1	Glacier	1173191	51.30N	117.52W	n/a	n/a	n/a	1323	1.8	1401
C2	Mistaya	05DA007	51.88N	116.69W	248	23.4	803	n/a	n/a	n/a
C2	Cataract	05BL022	50.29N	114.59W	166	0.0	256	n/a	n/a	n/a
C2	Banff	3050519	51.20N	115.55W	n/a	n/a	n/a	1397	2.8	510
N1	Berget	156.10.0	66.46N	13.88E	211	34.9	3577	n/a	n/a	n/a
N1	Forsbakk	156.15.0	66.29N	13.81E	56	0.0	2811	n/a	n/a	n/a
N1	Mo i Rana	79480	66.31N	14.15E	n/a	n/a	n/a	41	3.6	1443
N2	Bjørnfoss	156.13.0	66.47N	14.32E	306	20.4	2081	n/a	n/a	n/a
N2	Junkerdalselv	163.5.0	66.82N	15.41E	422	0.6	1059	n/a	n/a	n/a
N2	Dunderlandsdalen	79740	66.51N	14.91E	n/a	n/a	n/a	155	n/a	1366
N2	Glomfjord	80700	66.81N	13.98E	n/a	n/a	n/a	39	5.3	n/a

In Norway (Figures 1 and 2; Table 1), data group N1 included Q data from the glacially influenced Berget and non-glacial Forsbakk Rivers, and T and P data from the Mo i Rana climate station. Data group N2 included Q data from the glacial Bjørnfoss and non-glacial Junkerdalselv Rivers, and T and P data from, respectively, Glomfjord and Dunderlandsdalen climate stations. The two data groups are located close to each other and possess very similar overall hydroclimatic characteristics, although data group N2 is somewhat farther inland and its rivers are somewhat larger. Both are located in the region of Nordland, lying at the Arctic Circle but with a climate strongly moderated by immediate proximity to the Norwegian Sea. In low-elevation areas, the climate is predominantly maritime and characterized by mild summers and cold winters; continentality increases with elevation and distance from the coast. On average, precipitation is greater during fall and winter. Both N1 and N2 lie on the western flanks of the Scandinavian Mountains and drain immediately to the sea. Similarly to their Canadian counterparts, discharges in the Norwegian study rivers show a strong seasonal peak in summer, corresponding to snow and (for glacial rivers) ice melt, with low flows in winter.

Figure 2. Example of filtered daily data, over period of record used in the analysis, for data group N2. Data are presented as a raster hydrograph (e.g. Koehler 2013) or its meteorological equivalent. Horizontal-axis tick marks denote the start of every calendar month, beginning with 1 October (day of the hydrologic year, DOHY = 1); labels are provided for every third month. The figure shows ln(Q) and P^1/2, where Q is streamflow (m³/s) and P is precipitation (mm). These transforms are used only to facilitate effective visual display, and analyses were conducted on non-transformed data. Temperature, T, is given in °C. Some gaps or uneven record lengths are present in the T and P data. Overall hydroclimatic regimes for the other three catchment pairs considered are generally similar (see text for details).

ENSO and AO indices were obtained from the US National Oceanographic and Atmospheric Administration (NOAA) Climate Prediction Center. Recall that the AO is very closely related to the North Atlantic Oscillation (NAO); there is some consensus that the latter is a more regional expression of the former. Several ENSO indices are available. The Oceanic Niño Index (ONI), which is used by NOAA for operational classification of ENSO state (e.g. Kousky and Higgins 2007), was chosen here. Mean wintertime December-January-February (DJF) values were used, as is common practice, and were assigned to the corresponding hydrologic (i.e. JF) year.

Data processing

Hydrologic and climatic station data were organized into hydrologic years (HY), spanning 1 October of one calendar year to 30 September of the next, and assigned to the second calendar year. Day of the hydrologic year (DOHY) was consecutively numbered from the start of the hydrologic year. Thus, DOHY 18 of HY 2004 was 18 October 2003, for example; DOHY 168 of HY 1984 was 17 March 1984. For a given data group, the glacial and non-glacial flow records were then truncated to their common period of record to maximize comparability. Temperature and precipitation datasets within that data group were subsequently truncated to the same period. Data gaps were not interpolated.

A 21-point (n = 10) binomial filter (e.g. Mitchell et al. 1966), giving a 50% filter response at a timescale of 3n^1/2 ~ 9 days, was applied to each daily time series prior to analysis. Although results were not dramatically different from those obtained using shorter filter lengths or unsmoothed data (not shown here in the interest of conciseness), it was found that clearer analytical outcomes were obtained when this moderate low-pass filter was applied to the data prior to analysis. The resulting temporal resolution is roughly comparable to the 5-day sequent means used by Déry et al. (2009), for example, in their seasonally high-resolution analysis of hydroclimatic change. The filter was applied to daily data independently for each hydrologic year, starting at DOHY 1 and ending at DOHY 365; the first and last n days of data for a given HY were left unfiltered. The filter was not run across hydrologic years, so that year-to-year serial correlation was not artificially induced. All data within any window position of the 21-point filter which included a data gap were in turn transformed to data gaps, a conservative approach.

To be on the safe side, the analyses span the full hydrologic year. That said, interest here lies primarily with the late-summer to early-autumn period, when differences between glacial and non-glacial hydrologic regimes are usually most pronounced for mid- to high-latitude rivers in the northern hemisphere (e.g. Stahl and Moore 2006). The temporally high-resolution data employed here allow us to examine this narrow seasonal window in detail. There are also more fundamental motivations for using relatively high-resolution data. Often, the focus in statistical hydroclimatology lies with summary measures like seasonal and annual mean discharges. These metrics speak well to overall water supply issues, and the summation process involved in their estimation increases the signal-to-noise ratio. However, coarse-graining a time series (by temporally binning into lower-resolution averages, for example) can, in effect, simplify and linearize the system response (Penland 1996; Newman 2007; Hsieh 2009; Micovic and Quick 2009). This is convenient for some purposes but could in principle result in the loss of important nonlinear hydroclimatic signals, which are (as discussed above) of central interest to us, further militating for the use of non-aggregated daily data.

Methods

Spearman rank correlation

Teleconnections were assessed by correlating the hydroclimatic variate of interest, x(DOHY, HY), against a corresponding climate index, v(HY), independently for each DOHY using Spearman rank correlation: (1)

where y = 1, Y indexes the year of record considered, r^x and r^v are the ranked values of x and v respectively, and the overbar denotes the corresponding means. For example, the annual time series of Mistaya streamflow for DOHY = 1 was correlated against the annual time series of DJF ONI values; this was repeated separately for each DOHY = 2, 365. Corresponding p-values were also obtained. ONI was considered for the Canadian basins, and the AO index for the Norwegian watersheds.

For trend analysis, the same method was applied with the difference that the independent variable was time (specifically, HY) rather than the climate index, v(HY). As a trend detection method, Spearman correlation performs very similarly to the (also common) Mann-Kendall test (Dahmen and Hall 1990; Yue et al. 2002). The main focus of the trend analysis was to identify gradual monotonic shifts in the mean over the period of historical record common to the glacial-nival catchment pair, and assess similarities and differences in how the two distinct terrestrial hydrologic regimes respond to the same long-term forcing changes, whatever their sources may be. Thus, such additional steps as PDO signal removal or correction for serial correlation, for example, were not employed. Discussion of these complex topics is beyond the scope of this manuscript; for details, see for example Fleming and Weber (2012) and Fleming and Sauchyn (2013).

Spearman correlation is a robust, rank-based, distribution-free method and a standard tool of statistical hydroclimatology. Unlike the linear (Pearson product-moment) correlation coefficient, rank correlation is robust to mild nonlinearity. Like linear correlation, however, it has the limitation that it cannot be used to quantify the degree of non-monotonic association. For this reason, standard methods like correlation and composite analysis cannot be as effectively used for examining the parabolic teleconnections of particular interest to us here (Wu et al. 2005). In contrast, the AIC-based methods described below should readily accommodate non-monotonic teleconnections. The Spearman method is implemented for purposes of comparison.

AIC-based polynomial selection

As many hydrologists and climatologists may have limited prior exposure to AIC-based analysis and modelling, a brief summary as relevant to this study is provided below. Comprehensive reviews of the theory and application of AIC-based modelling in environmental science can be found in Anderson et al. (2000), Burnham and Anderson (2002), and Hobbs and Hilborn (2006). More mathematically formal treatments of AIC can also be found in the statistical literature.

The AIC has its roots in information theory. The most fundamental and broadly familiar aspect of information theory is the Shannon entropy:(2)

where the probabilities, p, are frequencies of occurrence for each of K distinct states in which the system of interest, captured by the discrete variate x, resides. Conceptually, it measures the quantity of information in a datastream, in terms of the degree to which uncertainty about the system is removed by monitoring the signal. The maximum-information case occurs if each state is equally probable, and therefore one has, in principle, no idea what the value of the signal will be without monitoring it. The minimum-information case occurs if only one state is present, and one could in principle know the value of the signal without monitoring it. The Kullback-Leibler (K-L) information (also K-L “distance”) possesses a similar form to the Shannon entropy:(3)

but there is additionally an approximating probability distribution, π. This is essentially the statistical distribution of the predictions from the model and, for a perfect model, π and p would be identical. The K-L information quantifies how much information is lost when approximating truth using a model. The K-L distance may be extended to the case of a continuous variate as:(4)

where f is full reality or truth, g is the model with parameters θ and I is the information lost when g is used to approximate f. Combining these information theoretic considerations with maximum likelihood concepts, Akaike (1973, 1974) developed an information criterion for model selection:(5)

where is the likelihood, represents the estimated model and its parameters, x are the data and λ is the number of parameters in the model. Assuming homoscedastic Gaussian model residuals, Equation (5) simplifies to:(6)

where Y is sample size and is RSS/Y. RSS is the residual sum of squares:(7)

where x_y and are the observed and model-simulated values for year, y. Note that λ is the number of parameters explicitly appearing in the model plus one, as in this context the fit metric (RSS) is viewed as an additional estimated parameter. For the small-sample case (Y/λ < ~40), a modified criterion is used:(8)

Note that AIC_c approaches AIC as the ratio of the sample size to the number of model parameters increases. In all cases here the small-sample approximation was used, which is hereafter referred to simply as AIC for notational simplicity.

The AIC is a model selection tool. From several models, g_m=1,M, each having different parameter values or model structures, the one having the lowest AIC value is considered best. By “best”, it is meant that model m provides the optimal combination of performance and parsimony: AIC increases with both model error (RSS) and model complexity (λ). That is, AIC-based model selection amounts to a rigorous mathematical formalization of Occam’s razor, and by construction avoids model over-fitting.

Importantly, estimated AIC values can also be extended to provide information about the conditional probability that a particular model is true, given the data. Amongst the various a priori candidate models that have been defined and supported as potential explanations for a particular dataset on the basis of physical reasoning, past experience or other considerations, model ranking can be accomplished by rescaling the calculated AIC values such that the best-performing model has a value of 0:(9)

The likelihood of a particular model given the data is proportional to exp (–Δ_m/2). These values may then be normalized to provide the Akaike weights:(10)

The value of is the approximate probability that model m is the Kullback-Leibler best model in the set of models considered. If this set of models adequately spans the range of physically plausible models for the situation at hand, then can in turn be treated as an approximate estimate of the probability that model m is true given the data. Additionally, the ratios of the Akaike weights to each other (/, i ≠ j), termed evidence ratios, provide a useful means for assessing the relative likelihood of different models.

AIC values are fairly routinely used in stepwise multiple linear regression modeling, including some hydrological applications. The probability information that can additionally be derived from AIC values has enjoyed far less widespread use but has nevertheless seen substantial application in some fields, such as ecology. However, it appears that neither AIC values nor (in particular) the model probability estimates provided by Akaike weights have been applied to routine tasks of statistical hydroclimatology, such as the detection of climate change signals or teleconnections, as an alternative to conventional methods like Spearman rank correlation as discussed above. The foregoing information theoretic approach was adapted to the needs of this study as follows.

For teleconnection analysis, the aim is to determine whether there is some relatively simple relationship between a hydroclimatic variable of interest (e.g. streamflow) and a climate index (e.g. the ONI). To this end, polynomial fits of varying degrees were used. Specifically, M = 3 models are physically plausible and consistent with prior work on teleconnections between ENSO or AO and Northern-Hemisphere temperature and precipitation:

where β are coefficients, and x(DOHY, HY) and v(HY) are as before. Note that a fundamental advantage of this approach over conventional null-hypothesis significance testing is that more than two models can be considered simultaneously, and that this can be done without making assumptions about which model must bear the burden of proof, except insofar as more complicated models must more than make up for relative lack of parsimony through improved simulation accuracy. The first model is a polynomial of degree 0 and simply amounts to the mean of the hydroclimatic variable, i.e. no relationship to the climate index (no-effect model). Model m = 2, a polynomial of degree 1, is equivalent to linear regression (linear effect model). The final model, a 2^nd-degree polynomial, is a quadratic function which can take on various graphical forms over a given v(HY) range, including the parabolic shape described by Wu et al. (2005) and Hsieh et al. (2006) for climate teleconnections (nonlinear effect model).

Wu et al. (2005) and Hsieh et al. (2006) reported finding no evidence for relationships corresponding to polynomials of degree higher than 2. As such, it is reasonable to assume that the three models above fully span the potential range of plausible univariate models relating a local hydroclimatic variable to a circulation index: none, linear and nonlinear of polynomial degree 2. In that case, the corresponding Akaike weights can be interpreted as approximate estimates of the probability that each model is true, given the data. The probability that any relationship exists between a hydroclimatic variate and a climate index is estimated here as p(^o = 1 or 2) = p(^o = 1) + p(^o = 2), where ^o is polynomial degree. Further, a nonlinearity ratio is defined here as the Akaike evidence ratio between Models 3 and 2, p(^o = 2)/p(^o = 1). When both p(^o = 1 or 2) and p(^o = 2)/p(^o = 1) are large, there is substantial evidence specifically for a nonlinear teleconnection. Note that while a 2^nd-degree polynomial with suitable coefficients can essentially reproduce the form of a 1^st-degree polynomial, the corresponding AIC is higher because this model is less parsimonious. Thus, even in the very specific case of a genuinely and strictly linear underlying relationship, the AIC-based polynomial selection method still distinguishes appropriately between Models 2 and 3.

The foregoing AIC-based methodology is altered for long-term trend analysis by replacing the independent variable in the polynomial models with time (specifically, HY). The AIC-based trend analysis methodology is able to detect both linear and nonlinear trends. By definition, however, interest lies solely with monotonic trend relationships, whether linear or nonlinear. This is a potential issue for the AIC-based polynomial selection method, as monotonicity is not guaranteed by a standard, unconstrained degree-2 polynomial fit (indeed, that is precisely why the method is used here for analysis of strongly nonlinear teleconnections as discussed above). The 2^nd-degree polynomial fitting exercise is therefore rephrased here as a constrained optimization problem:(11)

subject to:(12)

where are again the predicted values from the polynomial model, sgn is the sign function, and a(y) simply denotes any non-constant function of time. That is, while the local-in-time slope of the best-fit polynomial is permitted to vary (unlike the linear case), it is forced to be of constant sign for all y (and thus monotonically increasing or decreasing over the entire period of record). This constrained optimization was accomplished using a Nelder-Mead simplex direct search algorithm, with the addition of a penalty to term to the RSS cost function when the constraint is violated. Exploratory testing on Mistaya and Cataract data revealed no evidence for nonlinear monotonic trends. However, it was also found that fits could be sensitive to the polynomial parameter guesses employed for initializing the nonlinear optimization. As a check, therefore, the issue was further explored by removing the monotonicity constraint and simply performing a standard 2^nd-degree polynomial fit, which in principle should yield very liberal outcomes. Even in this case, however, little evidence was found for nonlinear relationships.

The foregoing does not imply that long-term monotonic trends in hydroclimate are inherently linear, but instead that they can be reasonably treated as linear for the record lengths available in this study. That is, a linear fit is sufficient for the purposes of this assessment, consistent with other recent findings in hydroclimatic trend detection (e.g. Fleming and Weber 2012). Thus, when using AIC-based polynomial selection for monotonic trend analysis, only two models were considered: 0^th- and 1^st-degree polynomials, i.e. the no-trend and linear-trend cases.

Results and discussion

A conservative interpretive framework was adopted, focusing on identifying key features in the resulting suite of 11,680 statistical analyses. In particular, a given hydroclimatic relationship is claimed only if the following criteria are generally satisfied: (1) consistency of that relationship across multiple climate or hydrometric stations, as applicable, within a given region, (2) overall consistency of the relationship with prior work, (3) physical consistency of dynamics across hydroclimatic variables for a given data group or region, (4) for the AIC-based teleconnection analyses, consistency of teleconnection form across multiple sites within a given region for a particular hydroclimatic variable, between physically related variables (such as P forcing and Q response) for a given hydrologic regime type in a given region, and/or with prior results for each region by Wu et al. (2005) and Hsieh et al. (2006) and (5) consistency, at least loosely, between Spearman rank correlation p-values, vs. AIC-based p(teleconnection) for the high-frequency analyses or p(trend) for the low-frequency analyses. The discussion below focusses on two main outcomes for summer-autumn flows: differential (glacial/non-glacial) linear trends, and differential (glacial/non-glacial) nonlinear teleconnections to ENSO in Canada and the AO in Norway.

Differential trends in late-season ecological flows

All four glacial rivers considered, in both Canada and Norway, exhibited strong and consistent flow declines in summer which were not consistently apparent in the corresponding non-glacial rivers. An example is provided in Figure 3 for data group N2. Glacier melt-influenced Bjørnfoss streamflows exhibit a consistent and strong decline over approximately mid-June through September (approx. DOHY 260–330) (Spearman p < 0.01; p(trend) > 0.90 and often near unity). In contrast, the largely non-glaciated Junkerdalselv catchment with which Bjørnfoss is paired shows no evidence for a decline over this period (Figure 3). Similarly, the glaciated Berget watershed from data group N1 demonstrates an analogous though less consistent or powerful decline in approximately late June (~DOHY 260–280) (Spearman p < 0.05; p(trend) ~ 0.70), whereas its pair (the non-glaciated Forsbakk catchment) again shows comparatively little evidence for a decline over this period. That is, there is a direct correspondence between declining late-season flow trends in the Norwegian study catchments and the presence of glacial cover.

Figure 3. Example of monotonic trend analysis results (data group N2). Left column: Spearman rank correlation (correlation coefficient, and p-value). Right column: Akaike Information Criterion (AIC)-based polynomial selection (linear trend rate, and approximate probability that trend is present). Heavy box outlines late-summer to early-autumn period.

The general worldwide pattern of glacial recession, beginning roughly at the nineteenth-century end of the Little Ice Age and perhaps accelerating through the late twentieth century due to anthropogenic climate change, has been identified as having different downstream hydrologic implications in different regions, apparently corresponding to different overall stages of glacial retreat (e.g. Fleming and Clarke 2003; Stahl and Moore 2006; Casassa et al. 2009; Moore et al. 2009; Baraer et al. 2012). In many cases, however, the effect has been a decline in glacier-fed streamflow, presumably due to a smaller contributing glacier area; this would be consistent with the current study’s findings for the glacierized Nordland catchments. It may also be loosely consistent with prior observations of declining summer flows in this study area by Stahl et al. (2010) and Førland et al. (2000).

It is interesting to compare glacier-driven summer flow declines against net glacier mass balance data. Engabreen is considered here, an outlet glacier of the Svartisen complex; this complex feeds both of the Norwegian glacial rivers considered here. Long-term balance data for Engabreen are provided by the World Glacier Monitoring Service (WGMS) Fluctuations of Glaciers (FoG) reports. A compilation of the available FoG data is provided in Figure 4. Annual balances have generally been positive over the period of historical record for Engabreen, unlike glaciers in many other parts of the world. Nevertheless, values have been more consistently neutral or negative over the last decade or more, and, overall, Engabreen has demonstrated a more globally typical negative mass balance trend over the available record, averaging about –23 mm/yr. The observation of declining summer flow trends for Berget and Bjørnfoss, but not for nearby non-glacierized Norwegian catchments, seems generally consistent with this declining mass balance trend.

Figure 4. Annual (net) glacier mass balance for Engabreen (Scandinavian Mountains, Nordland, Norway) and Peyto (Rocky Mountains, Alberta, Canada). Dashed lines provide simple linear regression trends for reference.

A similar but not identical pattern was observed in the Canadian results: strong late-summer flow declines were inferred, though the correspondence with presence of glacial cover is not as close as in Norway. A declining pattern was observed over August through to perhaps early September (approximately DOHY 305–340). The effect is significant to highly significant for Mistaya, Illecillewaet and Kootenay (Spearman p < 0.05 and in some cases < 0.01; p(trend) > 0.80), and nominally present though statistically insignificant for Cataract (Spearman p < 0.10, but p(trend) < 0.50). This trend may be mainly caused by a statistically significant decline (Spearman p < 0.05, and p(trend) often > 0.90) in precipitation throughout much of the year, including late summer. The precipitation declines may lead to both a reduction in snowpack, resulting in an earlier end to the melt freshet, and a loss of late-summer rainfall inputs. These discharge and (to a lesser degree) precipitation trends are consistent with prior work (Whitfield 2001; Zhang et al. 2000; Zhang et al. 2001), particularly if focus is placed on results from the most recent periods of record these workers considered, which are the most comparable to those in this study. The increasing T and decreasing P trends noted above also appear consistent with snowpack declines for southeast British Columbia reported by Rodenhuis et al. (2007). However, glaciers may also play a role in these late-summer baseflow declines in Canada. The statistical significance of post-freshet flow reductions in the glacial Illecillewaet basin persists later (into autumn) than for the paired, non-glacial Kootenay basin. Similarly, the late-summer flow reduction in the heavily glacierized Mistaya River is clearly statistically significant whereas that for the paired non-glacial Cataract Creek is not. That is, while late-season negative flow trends in these catchments are not uniquely associated with glacial cover as they were in Norway, the Canadian trends appear to be stronger, more seasonally persistent and more consistent basin to basin for glacial rivers relative to non-glacial rivers. This inference also generally agrees with prior hydrologic findings in and near the study area (Stahl and Moore 2006; Marshall et al. 2011; Fleming and Weber 2012). Additionally, the greater strength and consistency of late-season flow declines in glacial rivers in Canada is loosely consistent with generally negative, and downward-trending (–14 mm/yr), mass balances apparent in the FoG data for nearby Peyto glacier (Figure 4).

Differential teleconnections in late-summer ecological flows

Most of the rivers showed significant evidence for late-season teleconnections. However, glacial rivers were found to additionally exhibit strongly nonlinear teleconnections not present in non-glacial rivers.

An example is provided in Figure 5 for Canadian data group C1. The discussion here begins with monotonic associations, as identified by low Spearman p-values, or the combination of a high AIC-based probability of teleconnectivity and a low corresponding nonlinearity ratio. A negative correlation between ONI and Q is observed for both the glacial Illecillewaet and its non-glacial pair, the Kootenay, in late summer extending approximately into early autumn (Spearman p < 0.01). This mid-August to early-September flow decrease under El Niño conditions may arise from an observed negative winter precipitation correlation with ONI (Spearman p < 0.05 and sometimes < 0.01), which should produce anomalously low snowpack accumulation the preceding winter. Thus, seasonal snowpack disappears sooner so that the switch to groundwater-driven late-summer baseflow is earlier, giving highly significant late-season streamflow decreases. The converse occurs under cold-phase conditions. Overall, these findings are broadly consistent with prior work in the southern Canadian Rocky Mountains and nearby areas, although there is substantial variability between studies and rivers (e.g. Redmond and Koch 1991; Shabbar and Khandakar 1996; Shabbar et al. 1997; Lafrenière and Sharp 2003; Woo and Thorne 2003; Gobena and Gan 2006; Gan et al. 2007; Rodenhuis et al. 2007; Fleming and Whitfield 2010; Gobena et al. 2013).

Figure 5. Example of teleconnection analysis results (data group C1, El Niño-Southern Oscillation, ENSO). Left column: Spearman rank correlation (correlation coefficient, and p-value). Right column: Akaike Information Criterion (AIC)-based polynomial selection (approximate probability of any teleconnection, and evidence ratio for a nonlinear relationship). Heavy box outlines late-summer to early-autumn period.

However, the AIC-based results reveal an additional layer of complexity (see again Figure 5). Clear evidence is found for a strongly nonlinear streamflow teleconnection during late summer (roughly DOHY 300–310) for the glacial Illecillewaet River with p(teleconnection) ~1 and a nonlinearity ratio approaching 20, but not for the non-glacial Kootenay River, which has p(teleconnection) again ~1 but a nonlinearity ratio of ~0. The Illecillewaet’s August nonlinear Q relationship to ENSO is physically consistent with a contemporaneous (late summer) nonlinear T relationship to ENSO visible in Figure 5. Specifically, there is a quadratic relationship between ENSO state and T, such that strong El Niño and strong La Niña episodes both tend to give elevated temperatures (see scatterplots and polynomial fits in Figure 6). In August, when seasonal snowpack has been depleted and glacial ice is exposed, this increased T generates greater glacial melt production and streamflow. No corresponding flow response is seen in the non-glacial Kootenay River, because seasonal snowpack has been depleted by this point in the year and there is no glacial ice potentially available for melting; only the monotonic teleconnection apparently associated with winter snowpack modulation is apparent (see above). Note that the form of the ENSO-temperature relationship is consistent with that found by Wu et al. (2005) and Hsieh et al. (2006), though they only considered wintertime effects. There does not appear to be a corresponding late-summer precipitation teleconnection (Figures 5 and 6).

Figure 6. Example scatterplots for teleconnections between data group C1 and an El Niño-Southern Oscillation (ENSO) index (the Oceanic Niño Index, ONI). Results for early August are shown. M1, M2, and M3 correspond to 0^th, 1^st and 2^nd degree polynomial models; these in turn reflect no-effect, linear teleconnection and nonlinear teleconnection outcomes, respectively. Corresponding correlation skills, R, for each model are also provided; values are generally modest, reflecting the fact that any given climate mode explains only some limited portion of year-to-year variability. The most likely model, using probability estimates based on Akaike weights formally balancing model performance against model parsimony, is noted in each case. For this particular example, nonlinear teleconnections are identified as the preferred model for air temperature and glacial river flows, whereas a non-glacial streamflow teleconnection is present but linear, and precipitation exhibits no statistically defensible relationship to the ONI.

Similar ENSO relationships were obtained for data group C2; the nonlinearity ratios for both Mistaya Q and Banff T were in fact considerably higher than seen in C1, further reinforcing the foregoing interpretation. Further, the late-summer, strongly nonlinear, glacier-mediated hydrologic teleconnection to ENSO in Canada appears (in terms of watershed processes) to be identical to late-summer, strongly nonlinear, glacier-mediated hydrologic teleconnections to the AO in Norway. In late summer, a nonlinear flow relationship with the AO index was found for the Norwegian glacial rivers, which appears to be physically consistent with a contemporaneous nonlinear temperature relationship to AO. Specifically, there is a parabolic association between AO state and T in this region, such that strong positive-AO and strong negative-AO years both tend to be warmer. There is no evidence for nonlinear contemporaneous P responses. In August, when seasonal snowpack has been depleted and glacial ice is exposed, increased T yields greater glacial melt production and increased Q in both of the glacial rivers. In contrast, no corresponding nonlinear flow response is seen in the unambiguously non-glacial catchment (Forsbakk), because seasonal snowpack has been depleted and no glacial ice is potentially available for melting. A complication is that the other nominally non-glacial catchment (Junkerdalselv, with 0.6% glacial cover) shows a strongly nonlinear response similar to the glacial rivers. This discrepancy might reflect late-season, high-elevation snowmelt in the slightly more continental Junkerdalselv basin, or it might suggest that even very low levels of watershed glacial cover may be sufficient to induce a distinct hydroclimatic signal in this particular respect. It perhaps bears mentioning that there might in some sense be a partial ambiguity between the presence of glaciers and the presence of cold, wet conditions that tend to produce and sustain glaciers. That is, the same type of “snow-friendly” environment which yields catchment glacierization might also give seasonal snowpack lasting into late summer in most years, even for catchments with very limited glacial cover per se – and hence a late-summer AO sensitivity that somewhat lower, less snowy catchments do not exhibit. Note that the form of the AO-T relationship is fully consistent with that found by Hsieh et al. (2006), though they only considered wintertime effects. Note also that while glacier-mediated selective streamflow teleconnectivity to the AO has been previously identified elsewhere (Fleming et al. 2006), the teleconnections apparent here are strongly nonlinear.

In principle, such seasonally short-lived and strongly nonlinear teleconnections might be missed without the use of temporally high-resolution data (see discussion and references in Data section above). Confirmation was provided by re-running the ENSO analyses for the Canadian watersheds using unfiltered monthly mean time series in place of daily data. Nonlinearity in the short-duration, late-summer Illecillewaet glacial melt response could no longer be detected, and the corresponding Mistaya teleconnection disappeared altogether, under the much coarser sampling interval. To summarize, while temporally (monthly, seasonally or annually) integrated time series should provide a clearer picture of broad water supply conditions due to additive noise cancellation, the higher-resolution data employed here better capture variability and change in detailed catchment processes and dynamics, including nonlinear effects.

Other effects

Though the main interest in this study lies with late-season differences in the hydroclimatic dynamics of glacial vs. non-glacial rivers, other phenomena were also identified. Two of the more noteworthy effects are summarized here. First, all the Norwegian rivers showed strong positive flow trends in winter and early spring (roughly December through mid-April; see Figure 3). These Q increases are consistent with contemporaneous T and P increases, most clearly visible for data group N2. Note that high latitudes notwithstanding, manual data inspection suggested winter temperatures can rise above freezing here, reflecting moderating maritime influences. These wintertime shifts to warmer, wetter conditions are broadly consistent with prior literature on the region (Korhonen and Kuusisto 2010; Wilson et al. 2010 [particularly over 1961–2000, the most recent period they consider and most closely matched to the period of record used here]; Førland et al. 2000; Nesje et al. 2000; Dyrrdal and Vikhamar-Schuler 2009; Førland et al. 2009; Stahl et al. 2010; Lehmann et al. 2011), and the flow changes in particular are obvious even from simple visual inspection of the data (see Figure 2).

Second, nonlinear winter-spring T and P teleconnections were inferred in Norway (to the AO) and Canada (to ENSO), but consistency with Wu et al. (2005) and Hsieh et al. (2006) is partial and geographically and seasonally variable, and implications for Q variability appear complex. In summary, winter-spring T in the Nordland basins exhibits only a monotonic, approximately linear (increasing) relationship with the AO index, consistent with prior work on NAO responses in the region (e.g. Chen and Hellström 1999; Nesje et al. 2000; Kingston et al. 2006; Arguez et al. 2009; Kingston et al. 2009; Wrzesiński and Paluszkiewicz 2011; Dahlke et al. 2012). In contrast, quadratic forms were favoured by AIC polynomial selection for winter-spring P teleconnections. However, manual inspection of the results suggested that the teleconnection form may vary within the season, switching from a nonlinearly but monotonically increasing P relationship with the AO index (again generally consistent with prior work using techniques assuming monotonic relationships), to a parabolic relationship in which P increases with both very positive and very negative AO conditions (fully consistent with Hsieh et al. 2006). The net outcome appears to be that relatively warm, wet (cool, dry) mid-winter weather associated with positive- (negative-) phase AO conditions induces a winter Q increase (decrease); this response appears to shift to a modest but more distinctly parabolic Q relationship in spring. In Canada (see again Figure 5), highly nonlinear ENSO teleconnections were also found in winter T and P, which both tend to increase during extreme El Niño and extreme La Niña events. The results are strongly consistent with Wu et al. (2005) and Hsieh et al. (2006). However, there is little indication of this nonlinear climate signal in winter Q data. Though further south than the Nordland basins, the Rocky Mountain data groups lie at higher elevation and are farther inland from moderating maritime climatic influences. During mid-winter, temperatures seldom venture above freezing, and there seems little reason to expect that mid-winter ENSO signals in P and T have a significant contemporaneous effect on streamflow, unlike in the Norwegian catchments. Interestingly, there is no clear evidence for a nonlinear Q effect during the freshet either, as might be anticipated given the general snowpack implications of winter P anomalies. Reasons for this negative result are unclear, but a potential explanation is that the temporal integrating effect of seasonal snowpack accumulation may damp more subtle, nonlinear components of the winter climate input signal in this particular case. Overall, the details of nonlinear wintertime teleconnections are complex, and clearly demand further investigation.

Conclusions

Glacial ice plays a major role in the headwaters to many of the world’s major river systems, and hints of hydrologically relevant and strongly nonlinear climate variations have appeared in the prior literature. It was therefore assessed in this study how the presence or absence of watershed glacial cover affects streamflow responses to both low- and high-frequency climatic variations, with a particular emphasis on the consideration of nonlinear effects. An information theory-based approach was adapted to this purpose and applied to seasonally high-resolution data from paired catchments in the southern Canadian Rocky Mountains and arctic coastal Norway.

Two key hydrologic responses were as follows. First, in both the Norwegian and Canadian study watersheds, declining late-summer flows were identified and found to be stronger, more seasonally persistent and more consistent for glacier-fed rivers than for non-glacial basins. Although the results may reflect a number of region- and catchment-specific factors, progressive glacial area loss appears to be the primary contributor in both regions. Second, glacier-controlled, highly (non-monotonically) nonlinear streamflow teleconnections were identified in both Norway and Canada during late summer to early autumn. Parabolic associations between local air temperature, and the Arctic Oscillation (Norway) or El Niño-Southern Oscillation (Canada), yielded corresponding parabolic teleconnections in glacial melt production. The non-glacial rivers, lacking a glacial ice reservoir potentially available for melting, gave little or no evidence for highly nonlinear AO or ENSO teleconnections in late summer or early autumn. Taken together, these two outcomes reveal a decoupling between climatic forcing and hydrologic response, reflecting the modulating impacts of local-scale terrestrial hydrologic characteristics. That is, the presence and nature of water resource responses to a given type of climatic forcing can vary from river to river within an otherwise approximately uniform region, depending on the local presence or absence of upstream glacial cover.

Little quantitative support was found for nonlinearity in long-term trends, but as noted above there was considerable evidence for moderately to highly nonlinear teleconnections. These outcomes may help guide methodological choices in future studies. Additionally, the AIC-based polynomial-selection approach to trend and teleconnection detection introduced here proved effective, offering advantages over conventional null-hypothesis significance testing methods by simultaneously accommodating multiple candidate models, permitting easy identification of strongly nonlinear relationships, avoiding assumptions about assignment of null and alternative hypotheses and corresponding burdens of proof, and directly returning an estimate of the probability that the effect of interest is present. Finally, seasonally high-resolution data facilitates the study of nonlinear or transient effects, though this may come at the cost of compromised signal-to-noise ratios.

Acknowledgements

The authors thank Alex Cannon (Pacific Climate Impacts Consortium) for bringing Anderson et al. (2000) and Wu et al. (2005) to their attention, Lucy Vincent and Éva Mekis (Environment Canada) for providing daily AHCCD time series, Tommy Diep and Judy Kwan (Environment Canada) for geographic information system (GIS) support and Trine Fjeldstad (Norwegian Water Resources and Energy Directorate) for providing Norwegian streamflow data. Dan Moore (University of British Columbia) and Gia Destouni (Stockholm University) provided useful suggestions on an earlier version of this manuscript, as did two anonymous reviewers, an anonymous CWRJ associate editor and the CWRJ co-editor, Paul Whitfield.