Long-term reliability of the Athabasca River (Alberta, Canada) as the water source for oil sands mining

Significance We show that current and projected surface water allocations from the Athabasca River, Alberta, Canada, for the exploitation of the Alberta oil sands are based upon an untenable assumption of the representativeness of the short instrumental gauge record. Our trend analysis of the instrumental data shows declining regional flows. Our tree-ring reconstruction shows periods of severe and prolonged low flows not captured by the instrumental record. Exploitation of the Alberta oil sands, the world’s third-largest crude oil reserve, requires fresh water from the Athabasca River, an allocation of 4.4% of the mean annual flow. This allocation takes into account seasonal fluctuations but not long-term climatic variability and change. This paper examines the decadal-scale variability in river discharge in the Athabasca River Basin (ARB) with (i) a generalized least-squares (GLS) regression analysis of the trend and variability in gauged flow and (ii) a 900-y tree-ring reconstruction of the water-year flow of the Athabasca River at Athabasca, Alberta. The GLS analysis removes confounding transient trends related to the Pacific Decadal Oscillation (PDO) and Pacific North American mode (PNA). It shows long-term declining flows throughout the ARB. The tree-ring record reveals a larger range of flows and severity of hydrologic deficits than those captured by the instrumental records that are the basis for surface water allocation. It includes periods of sustained low flow of multiple decades in duration, suggesting the influence of the PDO and PNA teleconnections. These results together demonstrate that low-frequency variability must be considered in ARB water allocation, which has not been the case. We show that the current and projected surface water allocations from the Athabasca River for the exploitation of the Alberta oil sands are based on an untenable assumption of the representativeness of the short instrumental record.

Exploitation of the Alberta oil sands, the world's third-largest crude oil reserve, requires fresh water from the Athabasca River, an allocation of 4.4% of the mean annual flow. This allocation takes into account seasonal fluctuations but not long-term climatic variability and change. This paper examines the decadal-scale variability in river discharge in the Athabasca River Basin (ARB) with (i) a generalized least-squares (GLS) regression analysis of the trend and variability in gauged flow and (ii) a 900-y tree-ring reconstruction of the water-year flow of the Athabasca River at Athabasca, Alberta. The GLS analysis removes confounding transient trends related to the Pacific Decadal Oscillation (PDO) and Pacific North American mode (PNA). It shows long-term declining flows throughout the ARB. The tree-ring record reveals a larger range of flows and severity of hydrologic deficits than those captured by the instrumental records that are the basis for surface water allocation. It includes periods of sustained low flow of multiple decades in duration, suggesting the influence of the PDO and PNA teleconnections. These results together demonstrate that low-frequency variability must be considered in ARB water allocation, which has not been the case. We show that the current and projected surface water allocations from the Athabasca River for the exploitation of the Alberta oil sands are based on an untenable assumption of the representativeness of the short instrumental record. paleohydrology | statistical hydrology | oil sands | Alberta | climate variability O ver the past several decades, the province of Alberta has had Canada's fastest growing economy, driven largely by the production of fossil fuels. Climatic change, periodic drought, and expanding human activities impact the province's water resources, creating the potential for an impending water crisis (1). The Athabasca River (Fig. 1) is the only major river in Alberta with completely unregulated flows. It is the source of surface water for the exploitation of the Alberta oil sands, the world's third-largest proven crude oil reserve at roughly 168 billion barrels. The oil and gas industry accounted for 74.5% of total surface water allocations in the Athabasca River Basin (ARB) in 2010 ( Fig. 2) (2). An almost doubling of ARB water allocations since 2000, or 13 times the provincial average, is attributable to expanding oil sands production, which began in 1967 (Fig. 2).
According to the Canadian Association of Petroleum Producers, in 2012, surface mining of the oil sands and in situ extraction (drilling) required 3.1 and 0.4 barrels of fresh water, respectively, to produce a barrel of crude oil (3). This amounted to 187 million cubic meters of fresh water use in 2012, or the equivalent of the residential water use of 1.7 million Canadians (4). Within the next decade, cumulative water use for oil sands production is projected to peak at about 505 million cubic meters per year or a rate of 16 m 3 ·s −1 (5). The current (2010 data) total water allocation represents only 4.4% of the mean annual Athabasca River flow (2) (Fig. 3); however, allocation and use as a proportion of average water levels does not account for the large variability in flow between seasons and years (interannual coefficient of variation = 22%). During 1952-2013, the Athabasca River gauge at Athabasca recorded mean seasonal flows of 100 m 3 ·s −1 in winter (December−February) and 911 m 3 ·s −1 in summer (June−August), with extreme monthly flows of 48 m 3 ·s −1 in December 2000 and 2,280 m 3 ·s −1 in June 1954. The mean, maximum, and minimum annual flows were 422 m 3 ·s −1 , 702 m 3 ·s −1 , and 245 m 3 ·s −1 , respectively (Water Survey of Canada, wateroffice.ec.gc.ca/search/ search_e.html?sType=h2oArc).
Assessing the sustainability of ARB surface water availability is difficult because most regional streamflow gauges have been operational for a few years to a few decades, with long intervals of missing values, including during the 1930s and 1940s when drought was prevalent throughout the North American western interior. Various critics (e.g., ref. 4) emphasize the intraannual variability and express concern about water withdrawals in the low-flow winter season. Less attention has been given to flow variability at interannual to decadal scales associated with climate oscillations such as the Pacific Decadal Oscillation (PDO) and Pacific North American mode (PNA), which are known to have significant impacts on runoff from the Rocky Mountains, including in the ARB (6)(7)(8)(9), and the potential consequences of a long period of predominantly low flow. We investigate the response of ARB flows to the low-frequency components of the PDO and PNA. This paper is the first, to our knowledge, to explicitly model the variability linked to the PDO and PNA, when testing for trends in discharge, thereby removing the tendencies associated with these climatic oscillations, which can confound trend detection in the relatively short instrumental records (6,7). However, the available instrumental hydrologic records provide a limited sample of these low-frequency fluctuations; exploring longer records of this variability is critical. Therefore, we also reconstructed the Athabasca River annual flow from the growth rings in moisture-sensitive conifers at a network of sites in the upper reaches of the ARB and in an adjacent watershed (Fig. S1). We compare 900 y of inferred

