Examining controls on peak annual streamflow and floods in the Fraser River Basin of British Columbia

The Fraser River Basin (FRB) of British Columbia is one of the largest and most important watersheds in western North America, and home to a rich diversity of biological species and economic assets that depend implicitly upon its extensive riverine habitats. The hydrology of the FRB is dominated by snow accumulation and melt processes, leading to a prominent annual peak streamflow invariably occurring in May–July. Nevertheless, while annual peak daily streamflow (APF) during the spring freshet in the FRB is historically well correlated with basin-averaged, 1 April snow water equivalent (SWE), there are numerous occurrences of anomalously large APF in below- or near-normal SWE years, some of which have resulted in damaging floods in the region. An imperfect understanding of which other climatic factors contribute to these anomalously large APFs hinders robust projections of their magnitude and frequency.

We employ the Variable Infiltration Capacity (VIC) process-based hydrological model driven by gridded observations to investigate the key controlling factors of anomalous APF events in the FRB and four of its subbasins that contribute nearly 70 % of the annual flow at Fraser-Hope. The relative influence of a set of predictors characterizing the interannual variability of rainfall, snowfall, snowpack (characterized by the annual maximum value, SWE_max), soil moisture and temperature on simulated APF at Hope (the main outlet of the FRB) and at the subbasin outlets is examined within a regression framework. The influence of large-scale climate modes of variability (the Pacific Decadal Oscillation (PDO) and the El Niño–Southern Oscillation – ENSO) on APF magnitude is also assessed, and placed in context with these more localized controls. The results indicate that next to SWE_max (univariate Spearman correlation with APF of $\widehat{\mathit{\rho }}$ = 0.64; 0.70 (observations; VIC simulation)), the snowmelt rate ($\widehat{\mathit{\rho }}$ = 0.43 in VIC), the ENSO and PDO indices ($\widehat{\mathit{\rho }}$ = −0.40; −0.41) and ($\widehat{\mathit{\rho }}$ = −0.35; −0.38), respectively, and rate of warming subsequent to the date of SWE_max ($\widehat{\mathit{\rho }}$ = 0.26; 0.38), are the most influential predictors of APF magnitude in the FRB and its subbasins. The identification of these controls on annual peak flows in the region may be of use in understanding seasonal predictions or future projected streamflow changes.

1 Introduction

1.1 Study domain and motivation

The response of nival watersheds across the Northern Hemisphere to ongoing warming of the climate system is a topic of intense interest to those interested in the water–climate nexus. Since the early reviews of Barnett et al. (2005) and Milly et al. (2005), which offered a global perspective, more recent studies have investigated the response on continental and regional scales. In mountainous regions that supply fresh water to downstream populated areas, reduced snowpack in recent decades has led to concerns about water availability for human consumption, agriculture and fisheries (Stewart, 2009; Jiménez Cisneros et al., 2014), while changes in the timing of the freshet and a widespread decrease in the snow-to-rain ratio (Danco et al., 2016) have implications for flood magnitude and frequency (Matti et al., 2016). An understanding of the mechanisms behind these changes has been greatly aided by the use of process-based hydrologic models (Park and Markus, 2014; Duan et al., 2017), which permit the analysis of a host of variables that respond to historical climate forcings via physically consistent relationships. We apply such a model to the Fraser River Basin (FRB) of British Columbia (BC), Canada, a large, representative nival watershed, to study the dominant climatic controls on both observed and modelled hydrologic change.

The FRB is one of the largest watersheds draining the western slopes of the North American Cordillera, and home to both densely populated urban centres and a diversity of ecosystems closely linked to the economic prosperity of the region (Fraser Basin Council, 2010). The basin is also a focal point for salmon migration and spawning, constituting a key component of the highly valued commercial and subsistence fisheries of the region. The FRB lies between the Coast Mountains and the Continental Divide, with the headwaters of the Fraser in the northeast (∼ 53^∘ N, 118^∘ W) and its outlet at the Pacific Ocean in the southwest (Fig. 1). Its vast 240 000 km² area encompasses a range of climatic zones, from the snowy mountains of the eastern Rockies to dry interior plateaus and wet fertile valleys nearest the Pacific west coast. The hydrological diversity of the basin is well described in Eaton and Moore (2010), who reviewed seasonal streamflow regimes at the catchment scale across BC, and also in Moore (1991), who focused specifically on the FRB. For the most part, streamflows are unregulated in the FRB, with the exception of Kenney Dam on the Nechako River in the northwestern sector of the FRB (operational since the early 1950s) and several dams associated with the Bridge River power project in the Interior Plateau (completed in 1960). Regulated subbasins represent less than 10 % of the total area of the FRB (Bawden et al., 2015).

Streamflow in the interior catchments of the FRB is dominated by the snowmelt-fed spring freshet in April–July, leading to the usual characterization of the FRB as a nival basin. Hydrographs of major rivers in the FRB do not vary their form significantly from year to year due to the large amount of storage in the multi-year, high-elevation snowpack. Indeed, the longest record of gauged flow at the southwestern outlet of the basin, at Fraser-Hope, has never recorded an annual peak daily flow (APF) outside of the April–July period. Nevertheless, several catchments in the lower FRB, and to the west of the basin in the Coast Mountains and Western Cascades, exhibit annual streamflow peaks coinciding with the maximum of Pacific frontal rainstorm occurrence in October–December (Eaton and Moore, 2010; Padilla et al., 2015).

Figure 1 The Fraser River Basin of British Columbia, Canada, and four of its subbasins examined in this study. Colours correspond to elevations as used by the VIC model, at 1∕16^∘ horizontal resolution. Locations of streamflow gauge stations at the outlet to each subbasin are shown with black dots, while the major outlet for the FRB at Hope is shown with a red dot. Details of each subbasin are provided in Table 1.

1.2 Previous studies of climate and streamflow change in the FRB

The focus of this work is upon the main drivers of the interannual (and to some extent, interdecadal) variability in annual peak daily flow (hereafter APF), and not upon secular trends in APF in the FRB and its subbasins. Nevertheless, as there is a growing body of work on such trends, which have implications for non-stationarity of the climate, we briefly review this research here, limiting the discussion to observations concerning APF magnitude. Significant long-term trends in temperature and precipitation have been detected over a network of meteorological stations in the FRB and its surroundings; a summary of these trends is available in a recent report from the BC Ministry of the Environment (BC MOE, 2016) (in what follows, we average results for three regions used in the report which cover most of the FRB area: Sub-Boreal Interior, Central Interior, and Southern Interior). Annual (summer; winter) mean temperatures have increased at the rate of +0.9 ^∘C century⁻¹ (+0.7; +1.4 ^∘C century⁻¹) over the 1900–2013 period, while minimum winter temperature has increased at the more rapid rate of 2.2 ^∘C century⁻¹. The latter coincides with reduced winter snow depth (−10 % decade⁻¹) and snow water equivalent (−6 % decade⁻¹; both over 1950–2014), mirroring a much larger-scale phenomenon occurring over Canada and the US (Gan et al., 2013; Knowles, 2015) and much of the Northern Hemisphere (Rupp et al., 2013; Jeong et al., 2017). On the other hand, precipitation amounts have increased in the FRB over 1900–2013, by about +2 % decade⁻¹ in spring, summer, and fall, with no significant change in winter precipitation (BC MOE, 2016). However, sparse station coverage in the early decades of the twentieth century casts doubt on the robustness of this result. Restricting attention to the period 1950–2014, we used the regional data to compute precipitation trends of +3.3 % decade⁻¹ in spring, +2.3 % decade⁻¹ in fall and −4.1 % decade⁻¹ in winter, with no significant trend detected in summer. Over the last century (1912–2012), an earlier onset of the half-annual streamflow volume date was detected in one subbasin of the FRB (Stellako, 9-day advance) as well as at the major outlet to the basin at Hope (6-day advance). The latter conclusion was reinforced by Kang et al. (2016), who noted a statistically significant advance, by ∼ 10 days, of the APF at Fraser-Hope over the 1949–2006 period.

Recent decadal trends in the magnitude of annual and/or seasonal mean streamflow at gauge stations in BC and western Canada have been investigated by several authors (Pike et al., 2010; Bawden et al., 2015; DeBeer et al., 2016; BC MOE, 2016), with no significant trends detected at stations in the FRB. Hernández-Henríquez et al. (2017) conducted analyses on the FRB specifically, using a network of long-term measurements at 139 streamflow gauge stations available from the online Hydrometric Database (HYDAT; Water Survey of Canada, 2016). At Fraser-Hope specifically, Déry et al. (2012) found an increasing trend in interannual streamflow variability between 1960 and 2005, in both annual and seasonal means, suggesting an increase in the frequency and/or intensity of low- and/or high-flow conditions, but with no detectable trend in annual mean streamflow at Hope.

With respect to trends in the APF specifically, Cunderlik and Ouarda (2009) found no significant trend at most Reference Hydrometric Basin Network (RHBN) stations in the FRB over the 1974–2003 period, but noted a weak trend (∼ 1–5 %) toward decreasing magnitude and earlier APF occurrence at two stations. A Canada-wide nonstationary analysis by Tan and Gan (2015), suggested that APF increased over recent decades at two stations in the FRB: the Stellako R. at Glenannan and Chilliwack R. at Chilliwack Lake. More aligned with our interests here, Burn and Hag Elnur (2002) found a weak negative correlation between APF and annual mean temperature at Quesnel in the FRB over the 1950–1997 period. We look for similar relationships between a range of climatic variables and APF at the scale of both the FRB and its major subbasins in Sect. 3 below.

Table 1 Water Survey of Canada (WSC) hydrometric stations and corresponding basin and subbasin characteristics.

Finally, it is important to recognize the influence of large-scale climate teleconnections on streamflow in the FRB. Previous researchers have investigated the influence of modes of large-scale climate variability on various river basins in western North America. Specifically, relationships between total and/or peak annual river discharge and the Pacific Decadal Oscillation (PDO), El Niño–Southern Oscillation (ENSO), and the Pacific North America index have been examined (Shabbar et al., 1997; Rood et al., 2005; Gobena and Gan, 2006; Bonsal and Shabbar, 2008; Whitfield et al., 2010; Gurrapu et al., 2016), while impacts of large-scale teleconnections on snowpacks in BC were examined by Moore and McKendry (1996) and Hsieh and Tang (2001). El Niño and La Niña typically occur every 3 to 5 years, often separated by 1 to 2 years of neutral conditions. The lower frequency PDO is a superposition of several interacting climate processes, including ENSO, mid-latitude ocean currents and atmospheric influences on mid-latitude sea surface temperatures (Alexander, 2010). El Niño events are more likely to occur during the positive phase of the PDO, while La Niña events are more common during the negative phase.