Significance
We show that current and projected surface water allocations from the Athabasca River, Alberta, Canada, for the exploitation of the Alberta oil sands are based upon an untenable assumption of the representativeness of the short instrumental gauge record. Our trend analysis of the instrumental data shows declining regional flows. Our tree-ring reconstruction shows periods of severe and prolonged low flows not captured by the instrumental record.
proxy flows to the streamflow variability recorded in recent decades, identifying extended periods of paleo low flow that exceeded the worst-case scenario in the instrumental records. We consider (i) the extent to which gauged flows (e.g., Fig. 3), which are the basis for surface water allocation, fail to capture the full range of hydroclimatic variability and extremes evident in a longer proxy record and (ii) the temporal evolution of low-frequency Athabasca River variability, related to the PDO and PNA, with its consequences for ARB surface water availability.

Results
We examined the nine longest and most complete ARB instrumental streamflow records for significant trends, while explicitly modeling and removing the possible effects of the PDO, the PNA, and the El Niño−Southern Oscillation (ENSO) using generalized least-squares regression (GLS) (Fig. 1, Table 1, and Tables S1 and S2). Low-frequency variability (i.e., data slightly smoothed by a 5-y binomial filter) was analyzed given the associated severe socioeconomic and ecological impacts of prolonged low flows. We found declining flows throughout the ARB: in the upper reaches (Athabasca River at Hinton and McLeod River), the middle reaches (Athabasca River at Athabasca), and the lower reaches (Clearwater River, Athabasca River at Fort McMurray, and the McMurray−Athabasca segment). The PDO or PNA are significant predictors of annual river flow variability for all nine records, whereas the trend and ENSO (the SOI term) enter in most, but not all, regression models.
We produced a tree-ring reconstruction of the water-year flow of the Athabasca River at Athabasca for AD 1111-2010. Table 2, Tables S3 and S4, and Fig. S2 give the statistical details for the five tree-ring models that comprise the nested reconstruction. The significant F statistics, positive reduction of error (RE), and similar values of standard error (SE) and root-mean-square error of validation (RMSE v ) indicate robust models with considerable predictive power. Values of R 2 adj range from 47% to 80%, with the exception of model 1, which extends back to AD 1111 based on data from three sites. Nonetheless, model 1 is statistically significant. For the five nested models, Fig. 4 shows the water-year reconstructions versus the gauged flow at Athabasca for the common period 1952-2010. For model 5, based on abundant tree-ring data since 1868, the two curves are nearly indistinguishable. As the models encompass more time, and thus tree-ring data from fewer sites, the corresponding modeled flow deviates from the recorded flow. The largest deviations result from the underestimation of high annual recorded flows. However, the tree-ring models consistently capture low flows. Water deficits, rather than excess water, are the major supply issue in the ARB.
The full Athabasca River reconstruction spanning AD 1111-2010 is plotted in Fig. 5. It indicates that AD 1937-1949 was the longest interval of low flow since the installation of hydrometric gauges on the Athabasca River, although the instrumental data are missing during this period (Fig. 3). Considering only the past several centuries, this modern drought of record is exceeded in severity and duration by periods of sustained low flow in AD 1790-1806 and AD 1888-1896. To depict the relative severity of these paleodroughts, we identify analogous years in the hydrometric record according to annual discharge anomaly (percent) ( Table 3). Thus, for example, the drought of AD 1888-1896 can be portrayed as the low flow years of the 2000s plus 1960 and 2 y for which there is no modern analog; discharges in these years were much lower than the minimum recorded annual flow. The drought of AD 1790-1806 also includes a year that has no modern analog, but the most notable property of this drought is its duration: 17 consecutive years of below-average water levels. These two paleodroughts are highlighted because they are relatively recent; however, neither one is a worst-case scenario. Droughts of greater intensity and/or duration occurred in preceding centuries, including the megadroughts of AD 1143-1198 and AD 1311-1409. The reconstructed flows before AD 1482 are based on model 1, with the least explained variance, lowest RE, and highest SE (Table 2). Thus, we have the least confidence here; the large range of reconstructed flows may be an artifact of the lower sample depth. However, the sequence of wet and dry years is valid, and the prolonged periods of low flow indicate the persistence of long dry spells. Fig. 5 is a single best streamflow reconstruction based on the limited number of tree-ring chronologies selected by the regression algorithm from the large pool of potential predictors: 43 ringwidth chronologies, (early-wood width and annual increments at   A wavelets evolutionary spectrum of the full reconstruction (AD 1111-2010) shows that the modes of variability have not been consistent over the past millennium, but rather differed among the putative Medieval Climate Anomaly (MCA; ∼AD 1111-1550) (11), the putative Little Ice Age (LIA; ∼AD 1551-1850) (11), and the modern period (AD 1851-2010) (Fig. 6). During the MCA, there are low-frequency modes of variability in the 16-to 20-, ∼32-, ∼55-, and ∼220-y bands. During the LIA, there is significant low-frequency variability in the 16-to 20-y band only. Variability in the higher-frequency 2-to 5-y band, representative of ENSO, is present throughout. Significant power at the ∼32-y band reappears during the modern period. These patterns are confirmed using a multitaper method (MTM) with its superior spectral estimation (12) (Fig. S4). MTM analysis confirms the presence of significant peaks (95% level) in the 2-to 5-, 16-to 20-, 54-, and 220-y frequency bands during the MCA (Fig. S4A). It also confirms the absence of the lower-frequency variability during the LIA, which instead has significant variability only at the 2-to 5-, 10-, and 20-y frequencies (Fig. S4B).

Discussion
Our GLS analysis of the instrumental data reveals generally declining flows throughout the ARB (Table 1), confirming concerns of Schindler and Donahue (1) and Peters et al. (8). Other recent studies of historical streamflow trends in the ARB have produced inconsistent results, reflecting the use of various gauge records of varying lengths and, to a lesser extent, different statistical methods; although most used variants of the low-powered Mann-Kendall (MK) nonparametric trend test. Results from Bawden et al. (13) show strong decreasing trends in annual, warm season (March to October) and summer flows over most of the ARB. They used 19 hydrometric records, most of them short (<50 y). Similarly, Rasouli et al. (14) analyzed trends only since 1960, discovering marked declines in the streamflow input to Lake Athabasca. In contrast, Rood et al. (9) and Peters et al. (8), like us, examined streamflow data back to 1913, spanning major data gaps or interpolating shorter ones. Rood et al. (9), using parametric linear regression, found declining trends in summer and in the upper reaches, but found no significant trends in the middle and lower reaches using the centennial length records (  (9) and Peters et al. (8) emphasized the teleconnections between the regional hydroclimatic regime and North Pacific climate oscillations (6,7). They referred to the associated pentadecadal hydrologic variability as the rationale for using only the longest ARB records. However, unlike our study, previous researchers have not modeled and removed the confounding transient trends of the low-frequency components of the PDO and PNA before performing the trend test. Our GLS analysis also correctly handles residual autocorrelation and thus allows a high-powered parametric test. With a higher ratio of trend signal to noise, from the explicit modeling of hydrologic variability related to the climate oscillations, our approach has an improved chance of trend detection. The use of more-tractable hydrometric data, filtered with a 5-y binomial smoother, is not a major limitation, as negative trends in these filtered data represent sustained periods of declining flows with potentially serious consequences. Ecosystems and communities in the ARB can cope with a single severe low-flow year or two, but a prolonged period of lower flows is much more challenging.  The collection and analysis of tree rings from dry sites in the montane forest of western Alberta (e.g., refs. 10 and 15-19) has shown that the rate of tree growth at these sites is limited by the availability of soil moisture. However, there has been no previous attempt to produce a long streamflow reconstruction for the ARB using a network of tree-ring chronologies from the headwaters of the Athabasca River basin. Most of the previous research into the paleohydrology of the ARB has been conducted near the mouth of the Athabasca River at the Peace−Athabasca Delta. There have been two tree-ring reconstructions of the past levels of Lake Athabasca (20,21); both were relatively short (<200 y). A major study of the paleolimnology of the Peace−Athabasca Delta by Wolfe et al. (22)(23)(24) provides important context for the interpretation of recent hydrological changes, although at coarser resolution than that enabled by annual tree rings. Rasouli et al. (14) interpreted trends in recorded streamflow in the context of the Wolfe et al. (23,24) millennial-length lake level time series, concluding that Lake Athabasca levels could drop by 2-3 m by the end of the 21st century.
Our 900-y reconstruction of water-year flow revealed that the long-term variability of the Athabasca River is considerably greater than recorded by the gauge near the community of Athabasca (Figs. 3 and 5). There are reconstructed pregauge periods of consistent low flows that far exceed the duration of the worst recorded hydrologic droughts. A prolonged interval of low flow occurs in the tree-ring reconstruction from AD 1936-1949, during the period of missing data in the gauged flow of the Athabasca River. This low-flow period is seen throughout southern Alberta in various instrumental records. This modern drought, and more severe hydrologic droughts from the paleorecord, such as in AD 1790-1806 and AD 1888-1896, should be considered in the allocation and use of ARB surface water. The paleodroughts include years when discharges were so low that there is no modern analog in the gauge records. Further evidence of extremely low flows in the AD 1790s comes from Hudson Bay Company journals from Fort Edmonton on the North Saskatchewan River (25). "Amazing" is used to describe the warmness of the winter and shallowness of the water. In the spring of 1796, low river flows caused damage to light canoes and prevented the transport of furs, timber, and supplies. Reports of extensive grass and forest fire, and diminished populations of bison and beaver, are further indications of unusually dry conditions. The adjacent Athabasca and North Saskatchewan Rivers have a similar hydroclimate and emanate from the same ice field. Thus, the ARB paleohydrology presented here shares many characteristics with an initial reconstruction of the flow of the North Saskatchewan River (26) and a subsequent reconstruction (16) using more recent tree-ring data and from more sites. Geographically extensive hydroclimatic events include centennial-length lowflow periods during the megadroughts of AD 1143-1198 and AD 1311-1409, which are well documented in proxy records from other regions of western North America (27).
Our reconstruction provides context that the short instrumental record lacks (Figs. 3 and 5). Never during the last 62 y has the average annual flow of the Athabasca River at Athabasca been less than 200 m 3 ·s −1 . Our 900-y reconstruction shows 36 occurrences of flow below 200 m 3 ·s −1 . Within the instrumental period of 1952-2013, winter (DJF) flow was, on average, 23% of water-year flow. If we assume this ratio roughly holds over a longer time period, then a year with an annual flow of less than 200 m 3 ·s −1 , would have had a winter flow of less than 46 m 3 ·s −1 . The Government of Alberta (5) estimates that oil sands water use will reach 16 m 3 ·s −1 in the next decade. This water is withdrawn from the Athabasca River in the Fort McMurray region, not at Athabasca. Flow at Athabasca is 68% of that downstream at Fort McMurray (1958−2012 average). If a low flow year of less than 200 m 3 ·s −1 occurs at Athabasca, this is roughly equivalent to 68 m 3 ·s −1 in winter at Fort McMurray. In this case, water permanently withdrawn for oil sands extraction is an appreciable amount, particularly if these low flows repeatedly occur for a decade or more, as our reconstruction shows can be the case.
Our spectral analysis demonstrates the low-frequency variability in ARB flows linked to teleconnections with large-scale climate oscillations ( Fig. 6 and Fig. S4). The GLS trend analysis of instrumental data confirms that the PDO and the PNA, with their similar low-frequency components (28), have a significant impact on historical ARB flows (8). In the 900-y proxy flow record, the MCA spectrum contains the modes of variability typical of the PDO and the associated North Pacific Index. These frequencies (i.e., pentadecadal and bidecadal) are also seen in spectra of modern The regression model equations are given in Table S4. Yr1, first year of reconstruction; P, P value of F statistic: ** < 0.0001, * < 0.001.  Table 2 and Table S4. instrumental records along with the ENSO band at 2-5 y (29,30). The LIA spectrum is very different, without a strong pentadecadal band. This suggests that the low-frequency variability of the PDO and PNA has varied between the MCA, LIA, and modern periods, or that the teleconnections between these climate oscillations and ARB hydroclimate naturally varies. Either possibility would imply nonstationarity in a major driver of ARB hydroclimate over the past millennium, as speculated by some researchers (e.g., ref. 31).

Conclusions
Water resource engineering, and decisions about water allocation, supply, and storage, have relied on the assumption of stationaritythat hydroclimatic regimes fluctuate within a consistent range of variability captured by relatively short instrumental records (32). This assumption is undermined by research on the hydroclimate of the past millennium (e.g., refs. [33][34][35], including our 900-y tree-ring reconstruction of the Athabasca River, which encompasses periods of sustained low flow spanning multiple decades, unlike any hydrologic drought recorded instrumentally. When these prolonged droughts reoccur, it will be in a much warmer climate, presenting significant challenges for ARB surface water management (36). Our GLS trend analysis, with its superior trend detection ability, shows declining flows in much of the ARB, consistent with regional climate warming and the loss of glacier ice and snowpack at high elevations in the Rocky Mountains (33). Our study provides evidence of multidecadal variability of ARB discharges, at the low frequency of the large-scale climate oscillations (i.e., the PDO and PNA). At the very least, risk assessment in the oil sands industry should consider the repeated decadal droughts that are a common feature of the regional hydroclimate but have not occurred since the industry was established. The industry, and government agencies that regulate the use of surface water, may also want to consider the implications of a future period of sustained low flow in a warmer climate. The long-term trends and variability in river flow examined in this paper extend beyond the ARB, and thus the lessons for water resource management are transferrable to watersheds across Canada's western interior.

Methods
We used the nine longest gauge records and segments from the ARB (i.e., of at least 55 y in length, approximately one full PDO cycle) for a GLS multiple linear regression trend analysis, modeling the low-frequency variability associated with the PDO and PNA and also the higher-frequency variability driven by ENSO, using the South Oscillation Index (SOI) (see Table S2 for more details) (6). We spliced some records from sequential gauges following Rood et al. (9), but performed only minor gap filling. These nine gauge records were the predictands in the GLS regression models ( Table 1). The predictors were the winter (November−March average) PDO, winter (December−February average) PNA and (previous year's June−November average) SOI indices, and a monotonic trend. Because streamflow is lagged and smoothed relative to precipitation by surface and subsurface hydrology, and these large-scale climatic phenomena act at interannual time scales, the stream records were lagged from the climate indices by 0 y, ±1 y, and +2 y, and a 5-y binomial smoother was applied to the hydrologic and climatic data. The statistical model used in this study is Q t = μ + λT t + β 1 x 1,t + . . . + βkX k,t + « t , t = 1,. . ., L, where {Q t } is mean daily streamflow; index t runs over L years; μ is the mean streamflow over these years; T t is a linear trend with coefficient λ; {x i,t , t = 1,. . .,L} is the i th explanatory variable; k is the number of explanatory variables; β i is the ith explanatory variable's coefficient; and {« t } is the residual time series, which is an ARMA(p, q) process. An optimum minimal subset of significant predictors and an optimum minimal ARMA(p,q) residual model were chosen using the corrected Akaike Information Criterion (AICc) statistic applied to all predictor subsets of size ≤6, and for all p + q ≤ 6. The nonzero significance of the trend coefficient λ was tested by the Neyman−Pearson statistic (RP) using the null model of the optimum set of explanatory variables (minus the trend variable if included in the optimum set) versus the alternative model of the optimum set of explanatory variables together with the linear trend (if not already added).
The Athabasca River at Athabasca gauge provides the longest continuous hydrometric record from 1952 to 2013, although there are earlier data for 1913-1930 (Fig. 3). Water-year (October 1 to September 30) flow over 1952-2013 was used for the calibration and validation of the tree-ring models. The tree-ring data are from 16 sites within and near the Athabasca River Basin in north-central Alberta, where, during 2007-2013, increment cores and cross sections were collected in old-growth coniferous forest. Mapped in Fig. S1 and listed in Table S3 are the tree-ring sites in the Athabasca River's upper reaches and in the adjacent North Saskatchewan River basin headwaters. Nearly all of the tree-ring width data were generated using the WinDendro tree-ring image analysis system. For most sites, we also measured the earlyand late-wood widths. Cross-dating of the tree-ring series was verified with COFECHA (37). We standardized the tree-ring series using ARSTAN (38), removing the nonclimatic growth trends in the tree-ring series through Table 3. Two severe sustained paleodroughts in the ARB in recent centuries and analogous years from the instrumental hydrometric record Year %Δ* Modern analog † Fig. 6. Wavelet evolutionary spectrum of the Athabasca River tree-ring reconstruction for AD 1111-2010. Significant (P < 0.05) power exists at the frequencies and time intervals depicted in red and enclosed by a black line.
conservative detrending based on the use of a negative exponential curve or a cubic smoothing spline (a low-pass digital filter with a 50% frequency response cutoff). The standardized ring-width series of various lengths were averaged for each site, using a mean value function minimizing the effects of outliers (39). Table S3 gives statistical data for the 16 standard index chronologies. The number of time series per site varies from 19 to 148. Significant interseries correlation, as high as 0.836, and high mean sensitivity, mostly greater than 0.3, indicate a strong common response to interannual hydroclimatic variability. At all sites, first-order autocorrelation over the common period (1868−2004) exceeds the 95% significance level of 0.134, matching the similar significant persistence in stream flow (0.265). Table S5 shows the many significant (P < 0.05) correlations of the ring-width index data with streamflow and annual and seasonal precipitation measured at Jasper, Alberta, from 1971 to 2011, demonstrating that these tree-ring chronologies are suitable for reconstructing Athabasca River flow. These correlations are across all seasons, reflecting the role of snowmelt water and spring rain in recharging soil moisture and the role of summer rain in sustaining it.
The tree-ring reconstruction of water-year flow followed standard methods, using stepwise multiple linear regression with forward selection. The pool of potential predictors consisted of the standard index chronologies from the 16 sites. There were one to three ring-width index chronologies per site, depending on the availability of measurements of early and late wood. Lagging the predictor chronologies, by up to 2 y behind and ahead of the streamflow data, accounts for an offset between climatic conditions in a given year and the response of tree growth versus water levels. The regression models were validated using the leave-one-out method. The optimal number of predictors per model ranged from 7 to 10, according to various criteria: R 2 adj , RMSEv, RE, and SE.