During the negative (cold) phase of the PDO, and La Niña periods, winters in western Canada are typically cooler and wetter than average, with a larger snowpack at high elevations leading to higher annual discharge than average. Roughly opposite behaviour occurs during the positive (warm) PDO phase and occurrences of El Niño. Extending earlier work by Woo and Thorne (2003) and Thorne and Woo (2011), Gurrapu et al. (2016) determined that the PDO index is significantly anti-correlated (p < 0.05) with APF at eight hydrometeorological stations in the FRB, while APF is correlated with the Southern Oscillation Index (SOI), which tracks ENSO phase, at 11 stations (note that SOI is anti-correlated with the NINO3.4 index, and so also the PDO index). Specifically, both Thorne and Woo (2011) and Gurrapu et al. (2016) found that the observed APF was generally higher during the negative PDO phase and during La Niña years. For this reason, decadal trends in climatic variables, including streamflow, are sensitive to the phase of the PDO, which argues for a cautious interpretation of the results cited above.

The paper is structured as follows. Data sources, the VIC model and methods of analysis are described in Section 2. The main results of the study are gathered in Sect. 3, which begins with the regression analysis of observations before presenting insights from the VIC simulation results. Section 4 presents a few case studies that illuminate the precursors of high streamflow in particular years, and also reinforce the regression-based results. We conclude in Sect. 5 with a short discussion of outstanding issues and conclusions.

2 Data and methods

2.1 Study domain and observational data

The location of the FRB within western Canada and BC is shown in Fig. 1. The observational data set used for analysis and for driving the hydrological model (see Sect. 2.2 below) is taken from the gridded data set of surface temperature and precipitation at daily temporal and 1∕16^∘ spatial resolution described by Schnorbus et al. (2011), hereafter referred to as PCIC-OBS. The original station data, which span the period from January 1950 to December 2006, are interpolated to the grid and corrected for elevation using the Climate Data for Western North America (Climate WNA) package (https://sites.ualberta.ca/~ahamann/data/climatewna.html; Hamann et al., 2013), based on very high resolution climatologies for the region that are developed with the Parameter-elevation Regressions on Independent Slopes Model (PRISM; Daly et al., 2008). Precipitation is not partitioned into rain and snow separately in PCIC-OBS. Since we are interested in investigating the relationship between rainfall and APF specifically, we need to estimate the rainfall-only portion, R. We do this using the empirical fit of Dai (2008), which relates the fractional rain frequency to surface temperature via a hyperbolic tangent function having four fitted parameters. We chose the parameter values corresponding to land-only, annual mean precipitation (seasonal fitted coefficients were also given by Dai (2008), but do not differ much from the annual values).

Daily streamflow data were obtained from the Water Survey of Canada (WSC) Hydrometric Database (HYDAT; Water Survey of Canada, 2016) for five hydrometric stations located within the FRB, as summarized in Table 1. The main outlet at Hope, which integrates the flows from all upstream locations, receives the most attention in the paper, but four subbasin outlets are also considered. Three of these were selected on the basis of their leading contributions to the observed mean annual discharge at Fraser-Hope: Upper Fraser (29 %), Thompson–Nicola (28 %), and Quesnel (9 %) (Kang et al., 2016). These subbasins are located in the eastern FRB, cover 45% of the total area, and represent nival environments. The smaller Chilko basin, by contrast, lies on the southwestern edge of the FRB and intercepts a significant amount of rain falling on the east-facing side of the Coast Mountains, making it a hybrid (nival–pluvial) catchment. Hence, this subbasin was included in an attempt to probe the sensitivity of streamflow to rainfall. Manual snow survey (MSS) SWE measurements taken at the beginning of each non-summer month (eight times per year) were obtained from the BC Snow Survey Network Program, distributed by the BC River Forecast Centre (http://bcrfc.env.gov.bc.ca/data/). The data do not permit the exact determination of the annual maximum SWE, so we extracted the 1 April SWE from each year for analysis. Data from 19 locations that are at least 81 % complete spanning the period 1956–2014 were averaged to obtain an annual 1 April SWE time series for the entire FRB. The MSS stations are mainly located at high-elevation sites in the Coast and Rocky Mountain ranges; exact locations are shown in Fig. 1 of Najafi et al. (2017). Another source of snow cover data, from automated snow pillows at high-elevation sites, was also examined. However, as only a few of these sites have records longer than 20 years (within the 1956–2006 period of PCIC-OBS), we decided to exclude them from the analysis.

Table 2 Predictor variables and descriptions including their origin. Acronyms are as defined in the text, with the exception of ESRL/GCOS, Earth System Resource Laboratory/Global Climate Observing System (National Oceanic and Atmospheric Administration, NOAA), and JISAO, Joint Institute for the Study of the Atmosphere and Ocean (University of Washington/NOAA).

We use the PDO and NINO3.4 indices to characterize the relevant large-scale climate modes affecting the region. The PDO index is derived as the leading principal component of monthly sea surface temperature (SST) anomalies in the North Pacific Ocean northward of 20^∘ N (http://research.jisao.washington.edu/pdo/PDO.latest.txt). For the corresponding predictor variable, we use the mean PDO index from the preceding November to March, following Gurrapu et al. (2016). The NINO3.4 index is calculated from the Hadley Centre SST and sea ice gridded data set, HadISST1, as the area average of SST from 5^∘ S–5^∘ N and 170–120^∘ W, available from http://www.esrl.noaa.gov/psd/gcos_wgsp/Timeseries/Nino34/. We use the mean value from June to November of the preceding year as a predictor. Note that this differs slightly from Gurrapu et al. (2016), who instead made use of the SOI over the same months (NINO3.4 and the SOI are strongly negatively correlated). A list of all variables analyzed may be found in Table 2.

2.2 Hydrological model and boundary forcing

The gridded temperature and precipitation data described above were used to drive the Variable Infiltration Capacity (VIC) hydrological model (Liang et al., 1994) over a large portion of BC, including the FRB, from 1950 to 2006. The VIC model is applied at a horizontal resolution of 1∕16^∘ (∼ 5–6 km, depending on latitude), and solves the one-dimensional water and energy balance equations within each grid cell at a daily time step (with the exception of the snow sub-model, which runs at hourly resolution). The VIC implementation used in this study incorporates five elevation bands corresponding to 200 m vertical resolution, with the number of elevation bands in any one 1∕16^∘ grid cell depending on the topography within that cell. Each VIC grid cell can be assigned up to eight major vegetation classes, with a fractional cell area assigned to each, and these land cover fractions are identical for each elevation band within a VIC grid cell. Three of the input PCIC-OBS fields, maximum and minimum temperature and precipitation, are vertically interpolated to provide values for each elevation band (the other input field, surface wind speed, is not interpolated). In VIC, precipitation type (rain or snow) at the grid scale is determined by fixed temperature thresholds, with snow falling when T < 0 ^∘C, rain when T > 6 ^∘C, and a linearly interpolated mix of the two at intermediate temperatures. While this differs somewhat from the partitioning applied to the PCIC-OBS precipitation (Sect. 2.1) used in the regression analysis of Sect. 3 for both observations and models, the treatments are sufficiently similar that we do not expect a significant bias to arise. Snowpack is represented by a two-layer scheme – a thin surface overlying a thick deeper layer – subject to mass and energy balance, similar to other cold land process models. The snow albedo parameterization in VIC is based entirely on snow age, with a different albedo decay rate used during accumulation and snowmelt seasons. Interested readers are directed to Andreadis et al. (2009) for further details.

Glaciers are not explicitly parameterized in VIC; however, the version of VIC used here produces a perennial, accumulating snowpack at several high-elevation grid cells, many of which coincide with the locations of glaciers in the real system (also noted by Islam and Déry, 2017, who used a version of VIC with lower horizontal resolution but finer vertical resolution). These anomalous “glacier” cells cause the annual maximum of the basin average SWE, SWE_max, to increase approximately linearly with time t: i.e., SWE_max = at + SWE₀, where SWE₀ is the value of SWE at an affected cell at the start of the simulation (t = 0) and a is a fitted trend. The linear increase in SWE_max is clearly unrealistic and is not seen, e.g., in the corresponding observed 1 April SWE time series averaged over the FRB, which has a trend indistinguishable from zero. In order to reproduce this basin-averaged trend of approximately zero in the VIC model, we masked out the affected cells according to the criterion a > 60 mm yr⁻¹, with the slope a determined from a linear least-squares fit over the 1950–2006 period of the calibrated VIC simulation (see below). In addition, we confirmed a posteriori that the removal of these anomalous cells had little effect on the correlation structure of the basin average SWE_max with any of the variables examined in Sect. 3 below.

Soil parameters in VIC are defined for each grid cell. Soil classification and parameterization are based primarily on physical data from the Soils Program in the Global Soil Data Products CD-ROM (GSD, 2000), which originate from a global pedon database produced by the International Soil Reference and Information Center (ISRIC) (Batjes, 1995) and the FAO-UNESCO Digital Soil Map of the World (DSMW) (FAO, 1995). Although VIC contains a parameterization for seasonally frozen soils (Cherkauer and Lettenmaier, 1999), it increases the computation time significantly and so was not activated (the simulation we analyzed was originally produced for a different purpose). Surface and subsurface runoff are generated for each grid cell, and subsequently directed into a surface routing network (Lohmann et al., 1996, 1998; Schnorbus et al., 2011). Simulated streamflow can be extracted at grid cells representing outlets of the FRB or any one of its subbasins, which can be evaluated against stream gauge measurements.

VIC has been calibrated and evaluated in the FRB and its subbasins (Schnorbus et al., 2010; Shrestha et al., 2012, 2014) and also in other nearby hydrological basins (Schnorbus et al., 2011). Recently, Islam and Déry (2017) used a lower horizontal resolution version of VIC (1∕4^∘) to study its sensitivity to several different gridded input data sets, including PCIC-OBS. They found that while driving VIC with PCIC-OBS tended to overestimate SWE, the resultant hydrographs for the FRB and its basins were in better agreement with observations than when competing data products were used as driving data sets. VIC model output variables, including SWE, total column soil moisture (both over the period 1950–2006), and routed streamflow at the location of WSC stations (1955–2004), were obtained from the Pacific Climate Impacts Consortium Data Portal (Pacific Climate Impacts Consortium, 2014).

2.3 Selection of predictors and analysis methods

We use regression models to study the relationship between APF at Fraser-Hope and climate variables, termed “predictors”, that may influence APF. Predictors were selected based on physical intuition, inspection of the relevant literature and initial exploratory data analysis. Prior studies that were particularly helpful in this regard were Gurrapu et al. (2016), who examined the influence of the PDO on streamflow in western Canada, Jenicek et al. (2016), who examined the influence of snow accumulation and other variables on summer low flows, Coles et al. (2017), who studied snowmelt-runoff generation on Canadian prairie hillslopes, and Wever et al. (2017), who conducted model simulations of the joint effect of snowmelt and soil moisture on streamflow in a Swiss Alpine catchment. The predictor variables chosen are listed in Table 2.

Both univariate and multivariate linear regression models were constructed to explore predictor–predictand relationships. The APF was determined as the maximum annual value of the running 3-day mean discharge at Fraser-Hope. Most of the predictors are represented as basin averages, but two, the ENSO and PDO indices, describe large-scale, multi-year climate modes of variability. The influence of large-scale climate can be manifested as nonstationary behaviour (e.g., trends) in the predictand and/or predictors. Since regression models are sensitive to both trends and autocorrelation in the underlying time series, all univariate regressions were checked for the presence of both (e.g., see Shumway and Stoffer, 2017, and specifically their online supplement, https://onlinecourses.science.psu.edu/stat510/node/53).

Specifically, for the predictand (APF) Y_t and each predictor X_t, we first fit the models (by linear least-squares regression):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{\mathit{Y}}_t={\mathit{\alpha }}_{\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{1}}t+{\mathit{y}}_t,{\mathit{X}}_t={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}t+{\mathit{x}}_t,\end{array}$$

where y_t and x_t are the corresponding detrended annual time series. If, as determined from the fit, there was no statistically significant temporal trend (with p-value < 0.05) in either Y_t or X_t, then the latter were used in the subsequent analysis in place of y_t and x_t. Trends were detected in certain time series, as reviewed in Sect. 3.1 below. Regression models including variables with significant trends were generally found to have inflated correlation coefficients (e.g., Pearson's R² or Spearman's $\widehat{\mathit{\rho }}$) compared to those using the same variables after detrending. The analysis proceeds using the fitted linear regression model to the detrended series y_t and x_t, i.e., ${\stackrel{\mathrm{̃}}{\mathit{y}}}_t$ = γ₀ + γ₁x_t. If a statistically significant relationship was found (again with p-value < 0.05), then the associated correlation coefficient was taken to be a conservative measure of the relationship between the predictand and predictor. We say “conservative” because it is possible that the residuals ε_t = y_t − ${\stackrel{\mathrm{̃}}{\mathit{y}}}_t$ may yet possess an autoregressive structure, i.e., ε_t = φ₁ ε_t−1 + φ₂ ε_t−2 + … φ_i ε_t−i + w_t, where φ₁, φ₁, …, φ_i are constants (φ_i ≠ 0) and w_t represents white noise (Shumway and Stoffer, 2017). If ε_t > 2/${\sqrt{N}}_{\mathrm{yr}}$ at any lag i > 1, where N_yr is the length of the time series, then further analysis would be required to fully specify the regression model. However, including autoregressive terms as predictors in the model will only increase the explained variance (e.g., Pearson's R²), so we can be confident that the correlation coefficient after detrending, but without including autoregressive terms, underestimates the correlation coefficient. In practice, however, autocorrelation was detected only in a few cases (using the function acf2 in the R library astsa from the R Statistical Computing package; R Core Team, 2017), and principally amongst predictors themselves, not in the regressions of APF on predictors.

Ultimately, the nonparametric Spearman rank correlation, with sample estimator $\widehat{\mathit{\rho }}$, was used to characterize the interannual regression results as it makes no assumptions regarding the distribution of the climatic and hydrologic data. A correlation matrix of $\widehat{\mathit{\rho }}$ was constructed by repeated univariate regression over all variables (including detrending where necessary), and the results summarized in the correlograms presented in Sect. 3. Another nonparametric measure, the Theil–Sen slope, is used to calculate temporal trends of the predictand and predictors in Sect. 3.1. The Theil–Sen slope is a median of slopes calculated for each pair of data points and is less sensitive to outliers than the standard least-squares regression line (Sen, 1968). Multiple linear regression (MLR) analyses were also conducted by including all possible combinations of the predictors listed in Table 2. We retained the MLR model that featured 1) the most variables (N) with partial p-values less than 0.05 (N_sig), and (2) a Pearson-adjusted coefficient of determination $R_{\mathrm{adj}}^{\mathrm{2}}$ larger than any MLR with N < N_sig predictors. We also tried adding predictors in a stepwise manner, obtaining similar results.

Table 3 Trends in observed and VIC-simulated variables over the period of each record. Results are only shown for those variables that contain a statistically significant trend at the p < 0.05 level. The trend was calculated using the Theil–Sen median slope estimator (Sect. 2.3). Slightly different results are obtained for observed inputs to the VIC model due to the difference in start and end years of the time series. The residuals of each fit were checked for autocorrelation, with none detected. NA indicates this info is not available in the observations we used.

In order to discern differences between data samples drawn from potentially distinct distributions, we apply a permutation or resampling test at several points in the analysis. For two input data sets, X and Y of length n_X and n_Y, respectively, the test statistic S ≡ mean(X)/mean(Y) is calculated. The two sets are then combined, i.e., Z ≡ X ∪ Y, and a large number (we chose N = 10⁴) of paired, random samples of lengths n_X, n_Y are drawn from Z. For each sample (x_i, y_i), the same statistic S_i is recalculated, and the number of exceedances, N_exc(S_i > S), summed over all i = 1, …, N samples. If the ratio p = N_exc∕N is sufficiently near either 0 or 1, then the original test statistic S is taken to indicate a significant difference between the samples at the p or 1 − p level. We adopt the commonly used threshold value of min(p, 1 − p) = 0.05.

Finally, while most predictors are well characterized by their annual maximum or seasonally averaged values, daily rainfall needs special attention due to its high temporal variance even when averaged over the entire FRB. To better identify significant, multi-day rainfall episodes, we computed the current rainfall index (CRI), following Fedora and Beschta (1989) and Smakhtin and Masse (2000):

$$\begin{array}{ll}{\displaystyle }{\mathrm{CRI}}_t& {\displaystyle }=K<\times {\mathrm{CRI}}_{t-\mathrm{1}}+R_t=R_t+K{R}_{t-\mathrm{1}}+K^{\mathrm{2}}{R}_{t-\mathrm{2}}\\ \text{(2)}& {\displaystyle }& {\displaystyle }+\mathrm{\dots },K<\mathrm{1}.\end{array}$$

The CRI incorporates both the most recent daily rainfall amount (R_t) and an exponential decay in the weight given to previous rainfall events. K is the daily recession coefficient and reflects the storage capacity of the basin which depends somewhat on its physical properties (e.g., topography and area). We set K = 0.9 following a previous application in Pacific Northwest basins (Fedora and Beschta, 1989). The CRI proves useful for investigating intra-annual correlations between R and daily streamflow, which is covered in Sect. 3.3.

3 Results

3.1 Trends in observed time series

We begin by searching for temporal trends in the input time series. These are of intrinsic interest, but also for interpreting the results of the univariate and multivariate regression models later in this section. In addition to APF, we computed linear trends using the Theil–Sen slope estimator for the eight observed variables listed in Table 2 over the common period of record 1956–2006, with the results summarized in Table 3. Of the nine variables, four display statistically significant trends at the p < 0.05 level: APF, freezing degree days (absolute value of the sum of negative daily mean T < 0 ^∘C from 1 October to 31 March, hereafter FDD), April–June mean temperature (T_amj), and the PDO index. The correlation coefficient R² is small for all the fits (R² = 0.10–0.16), indicating that the trends are modest compared to the scatter in the annual data. The decreasing trend in FDD and increasing trend in April–June mean temperature (0.26 ^∘C decade⁻¹) are qualitatively consistent with the regional temperature trends summarized in Sect. 1.1. The trend in APF of −37 m³ s⁻¹ yr⁻¹, or −4.3 % decade⁻¹, is of the same sign but lower magnitude than that estimated for annual mean flow at Hope by BC MOE (2016) (−5.7 % decade⁻¹ between 1958 and 2012). However, it is worth noting that over the entire 1912–2014 period of the gauge record at Fraser-Hope, there is no significant trend detected in APF, even at the less conservative p-level of 0.1. The same is true of the PDO index over its much longer period of record (1900–2015). Nevertheless, we detrended the 1956–2006 observed time series of APF before computing the correlograms shown later in this section.

Figure 2 (a) Regression of observed annual peak daily streamflow (APF) magnitude at Fraser-Hope hydrometric station (vertical axis) against the NINO3.4 index over the period 1912–2014. Year labels are plotted as individual points. (b) Same as in (a) but for the PDO index. Years in boldface font are discussed in the text of Sect. 3.2. (c) Percentiles of APF in years when the NINO3.4 index is in the negative (La Niña) phase (filled circles) versus the positive (El Niño) phase (open circles). The solid curve shows the median result of resampling the combined data set of both ENSO phases 10⁴ times. The significance of differences between APFs in either phase and corresponding samples of the combined data set is assessed using a permutation test (Sect. 2.3), resulting in the indicated p-values. (d) Same as in (c) but for the PDO index.

3.2 Influence of large-scale climate modes on observed streamflow at Fraser-Hope

Before investigating relationships between climate and APF at the basin scale, we examine the influence of the PDO and ENSO on the observed APF at the Fraser-Hope stream gauge station. This station was not included in the Gurrapu et al. (2016) study, but we derive results similar to theirs at other stations, as shown in Fig. 2. The figure shows that higher APF is associated with negative (cold) phases of the NINO3.4 (Fig. 2a) and PDO (Fig. 2b) indices over the 103-year record (1912–2014) at Fraser-Hope. While the statistical relationships between the large-scale climate modes and the APF are robust (Spearman correlation of $\widehat{\mathit{\rho }}$ = −0.40 for NINO3.4, $\widehat{\mathit{\rho }}$ = −0.35 for PDO, both with p < 10⁻³), Fig. 2a and b show that the largest APF in the record (1948; shown in boldface font on the plot) occurred during a neutral PDO phase and weak La Niña conditions, as did the fifth- and sixth-largest (1964 and 1997, also in bold). The second- to fourth-largest APFs (1972, 1950, 2012) occurred during strong negative PDO phases but weak La Niña conditions. This suggests that while APFs occurring during a negative PDO or La Niña phase may be larger than average, the very largest APFs may be influenced by more local drivers, of either climatic or non-climatic nature (e.g., elevation and aspect). Indeed, our results indicate that the PDO and NINO3.4 indices explain only ∼ ${\widehat{\mathit{\rho }}}^{\mathrm{2}}$ ∼ 0.12–0.16 or 12–16 % of the variance (with considerable interdependence between the two modes). Correlations of similar magnitude and sign were found at many more stations within the FRB by Thorne and Woo (2011).

Figure 3 (a) Regression of observed APF against basin-averaged 1 April SWE over the Fraser Basin, within the period of overlapping records, 1956–2014. The solid blue line is the least-squares linear regression fit, while the dashed curves show the 95 % confidence interval. Adjusted Pearson and Spearman rank correlation coefficients, along with their associated p-values, are indicated at upper left. (b) Correlogram of all observed variables, averaged over the Fraser Basin over the common period of records for all variables, 1956–2006. Numerical values are Spearman rank correlation coefficients, and coloured squares indicate significant correlations at the p < 0.1 level.

Another way of visualizing streamflow sensitivity to ENSO and PDO phase is exhibited in the percentile plots of Fig. 2c and d. These plots can reveal differences in sample distributions more clearly than a histogram or density plot. For example, for the PDO the time series of APF was first divided into years with positive (PDO_pos) and negative (PDO_neg) PDO index. Then a resampling (permutation) test was applied to check whether the test statistic S = [mean(APF)]${}_{\mathrm{PDO},\mathrm{neg}}/$[mean(APF)]_PDO,pos differed from that calculated from 10⁴ random samples of the combined (both positive and negative PDO phase) data set (Sect. 2.3). Figure 2d shows that the APF in years with negative (cold) PDO phase is significantly larger than in years with positive (warm) PDO phase, consistent with the relationship seen in the scatter plot (Fig. 2b). Similar results were found for ENSO, as seen in Fig. 2c.

3.3 Basin-scale relationships amongst observed variables

To initiate our investigation of basin-scale drivers for APFs, we revisited the observed data. Specifically, we regressed the observed APF at Fraser-Hope against FRB-wide averages of the observed predictors in Table 2. Figure 3b shows a correlogram of the univariate regression results for all observed variables over the common period of 1956–2006. Variables with detected trends over this period (Sect. 3.1) were detrended before constructing the correlogram, as described in Sect. 2.3. We summarize the results for individual predictors below, before looking at their joint influence on APF at the end of the subsection.

3.3.1 Snow water equivalent (SWE)

As expected given the mainly nival character of the FRB, we find that interannual variations in 1 April SWE and APF are strongly correlated, with a Spearman correlation of $\widehat{\mathit{\rho }}$ = 0.64 over the common period of 1956–2014 of the MSS and WSC time series (Fig. 3a). This figure shows the raw data, while Fig. 3b shows the correlogram of all variables over the shorter common period of 1956–2006, including detrending where necessary. These differences lead to a somewhat smaller $\widehat{\mathit{\rho }}$ = 0.51 between 1 April SWE and APF. Nevertheless, the conclusion that years with higher than average SWE tend to produce higher peak streamflow at Fraser-Hope remains robust.

3.3.2 Temperature

Several diagnostics of the effect of surface air temperature T on APF were examined. First, we computed FDD from the PCIC-OBS data averaged over the entire FRB. This helps determine whether exceptionally cold winters correspond to an unusually large snowpack, potentially producing an enhanced snowmelt contribution to the spring freshet. Second, we calculated the April–June mean temperature T_amj, which roughly coincides with the interval between the time of 1 April SWE and APF. Finally, we computed the warming rate over the freshet period, dT∕dt, by computing the linear least-squares slope of T between 1 April and the date of APF, or annual peak date (APD). Over the 1956–2006 period at Fraser-Hope, the APD ranges from 16 May to 23 July, with a median of 13 June. The corresponding warming rates in different years range from 0.093 to 0.26 ^∘C day⁻¹ (median = 0.16 ^∘C day⁻¹), implying a typical warming over the median freshet (under the linear approximation) of ∼ 11 ^∘C.

The correlations of these variables with APF are displayed in the correlogram of Fig. 3b. Here it is seen that amongst these variables, April–June mean temperature is anti-correlated with APF ($\widehat{\mathit{\rho }}$ = −0.42), while FDD and dT∕dt are positively correlated with APF [$\widehat{\mathit{\rho }}$ = 0.21 (not significant) and $\widehat{\mathit{\rho }}$ = 0.26, respectively]. All else being equal, an unusually cold winter (high FDD) would be expected to result in a larger snowpack and a larger snowmelt contribution to streamflow consistent with the positive correlation found between FDD and APF. By the same reasoning, an unusually cold spring (low T_amj) resulting in a delayed, but more rapid snowmelt during the freshet, would again be expected to increase that year's APF. In years with normal summer T, this situation would produce a larger than normal warming rate, consistent with the positive correlation found between dT∕dt and APF. Finally, we note from Fig. 3b that a positive correlation is seen between T_amj and dT∕dt ($\widehat{\mathit{\rho }}$ = 0.30), implying that high spring temperatures are associated with rapid warming.

3.3.3 Rainfall

As mentioned in Sect. 2.3, daily rainfall R shows a high temporal variance even when averaged over the entire FRB, which could introduce spurious noise into the regressions. For this reason, we looked at three integrated forms of R: (1) the summed rainfall over the winter months, R_Oct−Mar; (2) the sum of rainfall between 1 April and the day of APF, R_Spring; and (3) the sum of R over the 15 days prior to the APD, R_APF. These metrics offer differing probes of the influence of rainfall on the APF at different lags and temporal resolution. However, none of the rainfall measures displayed a significant relationship with observed APF (Fig. 3b), reinforcing the designation of the FRB as a primarily nival basin. Indeed, R_APF was completely unrelated to APF ($\widehat{\mathit{\rho }}$ = 0.02, not shown in Fig. 3b). Although severe flooding in several small catchments in the US Pacific Northwest has been attributed to rain-on-snow events (McCabe et al., 2007; Surfleet and Tollos, 2013), we see no evidence of this in the much larger FRB.

Unlike SWE, the basin-averaged, annual maximum CRI, CRI_max, generally occurs during summer and fall (calendar days 160 to 301), although multiple peaks in a given year are often observed (see ff., Fig. 9). This contrasts with the narrower distribution of APF dates occurring in spring–summer (days 136–204), suggesting a weak relationship may exist between the two on the seasonal timescale.

Figure 4 Results of lagged regression of FRB-averaged daily rainfall, characterized by the current rainfall index (CRI), on streamflow at Fraser-Hope. Only the subset of years wherein there exists a significant relationship between CRI and discharge are shown. The ENSO state, as specified in the legend according to colour (EN: El Niño, LN: La Niña), corresponds to the year preceding the rainfall and discharge, and was obtained from http://ggweather.com/enso/oni.htm, based on the same definition used in the analysis of Sect. 3 (i.e., NINO3.4). The maximum correlation, max(${\widehat{\mathit{\rho }}}_{XY}$), is shown on the vertical axis, with the corresponding lag, τ, between rainfall and discharge on the horizontal axis. A positive lag means that streamflow lags CRI. Circles: observed precipitation versus observed discharge over the 1950–2006 period. Triangles: observed precipitation versus VIC-simulated discharge. All correlations are significant at the p < 0.05 level. Two outlier points with lags > 60 days (with none exceeding max(${\widehat{\mathit{\rho }}}_{XY}$) = 0.42) are omitted from the plot.

Focusing exclusively on annual maxima or summed values of rainfall and APF might cause us to miss important relationships between rain and discharge on shorter timescales. To investigate these types of linkages, we interrogated the respective daily time series. Specifically, we computed the annual cross-correlation function between deseasonalized streamflow (X) and CRI (Y) anomalies for each year i:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{\widehat{\mathit{\rho }}}_{XY,i}\left(\mathit{\tau }\right)={\displaystyle \frac{\mathrm{mean}\left[{\mathit{X}}_i^{\prime }\left(t\right){\mathit{Y}}_i^{\prime }(t+\mathit{\tau })\right]}{{\widehat{\mathit{\sigma }}}_{X,i}{\widehat{\mathit{\sigma }}}_{Y,i}}},\end{array}$$

where ${\mathit{X}}_i^{\prime }\left(t\right)$ is the anomaly time series of X (i.e., difference between X and its multi-year mean annual cycle) in year i, ${\widehat{\mathit{\sigma }}}_{X,i}$ is the standard deviation of ${\mathit{X}}_i^{\prime }$ (and similarly for Y), and τ is the applied lag (days) between ${\mathit{X}}_i^{\prime }$ and ${\mathit{Y}}_i^{\prime }$, ranging over ±183 days. If max(${\widehat{\mathit{\rho }}}_{XY,i}$) > δ at τ > 0, where δ = 2$/\sqrt{N}$ = 2$/\sqrt{\mathrm{365}}$, then CRI both leads and positively correlates with streamflow in that year, suggesting a causal relationship. The results are displayed in Fig. 4. Each point in the scatterplot represents max(${\widehat{\mathit{\rho }}}_{XY,i}$) calculated over (365 − τ) pairs taken from the two daily time series for the ith calendar year. Figure 4 shows that about half (26 of 50) of the years exhibit a significant relationship at some time during the year, with max(${\widehat{\mathit{\rho }}}_{XY,i}$) ranging from 0.24 to 0.72 with lags of τ = 0–51 days (Fig. 4, open circles). Two clear outliers, at τ = 93 and 114 days, were disregarded (not shown). For each point plotted in Fig. 4, the corresponding ENSO state for that year, taken from historical data, is indicated; 3 of the 5 (and 4 of the 10) largest correlations were associated with strong El Niño years, each with a lag of 2 days or less between rainfall and streamflow (Fig. 4). Rainfall events in the upper reaches of the FRB can have a delayed effect on streamflow at Hope of up to ∼ 1 week (BC River Forecast Centre, 2012), while in years with τ = 1 to 8 weeks, soil moisture that is near field capacity might be responsible for successive overland flow following rainfall events. However, Fig. 4 shows that these situations tend to occur during weak or moderate La Niña conditions.

Finally, we note the possibility that significant multi-day rainfall events coinciding with frozen or saturated soils might lead to occasional flooding at local scales within the FRB. While suitable soil moisture observations are not readily available at this time (see Sect. 5), we employ the VIC model to examine the role of soil moisture in both the FRB and its sub-basins using the regression framework in Sect. 3.4.

3.3.4 Relationships amongst observed predictors

Finally, we investigated relationships amongst observational variables unrelated to streamflow revealed in the correlogram of Fig. 3b. April–June mean temperature is positively correlated with spring warming rate and negatively correlated with 1 April SWE, both intuitively reasonable results, as is the positive correlation of FDD with SWE. Strong inverse relationships also exist between spring rainfall and (i) T_amj ($\widehat{\mathit{\rho }}$ = −0.48), and (ii) warming rate ($\widehat{\mathit{\rho }}$ = −0.69). Both of these relations are reasonable given that the rainfall fraction diminishes approximately linearly with decreasing temperature toward 0 ^∘C, according to the Dai (2008) parameterization used here. Finally, cold season (October–March) rainfall is anti-correlated with FDD ($\widehat{\mathit{\rho }}$ = −0.32) (Fig. 3b).

3.3.5 Multivariate regression

The above results were derived using univariate linear regression and correlations computed using the Spearman rank method. However, we also would like to attribute interannual variance in the APF to the combined variances of the predictors. As described in Sect. 2.3, we constructed a number of MLR relations including all possible combinations of the above variables that showed a significant Spearman $\widehat{\mathit{\rho }}$ when regressed individually against APF.

The MLR procedure yields the following multilinear regression fits:

$$\begin{array}{ll}{\displaystyle }{\widehat{\mathrm{APF}}}_{\mathrm{obs}}\left({\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}\right)& {\displaystyle }=\mathrm{4535}+\mathrm{11.40}{\mathrm{SWE}}_{\mathrm{Apr}\mathrm{1}}\\ {\displaystyle }& {\displaystyle }+\mathrm{17}\mathrm{068}(\mathrm{d}T/\mathrm{d}t)-\mathrm{603.2}{T}_{\mathrm{amj}}\\ \text{(4)}& {\displaystyle }& {\displaystyle }\left(R_{\mathrm{adj}}^{\mathrm{2}}=\mathrm{0.63},p<{\mathrm{10}}^{-\mathrm{9}}\right),\end{array}$$

$$\begin{array}{ll}{\displaystyle }{\widehat{\mathrm{APF}}}_{\mathrm{obs}}\left({\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}\right)& {\displaystyle }=\mathrm{4106}+\mathrm{10.48}{\mathrm{SWE}}_{\mathrm{Apr}\mathrm{1}}\\ {\displaystyle }& {\displaystyle }+\mathrm{17}\mathrm{815}(\mathrm{d}T/\mathrm{d}t)-\mathrm{494.3}{T}_{\mathrm{amj}}\\ {\displaystyle }& {\displaystyle }-\mathrm{402.3}\mathrm{NINO}\mathrm{3.4}\\ \text{(5)}& {\displaystyle }& {\displaystyle }\left(R_{\mathrm{adj}}^{\mathrm{2}}=\mathrm{0.65},p<{\mathrm{10}}^{-\mathrm{9}}\right),\end{array}$$

where the units of the predictors on the right-hand sides of Eqs. (4) and (5) are provided in Table 2. Equation (4) includes only local, basin-averaged predictors, while Eq. (5) includes the non-local influences of ENSO and PDO. Interestingly, only NINO3.4, not the PDO index, is a significant predictor of APF in the MLR, despite the comparable importance of both indices in the univariate regressions (Sect. 3.2). The residuals from both fits are indistinguishable from a Gaussian distribution at the 5 % significance level, according to the Kolmogorov–Smirnov statistic. Additional parameters of the fits are summarized in Table 4. The predictors on the right-hand side of Eqs. (4) and (5) are ordered from left to right by decreasing partial F value: e.g., SWE_Apr1 contributes the majority of the interannual variance in APF, followed by dT∕dt, T_amj and NINO3.4 (Table 4). Together, these variables explain 65 % of the variance, 2 % more than if the influence of ENSO is ignored. As mentioned above, rainfall variables do not contribute significantly to the interannual variability of the observed APF at Fraser-Hope.

3.4 Observed versus VIC-simulated streamflow and SWE

Before exploring relationships between APF and the wider array of variables available in the VIC hydrological model, as an evaluation exercise we compare several features of the VIC simulations with available observations in the FRB. In Fig. 5a, we compare the 10 largest APFs in observations and VIC at Fraser-Hope station over the simulation period of 1955–2004. The APFs have been ranked by their observed magnitude, with the highest flow years on the left-hand side of the bar graph. In most years, the VIC-simulated APFs are close to the observed values. Figure 5b compares the daily climatology and interannual variability of VIC and observations over the same period. VIC tends to overestimate the magnitude of the APF by ∼ 8 % (multi-year mean) to 16 % (multi-year maximum), and also simulates an APD ∼ 5 days later than in the observations. VIC underestimates interannual streamflow variability over most of the year, except over the period of peak flow from June to mid-August, when it displays higher variability. However, given that the interannual coefficient of variation in observed APF (CV = σ/mean[APF]) is 18 %, that the VIC-simulated CV = 20 %, and their degree of overlap, we conclude that the hydrographs are not significantly different.

Figure 5 (a) Annual peak daily streamflow (APF) magnitude in observations (red) and in the VIC simulation (blue) at Fraser-Hope, ranked by observed values over 1955–2004. (b) Annual cycle of daily streamflow in observations (red) and VIC simulation (blue) over 1955–2004 at Fraser-Hope. (c) Quantile plots of VIC-simulated versus observed APF and APD at Fraser-Hope. Differences between the respective distributions are assessed using a permutation test with 10⁴ resamples; the resulting p-values are indicated in each sub-panel. (d) 1 April observed SWE (red) compared with annual cycle of SWE from VIC, both averaged over the FRB for 1956–2006. Heavy curves in (b) and (d) show the multi-year mean, while shaded areas show the interannual range.

Figure 5c shows quantile plots of APF and APD for VIC compared to observations over the 1955–2004 period. Here, the permutation test was applied using the test statistic S = [mean(APF)]_VIC/[mean(APF)]_OBS (Sect. 2.3). Figure 5c shows that for APF, the VIC-simulated APF is indistinguishable from the stream gauge observations at the 5 % significance level (i.e., p = 0.07 > 0.05; p = 0.24 if the two outliers with VIC-simulated APF > 13 000 m³ s⁻¹ are removed), while the late bias in VIC-simulated APD is significant at the ∼ 1 % level. As our main concern in this work is with identifying key predictors of the APF, with less emphasis on APD, we conclude that VIC simulates APF reasonably well compared to observations.

Table 4 Fitted parameters from multilinear regression of the form: Y = a₀ + a₁X₁ + a₂X₂ + a₃X₃ + … For each basin or subbasin, the first row in columns 2 to 6 indicates the relevant predictors, while the second and third rows list the fit coefficients and partial F-values, respectively (all partial p-values are less than 0.05). The total F-value and adjusted Pearson correlation for each fit is provided in columns 7 and 8, respectively. For the subbasins, only VIC results are given, due to the small number of 1 April observed SWE values available. Results given for OBS include basin-averaged predictors over the FRB, while those for OBS_ext include indices of large-scale climate variability.

Figure 6 Correlogram of all observed and VIC-simulated variables, averaged over the Fraser Basin within the common period of records, 1955–2004. All values are Spearman rank correlation coefficients, and coloured squares indicate significant correlations at the p < 0.1 level. A hierarchical clustering algorithm has been applied to group like-with-like correlations.

Due to the dominant influence of annual snowpack on APF, it is also desirable to compare the VIC-simulated SWE with observed 1 April SWE measurements. In the analysis of the VIC simulation, we use SWE_max as a predictor instead of SWE_Apr1 since it is expected to be more closely linked to APF on physical grounds. In Fig. 5d, the annual cycle of the corrected, VIC-simulated SWE (Sect. 2.2) is compared with the observed 1 April SWE, with both quantities represented by their basin averages. While there is significant overlap in the interannual ranges of simulated and observed SWE, VIC appears to systematically underestimate the snowpack amount. However, it is important to recognize the sampling bias inherent in the observed, basin-average SWE: since the MSS stations tend to be located at high elevation (range: 750–2200 m), the observed estimates are likely biased high compared to a true basin average, as represented in VIC. In addition, we note that the apparent VIC biases in peak streamflow and SWE over the FRB are qualitatively consistent with VIC simulations over other nival basins in the Pacific Northwest (Salathé Jr. et al., 2014).

3.5 Relationships amongst VIC-simulated variables

In this section we apply the methodology of Sect. 3.3 to probe relationships amongst the wider array of variables available in the VIC hydrological model over the FRB. We begin with an analysis of the routed model streamflow at Fraser-Hope, as a representative outlet of the entire FRB, before examining VIC-modelled discharge at the outlets of four subbasins within the FRB in Sect. 3.6.

In addition to the PCIC-OBS variables used in the regression analysis of Sect. 3.3 – which were used as daily forcings for the VIC simulation – we include a number of variables available in the model but not in observations, namely annual maximum of SWE averaged over the FRB, SWE_max; calendar date of SWE_max; snowmelt rate, d(SWE)∕dt; melt season length, SWE_len; and antecedent total column soil moisture (detailed definitions are given in Table 2). The regression results are again summarized as a correlogram (Fig. 6) and in Table 4. The annual snowmelt rate, d(SWE)∕dt, was calculated as the best-fit, least-squares linear slope of SWE between the dates of SWE_max and APF. Formally speaking, d(SWE)∕dt is the snow ablation rate, which in VIC includes snowmelt as the dominant contribution along with evaporation and rain-on-snow (other snow removal processes occurring in nature such as blowing snow or avalanches are not simulated). For convenience, however, we refer to d(SWE)∕dt as simply the snowmelt rate with the understanding that these other processes may make minor contributions to snowpack disintegration.

3.5.1 Influence of large-scale climate modes on streamflow at Fraser-Hope

Through the influence of the forcing variables, the PDO and ENSO may affect the VIC-simulated APF at Fraser-Hope. Indeed, we find a slightly more robust, but qualitatively similar, relationship between the large-scale climate indices and modelled APF than for observed APF: as seen in Fig. 6, the correlations are Spearman $\widehat{\mathit{\rho }}$ = −0.41 for NINO3.4 and $\widehat{\mathit{\rho }}$ = −0.38 for PDO (both at p < 10⁻²).

3.5.2 Snow water equivalent and snowmelt

As noted for the observed 1 April SWE, we find that SWE_max exercises a strong control on APF ($\widehat{\mathit{\rho }}$ = 0.70, p < 0.001; Fig. 6). The latter lags SWE_max by an average of 84 days (range 52–127 days) over the simulation period; examples for specific years are displayed in Fig. 9. Further, Table 5 shows that of the top 10 APF years in this period, 7 were in the top 10 of SWE_max. Yet, the highest APF corresponds to the 6th-largest SWE_max (1982), while the 4th-highest APF occurred in an average year for snow accumulation (1958). It is therefore of interest to investigate what conditions, independent of snow accumulation, led to the comparatively large APF in those years.

Snowmelt displays significant positive correlations with APF ($\widehat{\mathit{\rho }}$ = 0.43) and several other predictors, namely SWE_max ($\widehat{\mathit{\rho }}$ = 0.63), FDD ($\widehat{\mathit{\rho }}$ = 0.38), date of SWE_max ($\widehat{\mathit{\rho }}$ = 0.32), dT∕dt ($\widehat{\mathit{\rho }}$ = 0.25), and maximum annual soil moisture ($\widehat{\mathit{\rho }}$ = 0.64). In years with high SWE_max, snowmelt in spring recharges soil moisture (over unfrozen ground; see Sect. 3.5.3), sometimes to field capacity, which leads to increased runoff and the strong positive correlations between snowmelt, soil moisture and APF. Significant negative correlations are also found with melt season length ($\widehat{\mathit{\rho }}$ = −0.53) and the NINO3.4 and PDO indices ($\widehat{\mathit{\rho }}$ = −0.28 and −0.46, respectively). The inverse relationship with the large-scale indices should be considered in light of their similar relationship with SWE_max ($\widehat{\mathit{\rho }}$ = −0.47 and −0.62, for NINO3.4 and PDO respectively): higher snowpack accumulates during La Niña and negative PDO phases, with reduced snowmelt due to cooler temperatures.

3.5.3 Soil moisture

In our VIC simulation, the basin-averaged annual peak soil moisture, SM_max, occurs during the snowmelt period from mid-March to mid-June, most often near the end of May, well after the day of SWE_max and approximately 3 weeks in advance of the APF (for an illustration depicting individual years, see Fig. 9). The direct influence of both SWE_max ($\widehat{\mathit{\rho }}$ = 0.83) and snowmelt ($\widehat{\mathit{\rho }}$ = 0.64) on SM_max is detected at high confidence (p < 10⁻³), while SM_max is, in turn, strongly correlated with APF ($\widehat{\mathit{\rho }}$ = 0.65, p < 10⁻³) in the annual time series (Fig. 6). A cross-correlation analysis of the respective daily time series over the entire common period reveals that, on average, SWE leads SM by 68 days (maximum cross-correlation, R = 0.65), SM leads streamflow by 22 days (R = 0.85), and SWE leads streamflow by 105 days (R = 0.74). The higher cross-correlation between SM and streamflow indicates that the annual cycle of SM is a better predictor of the daily streamflow hydrograph than SWE, despite the above-mentioned superiority of SWE_max as an annual predictor of APF. This is a reasonable result, since during the freshet, daily SM integrates contributions from both snowmelt and precipitation. On the other hand, SM_max could be considered inferior to SWE_max as a predictor of APF with respect to its much shorter lead time. SM_max also exhibits significant relationships with FDD ($\widehat{\mathit{\rho }}$ = 0.33), T_amj ($\widehat{\mathit{\rho }}$ = −0.27) and the NINO3.4 and PDO indices ($\widehat{\mathit{\rho }}$ = −0.44 for both). The negative phases of ENSO and the PDO bring more rain and snow, which consequently enhances soil moisture.

Table 5 Top 10 APF years at Fraser-Hope in the VIC simulation over the period 1955–2004, along with corresponding ranks of other basin-averaged predictors. Years in boldface font are highlighted in Fig. 9.

Also of interest is a possible relationship between APF and SM preceding the snowmelt period, for example, in fall before snow accumulation begins. Numerous studies point to the influence of antecedent soil moisture on seasonal streamflow in nival catchments (Maurer and Lettenmaier, 2003; Berg and Mulroy, 2006; Williams et al., 2009; Harpold and Molotch, 2015). Using a suite of land surface models including VIC, Koster et al. (2010) and Mahanama et al. (2012) demonstrated improvement of March–July streamflow forecast skill over the western United States using antecedent (1 January) soil moisture in addition to snow amount as a predictor. The basic mechanism is that antecedent reduced storage capacity of wet or frozen soils leads to more of the snowmelt being routed to discharge during the spring freshet, and vice-versa for dry soils. To investigate the interannual sensitivity of APF to antecedent SM, we used monthly mean SM from the preceding August through November in turn as predictors of APF, but found no significant correlations. This insensitivity may be due to the lack of a frozen soil parameterization in the version of VIC used here. In this implementation, soil drainage over the cold season tends to be overestimated (Cherkauer and Lettenmaier, 2003), which could cause a decoupling between autumn SM and spring freshet flows.

3.5.4 Temperature and rainfall

As in the observations, FDD and spring warming rate dT∕dt are positively correlated with APF ($\widehat{\mathit{\rho }}$ = 0.24 and $\widehat{\mathit{\rho }}$ = 0.38, respectively; Fig. 6). The connection between dT∕dt and APF is stronger than in the observations but, interestingly, there is no significant correlation between T_amj and APF as exists in the observational data (Eqs. 4 and 5). Further, the results in Table 5 suggest that large APFs in average SWE years are characterized by large dT∕dt (e.g., 1982, 1958, 2002).

Figure 7 Correlogram of observed and VIC-simulated variables, averaged over each of the indicated subbasins within the 1955–2004 period, against APF at the subbasin outlet (Table 1). The cell values and colour scale indicate Spearman rank correlations, significant at the p < 0.1 level. The ordering of variable columns is the same as in Fig. 6, for ease of comparison.

Again as in the observations (Sect. 3.3.3), no significant correlations are found between any of the rainfall measures and VIC-simulated APF. However, as also found in the observations, there are indications of a rainfall–streamflow connection at other times of the year. Repeating the cross-correlation analysis between streamflow and CRI as in Sect. 3.3.3 reveals that 20 of 50 years exhibit a significant relationship at some time during the year, with max(${\widehat{\mathit{\rho }}}_{XY,i}$) ranging from 0.27 to 0.66 with lags of τ = 2–55 days. Two clear outliers, at τ = 91 and 142 days, were disregarded. Figure 4 (triangles) shows that most (nine) of the rainfall-influenced years occur during the El Niño phase, while four occur during La Niña and seven during neutral years. Twelve of the 20 rainfall-influenced years occur during the cool phase of the PDO. More often than not, however, CRI lags daily streamflow, suggesting little to no relationship. With regard to other relationships, we find that cold season rainfall is correlated with spring rainfall ($\widehat{\mathit{\rho }}$ = 0.32) and anti-correlated with FDD ($\widehat{\mathit{\rho }}$ = −0.32) (Fig. 6). As mentioned in Sect. 3.3.4, spring rainfall is negatively correlated with T_amj, dT∕dt, and snowmelt, but with slightly different $\widehat{\mathit{\rho }}$ values due to the different lengths of averaging period in the calculations compared to the observed case (as per the definitions in Table 2).

3.5.5 Multivariate regression

As in the case of the observed variables, we constructed MLR relations including all combinations of the variables in Fig. 6 that showed a significant Spearman $\widehat{\mathit{\rho }}$ and selected the optimal MLR based on the criteria specified in Sect. 3.3.5. One predictor that was excluded from the MLR was SM_max: its high correlation with SWE_max indicates that its independent explanatory power is limited (Sect. 3.5.3). The resulting relationship is

$$\begin{array}{ll}{\displaystyle }{\widehat{\mathrm{APF}}}_{\mathrm{VIC}}\left({\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}\right)& {\displaystyle }=\mathrm{3239}+\mathrm{37.14}{\mathrm{SWE}}_{\max }+\mathrm{26}\mathrm{842}(\mathrm{d}T/\mathrm{d}t)\\ {\displaystyle }& {\displaystyle }-\mathrm{57.69}{\mathrm{SWE}}_{\mathrm{len}}+\mathrm{1770}\mathrm{d}\left(\mathrm{SWE}\right)/\mathrm{d}t\\ \text{(6)}& {\displaystyle }& {\displaystyle }\left(R_{\mathrm{adj}}^{\mathrm{2}}=\mathrm{0.75},p<{\mathrm{10}}^{-\mathrm{9}}\right).\end{array}$$

In contrast to the MLR constructed for observed streamflow, Eq. (5), neither of the large-scale climate indices has a significant influence, nor is T_amj an important predictor. Nevertheless, the four variables on the right-hand side of Eq. (6) account for 75 % of the interannual variance in APF. Three of the four were identified in the univariate regressions as important, with the exception being the length of the melt season, SWE_len. The latter displays a strong anti-correlation with snowmelt ($\widehat{\mathit{\rho }}$ = −0.55, Fig. 6), suggesting that the two variables are not independent. However, since removing each of SWE_len and d(SWE)∕dt from the MLR in turn yields a significantly poorer fit ($R_{\mathrm{adj}}^{\mathrm{2}}$ = 0.71 for d(SWE)∕dt only and $R_{\mathrm{adj}}^{\mathrm{2}}$ = 0.69 for SWE_len only), it seems that both predictors have some explanatory value. The absence of both the PDO and NINO3.4 indices from the MLR could be a reflection of the limitations of the driving data, inasmuch as the processing of gridded station data to the regular grid may weaken the influence of the large-scale climate drivers (Sect. 2.1).

3.6 Relationships amongst VIC-simulated variables at the subbasin scale

The same analysis as conducted for the entire FRB was repeated for each of the four subbasins comprising nearly 70 % of the annual flow at Fraser-Hope (Table 1; Kang et al., 2016). The results of the univariate regressions are presented as a correlogram in Fig. 7, while those for the MLR are provided in Table 4.

Overall, the relationships in the subbasins mirror those seen in the FRB as a whole. The univariate analysis shows that in all four subbasins, SWE_max and SM_max are good predictors of APF, while snowmelt and FDD exhibit positive correlations with APF in three of the four subbasins, Chilko being the exception (Fig. 7). While SWE_max has a slightly stronger influence on APF than SM_max in the FRB, SM_max is more influential at the subbasin scale, explaining over 75 % of the variance in APF in Quesnel (compared to 40 % for SWE_max). The correlation between spring dT∕dt and APF is strong in Thompson–Nicola ($\widehat{\mathit{\rho }}$ = 0.50, p < 0.05) but weaker in Quesnel ($\widehat{\mathit{\rho }}$ = 0.25, p < 0.1) and not significant in the other two catchments. Three of the four subbasins exhibit a strong inverse relationship between APF and the NINO3.4 and PDO indices, again with the exception of Chilko. The weak dependence of streamflow on the PDO, SOI and Pacific-North American indices in this and other catchments in the western FRB was also noted by Thorne and Woo (2011), who attributed this insensitivity to the low magnitude of discharge in these tributaries. These authors speculated that the main trunk of the Fraser, by contrast, integrates the influence of regional climate forcings on upstream catchments, bolstering the teleconnections. The insensitivity of APF to the local, basin-averaged predictors might be related to the smaller size of the Chilko compared to the other three subbasins, which lowers the signal-to-noise ratio of the basin averages. This, along with the comparatively lower APFs at the basin outlet, makes the detection of significant correlations more challenging in Chilko compared with the other subbasins.

Moving now to the MLR analysis, the results in Table 4 demonstrate that in all four subbasins, SWE_max is the most skillful predictor of APF, again in agreement with the results for the FRB as a whole. Likewise, in three of the four catchments, dT∕dt emerges as the next most significant variable. The exception is the Upper Fraser basin, where SWE_len is the second most skillful predictor. In this subbasin, APF is inversely correlated with SWE_len in the MLR, while the latter is in turn anti-correlated with both dT∕dt ($\widehat{\mathit{\rho }}$ = −0.81, p < 0.05) and T_amj ($\widehat{\mathit{\rho }}$ = −0.25, p < 0.1) in the univariate regressions (Fig. 6). These results suggest that anomalously warm springtime temperatures (compared to the preceding winter) are at the heart of the APF–SWE_len relationship in the Upper Fraser, not unlike what is seen in the other subbasins.

In contrast to the FRB as a whole, rainfall does appear to have a weak influence on APF in two of the subbasins. However, in each subbasin, APF is sensitive to a different measure of rainfall. In Thompson–Nicola, the basin with the largest area and second-highest mean elevation, R_Oct−Mar is positively correlated with APF, while in the Upper Fraser, R_Spring displays a somewhat weaker positive relationship with APF. In the Quesnel, a smaller basin of lower mean elevation, the preceding September mean soil moisture is an effective predictor of APF. Interestingly, neither the NINO3.4 nor the PDO index is an effective predictor of APF in the MLR of any of the subbasins, despite the fact that strong inverse relationships are still seen in the univariate correlogram (Fig. 6).

3.7 Co-dependence of streamflow predictors

The above results make it clear that, of the predictors considered in the FRB and its subbasins, APF is primarily influenced by SWE_max and the rate of warming in spring, dT∕dt. To explore the co-dependence between predictors, we show in Fig. 8 a scatterplot of the relative anomaly of each of these two predictors, $\mathrm{\Delta }X/\stackrel{\mathrm{‾}}{X}$ = (X_i − $\stackrel{\mathrm{‾}}{X})/\stackrel{\mathrm{‾}}{X}$, where X_i is the value of predictor X at year i and $\stackrel{\mathrm{‾}}{X}$ is the long-term mean (1955–2004), with the corresponding streamflow relative anomaly $\mathrm{\Delta }Q/\stackrel{\mathrm{‾}}{Q}$ indicated by the point colour. This type of plot has been used previously to explore the elasticity (i.e., non-linearity) of streamflow as a function of covariates (e.g., Andréassian et al., 2016). However, here we employ it primarily as an additional illustration of the co-dependencies identified by the MLR analysis.

In the FRB, the roughly uniform pattern of scatter in the vertical and horizontal directions indicates a weak co-dependence between ΔSWE_max and Δ(dT∕dt), consistent with expectation – dT∕dt is computed only after the day of SWE_max each year – and with the results of Figs. 5 and 6. The corresponding $\mathrm{\Delta }Q/\stackrel{\mathrm{‾}}{Q}$ values display an overall gradient from bottom left to top right, i.e., from lower than average SWE_max and dT∕dt to higher than average SWE_max and dT∕dt, again consistent with the univariate and multivariate regression results. Qualitatively similar results are found in the Chilko, Thompson–Nicola and Quesnel subbasins.

In the Upper Fraser subbasin, the primary predictor remains SWE_max, but the secondary predictor is SWE_len. As a result, the scatterplot displays a weak co-dependence between relative changes in the two variables, with some clustering toward the diagonal. Furthermore, the overall gradient in $\mathrm{\Delta }Q/\stackrel{\mathrm{‾}}{Q}$ is from top left to bottom right, i.e., from lower than average SWE_max and higher than average SWE_len to higher than average SWE_max and lower than average SWE_len. That is, in years when an anomalously large snowpack (ΔSWE_max > 0) melts fairly quickly (ΔSWE_len < 0), streamflow tends to be larger than usual (ΔQ > 0).

Figure 8 Scatterplots showing the co-dependence of relative APF, $\mathrm{\Delta }Q/\stackrel{\mathrm{‾}}{Q}$ (colour scale) on the two principal predictors from the MLR analysis (x- and y-axes), expressed as relative anomalies, in each basin. Counterclockwise from top left: entire Fraser (FRB), Thompson–Nicola (THN), Upper Fraser (UFR), Chilko (CHK) and Quesnel (QNL). The primary predictor in all subbasins is SWE_max, while the secondary predictor is dT∕dt in all subbasins except UFR, where it is SWE_len. The Spearman rank correlation between the two predictors is shown at upper left in each panel; other than in UFR, none are significant at the 10 % (p < 0.1) level.

4 Case studies: high-flow years in the FRB

While the MLR approach captures the relationship between various predictors and APF in the FRB and its subbasins in a statistical manner, it can miss unusual weather factors that are important in specific years. Both historical data and anecdotal accounts of large historical floods in the region place a large emphasis on the precise seasonal development of warming temperatures and rainfall during the snowmelt period from early April to late June as being a critical factor in the development of damaging floods (BC River Forecast Centre, 2012; Septer, 2007). For this reason, we present in this section an “anatomy” of streamflow and the associated key hydrological variables in years of particularly high VIC-simulated discharge at the Fraser-Hope station. We focus on VIC-simulated variables since full time series of all variables of interest are available to construct both a climatology and annual cycles for specific years.

Table 5 shows the top 10 APFs occurring in the VIC simulation from 1955 to 2004, with corresponding rankings of other key predictors entering the regression analysis of Sect. 3. Seven of the top 10 SWE_max years are in this group. This highlights the dominant effect of spring snowpack on APF magnitude in the FRB. A similar result holds for observed basin-averaged 1 April SWE and APF at Fraser-Hope: 8 of the top 10 1 April SWE years are also amongst the top 10 observed APFs. Looking at results for the other predictors, 7 of the top 10 SM_max years (not surprising given the linkage between SWE_max and SM_max noted in Sect. 3.4.3), 4 of the top 10 spring warming rates and 3 of the top 10 FDDs are in the group of top 10 APFs. Only 1 of the top 10 spring rainfalls is in the group, which happens to be the highest simulated rainfall in 1999 (8th largest APF). These results are therefore consistent with the key predictors identified for APF in the regression analysis of Sect. 3.

Figure 9 shows the annual cycle of discharge at Fraser-Hope along with other variables of interest over the calendar year for 3 of the top 10 APF years: 1958, 1972 and 1999. The daily evolution of discharge Q at Fraser-Hope and basin-averaged SWE, T, CRI and SM are plotted, along with their 1955–2004 climatologies, as an aid to evaluating how anomalous a given year is. A 15-day smoothing filter was applied to T to better highlight sub-seasonal trends. Also shown in the bottom two subpanels of each panel are the difference between the slope of T calculated over a 61-day moving window in the year of interest and its climatology (the warming rate anomaly), and the difference of the snowmelt rate from its climatology calculated over a 15-day moving window. Both anomalies were set to zero if T < 0, since we are interested in behaviour during the melt season only. We consider each year in turn, as each exhibits a unique phenomenology that illuminates the development of the APF in that year. Also, in this section we provide descriptions of the two largest flooding events ever recorded in the FRB – the freshets of 1894 and 1948 – both of which precede the simulated period. It is instructive that there are similarities between the meteorological precursors and covariates in those years and the explicitly simulated years highlighted below.

Figure 9 Annual cycles of key variables from the VIC simulation, constructed from daily output. Results for three high-APF years are shown: (a) 1958; (b) 1972; and (c) 1999. In each panel, sub-panels show, from top to bottom: streamflow (black, left-hand axis); SWE (blue, right-hand axis); surface air temperature (red, left-hand axis); current rainfall index, CRI (green, right-hand axis); soil moisture (brown); warming rate anomaly; and melting rate anomaly. In the upper three sub-panels, solid curves show time series for the year indicated at top right, while dashed curves show the multi-year (1955–2004) mean climatology. The calculations underlying the lower two sub-panels are described in Sect. 4.

4.1 Warm spring in an average SWE year: 1958

APF in the year 1958 ranked 4th over the 1955–2004 period of the simulation (Table 5) and 13th in observed APF. Yet, as simulated by VIC, 1958 was a normal SWE year, with SWE_max achieved near the end of March. Rainfall and soil moisture were also near normal up to this point, while winter temperatures were above normal but still below freezing (Fig. 9a). At the beginning of the melt season in April, dT∕dt began to exceed its climatological value, prompting rapid snowmelt toward the end of that month. The anomalous snowmelt persisted until the end of May, about the time of the APF (APD), while the period of elevated warming rate continued into early June. The snowmelt pulse produced a coincident soil moisture anomaly, which is remarkable in that these quantities are averaged over the entire FRB, yet their behaviour is consistent with small-scale processes. Subsequent to the APD, the soil moisture anomaly became negative, apparently in response to the prolonged warming, which would lead to enhanced evaporation, and lower than normal rainfall over the summer. Figure 9a demonstrates that neither rainfall nor high soil moisture prior to freeze-up were pivotal factors in producing the anomalously large APF in this year. Indeed, rainfall amounts were well below normal between snowmelt initiation and the APD (49th out of 50 years; Table 5).

4.2 A high-SWE year: 1972

The second-highest simulated APF over the 1955–2004 period occurred in 1972, which also featured the second-highest SWE_max (Table 5). This year also exhibited the largest observed APF over the same period (Fig. 5a). While the preceding winter was colder than normal – evidently a factor in generating an anomalously large snowpack – the spring warming rate was above normal, resulting in a strong snowmelt anomaly lasting from May until just after the APD in early June (Fig. 9b). As in 1958, this snowmelt pulse generated a large soil moisture anomaly, which likely produced more overland flow during the freshet. This evidence suggests that the extreme APF in 1972 was primarily caused by a combination of elevated warming in spring and a heavy snowpack, which provided a large water surplus during the freshet. However, there is also an indication that a series of heavy rainfall episodes starting in early June, just days before the APD, affected the APF and particularly the duration of the high-flow period. An extended period of anomalously large rainfall is seen in the CRI between early June and late July, with a clear response in SM and in the persistence of the high-discharge anomaly until the end of August – both are clearly distinguishable from their counterparts in 1958 (Fig. 9a).

Finally, it is worth noting that 1982, the year with the highest simulated APF (not shown in Fig. 9), was characterized by a very similar evolution of predictors as in 1972, including an extended, coincident period of anomalously large rainfall. However, in that case, there were two additional elements: (1) the spring warming rate was the highest of any year in the record (Table 5), and (2) a strong warm anomaly ∼ +5 ^∘C occurred in mid-June just prior to the APD, which evidently generated enough additional snowmelt runoff to position that APF as the highest over the simulated period.

4.3 Influence of rainfall in a high-SWE year: 1999

The year 1999 holds the distinction of having the highest simulated spring rainfall in the record. It also had the fourth-largest VIC-simulated SWE_max, while ranking eighth in APF (fifth in observed APF; see Fig. 5a and Table 5). The temperature development was unremarkable, albeit a little colder than the climatology in the freshet months of April and June (Fig. 9c). The hydrograph reflects this increased role of rainfall by its evident synoptic timescale variability, in contrast to its smoother counterpart in the snowmelt-dominated years examined above (Fig. 9a and b). Indeed, the character of the hydrograph over the freshet period is such that isolating the single largest APF misses key features of the flow development, which involves not one but four peak flows of roughly equal magnitude spanning the period from mid-June to mid-August. Consequently, the use of annual predictors with a single APF as predictand, as adopted in our univariate and multivariate regressions, is likely to underweight a strong influence of rainfall in a particular year. Table 3 shows that the maximum cross-correlation of deseasonalized daily streamflow and CRI anomalies during 1999 is 0.61 at a lag of 35 days (CRI leading streamflow), while in 1958 the maximum is just 0.38 at a lag of 47 days. This implies that rainfall should be treated differently than SWE and other predictors to correctly capture its influence on APF. The evolution of SM is also different than in the prior cases considered, insofar as it displays step-like jumps in mid-March, early April, and mid-May followed by a plateau at a high level until early July (Fig. 9c). This SM anomaly persisted for the remainder of the year, except for a brief 2-week period in November. While in snowmelt-dominated years SWE approaches its climatological value in late summer (except in 1958, when it fell below that value), in 1999 a positive SWE anomaly persisted throughout the year.

4.4 The two highest flow years in the FRB: 1894 and 1948

A description of meteorological precursors to the 1894 and 1948 floods is available from the BC MOE (BC MOE, 2008) and also the historical review, based on newspaper accounts, of Septer (2007). According to limited archival data available at two interior locations in the FRB, snowfall amounts during the winter of 1893 and 1894 were near normal over the FRB as a whole, although the snowpack remaining from the previous summer was reportedly larger than usual. While the spring of 1894 was cold and wet in the FRB, temperatures rose rapidly in late May, becoming unseasonably warm near the end of that month. By the end of May, the Fraser overtopped its banks at many locations and existing dikes were breached, leading to the “Great Chilliwack Flood” in the densely populated Lower Fraser valley. The Fraser River at Mission (near Hope) peaked on 5 June, but floodwaters did not subside in many areas until early July. A hydraulic modelling exercise conducted in 2008 estimated a peak discharge at Hope in the range of 16 000 to 18 000 m³ s⁻¹, with a best guess of 17 000 m³ s⁻¹ (BC MOE, 2008).

The meteorological record for the FRB in 1948 is considerably more complete than for 1894. It shows that precipitation was well above normal for a long period in advance of the 1948 freshet, with excess rainfall in autumn 1947 and heavy snowfall in the subsequent winter and early spring. As a result, soils in many areas were likely near saturation at the time of freeze-up, an important factor for the partitioning of runoff from spring snowmelt. Also, the accumulated reservoir of snow was larger than normal. The early spring of 1948 was cooler than average to 20 May, followed by a rapid increase to unseasonable warmth in late May–June. Thus, sensible heat inputs to the existing high snowpack in the FRB were large, leading to a rapid and voluminous release of meltwater to lower elevations. The peak flow at Hope on 31 May was measured at 15 200 m³ s⁻¹, which remains the highest value in the instrumental record (Septer, 2007; BC MOE, 2008). Apparently, neither the 1894 nor the 1948 event was characterized by excessive rainfall coincident with the meltwater peaks.

5 Discussion and conclusions

In Sect. 3.4, we showed that the VIC model driven by gridded observations provides an adequate simulation of streamflow and maximum annual SWE in the FRB. The mutual resemblance of the correlograms (Figs. 3b and 6), which summarize the univariate linear regression fits to observed and VIC data, along with the similar forms of the respective MLR models (Eqs. 1–3), give one further confidence in the ability of the VIC model to simulate interactions between the key controls of streamflow acting in the real system. Furthermore, differences in these controls in the observed and modelled cases suggest fruitful avenues of further research. For example, the use of more complete snow survey or satellite products might permit the estimation of SWE_len or snowmelt rate, and confirmation of the influence of these terms appearing in the VIC-derived MLR, Eq. (6). Or conversely, improved driving data for the VIC model might reveal the influence of the PDO and/or ENSO explicitly in the MLR, or the inverse correlation with spring temperature and/or freezing degree days. The relationships derived between the basin-averaged predictors examined and streamflow at the major basin and subbasin outlets in the FRB should be relevant to other large, nival watersheds elsewhere on the globe, since the physical drivers (e.g., temperature seasonality and precipitation phase change) are generic and scale-independent. However, it is important to recognize that the influence of large-scale drivers (e.g., ENSO, PDO) will differ depending on the specific geographic setting.

Interestingly, in the historical climate record we find little prospect of using basin-averaged, rainfall-related indices as predictors of APF at Fraser-Hope. It is possible, however, that at subbasin scales, and/or in combination with an appropriate soil moisture indicator, some measure of rainfall might prove useful for predicting localized flooding. This has been a valuable strategy for APF/APD hindcasting in nival–pluvial and pluvial catchments (Neiman et al., 2011; Surfleet and Tullos, 2013). Moreover, given that the rain-to-snow ratio will increase in the FRB in response to further warming (see below), future studies of projected streamflow could uncover an explicit rainfall–streamflow connection at the basin scale.

We emphasize that the identification of the principal controls on APF at the interannual timescale within the MLR framework has its limitations. The relationships derived from historical inputs can only be considered broadly indicative of expected behaviour, and then only in a basin-averaged sense. Although based on only 50 years of simulated basin hydrology driven by observed data, the VIC results imply that a high SWE_max, say in the top quartile of historical values (i.e., 7 of the top 10 in Table 5), increases the chance of an upper quartile APF. The long interval between SWE_max and APF, of order 2–4 months, makes the former a valuable early warning indicator of possible flooding in the lower FRB (BC River Forecast Centre, 2012). However, it is also true that 3 of the top 10 SWE_max years (1976, ranked 3rd; 1991, 8th; and 1971, 9th) did not exhibit remarkable floods: in those years, the corresponding APFs ranked 20th, 30th and 23rd. And while SM_max has about the same success rate as a predictor of APF – 7 of the top 10 APF years are top 10 SM_max years – these are the same as for SWE_max, meaning that soil moisture during the freshet is of little additional predictive value. As mentioned in Sect. 3.5.3, the lack of a connection between fall SM and APF in the historical VIC simulation suggests that future model studies should include the effect of soil freezing to more realistically simulate the spring thaw. Further, the estimation of groundwater storage and its evolution using NASA's Gravity Recovery and Climate Experiment (GRACE), which has proved promising for flood forecasting in large, snowmelt-dominated basins (Reager et al., 2014; Wang and Russell, 2016), could find useful application in the FRB.

In Sect. 3, we demonstrated that the springtime warming rate, dT∕dt, was the next most skillful predictor of APF in the FRB. Moreover, dT∕dt appears to provide additional predictive value over that of SWE_max and SM_max, insofar as 2 of the top 10 dT∕dt years are not top 10 SWE_max and SM_max (1958 and 2002; Table 5). Thus, the combination of an upper quartile SWE_max with a high dT∕dt over the snowmelt period may presage flooding at Fraser-Hope station (9 out of the top 10 APFs were characterized by one or the other; 4 featured both). The definition of dT∕dt used here is dependent upon the APD, which is of course unknown at the time of SWE_max; however, in a predictive context, one could simply use the most recent daily T observation to compute the warming rate within a window of increasing duration. Better still, if hydrological forecasting on the daily to seasonal timescale is of paramount interest, employing a forecast methodology based on an ensemble of hydrological simulations and multiple initializations would constitute a superior approach (Bazile et al., 2017).

In the event that future simulations of the FRB were to furnish annual time series of the basin-averaged predictors appearing in Eq. (6) and/or Table 4 – for example from a coupled global climate model – could the MLR relationships derived in this paper be used for predictions of future APF? We could answer in the affirmative if both of the following conditions were to hold: (1) the dominant predictors of APF in the late twentieth century (defining a multi-decadal reference period) remain unchanged in the future period of interest, and (2) the APF predictions comprise multi-year or decadal-scale anomalies from the reference period: e.g., APF in 2080–2100 compared to 1980–2000. Assumption (1) is placed in some doubt due to the significant warming projected for the FRB in future (Shrestha et al., 2012), which will increase the rain-to-snow ratio over large areas of the basin already characterized by an annual mean temperature near 0 ^∘C. This suggests that rainfall may become a key predictor for high streamflow in the region, with perhaps a less pivotal role played by SWE_max. Regarding condition (2), the suppression of interannual variability in the MLR fits underestimates interannual variations in APF, meaning that only decadal-scale averages of the predicted APF would be meaningful.

Of what use, then, are the predictor sets derived from the MLR analysis? The application of nonstationary extreme value theory to APFs suggests at least one promising avenue of progress (Towler et al., 2010; Shrestha et al., 2017). In such models, the parameters of fitted distributions (e.g., location, scale and shape parameters of the generalized extreme value (GEV) distribution) are taken to be functions of climate covariates. Temperature and precipitation are often chosen for this purpose, due to their availability in observations and climate models and the presumption that they should influence streamflow at least under some circumstances. But as we have shown, in a large nival basin such as the FRB, rainfall has a weak influence on APF (although it would be expected to influence smaller magnitude flows), while the warming rate is more influential than temperature itself in determining APF magnitude. Hence, it may be worth introducing snowfall or SWE_max and/or spring dT∕dt as covariates of the GEV parameters in the fitting of nonstationary probability distributions of APFs. We leave this as an interesting topic for future investigation.