COMPARATIVE STATISTICAL HYDROCLIMATOLOGY OF GLACIAL AND NIVAL RIVERS IN SOUTHWEST Y U K O N AND NORTHWEST BRITISH COLUMBIA by SEAN WILLIAM FLEMING B.Sc, The University of British Columbia, 1994 M.S., Oregon State University, 1997 M.S., Oregon State University, 1998 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF T H E REQUIREMENTS FOR T H E DEGREE OF DOCTOR OF PHILOSOPHY in T H E F A C U L T Y OF G R A D U A T E STUDIES (EARTH A N D O C E A N SCIENCES) We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA October 2004 © Sean William Fleming, 2004 ABSTRACT A diverse suite of nonparametric statistical and time series analysis techniques was applied to historical streamflow data from five glacier-fed and four snowmelt-fed rivers in the southwest Canadian subarctic, in order to determine whether, and if so, how the two fluvial regimes respond differently to variability in climatic forcing. Four types of streamflow variability dynamics, corresponding to four timescales of climatic forcing, were considered. Methodological development was also performed as appropriate. Results were as follows. At the seasonal level, the annual hydrographs of glacial and nival rivers differed in virtually every aspect of their flow magnitude and timing. Glacial rivers exhibited higher water resource productivity throughout the year and an extended freshet, and glacial cover was of similar importance to basin scale in determining annual hydrologic cycle amplitude, probably reflecting negative mass balance trends. High-frequency interannual streamflow fluctuations were then corifirmed to be strongly attenuated by watershed glacierization, but the effect was statistically significant only for robust metrics of freshet flow and timing and overall yearly flow magnitude. The dampening mechanism responsible appears restricted to variability timescales of roughly two years or less, consistent with separate analyses of lower-frequency behaviour. Glacial and nival rivers in the study area were then demonstrated to exhibit selective teleconnectivity to the Arctic Oscillation and El Nino-Southern Oscillation, respectively, from a net annual water resource perspective; this reflects a broader decoupling of glacial and nival yearly discharge anomalies. Seasonal responses were more complex. At the longest timescale considered, presence or absence of watershed glacierization was shown to control whether rivers here grew larger or smaller, respectively, in response to historically-observed climatic changes. Progressive shifts in annual hydrograph form were also observed and linked to glacierization. The results were in all cases readily physically interpretable, underscored the dramatic impacts of glacierization upon water resource variability, and strongly confirmed and extended existing understanding of glacial-nival hydroclirnatologic contrasts. In particular, this is the first study to examine the comparative hydroclimatology of glacial and nival rivers at a broad range of variability timescales using a single dataset, affording a more comprehensive view of such watershed dynamical processes. ii TABLE OF CONTENTS ABSTRACT ii Table of Contents iii List of Tables vii List of Figures viii List of Symbols ix List of Abbreviations xi Preface xii Acknowledgements xiii CHAPTER 1 MOTIVATION AND APPROACH 1 CHAPTER 2 STUDY AREA, DATA, AND METHODOLOGY 5 2.1 Study Area and Data 5 2.2 Statistical Overview 12 2.2.1 Inferential Statistics and Hypothesis Testing 12 2.2.2 Distributional Characteristics 16 2.2.3 Nonparametric Statistical Methods 22 2.2.4 Notes on Method Selection, Development, and Implementation.. 25 CHAPTER 3 ANNUAL STREAMFLOW REGIMES 27 3.1 Introduction 27 3.2 Data Processing and Methodology : 28 iii 3.2.1 Data Processing • • • 28 3.2.2 Distributional Analysis: Timing 28 3.2.3 Distributional Analysis: Magnitude 30 3.2.4 Empirical Orthogonal Function Analysis 30 3.3 Results 34 3.3.1 Flow Timing 34 3.3.2 Flow Magnitude 36 3.3.3 Spatiotemporal Variability 37 3.4 Discussion 39 3.4.1 Hypsometry: Glaciers versus Snowpack 39 3.4.2 Generality and Relationships to Temporal Variability in Climate 41 3.5 Conclusions 43 CHAPTER 4 HIGH-FREQUENCY INTERANNUAL VARIABILITY 44 4.1 Introduction 44 4.2 Data Processing and Methodology 45 4.3 Results and Discussion 49 4.3.1 CVResults and Interpretation 49 4.3.2 Timescales of Variability Dampening 52 4.4 Conclusions 54 CHAPTER 5 CLIMATE MODE RESPONSES 56 5.1 Introduction 56 5.2 Data Processing and Methodology 58 5.2.1 Data 58 5.2.2 Methods 63 5.3 Results 68 5.3.1 EOF-Correlation Analysis 68 5.3.2 Hydrologic Composite Analysis 75 5.3.3 Composite Analysis of Surface Meteorological Effects 80 5.4 Discussion 85 5.5 Conclusions 87 iv CHAPTER 6 IMPACTS OF CLIMATIC CHANGE 88 6.1 Introduction 88 6.2 Data Processing and Methodology 90 6.3 Results 94 6.4 Interpretation 98 6.4.1 General 98 6.4.2 Baseflow 100 6.5 Implications and Limitations 101 6.6 Conclusions 103 CHAPTER 7 ECOLOGICAL IMPLICATIONS 104 7.1 Introduction 104 7.2 Ecological Background 105 7.2.1 Study Area 105 7.2.2 Data 106 7.3 Habitat Implications 108 7.3.1 Freshet Floods 108 7.3.2 Spawning Flows 108 7.3.3 Winter Habitat Availability 109 7.3.4 Attenuation of High-Frequency Interannual Variability 110 7.3.5 Selective Teleconnectivity 110 7.3.6 Differential Long-Term Trends I l l 7.4 Fish Species Richness 112 7.5 Conclusions 113 CHAPTER 8 SUMMARY AND SYNTHESIS 115 8.1 Review 115 8.2 Surnrnary of Results 118 8.2.1 Annual Flow Regimes 118 8.2.2 Interannual Variability Modification 119 8.2.3 Large-Scale Ocean-Atmosphere Circulation Patterns 119 v 8.2.4 Climate Change Impacts 120 8.2.5 Kjyal-RWthral Hydroecology 121 8.3 Synthesis I 2 2 8.3.1 General Synthesis I 2 2 8.3.2 Timescale Relationships 124 8.4 Recommendations I 2 ? References 130 APPENDIX 1 NOTES ON EOF ANALYSIS 155 A l . l General 155 A1.2 Gridding 160 A1.3 Rotation and Nonlinear Approaches 161 A l .4 EOF Analysis versus PCA 164 APPENDIX 2 MARKOV CHAINS, DESERIALIZATION, AND TRENDS 166 A2.1 Introduction 166 A2.1.1 Nonparametric Statistical Detection of Monotonic Trend 166 A2.1.2 Markov-Chain Processes and Deserialization 167 A2.1.3 Problem, Purpose, and Scope 169 A2.2 Method 170 A2.3 Results 171 A2.3.1 Mann-Kendall Test, AR( 1) Noise, and Deserialization 171 A2.3.2 Alternative Approaches 177 A2.3.3 Ancillary Results 179 A2.4 Discussion 180 A2.4.1 A Practical Example 180 A2.4.2 Serial Correlation from other Deterrninistic Processes 182 A2.4.3 Deserialization Assumptions and Real Physical Systems 183 A2.5 Conclusions 185 v i LIST OF TABLES 1 Hydrologic background Mormation 10 2 Leading-mode eigenvectors for three EOF analyses performed using different data pre-processing schemes 37 3 Hydrometric variables considered 46 4 Calculated CV values 50 5 Results of statistical comparison of CV values 50 6 Leading-mode eigenvectors and percentages of variance explained 71 7 (Symmetric) correlation matrix between LMUPC time series and climate indices 71 8 Stack-correlation matrix 74 9 Trend estimation results 95 10 Results of statistical comparison of best-fit trends 96 11 Fish species presence or absence 107 vii LIST OF FIGURES 1 Location map 6 2 1979 Alsek River annual hydrograph 6 3 Empirical CDF of Wann River daily mean flows 18 4 Empirical CDF of Wann River yearly mean flows 18 5 Empirical CDF of Alsek River annual minimum daily discharge 21 6 Empirical CDF of Dezadeash River annual interquartile range of daily flows 21 7 1986 Alsek River daily discharge 24 8 Empirical CDF of 1986 Alsek River daily discharge 25 9 Daily mean specific discharge averaged over period of record 29 10 Daily proportional cumulative flow volume, QC(T), over the average year . . . . 34 11 Empirical CDF of QC(T) 35 12 CDF of daily specific discharges over the course of an average year 36 13 Annual average regimes for air temperature and specific discharge 40 14 Power spectra 54 15 Standardized monthly streamflow over common record, mean annual cycle removed 59 16 Standardized ENSO monthly index 62 17 Standardized A O monthly index 62 18 Leading-mode unrotated principal component time series 69 19 Composite annual hydrographs of standardized hydrometric station data . . . 76 20 Composite annual regimes of mean temperature and total precipitation 81 21 Total annual flow volume 91 22 Case (1) - significance levels found for white noise time series '. 172 23 Case (2) - significance levels for AR( 1) time series 173 24 Case (3) - type-II error probabilities, linear trends plus white noise 175 25 Case (3) - slope M A E , time series consisting of trends and white noise 175 26 Case (4) - type-II error probabilities, linear trends plus AR(1) noise 177 27 Case (4) - slope errors, linear trends plus AR( 1) noise 177 28 Observed February median daily discharge, Yukon River 181 viii LIST OF SYMBOLS q river discharge (flow, streamflow) mean daily discharge over full record A basin area above streamgauge Rav mean annual runoff over period of record 0(1 Oy) of the order of 1 0 r TSS total suspended solids T water temperature i, j, k, r, n, y, v, r indices (of time, space, EOF mode, etc.) N, n, M, Q sample sizes in various contexts, with subscripts as appropriate X random variate; X - JC, for i=l,N xt i value of random variate, X a probability of a type-I error P probability of a type-H error p(.) probability that the argument, (.), is true; note lower-case p CDF(^f) cumulative distribution function of random variate, X PDF(JT) probability density function of random variate, X QC(T) daily proportional cumulative flow volume D Kolmogorov-Smirnov test statistic jl y^-mode eigenvalue indicator function with argument, x med(X) median value of random variate, X CV rank-based coefficient of variation R Pearson product-moment correlation coefficient z standard normal variate (or deviate) Rs Spearman rank correlation coefficient tsmw, S Mann-Whitney test statistic and sum inf / sup {XJ } inferior / superior; smallest / largest value of x, from the set i = 1, arg{k\ \y{k) = yQ | particular value of the argument, k, to y(k), such that y(k) = yo Qx(a) a • 100 t h quantile of random variate, X qmi„ annual minimum daily discharge qmed annual median daily discharge qmean annual mean daily discharge qmax annual maximum daily discharge qiQR annual interquartile range of daily discharges vr total annual flow volume tc annual hydrograph centroid tmax calendar date at which qmax occurs P(f) spectral power at frequency,/; note upper-case P dB decibels N(a) number of standard deviations above or below climatic normal p(l) autocorrelation function at lag, / +/'CT,>) composite value for positive / negative event year, corresponding to month of the year, T, and location, rj m linear regression slope ¥(tj) spatial data vector at time, u ke A*-mode eigenvector (empirical orthogonal fiinction) kQfti) A*-mode principal component time series Y spatiotemporal data matrix C data covariance or correlation matrix, as appropriate; upper-case C E eigenvector matrix A principal component matrix W Mann-Kendall test statistic sgn(n) sign function with argument, it c lag-1 serial correlation coefficient; note lower-case c e independent Gaussian (white) noise PWX deserialized (prewhitened) random variate (time series), X x LIST OF ABBREVIATIONS WRA Warm River at Atlin ARBR Alsek River above Bates River K R K L Kluane River at outlet of Kluane Lake WRAH White River at kilometre 1881.6 of the Alaska Highway T R K L Taldiini River at outlet of Kusawa Lake WRC Wheaton River at Carcross B C M Big Creek at the mouth MCRW M'Clintock River at Whitehorse DRHJ Dezadeash River near Haines Junction C L T Central Limit Theorem HCCD Historical Canadian Climate Database EOF empirical orthogonal function AR(1) first-order autoregressive ENSO El Nifio-Southern Oscillation PDO Pacific Decadal Oscillation A O Arctic Oscillation SLP sea level pressure SST sea surface temperature LMUPC leading-mode unrotated principal component DJFM December through March, inclusive E T evapotranspiration G C M general circulation or global climate model FISS Fisheries Information Summary System PCA principal component(s) analysis M A E mean absolute error xi PREFACE Don't believe anyone who tells you rivers have no time for stories. Linda Williamson, Listen to the River Listen! Wisdom is calling out. Reason is making herself heard. On the hilltops near the road and at the crossroads she stands. At the entrance to the city, beside the gates, she calls. Proverbs 8:1-3 The only good eigenvalue is a dead eigenvalue. Paraphrased from an e-mail to R. Haggerty from H Hollenbeck, 1997 ACKNOWLEDGEMENTS I wish to thank a number of people who helped make this dissertation possible. My advisor, Garry Clarke, offered tremendous support, advice, and intellectual freedom. I also wish to thank the additional members of my supervisory committee, Roger Beckie and Dan Moore, for their time. Dan deserves particular recognition for his repeated and valuable contributions of feedback and ideas throughout the completion of this research. My fellow EOS glaciology group members also provided useful specific or general information, advice, or constructive criticism at various points, as did Mike Church, Paul Whitfield, Ed Quilty, John Petkau, William Hsieh, and a number of journal reviewers and editors. The Canadian taxpayer also warrants recognition. This research was supported primarily by a National Science and Engineering Research Council (NSERC) Postgraduate Scholarship (PGS-B-242551-2001) and a University of British Columbia University Graduate Fellowship, and is a contribution to the Climate System History and Dynamics Program (CSHD), which is jointly sponsored by NSERC and the Meteorological Service of Canada. Additional support to the completion or presentation of this research was provided by the Thomas and Marguerite MacKay Memorial Scholarship, Walter W. Jeffrey Memorial Scholarship, UBC Department of Earth and Ocean Sciences, and Canadian Meteorological and Oceanographic Society. On a more personal note, I would like to thank my mother, Alice, for all that she has done; my wife, Kristine, for her understanding; Mackenzie and Peanut, for (sometimes) welcome distractions; and God, for providing interesting puzzles of natural history to solve, as well as the opportunity to do so. xiii CHAPTER 1 MOTIVATION AND APPROACH Interannual streamflow variability has long been a core issue in hydrologic science and engineering. Interest has only grown more acute in recent years due to increasing awareness of the fluvial impacts of ocean-atmosphere circulation phenomena, such as El Nifto events, and of historically-observed and model-forecast climatic warming. Watershed glacierization is emerging as a key link in the causal chain between climatic forcing and water resource responses. Glacier-fed (more properly, glacionival) rivers play a central hydrologic role in many regions, including much of North America, Europe, and the Indian subcontinent; they are typically spatially interspersed with nival rivers; and a growing body of evidence, discussed in detail in subsequent sections, suggests that the two types of fluvial system can respond quite differently to climatic variability. Glacierization may therefore introduce profound spatial heterogeneity, even within an otherwise uniform hydroclimatic region, into the already complex question of temporal streamflow variability. The potential impacts may extend well beyond traditional water supply, hydroelectric power, and flood hazard planning and management considerations. For example, streamflow time series are a useful alternative in empirical climate research to the spatial point measurements of meteorological data because they sample a much greater area and elevation range (e.g., Meko and Stockton, 1984; Cayan and Peterson, 1989; Burn and Soulis, 1991; Burn, 1994; Lettenmaier et al., 1994; Zhang et al, 2001). Moreover, meteorological driving forces interact with both each other and land surface and subsurface characteristics in very complex ways, so streamflow data can offer a far clearer picture of the net hydrological impacts of climatic variability than do individual meteorological time series. By the same token, however, accurate inference of climatic variability from hydrologic data requires a far stronger understanding than we currently possess of the filtering effects of watershed processes and properties, such as those potentially arising from the presence or absence of glacial ice. In addition, while distinctions between kryal (glacMy-dorninated) and rhithral (snowmelt-dominated) lotic ecosystems have been attributed solely to contrasts in the static, baseline conditions of glacierized and 1 nonglacierized catchments (e.g., Ward, 1994; Milner and Petts, 1994; Dorava and Scott, 1998; Dorava and Milner, 2000; Milner et al., 2001), it is conceivable that differences in the temporal variability of glacial and nival systems might also play a role. But most fundamentally, an empirically defensible comparative understanding of the impacts of watershed glacial cover upon interannual water resource variability is a key component - yet also one of the least-developed elements - of our knowledge of northern and alpine fluvial hydrology. The purpose of this dissertation, then, is to develop a more thorough understanding of similarities and differences in the interannual streamflow variabilities of glacial and nival river systems. A number of potential approaches to the question are available. I posit that it is preferable for such a study to possess the following properties, or perhaps more conservatively, that a study with these characteristics offers strong potential: • Comparative: Considerable attention has been devoted to the characteristics and dynamics of both glacial and nival rivers, and glacial-nival distinctions might plausibly be proposed on the basis of such knowledge. If our purpose is to specifically address similarities and differences between the two systems, however, direct comparative evaluation of forms and degrees of interannual streamflow variability is the more appropriate approach. • Empirical: Empirical science can be of great utility to investigations of complex earth systems. Such systems typically display behaviour which is either inherently stochastic, or deterrnined by complex interactions between coupled, possibly highly nonlinear, deterministic processes that are difficult to identify, describe, and parameterize effectively. Investigation of the patterns, nature, and controls of streamflow variability clearly represents a suitable venue for the application of statistical and time series analysis methods, which can in turn be reasonably expected to provide new physical insight, and additional focus and direction for the development of physically-based models. An empirical approach is primarily data- rather than model-driven, invokes a minimum of physical assumptions and preconceptions, and perhaps places a somewhat greater emphasis on exploration than detailed explanation. 2 • Water resource-oriented: Relative to a low-order headwater stream, a larger catchment integrates a more representative and useful range of watershed characteristics and processes (rather than being utterly dominated by glacial mass balance, for example), in general possesses greater biomass and biodiversity (e.g., Vannote et al, 1980; Milner and Petts, 1994; Roper and Scarnecchia, 2001), is more likely to be incorporated in extensive, long-term, and more-or-less continuous and permanent streamgauge networks (here, that operated by the Water Survey of Canada), and presents greater direct socioeconomic significance. To be of the broadest interest and utility, and also to make maximum use of available data, the study should thus be primarily concerned with downstream water resource issues and implications - that is, emphasizing rivers of substantial size monitored by the WSC. • Concerned with a range of timescales: Variability in climatic forcing occurs at a variety of timescales. In addition, scale issues play a central role in both surface and subsurface, and physical and chemical hydrology (e.g., Skoien et al., 2003). It is therefore appropriate to choose a variety of timescales of streamflow variability for analysis or, more specifically, a variety of potential differential glacial-nival streamflow behaviours corresponding to a range of timescales. In this study, a broad but carefully-selected suite of statistical techniques is applied to historical river discharge data to identify and describe glacial-nival contrasts in temporal water resource variability. Both univariate and multivariate methods are employed, with an emphasis on the use of techniques which are free of distributional assumptions. I use several decades of WSC streamgauge data from nine variously-glacierized watersheds in a climatically homogeneous region of southwest Yukon and northwest British Columbia, with catchment areas ranging from hundreds to tens of thousands of square kilometres, effectively emphasizing downstream water resource impacts. Four classes of phenomena, each representing a different timescale of variability, are considered. These are: • glacial-nival differences in the annual flow regime (seasonal timescales), which helps set the stage for subsequent analyses of interannual variability; 3 • glacial modification of the severity of year-to-year flow fluctuations (short interannual timescales); • contrasting responses of glacial and nival rivers to large-scale ocean-atmosphere circulation patterns (intermediate interannual timescales); and • differential glacial-nival streamflow trends under observed historical climatic shifts (long interannual timescales). As detailed later in this thesis, there exist a broad array of unanswered questions - some quite fiindamental - with regard to each of these types of differential glacial-nival hydroclimatic response, as well as with respect to how these phenomena may or may not relate to each other. To the best of my knowledge, this is the first comparative, empirical study of the water resource variability of glacial and nival systems to consider such a wide spectrum of timescales in a single region and using a single set of rivers and discharge data across analyses. The work is strongly unified in purpose and approach, but assessment of distinct phenomena at various timescales does imply a partitioning of the research. In particular, rather than attempting to force a strictly uniform data processing and analysis protocol upon all timescales of investigation, methods are selected for each topic of interest as deemed most appropriate. This leads to presentation of the hydroclimatological analyses and results in four, semi-free-standing sections (chapters three through six). These are preceded by a background chapter which describes the study area and data and provides a general methodological overview, and are followed by an evaluation of potential hydroecological impacts, a summary and synthesis chapter, and a list of references cited. Appendices provide additional background information, or incorporate material generated during the completion of this thesis which is relevant but considered out of place within the body of the dissertation. 4 CHAPTER 2 STUDY AREA, DATA, AND METHODOLOGY 2.1 Study Area and Data Southwest Yukon and northwesternmost British Columbia, Canada (Figure 1) possess a relatively unusual assortment of properties advantageous to a statistical analysis of the riverine impacts of catchment glacierization. While data availability is modest by the usual standards of empirical climate studies, the combination of reasonable spatial density and temporal duration of river discharge records, the spatial arrangement of hydrometric stations, and the overall glacial, climatological, and hydrological characteristics of the region render it a superb natural laboratory for studying the convergence and interactions of processes in the atmosphere, cryosphere, hydrosphere, and biosphere. A general description of the study area and data is provided below; certain aspects have also been discussed in Fleming and Clarke (2003) and Fleming (in press). Salient points will be reiterated or expanded upon as appropriate in subsequent sections. The study area has a subarctic, transitional maritime-continental climate. Although Yukon and B.C. climatic zones may be delineated in a number of ways, at a synoptic scale the region is largely internally homogeneous and distinct from surrounding areas (see Wahl et al, 1987 and Phillips, 1990 for summaries of Yukon weather and climate). The region is partially shielded from the Pacific by the St. Elias/Wrangell/Coast Range mountain complex, which houses the tallest peaks of Canada and the U.S., and is known locally as the Icefield Ranges. This leads to a markedly drier climate than adjacent, maritime parts of Alaska, but sufficient precipitation is available to support extensive taiga forests in valleys and glacierization at higher elevations; indeed, the largest icefields to be found outside Greenland and the Antarctic lie within and adjacent to the study area. The region is approximately bounded to the north and east by drier conditions and a consequent paucity of glaciers, to the south by a more temperate climate and data scarcity, and to the west by the climatic divide of the mountains. Winters are very cold, with essentially all precipitation 5 Figure 1: Location map. Filled and empty circles give locations of gauges monitoring glacierized and nonglacierized basins, respectively. WRA=Wann River atAtlin; WRAH=White River at kilometre 1881.6 of the Alaska Highway; KRKL=Kluane River at outlet of Kluane Lake; TRKL=Takhini River at outlet of Kusawa Lake; MCRW=M'Clintock River at Whitehorse; ARBR=Alsek River above Bates River; BCM=Big Creek at the mouth; DRHJ=Dezadeash River near Haines Junction; WRC=Wheaton River at Carcross. \ 01 E» IB •o I 800 700 600 -j 500 400 300 200 100 n n n l~l I n 1 2 3 4 5 6 7 8 9 10 11 12 month Figure 2: 1979 Alsek River annual hydrograph. Details vary between years and basins, but overall hydrograph form is typical of glacial and nival streams in the region, being dominated by the melt freshet, with very low winter base/lows sustained by aquifer and (where applicable) lake and wetland storage. 6 falling as snow, leading to freshet-dominated annual hydrographs (Figure 2). This is convenient for our purposes as it minimizes interpretational ambiguities that could arise from the addition of a substantial pluvial component to the overall river flow regime. Summers are mild, and daylengths exceed 18 hours in June and July. The region is strongly influenced by both Pacific and Arctic air masses and weather systems, and neither apparently dominates clearly or consistently during any season. Pacific influences are greater than in other parts of the Yukon. It is worthwhile to note, however, that the study area largely remains a frontier region to science, and only general properties of its hydroclimatology are understood. The Icefield Ranges form not only a climatic divide but also a continental drainage divide. However, this is a geologically very recent and possibly temporary development {Tempelman-Kluit, 1980; see also Bryan, 1972), and a number of rivers with headwaters in the continental interior penetrate this substantial but discontinuous topographic barrier to empty into the Pacific. Individual watershed boundaries can locally be topographically extremely subtle in spite of the considerable overall relief in the region, and both the current drainage patterns and their histories are unusually complex here. One of the study rivers (Kluane; see below) may have reversed flow directions and switched continental drainages in the recent postglacial Holocene {Bostock, 1969), and at least two of the study rivers (Dezadeash and Alsek) have experienced substantial, repeated drainage obstructions due to periodic surges of the Lowell Glacier, which formed large ice-dammed lakes within these watersheds; the most recent full blockage may have occurred around 1900 (Clague and Rampton, 1982). In addition, the terminus of the Kaskawulsh Glacier overlies a continental drainage divide, and the relative volumes of Kaskawulsh meltwater ultimately routed to the Kluane versus Alsek Rivers has been hypothesized to fluctuate from year to year (Bryan, 1972; Barnett, 1974; Johnson, 1986). While these complexities are presented mainly to provide context, one is of immediate practical concern to the current hydrologic analysis: a number of the study rivers are fed by (among other things) one or more alpine glaciers which, in effect, constitute lobes of the St. Elias ice field. Overall meltwater evacuation patterns may not follow ice surface slope or underlying bedrock topography in any straightforward way. This causes difficulties in quantitatively assessing levels of glacial influence within several of the glacierized catchments; similar challenges were encountered 7 by Bryan (1972) during his investigation of Slims River drainage v^thin the present study area, and by Braithewaite and Olesen (1988) in their Greenland hydrology work. While it is therefore unclear how conventional quantitative metrics of watershed glacial cover, such as percent-of-basin-glaciated, may be evaluated in a reliable and meaningful way for several of the catchments, qualitative rankings of the relative degrees of glacial influence are reasonably straightforward to assign between catchments on the basis of topographic maps from the region. Bedrock geology (Gordey and Makepeace, 2001) consists of extremely heterogeneous accreted terranes which show no discernable correspondence to presence or absence of watershed glacial cover, and so offers no potential to corrupt the hydrological comparisons performed here. The entire region was glaciated during the Last Glacial Maximum, leading to a reasonably uniform veneer of glacial, glaciofluvial, and glaciolacustrine sediments within river valleys {Fulton, 1995; Gordey and Makepeace, 2001). There is, however, some theoretical potential for currently-glaciated watersheds to possess more extensive and/or hydraulically conductive surficial deposits. While this hydrogeological contrast is purely conjectural and, even if present, does not seem likely to materially affect overall interannual flow variability, it could conceivably present impacts for average winter baseflows and aquatic habitat availability; this is discussed further in chapters three and seven. Although the current study is not limited to long-term glaciohydrologic shifts, such as those that may arise from global climate change, such longer timescales are indeed considered here (chapter six). Moreover, as we shall see, long-term shifts also apparently exert some control upon short-term behaviour, such as contrasts between the typical freshet magnitudes of glacial and nival systems (chapter three). It is important to note, then, that prior work here has identified regional trends in hydrology and meteorology, including a wanning trend, broadly consistent with (although not necessarily indicative of) anthropogenic climate change (Zhang et al, 2000; Whitfield and Cannon, 2000a,b; Zhang et al, 2001; Whitfield, 2001). Meteorological trends appear to be spatially homogeneous (Zhang et al, 2000; Whitfield and Cannon, 2000a,b; Whitfield, 2001). Note also that climate changes and associated hydrological impacts have been deduced from both paleoclimatological studies and G C M predictions to be most acute and, therefore, most readily identified and studied at 8 higher latitudes, such as those of the study area (e.g., Roots, 1989; Zaltsberg, 1990; Rouse et al, 1997; Zhang et al, 2000; IPCC, 2001b; Zhang et al, 2001). Water Survey of Canada river discharge data from the H Y D A T database were used (Environment Canada, 1999). Both daily and monthly records were utilized. Small data gaps were considered acceptable as the techniques applied here accommodate missing data points without the necessity of interpolation. A minimum record duration of 20 years was set in a compromise between the desire for longer datasets on the one hand, and for a larger number of datasets and adequate spatial coverage on the other. Only stations monitoring reasonably undisturbed watersheds were chosen (see Whitfield and Cannon, 2000a,b). Nine stations were ultimately selected, corresponding to five glacierized (Kluane, White, Alsek, Wann, Takhini) and four nonglacierized (Dezadeash, M'Clintock, Big Creek, Wheaton) basins; record lengths range from 22 to 47 years. Hydrometric station locations and details are provided in Figure 1 and Table 1. While the Dezadeash River is a tributary of the Alsek, the former provides a small proportion of latter's runoff, so it seems reasonable to treat hydrologic data from these rivers as independent samples for the purposes of this study; as we shall see, the two rivers do indeed display fundamentally different types of streamflow variability, increasing confidence that this approach is valid. The glacierized and nival basins are geographically interspersed, reducing concerns that systematic differences in the hydrologic responses of the two groups, if present, simply reflect spatial correlation or regional gradients in climate or in watershed properties other than glacierization. Drainage areas to the streamgauges are O(102 to 104) km2; the study catchments may therefore be categorized as large according to the classification scheme of Skvien et al. (2003). Such relatively sizable catchments place an emphasis upon downstream water resource characteristics, and largely eliminate the influence of aspect from the analysis as each watershed contains a multitude of hillslopes oriented in all directions. Depending on the analysis performed, either the full records from all nine rivers or only the mutually common period of record from eight rivers were used, as described in the following chapters. Similarly, data processing steps depend upon both the particular phenomenon under investigation and the specific methodology used for statistically assessing hydroclimatic behaviour at that timescale, and are therefore discussed in detail in 9 Table 1: Hydrologic background information. Station abbreviations are explained in Figure 1 caption. Glacier-fed rivers are ordered left to right with increasing degree of watershed glacierization (see text for details). Information is largely from Environment Canada (1999). ID is WSC streamgauge code; leading 0 is omitted here. Minimum and maximum elevations, from topographic maps, are approximate. A is catchment area above streamgauge; qm is mean daily discharge over full HYDAT record; mean annual runoff, Rm, was calculated from A and q&. Catchment lake and glacial covers as proportions of A were found using planimetric methods and rounded to nearest full percentage point. These glacial cover estimates are not reliable (see text) and are provided only for illustrative purposes. Total suspended solids (TSS) and temperature (T) data are spatially and temporally sporadic, and solely from April to October (typically June to August), but values are typical of freshet-time glacial and nival systems in general. Bering drainage is via the Yukon River. glacierized basins nival basins T R K L WRAH K R K L ARBR WRA WRC MCRW DRHJ B C M ID 9ac004 9cb001 9ca002 8ab001 9aa015 9aa012 9ab008 8aa003 9ah003 long (dec °) -136.12 -140.55 -139.05 -137.97 -134.21 -134.90 -134.46 -137.51 -137.02 lat (dec °) 60.61 61.99 61.43 60.12 59.43 60.08 60.61 60.75 62.57 mmelev (m) 750 750 750 450 750 750 750 450 450 maxelev (m) 8350 10500 10500 10500 7650 8050 6850 7600 6650 record start 1965 1975 1953 1975 1964 1966 1966 1953 1975 record end 1986 1999 1995 1999 1993 1999 1994 1999 1999 J(km z ) 4070 6240 4950 16200 269 875 1700 8500 1750 qav ( m V ) 52 110 76 218 7 8 10 43 7 (mm yr"1) 404 556 483 425 864 284 181 158 133 lake cover 5% 0% 8% 2% 4% 1% 1% 4% 0% glacial cover 3% 4% 6% 7% 9% 0% 0% 0% 0% TSS(mgUl) - 2942 - 348 - - 33 - 14 r ( ° C ) - 5.6 - 5.1 - - 14.4 - 8.9 drainage Bering Bering Bering Pacific Bering Bering Bering Pacific Bering 10 each subsequent section. Due to the local form of the annual hydrograph (Figure 2), it is reasonable to leave the time series parsed into calendar years as extracted from the H Y D A T database, rather than reorganizing into hydrologic years, a common but by no means universal processing step more appropriate to other hydroclimatic regions. A few of the study rivers have been the subject of, or otherwise included in, the small number of prior hydrologic studies that have been completed in the current study area. Salient results of previous work on these rivers, when available, will be discussed in subsequent chapters as appropriate. Note that no comprehensive review of the hydrology or hydroclimatology of rivers in southwest Yukon or northwest British Columbia appears to exist. The basic purpose of this thesis is to empirically determine, on the basis of historical discharge data, the impact of watershed glacial cover upon water resource responses to climatic forcing at a variety of timescales in a given geographic region. The overall approach used to attain this goal in all the analyses presented here is to delineate how the study rivers respond to climatic forcing at a given timescale, and then to determine whether and, if so, how the responses of glacial and nival rivers systematically differ. The underlying assumption is that no other watershed and/or data characteristics potentially capable of significantly controlling hydrologic responses of the types considered here are associated, either physically or by chance, with presence and degree of watershed glacial cover. As outlined above, the overall climatic homogeneity of the region, the general hydrogeologic characteristics of the study area, and the spatial interspersion of glacial and nival catchments provide us with some degree of confidence that this assumption is valid. In addition, two-sided Spearman rank correlation and Mann-Whitney tests (the nature and manner of implementation of these tests is discussed in subsequent chapters) were used to determine whether glacial cover is associated with: (1) lake cover, which is well-known to modify streamflow variability, at least at subannual timescales; (2) catchment area, which can result in a nonlinear scaling of many hydrological processes that cannot be adequately compensated for using simple linear transformations, like time series standardization or normalization to specific discharge (e.g., Eaton et al, 2002; Sivapalan et al, 2002; Skoien et al, 2003); and (3) record length. Impacts of dataset duration are also considered in a more 11 detailed fashion as appropriate in subsequent sections. No statistical evidence was found for any such associations. Regardless of what effects other watershed properties may or may not have upon interannual streamflow variability, they cannot, therefore, lead to the false positive identification of glacial impacts, although they could conceivably obscure such relationships. The apparent sole remaining association of potential interest in this regard is that between glacierization and hypsometry. This is, of course, a fundamental physical relationship insofar as the generally greater precipitation and colder temperatures of higher elevations contribute to basin glacierization; consequences of hypsometry will be discussed on a case-by-case basis as they arise. 2.2 Statistical Overview 2.2.1 Inferential Statistics and Hypothesis Testing Statistical methods can be roughly divided into four general types. (1) Descriptive statistical procedures simply consist of compact summarizations of data. Measures of central tendency or location, for example, collapse potentially very large volumes of information into a single, "typical" value, such as the arithmetic, harmonic, or geometric mean, median, or mode. (2) Inferential statistical methods use statistical theory to draw conclusions or make decisions on the basis of the available data. The procedure is often phrased in terms of a hypothesis test: the null hypothesis, denoted Ho, usually states the absence of a relationship or distinction, and this is tested against an alternative hypothesis, Hi, usually stating that a relationship or distinction is present. (3) Modelling and prediction procedures attempt to explicitly, although perhaps nonphysically, relate one or more dependent variables (predictands) to one or more independent variables (predictors). An example might be a regression model forecasting tomorrow's streamflow as a linear combination of, say, today's streamflow, air temperature, and precipitation. (4) A broad range of coupled statistical-physical modelling procedures exist. These include ensemble or Monte Carlo methods that use statistical distributions of physical parameters to drive process-based physical models, such as finite difference solutions to the groundwater flow equation, generating statistical distributions of output variables as the product; and coupled dynamical-statistical methods for El Nino forecasting, which combine a process-based physical ocean model with a 12 statistical atmospheric model. Clearly, this is not the only way to parse the wide field of statistical science, nor are these four approaches by any means mutually exclusive: inferential statistical tests are often phrased in terms of a descriptive statistic, for example, and modelling and prediction procedures are often used to test hypotheses. The primary means of investigation used here consists of a suite of inferential statistical methods. As noted above, these typically take the form of a hypothesis test in which one decides whether a null hypothesis of randomness or of no difference, effect, or association can be rejected in favour of an alternative hypothesis claiming some specific relationship. A typical Ho might be that two populations have the same mean, and this may be tested against the alternative that the means are different (a two-sided, or two-tailed, test) or that the mean of one dataset is larger than that of the other (a one-sided, or one-tailed, test). All inferential statistical tests, however exotic or mundane they may be, use a similar overall procedure: (1) define null and alternative hypotheses; (2) evaluate a test statistic, the magnitude of which generally increases with the strength of the relationship or distinction being tested. For example, the observed test statistic for the aforementioned means test will tend to increase with an increasing difference between the observed means of the two populations being compared; and (3) reject (or fail to reject) the null hypothesis on the basis of the likelihood that the observed value of the test statistic occurred by chance, the premise being that large (absolute) values of the test statistic are unlikely to happen by accident. The latter step is often completed by comparing the observed test statistic to critical values corresponding to specific probabilities of a type-I error (rejecting the null hypothesis when it is, in fact, true), denoted a. Typical values of a selected for hypothesis testing range from 0.01 to 0.10, depending on one's level of conservatism (e.g., Daniel, 1990). The exact probability,^, that the observed test statistic occurs by chance can also be directly determined in many cases, or bracketed by the use of tabulated test statistics, a and p are essentially synonymous; the former notation is used here to avoid ambiguity with the more general sy mbo l ,deno t ing the probability that the argument (.) is true, which is also employed in the thesis. The manner in which hypothesis tests are phrased is inherently conservative: our starting point is a presumption of dullness, and we only reject the null hypothesis of no interesting relationship if we can assure ourselves that there is a chance of only a 100% that we are 13 wrong, a being a small fraction. In addition, failure to reject the null hypothesis should not be interpreted as anything more than an indication that we cannot say much about the data in the context of the currently-selected Ho and Hi. Failure to reject the null hypothesis does not, per se, constitute either a proof of the null hypothesis or a disproof of the alternative hypothesis. Along these same lines, if a potential relationship is considered to be a signal and all other data components are noise, hypothesis tests essentially amount to a measure of the signal-to-noise ratio of the data, recast into a statistical decision-making framework; subtle signals may thus be missed. In general, there is a trade-off between the value of a selected and the statistical power of the test. Smaller values of a generally lead to larger values of ft, the probability of a type-II error: failing to reject Ho when, in fact, Hi is true, i.e., failing to identify a signal when one is indeed present. Statistical power is discussed in some further detail in appendix two. Overall, then, hypothesis tests only yield a statistically unequivocal conclusion when the null hypothesis is rejected. Most hydroclimatological statistical methodologies, such as the above means test, share a number of assumptions; three will be briefly discussed here. The first is that the frequentist approach to statistical inference is appropriate. The essence of this approach is that the observed frequency of occurrence is equivalent to the probability of occurrence; for example, if a particular event or result occurs once out of ten trials, then its probability of occurrence is taken to be one in ten. The competing approach is the Bayesian method, which draws a distinction between frequency and probability and aims to incorporate a priori information. Although philosophically appealing and of substantial utility in some circumstances, application of Bayesian methods can be difficult or ambiguous, and is at a relatively early stage in the earth sciences; in contrast, the frequentist approach used here is by far the more common method and has been proven to be of great utility. Second, inferential statistical hydroclimatologic methods generally assume stationarity. The stationarity condition consists in this case of time- and/or space-invariance of the probability density function of the random variate in question. The necessity for this assumption is well-summarized by Deutsch and Journel (1998) and Jury et al. (1991). The less stringent condition of weak or second-order stationarity, which requires the stochastic process to be of constant mean and variance, is satisfactory for most purposes in the earth sciences (e.g., Jury 14 et al, 1991; von Storch and Zwiers, 1999). This means that for time series, the mean and variance do not change systematically over time; for spatial data, that the mean and variance are not functions of location; and for spatiotemporal data, that the mean and variance are functions of neither space nor time. Note that stationarity in the variance is defined such that serial dependence does not constitute a violation, provided that the correlation or variance reaches a constant value at some lag, which may be considered a characteristic timescale or space scale for the process, and can be empirically evaluated as the integral time or space scale, range, or e-folding time or distance (e.g., Skeien et al, 2003). Stationarity in the variance and mean are sometimes referred to, respectively, as homoscedasticity (also spelled homoskedasticity) and, simply, stationarity. In theory, one might test for time series nonstationarity in the mean using techniques such as trend tests or split-sample tests for intervention. Stationarity in the variance could conceivably be evaluated using a straightforward split-sample F-test; several sophisticated (e.g., spectral/wavelet) methods have also been developed for this purpose, though these often assume specific statistical models (e.g., Cahill, 2002; Chen and Rao, 2002). An interesting potential issue in the general use of such methods, however, is that being phrased as hypothesis tests, and given the general characteristics of hypothesis testing described above, these approaches only yield strong results when there is powerful evidence to reject the null hypothesis, i.e., when the time series is strongly nonstationary. They cannot, therefore, be used to actually prove first- or second-order stationarity per se. As noted by von Storch and Zwiers (1999), it is often difficult to verify the assumption of even weak stationarity in practise, so it is generally accepted to be true in the absence of contradictory physical arguments. Indeed, considering the issue at a more fundamental level (see also Jury et al, 1991), Deutsch and Journel (1998) describe stationarity not as an assumption but, rather, as a decision which cannot be checked from data in most cases. Note finally that not all methodologies are equally sensitive to stationarity assumptions. Spectral analysis, for example, is highly susceptible - hence the removal of least-squares regression lines as a standard pre-processing step - whereas empirical orthogonal function analysis will simply report nonstationarity faithfully within the resulting modal maps and/or time series. Third, and of particular practical concern to us here, is that conventional, parametric statistical methods assume that the data are drawn from a normal, or Gaussian, parent 15 distribution. Failure of a dataset to satisfy the assumption of normality does not appear to lead to a universally liberal or conservative result. Rather, the farther the underlying distribution lies from the normal distribution, the less meaning the statistical test has. Minor departures can often be ignored, and in some cases - for example, if the data are drawn from a lognormal distribution - a simple transformation can render the data approximately normal. Techniques also exist for rendering an arbitrary distribution Gaussian, but these fitted transformations are complex and can induce a substantial loss of information (see Box and Cox, 1964 and Cahill, 2002). Discussions of the applicability of parametric methods to hydrologic data, and of statistical methods which avoid distributional assumptions, follow in the next two sections. 2.2.2 Distributional Characteristics Some of the distributional characteristics of the data analyzed in this thesis are briefly explored here. The distribution function (DF), more commonly known as the cumulative distribution function (CDF), is the probability that the random variate under consideration takes on a value less than or equal to some specified value. As such, it is the integral of the probability density function (PDF), which is the more common form for describing probabilistic distributions in theoretical work. The empirical CDF, also known as the empirical distribution function (EDF), of a dataset can be evaluated in various very slightly different ways; the choice generally does not affect the form of the CDF and grows asymptotically irrelevant with increasing dataset size. The following definition is most convenient as it corresponds in a simple way to memory space allocation, which is a useful computational characteristic, particularly for Monte Carlo analysis: CDF(X) = p(Xj xk (2) Equation (1) is a nondecreasing step function of xt taking on values of (in principle) 0 to 1, in increments of l/N or multiples thereof. While the random variate and, therefore, its true distribution may be either discrete or continuous, in the latter case the discrete nature of sampling yields the discrete empirical CDF described by the above equations. In practise, constructing an empirical CDF amounts to a data sorting and ranking procedure (for an example, see Jury et al., 1991), which is also used for both estimating quantiles (descriptive statistics used in various analyses presented in subsequent chapters) and randomization distributions (in Monte Carlo randomization significance testing; see below). The empirical CDF is a simple, stable, and assumption-free estimator of the statistical distribution of a dataset, unlike the PDF. The drawback to the CDF is that, since it is a monotonicalfy nondecreasing function, interesting features are not so visually distinctive as they are in the graphical presentation of a PDF. For example, a bimodal distribution shows dramatic twin peaks in a PDF, but appears simply as a flat spot in a CDF. This one limitation notwithstanding, empirical CDFs are used here to explore the general distributional characteristics of some of the datasets analyzed in subsequent chapters. Figures 3 and 4 illustrate empirical CDFs of 1993 daily and 1964-1993 yearly mean discharge data, respectively, from the Wann River. The daily data are presented as departures from the 1964-1993 mean annual cycle, and the yearly values shown consist of departures from the 1964-1993 hydroclimatic normal. Removal of the mean annual cycle from the daily data is important for an apples-to-apples comparison here because, as will be seen in chapter three, the seasonal cycle imparts its own strong signal to the daily discharge CDF but is, of course, absent from annually-discretized data. The daily data are strongly positively (right-) skewed and roughly lognormal, which is typical for daily discharges (although an appropriately-parameterized gamma distribution, for example, might also yield a good fit), whereas the yearly data show an approximately symmetric and roughly Gaussian distribution. Such shifts toward normality with increasing timescales of averaging have 17 0.75 + CDF 0.5 0.25 + ( . . .1 1...... ..... -10 10 20 ,3_-1« daily mean discharge (m s') 30 Figure 3: Empirical CDF of Warm River daily mean flows. Data are expressed as departures from 1964-1993 mean annual cycle. 0.75 C D F 0.5 0.25 yearly mean discharge (mV) Figure 4: Empirical CDF of Warm River yearly mean flows. Data are expressed as departures from 1964-1993 hydroclimatic normal. been previously recognized in hydroclirnatology and attributed (Lettenmaier, 1995; see also von Storch and Zwiers, 1999) to the Central,Limit Theorem (CLT). The C L T possesses more complexity than it is typically credited with (see Conover, 1971 for a relatively comprehensive description). The basic notion, however, is that if one additively combines 18 M random variables, XJ = l'M, to form a new random variable, Y, as would be done (for example) during averaging or integration, then Y will be assymptotically normally-distributed with increasing M , regardless of the forms of the distributions possessed by each of the M variables being averaged. The C L T is often cited as a justification for the assumption of normality and, as a consequence, the use of relatively simple, well-established parametric techniques. In many cases this is appropriate, as attested to by the extraordinary range of phenomena well-described by the classical bell curve. Unfortunately, the applicability of the C L T is far from universal. As noted, for example, by Lettenmaier (1995) and von Storch and Zwiers (1999), the C L T is often of little use in addressing finery-discretized hydroclimatological data because it only guarantees asymptotic convergence to a normal distribution, without providing any information about how much data or what timescale of averaging is required to attain a Gaussian form. In the present case, the result is that a statistical technique assuming normality is of little use in dealing with the strongly-skewed daily discharge data. Another question of timescale arises when we consider the CDF of yearly mean data. While far more Gaussian than the daily data insofar as the skewness has largely disappeared, one could argue from the lack of a highly convincing S-shape in Figure 4 that a uniform distribution describes the empirical CDF about as well. Perhaps the C L T applies in principle: if only we had a longer streamflow record and, therefore, more data to register the rare events that constitute the tails of a normal distribution, we would be more impressed. By the time we have enough annual samples to actually confirm this, however, large-scale climate system nonstationarities may have grossly altered the distributions of our hydrometric variables. Perhaps the strongest problem with the applicability of the CLT, however, is that many hydrologic processes are non-additive and, similarly, many hydrometric variables are derived from raw discharge data in a non-additive manner. Put simply, not everything is an average. Granted, there exist alternative theoretical rules of significant practical importance which describe the distributions to be expected from other types of processes. For example, the lognormal distribution is, in practise, at least as frequently encountered in hydrology as the normal distribution; this may be attributable to the fact that multiplicative processes lead 19 to lognormality in what amounts to a log-space application of the C L T (see Bras, 1990). Under certain conditions, annual maxima time series lead naturally to the Gumbel, or type-I extreme large-value, distribution (e.g., Bras, 1990). In addition, various subsurface hydrological processes which involve a range of transport timescales have recently been linked to the power-law and gamma distributions and, in turn, various types of fractal behaviour (e.g., Haggerty et al, 2000, 2002; Kirchner et al, 2000, 2001). However, these relationships do not help greatly with the question of using statistical techniques which assume a normal distribution, nor do there exist relationships describing the theoretical distribution for all hydrometric quantities. Some examples may help to illustrate the problem. The empirical CDF of annual minimum daily discharges of the Alsek River (Figure 5) is somewhat negatively (left-) skewed; this is in contrast to Gaussian symmetry, and is also unusual insofar as most hydrologic data that show asymmetric distributions are skewed to the right, as in the lognormal and gamma functions. More importantly, however, the clearly-defined near-zero CDF slope in the vicinity of the median is indicative of bimodality. The empirical CDF of annual interquartile range of daily flows in the Dezadeash River (Figure 6) appears roughly Gaussian at first blush, but a well-defined shallowing of the slope around the middle of the CDF is entirely uncharacteristic of a normal distribution and is, in fact, indicative of platykurtic behaviour. Indeed, it is now widely recognized that many hydrologic data types possess distributions which are not Gaussian and cannot be easily transformed to normality, and nonparametric techniques, which avoid distributional assumptions altogether, have thus become common in hydroclimatology. The foregoing exploration of the distributional characteristics of the data used in this thesis, although cursory, confirms that assumptions of normality are not broadly justified for the purposes of this investigation, particularly given that there is little theoretical basis in these examples for invoking the CLT. Application (or development, as necessary) of statistical methodologies which do not make any a priori assumptions as to the nature of the underlying statistical distributions is, therefore, warranted. 20 1 0.75 4-C D F 05 0.25 ~i—I—I—i—1—r 15 20 -T-+- -1—i—i—i—r~ 25 30 35 qmin (m3s1) 40 Figure 5: Empirical CDF of Alsek River annual minimum daily discharge. 1 -, -. • 1 • a a 0.75 -• • • CDF 0.5 • 0.25 - y • • 0 4 - A — i — i — I — i — i — i — I — i — ' — i — I — > — ' — i — 10 30 50 70 90 QIQR (mV1) Figure 6: Empirical CDF of Dezadeash River annual interquartile range of daily flows. 2.2.3 Nonparametric Statistical Methods Statistical methods that make no assumptions regarding data distributions are variously known as distribution-free, rank-based, or nonparametric. Note that we follow the practise, common within both the statistical and hydroclimatological communities, of using these three terms as synonyms; see Daniel (1990) for a more complete discussion of nomenclature. There are, in a very general sense, two kinds of such procedures, which will be referred to here as analytic and Monte Carlo. Analytic tests are essentially nonparametric equivalents to the usual litany of parametric statistical techniques described in standard statistics textbooks; similarly, then, they involve the calculation of a test statistic according to some standard, derived equation of typically simple form, followed by comparison to tabulated values. The Mann-Whitney median test for the difference between the location measures of two populations, for example, is a nonparametric equivalent to the parametric means test mentioned previously in this chapter. In one way or another, most or all of these tests make use of relative rankings of the data. Specific analytic nonparametric procedures used in this thesis include the Mann-Whitney, Spearman correlation, Mann-KendalL and Kendall partial tau tests. Full descriptions of these tests, including the mechanics of implementation and tables of the test statistics, can be found in Conover (1971), Daniel (1990), Kendall and Gibbons (1990), and Burn (1994); summaries or relevant details and quirks will also be provided as needed in subsequent chapters. An important additional advantage of this type of nonparametric test is that, in general, it can be used with ordinal (rank) data. This will prove to be of substantial utility in a number of instances throughout the thesis due to difficulties in reliably assigning precise quantitative measures of glacierization to a number of the study catchments (see first section, this chapter). The price of the increased robustness and flexibility of such methods is a potential decrease in statistical power; the evaluations that follow in subsequent chapters may, therefore, be conservative. Monte Carlo methods constitute a very general class of statistical techniques which simply use (generally computer-based) iteration to generate distributions in some way. Of particular interest to us here are Monte Carlo resampling procedures, which are in turn 22 divided into bootstrap and randomization (or permutation) tests. These are nonparametric statistical inference tools which avoid distributional assumptions by, in essence, using the data distributions rather than an assumed distribution; an additional, albeit related, attribute is an apparent insensitivity to nonrandom sampling. Good reviews are provided by Edgington (1995) and Manly (1997). Bootstrapping and randomization are generally identical in practise, except that resampling is performed with replacement in the former, and without replacement in the latter. It is not clear which technique, if either, is better, although some theoretical work appears to suggest that randomization can offer substantial power advantages over bootstrapping, at least under certain conditions (Corcoran and Mehta, 2002). Monte Carlo randomization is used in this thesis to perform nonparametric hypothesis tests on the significance of: empirical orthogonal function modes, an apparently new application for the technique; the linear Pearson product-moment correlation coefficient, an application first suggested by Pitman (1937); and particular patterns observed in long-term winter baseflow trends, a completely novel application. The latter example also illustrates the flexibility of resampling techniques. Conventional analytic procedures, whether parametric or nonparametric, consist of a set library of specific questions (are the means different? are the data serially correlated?) and pre-packaged approaches to answering them. In contrast, the hypothesis to be tested using randomization or bootstrapping is constructed by the user, and the distribution of the test statistic is developed directly from the data; as a result, the questions which can be asked and answered using such methods are limited less by procedural constraints than by the researcher's imagination. A few words regarding the relationship of descriptive statistics to distributional assumptions are also warranted. Such standard descriptive statistics as the mean and variance tie in very neatly with the normal distribution. In the absence of Gaussian distributions, or distributions which can be easily transformed to normality, such statistics lose much of their allure; moreover, they are less robust to outliers than quantile-based measures, such as the median and interquartile range. Note that outliers may reflect not only unusual physical conditions or spurious measurements, but also systematic non-normality in the data, such as that arising from skewness. Extensive use is therefore made in this thesis of quantile-based, and thus loosely nonparametric, descriptive statistical measures. However, the importance of doing so appears considerably less than that of applying nonparametric inferential methods, 23 particularly given that the mean and variance are not strictly defined solely for normal populations. Moreover, the choice of descriptive statistic must also be subject to physical common sense. Consider, for instance, 1986 daily discharge data for the Alsek River (hydrograph, Figure 7; empirical CDF, Figure 8). The annual cycle is not removed in this example. The data are clearly non-normal, and the annual median daily flow is considerably lower than the annual mean daily flow. In principle, the median is a more appropriate measure of central tendency in this case. In practise, however, we see that the median corresponds in time to early spring or late fall flows, whereas the mean provides a more intuitively reasonable summary measure of overall annual flow conditions, better incorporating the effects of the relatively short-lived but physically key summer freshet. The arithmetic mean thus remains a very useful descriptive statistic in this case, irrespective of the feet that the data being described are strongly non-Gaussian. discharge 1000 -f 400 + 600 800 + 200 + 0 mean =P 92 0 184 276 calendar day Figure 7: 1986Alsek River daily discharge. 24 CDF 0 i i I i i i i i i i i i i i i i i i i i 0 200 400 600 800 1000 discharge (mV1) Figure 8: Empirical CDF of1986 Alsek River daily discharge. 2.2.4 Notes on Method Selection, Development, and Implementation In general, emphasis was placed on the standard application of a suite of stock statistical methodologies, selected as appropriate for addressing particular physical questions. In a number of instances, however, it was felt that either the standard technique, or conventional statistical hydroclimatological wisdom, was in some way inadequate or inappropriate; or it was found that there appeared to be no existing, boxed statistical method available to properly answer a certain question. Consequently, a substantial amount of methodological development or refinement was also completed during the course of the dissertation research. Examples include blending the N-rule approach with Monte Carlo randomization modelling to obtain nonparametric estimates of empirical orthogonal function modal significance (chapters three and five); binomial and randomization tests upon hydroclimatological trend patterns (chapter six); and evaluation of the roles of Markov-chain processes and deserialization upon trend identification, using synthetic Monte Carlo models (appendix two). 25 Three approaches were used to perform data processing and statistical analysis: (1) calculation by hand, used for relatively short, one-time analyses; (2) manipulation or calculation in a spreadsheet program, usually for simpler tasks, performed on a personal computer; and (3) manipulation or calculation using Fortran software written for the purpose, generally for more complex tasks that could not be readily completed using other approaches, performed on a Sun workstation. The various Fortran codes made extensive use of subroutines from the Numerical Recipes (Press et al, 1992) and SLATEC (Jones, 1969) libraries. 26 CHAPTER 3 ANNUAL STREAMFLOW REGIMES 3.1 Introduction Differences between the typical annual hydrographs of rivers primarily fed by glacier-versus snow-melt, although noted previously by other workers (e.g., Meier and Tangborn, 1961; Meier, 1969), were first formally defined in a regime context by Church (1974). Seeking to broadly classify the flow regimes of northern Canadian rivers, Church (1974) observed that, in general, freshets in glacial rivers were much longer, and peaked much later, than those in nival streams, primarily as a result of the continued availability of glacial ice melt well past the date at which the annual snowpack is depleted. Such distinctions have since become widely-recognized, and our understanding of the characteristics and dynamics of both nival and glacial rivers has grown considerably. However, aside from some intriguing results presented by Fountain and Tangborn (1985; discussed as they arise in this chapter), little direct, systematic comparative analysis of the overall annual streamflow regimes of nival and glacial rivers appears to have been conducted following the seminal but schematic work of Church (1974). At least two considerations suggest that the issue should be revisited. First, it has been demonstrated that the streamflow responses of glacierized (currently glaciated) and nival catchments to climatic variability differ in a number of fundamental ways and at a variety of timescales (e.g., Fountain and Tangborn, 1985; Moore, 1992; Neal et al., 2002; Lafreniere and Sharp, 2003). It would seem, however, that a more thorough understanding of the average annual hydrographs of the two fluvial regimes should be a prerequisite to detailed investigations of how these characteristics vary differentially over time in response to climatic forcing. Second, there has been increasing interest over the last decade in the comparative ecology of glacial and nival rivers, but this work has been based on a relatively simplistic picture of sediment load, temperature, and spawning-season flow distinctions. A better understanding of average annual streamflow regimes is an important and appropriate 27 starting point for more thorough investigations of glacial-nival ecological distinctions, a topic we shall broach again in chapter seven. Differences between the average annual flow regimes of glacial and nival rivers are further explored here in a baseline study, elements of which are discussed by Fleming (in press), to serve primarily as fundamental context to the temporal variability analyses presented in subsequent chapters. The study consists of a comparison of glacial and nonglacial hydrologic regimes using Kolmogorov-Smirnov analyses of flow and timing characteristics (distributional analyses), followed by empirical orthogonal function analysis (spatiotemporal variability analysis). These approaches summarize overall flow regimes in a compact but effective manner. 3.2 Data Processing and Methodology 3.2.1 Data Processing Plots of average annual hydrographs for the study rivers over the duration of record are presented in Figure 9. Typical discharges and watershed scales of the study catchments vary by two orders of magnitude (Table 1), so throughout much of this study I adopt the common normalization of flow data by catchment area to facilitate interbasin comparisons. The resulting quantity (specific discharge) may be regarded as a measure of the water resource productivity of a basin. 3.2.2 Distributional Analysis: Timing It is well-known that watershed glacierization typically induces a longer freshet, extending later into summer and exhibiting a delayed peak. However, given the qualitatively apparent differences in discharge per unit catchment area between glacier- and snowmelt-fed rivers (Table 1 and Figure 9), it is worthwhile to examine this phenomenon more closely, explicitly distinguishing between flow magnitude and flow timing effects. This is accomplished by comparing daily proportional cumulative flow volumes over the course of the year, calculated as: 28 \q(t)dt & ( r > = 7 & . 2 j ^ observed where N is the length of the time series analyzed. The second is a method developed here as a nonparametric analog to Preisendorfer's (1988) N-rule, consisting of Monte Carlo randomization (Edgington, 1995; Manly, 1997) modelling of the eigenspectrum. The 7th mode is considered significant if its observed eigenvalue satisfies the following condition: j ^ observed * <*rE {j ^critical \ \p{j ^ * j ^critical ) = * "«]} (6) 32 where a is the probability of a type-I error. The randomization distribution of j mode eigenvalues, is generated by performing EOF analysis upon each of r = 1, NMC time series realizations, which are synthesized by resampling the original data without replacement, and holding the resulting eigenspectra to evaluate cumulative distribution functions for each mode; I is the indicator function. Equation (6) is evaluated by setting jl* = jkcrmcai in equation (7) once the randomization distribution has been determined. Note that any modal significance test is fundamentally one-tailed. No correction is made for the substantial serial dependence in these monthly data because we wish the procedure to identify any strong, regionally coherent signal as significant, regardless of its nature. As a simple pedagogical example, if 90% of the total dataset variance is explained by a single, spatially uniform AR(1) time series, such a signal would be of great physical interest to us irrespective of the fact that, in the strictest sense, it consists simply of red noise. Moreover, most of the autocorrelation in the non-deseasonalized time series considered in this section arises from the seasonal cycle (see also appendix two), which of course is our focus here. Third, restricting our attention in the following analyses to the leading mode (that which explains the greatest proportion of total dataset variance; i.e.,y" =1) serves as an informal but useful safeguard against overinterpretation of the EOF results. This is not uncommon in practise, as evidenced for example by the Arctic Oscillation and Pacific Decadal Oscillation index definitions (Thompson and Wallace, 1998; Mantua, 1997), and is heuristically justified by the observation that very often, the leading-mode spatial pattern and time series capture the most essential and interesting features of a spatiotemporal dataset. Finally, the proportion of variance explained by the leading mode, which in effect is tested using the aforementioned Monte Carlo randomization procedure, can also be qualitatively assessed for rneaningfulness. Higher percentages of total spatiotemporal variance explained (7) 33 are indicative o f a spatial pattern/time series combination that more effectively collapses overall spatiotemporal behaviour. Although fairly common, rotation procedures (including but not limited to the well-known varimax scheme) present a litany o f potential disadvantages and are contraindicated as a general processing tool unless the eigenvalues are degenerate (von Storch and Zwiers, 1999), which was not found to be the case here. Note also that Lins (1985; 1997) concluded that rotation did not materially alter (or improve) the interpretation o f streamflow EOFs. A more thorough discussion o f rotation is provided in appendix one. 3.3 Results 3.3.1 Flow Timing Glacier-fed rivers conserve a greater proportion o f their total annual flow budget for late summer, and indeed through to winter, than do nival catchments (Figure 10), irrespective o f the overall water volume available. The difference in temporal flow allocation primarily results from the continued availability o f glacial melt water to summer's end, well past the 1 -0.75 -QcO) 0.5 0.25 -0 1 92 183 274 365 calendar day Figure 10: Daily proportional cumulative flow volume, Qc(T) (dimensionless), over the average year. nonglacierized ——glacierized 34 date at which the finite annual snowpack is exhausted, and is consistent with a longer freshet. This melt availability arises primarily from the essentially unlimited (at an annual scale) glacial melt source, dubbed the urilimited-supply effect by Meier (1969), but likely also reflects progressive release from supra-, en-, and sub-glacial meltwater storage (Fountain and Tangborn, 1985; Fountain and Walder, 1998). A point o f particular note is that baseflow following the freshet (and continuing to the next year's melt season) is predominantly supported by aquifer storage and release. Apparently, a lesser degree o f groundwater depletion by the end o f summer and, therefore, allocation o f proportionately greater water volumes for winter thus also arise from the extended freshet in glacier-fed rivers. The C D F o f QC(T) gives the probability that, on any given day o f the year, the proportional cumulative flow volume from the start o f the year is less than or equal to some specified value. This empirical distribution is obtained by exchanging the abscissa and ordinate o f Figure 10 and performing an adjustment o f the new vertical axis based upon the number o f days in a year; the result is illustrated in Figure 11. Evaluation o f the resulting C D F s for Qc(T) Figure 11: Empirical CDF o/Qc(T). 35 glacierized and nonglacerized basins using a two-sided Kolmogorov-Smirnov test, preferable to the chi-square test in this application because the distributions are continuous (Daniel, 1990), confirms that the CDFs are different at a < 0.05. Although these distributions are means over space and time, they remain only samples of the true population distributions and a probabilistic comparison using the Kolmogorov-Smirnov test is clearly appropriate. 3.3.2 Flow Magnitude Results are shown in Figure 12. Unlike Figure 11, the temporal structure of the annual flow regime is not explicitly preserved in these CDFs, although the relatively simple nature of the annual hydrograph in this region permits interpretation of the CDFs in terms of seasonality. The difference between the distributions is highly significant (a < 0.01, two-sided Kolmogorov-Smirnov test). Glacier-fed rivers yield much higher flows per unit catchment area than nonglacierized watersheds. This again arises primarily from the far greater melt reservoir capacity of glacial ice relative to the annual snowpack; but unlike the case for flow 0.75 CDF 0.5 j 0.25 nonglacierized ^—glacierized 0.02 0.04 0.06 specific discharge (mV'km2) Figure 12: CDF of daily specific discharges over the course of an average year. CDF gives the probability that, on any given day of the year, daily mean discharge (normalized by drainage area) is less than or equal to some prescribed value. 36 tirning (see above), storage of meltwater on, within, and beneath the glacier does not play a role in this respect. Differences are greatest at high-percentile flows, corresponding to the melt season. Together with the results from the foregoing section, this indicates that the freshets of glacier-fed rivers are both longer and larger, due to both timing and magnitude effects. However, the entire distribution for glacierized catchments lies to the right of that for nival rivers. Of special interest here is that the lowest-percentile flows, which correspond to the winter discharge regime, are consistently and considerably higher (about double) for glacier-fed rivers in the study area. This is likely attributable to the impact of both a longer (as noted above) as well as larger freshet upon winter groundwater supplies. 3.3.3 Spatiotemporal Variability For all three EOF analyses, the leading mode was found to be significant by both randomization and eigenvalue separation criteria, and explained the vast majority of the total spatiotemporal variance. The three corresponding leading-mode PC time series are mutually associated at Pearson product-moment correlation coefficients of +0.94 or higher, and represent a generalized annual cycle with comparatively minor interannual variability. The leading-mode EOFs for the three analyses are presented in Table 2. Comparisons of the Table 2: Leading-mode eigenvectors for three EOF analyses performed using different data pre-processing schemes. Raw = discharge without normalization; specific = specific discharge; standard -standardized discharge. Results consist of map patterns, but tabular presentation is more convenient for our purposes. Proportion of the total spatiotemporal variance of each of the three datasets explained by the respective leading mode is also shown. run % glacierized basins nival basins var WRAH K R K L ARBR WRA WRC MCRW DRHJ B C M raw 95 0.49 0.23 0.83 0.03 0.03 0.02 0.09 0.02 specific 89 0.52 0.30 0.33 0.68 0.21 0.07 0.07 0.07 standard 74 0.38 0.32 0.38 0.39 0.38 0.33 0.36 0.26 37 eigenvector entries across individual basins and between data normalization schemes suggest that glacierization and drainage area are the key controls, and provide a clear picture of the relative influences of these two watershed properties upon the amplitude of the annual hydrologic cycle. Note that in the current study area, this is essentially synonymous with the freshet, which is in turn responsible for most of the total annual flow volume and thus generally indicative of overall river size. Specifically, under the strong normalization to zero mean and unit variance, the leading-mode eigenvector entries are quite uniform across the eight basins considered, and seasonal-and regional-scale spatiotemporal variance is largely explained by the generalized annual cycle alone. Under the weaker normalization by catchment area to specific discharge, eliminating the linear impacts of watershed scale but not those of glacierization, the combination of watershed glacial cover and a generalized annual hydrologic cycle explains the bulk of basin-to-basin seasonal-scale discharge variability, with glacierized basins exhibiting larger freshets and thus contributing more to the total variance, yielding consistently larger eigenvector entries relative to nival basins. Finally, without normalization, it appears that catchment scale, glacierization, and a generalized annual cycle are the dominant three factors, explaining in combination 95% of the total spatiotemporal variability in the raw discharge data. This is further investigated by applying rank partial correlation (Kendall, 1942; Kendall and Gibbons, 1990), which assesses multivariate association without assumptions of linearity or normality, to the leading-mode eigenvector entries obtained from the unnormalized data. A correlation of +0.56 was found between glacierization and EOF entry (holding the effects of catchment size fixed), and a correlation of +0.65 was found between catchment size and EOF entry (holding the effects of glacierization fixed). Both associations are significant at a < 0.05; a one-tailed test was used, which is appropriate in this case because our alternative hypothesis can be guided by strong a priori evidence that a correlation, if present, should be positive. The similarity of these observed partial tau statistics indicates that the influence of watershed glacierization upon the amplitude of the annual hydrological cycle is in the same direction, and of similar magnitude, to that of watershed area. This seems a powerful testament to the importance of glacial effects upon overall streamflow regimes within a given hydroclimatic region. 38 3.4 Discussion 3.4.1 Hypsometry: Glaciers versus Snowpack General hypsometric differences between glacierized and nival basins undoubtedly lead to greater winter precipitation in the former, which in turn promotes a larger and longer-lasting snowpack. In combination with later seasonal snow melting with increasing altitude, elevation effects upon snowpack thus produce a larger and longer freshet, potentially obfuscating the specific relationship between glacierization and runoff. However, given that the differences in the average annual hydrographs of glacierized and nival catchments are dramatic, that both groups possess similar minimum elevations, and that maximum elevations show only a general relationship to watershed glacial cover (Figure 9, Table 1), one might intuitively infer that glacierization is the dominant catchment property. Comparing the linear relationship between peak flow timing and elevation (reflecting snowpack effects) to the nonlinear relationship between peak flow timing and percent glacial cover for 32 rivers in northern Washington state, USA, Fountain and Tangborn (1985) also argued that differences between glacierized and nival freshet flows are primarily attributable to glacial, rather than hypsometrically-induced snowpack, effects. It is nevertheless worthwhile to further explore the issue, which is done here by comparing three composite annual time series: specific discharges averaged over glacierized and nonglacierized basins, and 1975-1995 mean monthly surface air temperature data averaged between the three temperature stations of the Historical Canadian Climate Database (HCCD) present within the study area (Burwash Landing, Haines Junction, and Whitehorse Airport; see Vincent and Gullet, 1999 and Mekis and Hogg, 1999). As with most meteorological information, these are low-elevation spatial point data and therefore are a poor absolute measure of temperatures throughout the region, but they provide a reliable portrayal of relative seasonal temperature variations. The time series are illustrated in Figure 13. Clearly, nival streams typically peak relatively early and are well into summer-fall recession by the time the year's maximum temperature is attained, whereas glacial streamflow and air temperature are very well-coupled and peak 39 in the same month (July, on average). Meier and Tangborn (1961) noted a similar concurrence between peak air temperatures and discharges for a glacial basin in Montana, and our conclusion is generally consistent with Meier's (1969) observation that peak glacial runoff production in a glacial basin in Washington occurred under a mid-summer combination o f minimum cloud cover, high insolation, and low snow/ice albedo. More sophisticated investigations were also conducted using lagged cross-correlation analyses o f the average annual time series o f Figure 13, as well as their first derivatives with respect to time so as to emphasize timing o f the annual peaks; cross-spectral analyses o f continuous monthly records; and comparisons o f the months during which peak flows and temperatures occurred over the period o f record, again using continuous monthly data. The ensuing volume o f results was not found to add substantially to the interpretation so, in the interest o f brevity, details are not presented here; they are nevertheless mentioned because they provided general quantitative support to the qualitative observations made above on the basis o f Figure 13. specific discharge month Figure 13: Average annual regimes for air temperature and specific discharges. 40 The key result is that the annual discharge regime of strictly nival streams in the study area is clearly snowpack-limited, whereas the annual flow regime of glacier-fed rivers is limited not by snowpack availability but by temperature (note that while supra-, en-, and sub-glacial meltwater storage offers a delay mechanism, the exact timing correspondence between glacial streamflow and air temperature in Figure 13 seems unmistakable). There are two potential explanations for such temperature-limited behaviour: annual snowpack in glacierized catchments may be consistently in excess of the annual energy available for melting, or large glacial ice reservoirs provide an effectively unlimited meltwater source on an annual basis. These explanations are essentially identical for our immediate purposes, since the former effect implies formation of permanent ice cover. In combination with the considerations listed above, and noting that such close coupling of temperature to flow timing requires an unlimited melt reservoir and therefore implies a coupling also to flow magnitude, it can be reasonably concluded that while hypsometrically induced precipitation differences between glacierized and nival catchments may help account for timing and magnitude differences in the annual hydrographs of the two regimes, these contrasts are primarily attributable specifically to the presence or absence of glacial ice. 3.4.2 Generality and Relationships to Temporal Variability in Climate The glacial-nival freshet timing and duration distinctions described here are a refinement of long-standing knowledge, and are believed to be general. Differences in freshet magnitude were predicted by Fountain and Tangborn (1985) under the assumption of negative glacial mass balance, but do not appear to have previously been empirically demonstrated in a comparative, regime context. These arise from physical considerations similar to those driving freshet duration differences and thus are also believed to be general. The winter baseflow differences identified here, apparently for the first time, are interpreted to arise primarily from greater and more extended melt-season aquifer recharge (or by the same token, less melt-season aquifer depletion), and in this sense would also appear to be geographically universal. Note, however, that local hydrogeological characteristics may also play a role: basin glacierization and aquifer properties may show some correlation due to the somewhat greater likelihood of encountering vertically and areally extensive, nonindurated deposits of relatively high-permeability, high-porosity till in currently-41 glaciated watersheds (see also Dorava and Scott, 1998), and Bradford and Grout (2000) attributed higher baseflow and greater winter habitat availability in the higher reaches of the Yukon basin to greater glacial till abundance. The analyses should therefore be replicated elsewhere to confirm the geographic generality of baseflow differences. This suite of regime distinctions may, however, be in part a product of our times. It was inferred in the previous section that annual glacial streamflow is temperature-limited, and that the greater specific discharges of glacierized basins are derived specifically from glacial melt. This is consistent with ablation in excess of accumulation, which is in turn fully consistent with climatic trend analyses (Whitfield and Cannon, 2000a,b; Zhang et al, 2000; Whitfield, 2001), glacial mass balance data (Arendt et al, 2002), and hydrologic trend analyses (chapter six) from the study area. These indicate regionally approximately uniform warming trends, a strong negative trend in net annual balance (indicating progressive loss of glacial ice mass), and positive and negative total annual flow volume trends in glacial and nival rivers, respectively; the latter has been interpreted to reflect the combination of increases in both glacial melt generation and evapotranspiration. The excess of ablation over accumulation implied by the current analysis is also consistent with Fountain and Tangborn's (1985) model suggesting negative mass balance as a basis for higher predicted specific runoff in glacier-fed rivers. That negative mass balance trends are the global norm for alpine glaciers due to historical climatic warming (Dyurgerov and Meier, 2000) further suggests the geographic universality of the characteristics identified here. However, this may also imply that such baseline streamflow regime differences are in fact specifically reflective of climatic warming conditions over the period of hydrologic record, a possibility apparently first noted by Meier (1969). Nevertheless, glacial-nival flow timing and magnitude distinctions would likely persist in muted form under long-term neutral or positive net mass balances due to both the three-way association between glacial cover, hypsometry, and snowpack, and the additional meltwater storage and progressive release offered by various glacial mechanisms. It is not entirely clear what effect inter-regional variation in mass balance cycle amplitude (for a given net balance) might have upon glacial-nival streamflow distinctions. 42 3.5 Conclusions Growing interest in the differential responses of glacial and nival rivers to climatic forcing, and in ecological distinctions between the two streamflow regimes, suggests the need for a better comparative understanding of how the annual hydrologic cycle differs with presence or absence of catchment glacial cover. In this chapter, timing and magnitude characteristics of the average annual hydrographs of five glacierized and four nival catchments in the southwestern Canadian subarctic were empirically identified and compared. The chief statistically significant hydrological conclusions using Kolmogorov-Smirnov and empirical orthogonal function analyses are: (1) watershed glacial cover results in freshets that are longer, larger, and peak later than those experienced by the nival regime; (2) the winter baseflows of glacial rivers are also much higher on a unit-catchment-area basis; and (3) basin scale and degree of watershed glacial cover are of comparable importance in determining the magnitude of the annual hydrologic cycle. These differences arise from the greater availability, both in volume and over time, of meltwater in glacial catchments, which in part reflects the consistently negative alpine glacial mass balances observed both in the present study area and globally under historical climatic warming. 43 CHAPTER 4 HIGH-FREQUENCY INTERANNUAL VARIABILTY 4.1 Introduction The most widely-demonstrated effect of watershed glacial cover upon interannual streamflow variability is an inverse relationship between the severity of year-to-year streamflow fluctuations and percent-of-basin-glaciated (Henshaw, 1933; Meier, 1969; Krimmel and Tangborn, 1974; Zuming, 1982; Fountain and Tangborn, 1985; Ferguson, 1985; Braithwaite and Olesen, 1988; Chen and Ohmura, 1990a; Moore, 1992). The following interpretation, presented by Meier (1969) and elaborated upon by Krimmel and Tangborn (1974), appears widely accepted. In high-precipitation years, nonglacierized catchments experience increased runoff due to the greater snow volume available for summer melting. Snow has a significantly higher albedo than ice, however, so a greater snowpack reduces the amount of glacial melt water released. This glacial water storage causes a dampening of flow increases during years of anomalously high precipitation. In low-precipitation years, when nonglacierized catchments experience a smaller snowpack and decreased runoff, low-albedo glacial ice becomes exposed more readily, and glacially-derived melt thus increases to compensate. Some studies suggest that increases in basin glacierization induce decreasing year-to-year runoff variability only until a critical value of ice cover is reached. Beyond this value, which appears to differ among study areas, streamflow variability begins to increase and ultimately reaches levels comparable to a nonglacierized catchment (Fountain and Tangborn, 1985; Braithwaite and Olesen, 1988; Chen and Ohmura, 1990a; Moore, 1992). Moore (1992) posited that this behaviour arises from an increasing sensitivity of basin runoff to fluctuations in glacial mass balance. Downstream from the glacier terminus, however, where the proportion of the basin glacierized is relatively small, the dominant effect is decreasing streamflow variability with small step increases in glacierization. Note also that all the aforementioned studies use the coefficient of variation (CV) of historical data as a measure of interannual variability; since year-to-year sfreamflow fluctuations typically far exceed longer-term changes, such as those 44 arising from monotonic trend, it seems reasonable to ascribe the variability captured by CV results primarily to relatively short-term hydrologic fluctuations of the type considered in the above physical model. Such glacial dampening of flow variability is of widespread interest to the planning and management of hydroelectric power and water supply resources (see Meier, 1969; Tangborn, 1984; Fountain and Tangborn, 1985; Moore, 1992); there may also be repercussions for the use of hydrometric data in climate analysis, and for kryal-rhithral ecological distinctions (see chapters one and seven). The importance of the phenomenon notwithstanding, a number of fundamental issues with respect to glacial dampening of streamflow variability remain largely or wholly unaddressed. Chief among these is a more precise analysis of how the interannual variabilities of various hydrograph components each relate to degree of watershed glacial cover. With the partial exception of Fountain and Tangborn (1985), previous work has been limited to investigation of a single streamflow metric: annual runoff, or some very similar quantity. It is also of interest, however, whether ice cover-runoff variability relationships apply to other salient hydrometric variables, such as annual minimum, median, and maximum daily discharges, peak flow timing, and amplitude of the annual freshet. Additional issues include applicability of these relationships to subarctic regions, which are climatologically, glaciologically, and hydrologically distinct from previous study areas; the greater appropriateness of rank-based statistical methods, in contrast to the parametric techniques used in prior work; and estimation of the time scales at which such dampening processes operate. All of these issues are addressed here; Fleming and Clarke (in review) summarize some of these results. 4.2 Data Processing and Methodology To maximize the utility of temporally limited data, a fully common period of record was not used, and record lengths range from 22 to 47 years. Note that CV analysis is likely relatively insensitive to dataset duration insofar as the results do not depend directly upon the temporal structure of the data but only upon their net relative variance; available empirical evidence appears to support this (Moore, 1992). Daily discharge data from each river were then 45 processed {Jones, 1969; Press et al, 1992) to obtain annual time series for the six magnitude and two timing metrics described in Table 3. Table 3: Hydrometric variables considered (continued on next page). All metrics are derived from daily streamflow data, cjj, and evaluated independently for each study river for each year, i, of a ~N-year record. Time index] - 1,M in these formulae is the sequential calendar (Julian) day of a given year; M = 365 or 366 depending on whether year in question is a leap year. Metric Definition Description and comments flow magnitude measures: Annual minimum daily mean discharge; measure of dominantly groundwater-derived winter baseflow; Qmin inf fe,} discbarge measurements under winter channel ice cover (m3/s) are potentially subject to larger observational errors than those under open-water conditions (e.g., Pelletier, 1990; Walker and Wang, 1997; Moore et al, 2002). > Annual median (50th percentile) daily mean discharge; qmed (m7s) 0,(0.50) tends to correspond here to latest winter/earliest spring, or latest fall/earliest winter flows; may be subject to the same caveats as qmi„ Qmean (m3/s) Annual mean daily discharge; the only nominally parametric variable used; qmed and qmean differ substantially due to non-normality of qj (Figures 8, 9) Annual maximum daily mean discharge; extreme event Qmax (m3/s) supfej measure; tends to reflect a combination of peak freshet conditions and a summer rainstorm, which may be spatially and temporally erratic in its occurrence Interquartile range of annual daily mean flows qiQR (m3/s) Qq(0J5)-Qq{0.25) (difference between 25 th and 75th percentile daily flows in a given year); measure of freshet magnitude and (see Figure 9) amplitude of the annual hydrologic cycle 46 VT (m3) M \qjdj Total annual flow volume; largely accumulated during the freshet; because of more sophisticated integration technique (Jones, 1969) and consequent differences in the way data gaps are handled, VT* M • qmean', robust measure of the in-channel portion of the hydrologic cycle, and thus of the total water volume available for human and ecological needs (see also chapter six) flow timing measures: tc (d) arg \qtdi J=i = 0.5 Centroid of the annual hydrograph; sequential calendar (Julian) day at which 50% of the year's total annual flow volume has occurred; a measure of freshet timing; see also Zhang et al. (2001) (d) Sequential calendar day at which qmax occurs; another ars\ J | {qj = qmax)} measure of freshet timing; generally subject to the same caveats as qmax Coefficients of variation for the 72 resulting time series were then calculated using an expression following very roughly from Henshaw (1933): CV = med me d{x) (8) where X = (xt, i = are the values of the annual metric under consideration (e.g., annual maximum daily flow) over a jV-year record, and med(X) is the median value of X over the full record duration. Both the above expression and most of the streamflow metrics (Table 3) are fully nonparametric, as are the statistical inference techniques described below. Relationships between degree of catchment glacial cover and calculated CVs were investigated for each hydrometric variable using Spearman rank correlation. The method is 47 free of distributional assumptions, robust to nonlinearity and outliers, and requires only ordinal measurement scales. The latter property is advantageous because, as discussed in detail in a previous chapter, accurate quantification of percent-of-basin-glaciated is challenging in the study area, whereas qualitative determination of relative degree of glacierization between basins is relatively straightforward. Ranks were assigned manually to both glacierization and CVs according to the protocol described, for example, by Daniel (1990); in the case of ties, each tied value was assigned the mean of the rank positions for which it was tied. The Spearman correlation coefficient, Rs, was evaluated as the conventional linear Pearson product-moment correlation coefficient, R, between the ranked data (e.g., Conover, 1971); the linear correlation coefficient between two populations, A'and Y, each of size N is in general: For each metric, the null hypothesis of no association was tested against the alternative hypothesis of an inverse relationship. Although this one-tailed test was deemed most appropriate for the current application, a two-tailed re-analysis indicated that conclusions drawn from this dataset are not strongly sensitive to such choices. As the number of samples is small (nine), Rs could not be readily transformed into a z-score; Rs was therefore used as the test statistic and manually compared against tabulated critical values (see Daniel, 1990). We make two procedural notes as follows. First, graphical plots (not shown for conciseness) were made of rank(glacial cover) versus rank(CP) to assess each metric for evidence of increasing CV at very high glacierizations, which would lead to a concave-up curve (see introduction). No such evidence was found, so Spearman correlation methods, which do not assume linearity but do require monotonicity, are appropriate. Second, for a given metric, the median CV value over all glacierized basins was compared to the median across all nival catchments using a Mann-Whitney test, as a check against the performance N (9) 48 of the Spearman correlations. This is potentially of some importance, as rank ties are ubiquitous in the Spearman analysis since four of the nine basins possess the same level of glacierization (i.e., none), and an abundance of ties can in some cases compromise the test. In contrast, level of glacierization is not ranked in the medians comparison and ties in CV are very rare, thus avoiding this complication. The Mann-Whitney procedure was implemented as described in detail by Daniel (1990). The overall procedure, for a particular hydrometric variable, consists of pooling and manually ranking all the CV values from both populations (i.e., glacial and nival), dealing with the rare tie in the manner described above for Spearman ranking; surnrning the pooled ranks thus determined for one of the populations (the choice of which population is arbitrary); calculating the observed Mann-Whitney test statistic: ft = 5 - ^ k ± ! ) (10) where H / is the number of data in the population selected for the rank summation and S is the aforementioned sum; and comparing the result to a table of critical values. To maintain consistency with the rank correlation analysis, one-tailed tests were completed; note that the highest a for which a tabulated critical value of the one-sided Mann-Whitney test statistic is readily available is 0.05 (Conover, 1971; Daniel, 1990). In the current analysis, conclusions drawn from the Mann-Whitney tests were in all cases fully consistent with those from the correlation tests. 4.3 Results and Discussion 4.3.1 CV Results and Interpretation Results of the CV analysis are given in Tables 4 and 5; these should be interpreted in light of the comments of Table 3. Although the median value of every hydrometric variable considered was lower for glacier-fed rivers than for nival catchments, the differences (by both Mann-Whitney and Spearman tests) were significant only for certain quantities. Specifically, interannual variabilities in robust measures of flow magnitude during the 49 Table 4: Calculated C V values (dimensionless). metric glacierized basins nival basins T R K L W R A H K R K L ARBR WRA WRC M C R W DRHJ B C M Qmin 12 9.9 33 13 12 11 13 20 63 Hmed 17 17 14 22 25 15 24 15 37 qmean 6.9 13 9.1 10 8.8 13 12 16 24 qmax 17 12 11 15 16 14 27 25 29 qiQR 18 21 20 16 18 20 28 32 38 VT 6.9 12 9.2 7.9 8.8 13 12 17 23 tc 0.9 2.5 1.8 1.0 1.5 3.1 . 4.2 3.1 6.6 tmax 3.0 9.0 2.4 6.2 6.0 7.1 3.6 7.7 16 Table 5: Results of statistical comparison of CV values. Median CV values across glacierized and nival basins are med(CV)G and medCCVV, respectively. Rank correlation coefficient is Rs and does not assume linearity. The most restrictive significance levels at which the null hypotheses of equal medians between glacial and nival groups, and of no association between CV and degree of glacial cover, are rejected using two-tailed tests free ofdistributional assumptions are On,w and cts, respectively. metric medians testing correlation med(CV)G med(CV)n CCmw Rs as qmin 12 17 >0.05 -0.28 >0.10 qmed 17 20 >0.05 0.08 >0.10 qmean 9.1 14 0.05 -0.54 0.05 qmax 15 26 >0.05 -0.55 0.10 qiQR 18 30 0.05 -0.73 0.025 VT 8.8 15 0.025 -0.71 0.025 tc 1.5 3.6 0.01 -0.77 0.025 tmax 6.0 7.4 >0.05 -0.37 >0.10 50 freshet, or which are otherwise largely reflective of freshet flows, exhibited statistically significant inverse relationships with degree of watershed glacial cover (qmean, OJQR, VT); these relationships were also found here to extend to freshet timing (fc). No correlation was unambiguously found with the CVs of qmax and tmax, probably because the highest daily discharge of the year tends to arise from a combination of freshet conditions and a summer rainstorm; the latter is spatially and temporally sporadic and may not bear any relationship to catchment glacierization.. However, Rs for qmax was only marginally insignificant at a< 0.05; application of a more liberal but still reasonable and common criterion (a < 0.10; e.g., Burn, 1994; Zhang et al, 2001; see also Daniel, 1990) would thus alter our conclusions for qmax. Interannual variabilities in baseflow (qmm), as well as qmed, which is something of a hybrid baseflow/typical annual flow metric (see chapter two and Table 3), showed no signs of a relationship to glacierization. Snow- and ice-melt contributions are negligible during non-freshet periods, clearly reducing systematic differences in variability between nival and glacial systems. However, glacial-nival variability contrasts in freshet discharge might imply variability contrasts in summer recharge to aquifer, lake, and wetland storage and, therefore, in winter baseflow conditions. The independent variability-dampening effects of these storage mechanisms may further level the playing field between glacial and nival streams during the baseflow period. It is conceivable that greater measurement error associated with channel ice cover conditions might yield less reliable and thus more scattered CV estimates, and thereby also contribute to a lack of significant relationships (see Table 3 notes). The moderating effects of glacierization upon interannual variability in the summary hydrograph measures considered in previous studies (qmean and vr, or metrics based thereon) thus apply also to the subarctic. Much previous work considered relatively small alpine glaciers in a temperate maritime Pacific climate (Henshaw, 1933; Krimmel and Tangborn, 1974; Fountain and Tangborn, 1985; Moore, 1992); these results were extended to central Asia (Ferguson, 1985), northwestern China (Zuming, 1982), and continental Europe (Chen and Ohmura, 1990a). Braithwaite and Olesen (1988) found that total annual runoff was less variable in a glacierized catchment in southern Greenland than in an adjacent nonglacierized watershed and provided an innovative statistical model of the phenomenon, but acknowledged that a two-sample analysis was insufficient to draw larger conclusions 51 regarding glacial variability dampening in the arctic or subarctic. Although we conclusively extend previous results to high latitudes and large ice fields, the analysis also demonstrates that such inverse relationships do not apply to all hydrograph components, a result requiring confirmation in other hydroclimatic regimes. Nevertheless, a comparison to the tentative analysis of monthly discharges in three Pacific Northwest rivers presented by Fountain and Tangborn (1985) suggests that the more detailed suite of relationships described here may indeed be geographically universal. The main limitations to the current study arise from data availability constraints. Consideration of a greater number of rivers sampling a broader range of glacial watershed covers could prove helpful. For example, a larger sample size might conceivably reveal a relationship between glacierization and the CV of qmax at a more restrictive confidence level; note that the observed Rs of 0.55 would be significant at a < 0.05 if the number of samples (rivers) was 11, rather than nine. In addition, while some of the study rivers effectively drain lobes of the massive St. Elias ice field, gauges are generally located well downstream from glacier termini and therefore sample relatively unglacierized watersheds. Consequently, failure to identify a concave-up relationship between glacial cover and CV for any metric does not imply that such a relationship does not occur here. 4.3.2 Timescales of Variability Dampening In a tentative but instructive complement to the coefficient of variation analysis, spectral methods are employed to evaluate the timescales at which glacial dampening of streamflow variability operate. Although most often used to identify periodicities, the power spectrum partitions the full temporal variability of a dataset into the portions of total variance explained by different timescales of variability, regardless of whether that variability is cyclic in nature. Means and least-squares linear trends were removed from daily runoff data for each of the nine stations over the full period of record; Bartlett windowing was applied; and spectra were found using the Lomb-Scargle method for unevenly-sampled (gappy) data (Press et al, 1992). Using simple left-endpoint numerical approximations to the integrals, sufficient for our purposes and in particular given the high sampling rate, and noting that the 52 frequency-domain sample interval, A/, is constant, these were then transformed to a logarithmic ratio (decibel) scale: = lOlog 10 f'+df ]>(/>/ f m a x \P(f)df ~101og 10 j*max / > 1 4 / P{f) = 101og 10 p(f) fmn (11) where the numerator and denominator give the spectral power over a particular band and over the entire frequency range sampled, respectively; fmin is the lowest frequency theoretically sampled, corresponding to the record length for each river (from 5.8x10'5 d*1 to 1.2X10"4 d"1); and fmax is the highest frequency directly sampled, corresponding to the Nyquist frequency (0.5 d'1). This normalization by total spectral power, a measure of the total variance of each dataset, facilitates graphical interbasin comparisons of periodogram form. Finally, 15-point Daniell estimates were produced. See Bras and Rodriguez-Iturbe (1993), von Storch and Zwiers (1999), and Fleming et al. (2002) for background and general explanations of the above processing steps. We neither generate confidence limits to the power spectra nor filter out periodicities, as we are interested in both periodic and aperiodic signal components. Results are illustrated in Figure 14; note that double-log (or equivalently, log-decibel) plots are standard when noise effects are a point of interest. Figure 14 demonstrates that all the study rivers possess inverse frequency-scaled (red) noise spectra, suggesting that glacial-nival spectral contrasts may be interpretable as differences in system memory, which might in turn be formalized using autoregressive or fractal models. However, relationships between serial correlation, red noise, and Markov-chain models may be more ambiguous and difficult (e.g., appendix two) than often recognized, and we limit ourselves here to a brief qualitative interpretation. With the notable exception of the one-year band if = 2.7xl0'3 d"1; obvious in arithmetic-space plots), which arises from the stronger freshets of glacierized watersheds (chapter three), glacial river spectra consistently exhibit lower relative power in comparison to nival rivers at timescales less than approximately two years if > about 1.3xl0"3 d'1). The result suggests that watershed 53 glacierization predominantly moderates short-term streamflow fluctuations. Glacial rivers may, therefore, be no less susceptible than nival streams to large-scale ocean-atmosphere circulation patterns (e.g., E l Nino-Southern Oscillation, E N S O ; Pacific Decadal Oscillation, PDO) , as such quasi-periodic phenomena have characteristic timescales well in excess o f two years. This is consistent with studies showing that glacial rivers respond to E N S O and P D O fluctuations in a manner at least as strong (but not the same) as nival watersheds (Neal et al, 2002; Lafreniere and Sharp, 2003), a subject considered in much deeper detail in the following chapter. Figure 14: Power spectra. Dark grey: nival streams; light grey: glacial streams. Two panels are identical in content, differences being only in display; glacial rivers are in the background in the left panel, and in the foreground in the right panel. 4.4 Conclusions Watershed glacial cover has been known for some time to modify streamflow variability. A t locations well downstream from the glacier terminus, where the percent o f basin glaciated is relatively small and human needs and biodiversity and biomass are greatest (e.g., Vannote et al, 1980; Roper and Scarnecchia, 2001), this modification takes the form o f a general 54 attenuation of year-to-year runoff fluctuations. Understanding of streamflow variability attenuation by glaciers nevertheless remains rudimentary. This study extends current knowledge by demonstrating that: (1) rank-based coefficients of variation (CV) in robust measures of freshet-related flow magnitude and timing (total yearly flow volume, annual mean daily discharge, and freshet magnitude and timing) scale inversely with watershed glacierization; (2) CVs in extreme events measures (magnitude and timing of annual peak daily discharge) do not exhibit an umabiguous relationship to degree of glacial cover; and (3) CVs in winter baseflow measures (annual minimum and median daily flows) also show no glacial moderating influences. The CV results are readily physically interpretable. Additionally, spectral analyses of daily discharge records from the nine study rivers tentatively suggest that watershed glacial cover most strongly attenuates streamflow variability occurring at timescales of about two years and less. 55 CHAPTER 5 CLIMATE MODE RESPONSES 5.1 Introduction Interannual streamflow variability in western North America has been reliably linked to large-scale ocean-atmosphere circulation phenomena, such as ENSO events and PDO regime shifts (e.g., Cayan and Peterson, 1989; Redmond and Koch, 1991; Kahya and Dracup, 1993; Dracup and Kahya, 1994; Lall and Mann, 1995; Mantua et al, 1997; Neal et al, 2002; Brito-Castillo et al, 2003; Woo and Thome, 2003). Similar or analogous associations between various large-scale circulation indices and streamflow variability have also been identified in other parts of the world (e.g., Trenberth, 1977; Jury, 1998; Franks and Kuczera, 2002; Cullen et al, 2002; Fowler et al, 2003). Aside from helping to iUuminate the nature of water resource fluctuations, these relationships hold some promise as a streamflow forecasting tool (e.g., Simpson et al, 1993; Piechota and Dracup, 1999; Hamlet and Lettenmaier, 1999a; Coulibaly et al, 2000; Hsieh et al, 2003), particularly given increasing ENSO prediction capabilities (e.g., Latif et al, 1998; Tang and Hsieh, 2002). Relative to pluvial systems, freshet-dominated runoff environments would seem to be especially amenable to circulation index-based forecasting approaches, as the snowpack effects of wintertime variability in large-scale circulation patterns generally lead the summer melt pulse by several months. However, Neal et al. (2002) found that glacier- and snowmelt-fed rivers in the Alaskan panhandle responded differently to the 1976-77 PDO regime shift, and Lafreniere and Sharp (2003) observed that a glacial river in Alberta was affected by the 1997-98 El Nino in a manner opposite to that in an adjacent nival basin. In addition, a larger-scale study by Woo and Thome (2003) hypothesized that the generally clear ENSO responses of western Canadian snowmelt-fed rivers may not extend to glacial systems. The possibility that glacial and nival streamflows consistently respond differently to interannual variability in circulation patterns thus presents fairly clear implications for a wide range of water resource issues, as well as the potential obfuscating impacts upon empirical climate studies which 56 utilize streamflow time series as a surrogate for meteorological point measurements. In addition, the streamflow impacts of such organized modes of large-scale climatic variability have been empirically and/or physically linked to sediment transport variability and fluvial geomorphic events (e.g., Smolders et al, 2002; Viles and Goudie, 2003; Aalto et al, 2003) and to fluctuations in lotic and lentic fish populations (Gunn, 2002; Smolders et al, 2002). Differential glacial-nival streamflow responses to circulation events may thus imply parallel contrasts in fluvial geomorphic and hydroecological responses, contributing an additional layer of complexity to baseline contrasts in the geomorphology, habitat, and ecology of glacial and nival systems (e.g., Ward, 1994; Smith et al, 2001; Dorava and Milner, 2001; see also chapter seven). This may be of particular practical significance to the management of severely-stressed salmon stocks in the variably-glacierized watersheds of the Pacific Northwest, although efforts to understand such impacts are greatly complicated by the need to disentangle freshwater and marine environmental influences upon anadromous fish (e.g., Mote et al, 2003). Finally, understanding the differential forcing effects of large-scale circulation patterns upon glacial and nival rivers may ultimately be crucial to a detailed assessment of the water resource impacts of potential global climate changes, insofar as such changes may be intimately tied to circulation behaviour (e.g., Trenberth and Hoar, 1996; Latif et al, 1998; Wallace and Thompson, 2002; Tsonis et al, 2003; see also Eltahir and Wang, 1999). Although intriguing, prior empirical delineations of whether and how glaciers modify streamflow responses to large-scale ocean-atmosphere circulation patterns remain tentative. In this study, we more thoroughly assess the potential for such effects by applying empirical orthogonal function (EOF), correlation, and composite analysis techniques to historical monthly discharge records from eight of our study basins, and comparing results between the two groups. The potential effects of two circulation patterns, ENSO and the Arctic Oscillation (AO), are considered. 57 5.2 Data Processing and Methodology 5.2.1 Data Eight of the nine study basins, four glacierized and four nival, were selected for analysis. As noted in chapter three, while the time series need not be continuous, EOF analysis (application of which is discussed further below) does require, as performed here, data for all locations at a given time, necessitating truncation of records to a common interval. As in chapter three, then, the Takhini River was omitted as the station was closed in 1986, which would have led to an unacceptably short common period of record. Monthly mean streamflow data over the period of record common to all selected stations span the interval 1975 to 1993 with gaps, yielding time series of uniform length N - 197. These are the same raw data used in the EOF analysis of chapter three, but different processing is used to emphasize interannual varability. Specifically, the time series were seasonally standardized (Figure 15) to zero mean and unit standard deviation (e.g., Salas, 1980) to facilitate interbasin comparisons and remove the seasonal cycle, which dominates the raw data: "q,^ qv,T„ = X"*»~X*» , with / = 12(v -1975) + r (12) ST,tJ where the mean and standard deviation correspond here to a given month (r = 1,12) over all years, v, on record; the i subscript simply denotes cumulative months from January 1975; and n = 1,8 rivers. The resulting standardized time series thus provide positive and negative anomalies relative to the 1975-1993 hydroclimatic normal, expressed as the number of standard deviations away from that zero mean, following removal of the average 1975-1993 annual cycle. While most of the hydroclimatic information content of these data - as well as, perhaps more importantly, the dominant water resource impacts - are associated with the summer flow period, rather than the groundwater-supported winter baseflow regime, it nevertheless seems prudent to retain the full time series for analysis. Given the purpose of the study, station data are not interpolated to a grid, a reasonably common step in climatological applications (see appendix one for a more detailed discussion). 58 Jan-75 Jan-81 Jan-87 Jan-93 Jan-75 Jan-81 Jan-87 Jan-93 Jan-75 Jan-81 Jan-87 Jan-93 Jan-75 Jan-81 Jan-87 Jan-93 Jan-75 Jan-81 Jan-87 Jan-93 Figure 15: Standardized monthly streamflow over common record, mean annual cycle removed. Left column: glacial rivers (decreasing glacierization from top). Right column: nival streams. Dimensions are number of standard deviations above or below 1975-1993 monthly normals for each river. 59 In part on the basis of the regional climatic characteristics sumrnarized in chapter two, and because the Southern Oscillation and Northern and Southern hemisphere annular modes may be the dominant patterns of variability in the global atmosphere (Wallace and Thompson, 2002), two modes of large-scale climatic variability were considered: ENSO and the Northern Annular Mode. ENSO is a coupled ocean-atmosphere interaction characterized by zonal (east-west) seesaws in equatorial Pacific sea level pressure (SLP) and sea surface temperature (SST); recent reviews and syntheses include Wallace et al. (1998) and Harrison and Larkin (1998). ENSO impacts propagate to distant regions by directly modifying general atmospheric circulation and, to a lesser degree, via coastal waveguide effects (Horel and Wallace, 1981; Mysak, 1986). We use the most common Southern Oscillation Index (SOI) definition, the difference in standardized SLP between Papeete, Tahiti and Darwin, Australia. A negative SOI value is indicative of El Niflo (warm phase) conditions, with slackened or reversed (normally, easterly, i.e., west-directed) trade winds and a pool of anomalously warm water in the central and eastern equatorial Pacific. La Nina (cold-phase) events are indicated by positive SOI values. El Nifios and La Ninas are not perfect opposites, and both their tropical Pacific SST and SLP signatures and (in particular) their global teleconnections show distinct asymmetries in both magnitude and location (e.g., Hoerling et al, 1997). Simple decomposition of the SOI into sinusoidal basis functions yields a series of spectral peaks with periods from 2 to 9 years; the most significant timescales of variability have been identified as about 2 and 4 years (Latif et al, 1998). This does not imply, however, that the SOI is a linear superposition of independent sinusoids with differing frequency as per conventional Fourier analysis, so it cannot be simply forecast like (for example) the tides. Rather, ENSO is a quasi-periodic process, which may be interpretable as a fundamentally perfect oscillator whose periodicity and, therefore, full predictability is marred by the effects of random noise, nonlinear interactions with the annual cycle, and/or decadal variations in the mean state (see Latif et al, 1998). Nevertheless, significant advances have been made in ENSO prediction. The Arctic Oscillation (AO), synonymous with the Northern Annular Mode (NAM), consists of a hemispheric-scale, but primarily high-latitude, meridional (north-south) seesaw in atmospheric mass between the Arctic and midlatitudes (Thompson and Wallace, 1998, 2000, 2001; Wallace and Thompson, 2002). An analogous oscillation occurs in the southern 60 hemisphere (the Southern Annular Mode, SAM), but seems of little direct relevance to the current study. The AO, of which the North Atlantic Oscillation (NAO) may be a regional expression, appears to be associated with fluctuations in the strength of the winter stratospheric polar night jet and may arise from such high-altitude circulation variability. Its index consists of the leading-mode unrotated principal component (hereafter referred to as LMUPC) time series of SLP poleward of 20°N. A positive AO index value is broadly indicative of negative and positive SLP anomalies in the Arctic and midlatitudes, respectively, and relatively strong 55°N (surface) westerlies; a negative index value indicates the converse (Thompson and Wallace, 1998; Wallace and Thompson, 2002). A preliminary spectral analysis of the 1900-2000 annual index yielded a series of very roughly equal peaks with periods ranging from about 2 to 7 years, with an additional, perhaps questionable, local maximum at about 80 years. The oscillation is thus not straightforwardly periodic; its overall predictability has yet to be assessed. AO and ENSO indices were downloaded as monthly time series from the University of Washington's Climate Data Archive (JISAO, 2002) and are illustrated, following standardization, in Figures 16 and 17. The ENSO and A O indices are not linearly correlated over the study period (R = 0.03 to 0.09 using monthly and annual data). The PDO (Mantua et al., 1997 and references therein) is also a point of interest; unfortunately, the available common period of streamflow record does not adequately span the 1976-77 PDO regime shift for meaningful analysis of local fluvial impacts. The Pacific-North America (PNA) pattern (Wallace and Gutzler, 1981) is physically associated, and statistically correlated, with ENSO (e.g., Shabbar et al, 1997; Woo and Thome, 2003), so in the interest of compactness it is not separately considered here. The surface meteorological and, to a far lesser degree, streamflow impacts of interannual variability in ENSO and A O indices have been tentatively defined within the study area by regression and composite mapping at very large spatial scales (Cayan and Peterson, 1989; Thompson and Wallace, 1998, 2000, 2001; Shabbar and Khandekar, 1996; Shabbar et al, 1997; note that ENSO mass balance impacts were considered by Walters and Meier, 1989, Hodge et al, 1998, and Bitz and Battisti, 1999 for a few British Columbian and Alaskan glaciers, but these are too far from the study area to be representative). Both Pacific and 61 4 N[a) 0] -2\ 2\ .4 4 — Jan-75 Jan-81 Jan-87 Jan-93 Figure 16: Standardized ENSO monthly index. N(o) is the number of standard deviations from the 1975-1993 climatic normal. Arctic circulation patterns appear to be locally expressed. However, the correlations are low and the study area seems to be at the local spatial limits of influence of such circulation phenomena, apparently forming something of a dead zone. Given the large mapping scales and the relative scarcity of meteorological and hydrometric station data within the study region, it is not entirely clear from previous studies whether this is in fact the case. If so, however, such weak and mixed responses are in one sense a disadvantage: glacial-nival distinctions in streamflow responses might be more clearly elucidated by considering a -4 4 - r -Jan-75 Jan-81 Jan-87 Jan-93 Figure 17: Standardized AO monthly index. Dimensions are as in Figure 16. 62 study area where the effects of circulation patterns are very strong, and perhaps in which one set of patterns (Arctic or Pacific, but not both) is clearly dominant so that interpretation is not complicated by interference between (apparently independent) processes. As we shall see, however, this subtlety and complexity ultimately provides us an opportunity to investigate a richer diversity of behaviour. 5.2.2 Methods EOF-Correlation Analysis Simple correlation mapping of the ENSO and A O indices upon the eight hydrometric stations was not found to yield clearly interpretable results. This is not at all surprising given the aforementioned apparent weakness and potential complexity of circulation pattern effects in the study area, in addition to the modest streamflow dataset length, and appears consistent with the preliminary station correlation analyses of Woo and Thome (2003). On the reasonable premise that the regionally coherent - though not necessarily regionally uniform - component of spatiotemporal flow variability represents our best bet for identifying and assessing the local water resource impacts of hemispheric- or planetary-scale circulation patterns, we shift the focus of our correlation efforts from individual stations to regional fluvial signals. This is accomplished using empirical orthogonal function (EOF) analysis. EOF analysis was performed in a fashion similar to that of chapter three. Procedural steps which were different from those of chapter three, or which require clarification or reiteration, are as follows. Analyses were performed separately using three datasets, in this case consisting of time series from: all eight rivers (hereafter denoted as the regional results); only the four glacier-fed rivers (glacial); and only the four snowmelt-fed rivers (nival). As a consequence of the processing steps described previously, the data correlation (rather than covariance) matrix was used, which is appropriate for our purposes; note also that removal of the mean annual cycle emphasizes entirely different aspects of the data than were considered in the EOF analysis of chapter three. As discussed previously (and in more detail in appendix one), rotation procedures present a litany of potential disadvantages and 63 are contraindicated as a general processing tool unless the eigenvalues are degenerate; this was again not found to be the case here for the leading modes. Finally, the same four methods employed in chapter three for assessing modal significance were again applied. These are: restriction of our attention to the leading mode; qualitative consideration of the proportion of variance explained by the leading mode; a quantitative modal separation test on the basis of estimated EOF sampling error; and a quantitative nonparametric Monte Carlo randomization test on the eigenspectrum. The latter was again performed without explicit adjustments for serial correlation, for essentially the same reasons specified in chapter three. The three resulting monthly LMUPC time series were then compared to each other, and against the ENSO and A O indices. ENSO and A O events and their surface meteorological effects are by no means restricted to winter months. However, the A O appears predominantly a wintertime phenomenon; the canonical ENSO event is typically strongest from about early fall of the onset year to late spring of the following year, peaking in the intervening winter; and ENSO temperature and precipitation impacts are locally strongest during winter (e.g., Shabbar and Khandekar, 1996; Harrison and Larkin, 1998; Thompson and Wallace, 2000). Inspection of monthly SOI and A O index time series suggests that their mean state over the period December through March provides a generally reliable indication of whether or not an ENSO or A O event is occurring, and appears slightly more sensitive in this regard than consideration of a full calendar or hydrologic year of monthly data. It is therefore reasonable to replace the monthly indices with their annual DJFM averages (assigned to the JFM year in this study; see subsequent composite analysis discussion for further details regarding year assignment), and to compare these with annually-averaged monthly L M U P C time series. This general type of approach is not uncommon (e.g., Horel and Wallace, 1981; Shabbar et al, 1997), and is of particular value when considering potential streamflow consequences. Relationships between the local meteorological effects of circulation index variability and their streamflow impacts can follow a complex system of lags. For example, freshet flow depends (at a minimum) upon both winter snowpack (accumulated over lags of about 1 to 5 months prior the start of the freshet and roughly 6 to 10 months prior to its end), and spring and summer temperatures (perhaps a 0 to 1 month lag). The advantage of using annual time series, then, is that it avoids the necessity of explicitly considering such lagged cross-correlations, rendering the analysis simpler and 64 perhaps more reliable. An annual mean LMUPC data point was calculated only for years in which a full 12 months of data were available; the resulting time series are thus of length NA = 10, spanning the interval 1979-1993 inclusive, with gaps in 1980 and 1983-1986. Two additional procedural notes should be made. First, this annual averaging of monthly LMUPC time series facilitates a substantial reduction in EOF sampling error (see equation (5)) relative to that obtained through the application of EOF analysis directly to annual average discharge data, which of course would immediately yield annual LMUPC time series; the issue is particularly important for such relatively short records. Second, while annual time series by definition do not contain the seasonal cycle, and it might thus seem that the LMUPC time series generated using the non-deseasonalized time series of chapter three could be used as the basis of the annual averaging process, the foregoing procedure is in fact far preferable. The mission, as it were, of a leading EOF mode is to capture as much of the total dataset variance as possible (subject to the constraint that all modes are orthogonal and uncorrelated), and the vast majority of the spatiotemoral variance in non-deseasonalized streamflow arises from the annual cycle. The LMUPC time series of chapter three, therefore, have in essence been fitted to the seasonal cycle, and interannual variability may not be reliably represented. Similarly, the calculated proportions of variance explained, the Monte Carlo randomization test, and the test for eigenvalue separation would be meaningless from an interannual perspective. Such problems are not encountered using the procedure outlined above. Pearson product-moment correlation coefficients between the various time series were calculated, and two-sided tests for significance were performed using a Monte Carlo randomization approach, closely analogous to that applied above and in chapter three for assessing EOF modal significance. The method, which avoids distributional assumptions, was first suggested in a correlation context by Pitman (1937); see also Edgington (1995) and Manly (1997). As implemented here, the absolute value of the correlation coefficient is taken to be the test statistic. R0bs is the value obtained from the original time series. For a given pairwise comparison, the sequence of observed data values is randomly reordered, individually for each time series, and the absolute value of the correlation coefficient between the two reshuffled records is calculated. This is independently performed NMC 65 times, resulting in NMC correlation coefficients; their empirical CDF provides the distribution of the test statistic, known as the randomization distribution, and its (I-a) • 100th percentile value constitutes the critical value, Rcrih against which the observed test statistic is compared: if R0bs > Rent, then the null hypothesis of no association is rejected against the alternative hypothesis of linear correlation at confidence level, a. While this form may be suggestive of a one-tailed test, the use of absolute values in fact yields a two-tailed test. One could equally well employ a somewhat more conventional, but more cumbersome, approach to two-tailed significance testing by abstaining from absolute values and comparing Robs, which would then simply be the observed correlation coefficient exactly as calculated using equation (9), against two Rcru values corresponding to the (1 -a/2) • 100th and (a/2) • 100th percentiles of the randomization distribution; note that it is unclear whether symmetry and centredness of the empmcalJy-determined randomization distribution is guaranteed (unlike z-scores, for example), so that \Rcru(ai might not equal \RCrtt(\-ai-, and it would thus be prudent to examine both Rcrtt values if such an approach was utilized. In contrast to the foregoing EOF application of randomization modelling, potential serial dependence concerns cannot be swept aside so readily with respect to the Pearson correlation analysis. Nonetheless, simple tests on the significance of the lag-1 serial correlation coefficients of the various annual time series using the method of Anderson (1942) as interpreted by Salas et al. (1980) did not indicate strong serial dependence. Moreover, integral time scales: were generally found to be roughly equal to or less than unity, where / is the lag, p(l) is the autocorrelation function, and A/ is the nominal sampling interval of one year. The empirical autocorrelation functions underlying both the lag-1 correlation significance tests and the decorrelation time estimates, and calculated blindly using all subsequent data points in the time series, are ad hoc due to the uneven sampling interval arising from missing data. By the same token, however, data gaps also directly reduce whatever serial dependence may or (13) 66 may not have originally been present by breaking the Markov chain, so to speak; indeed, although the approach may require closer scrutiny (see appendix two), time series decimation has been proposed as a method of eliminating the effects of autocorrelation upon statistical tests (e.g., von Storch, 1995). Nevertheless, recompletion of these tests in a modified manner which makes use only of truly time-adjacent data points still did not strongly suggest serial correlation. This is consistent with the general comments of Lettenmaier et al. (1994) and the discussion of appendix two, which suggest that autoregressive noise is typically not strongly present in annually (but not more finely) discretized river discharge data. In any event, it is clearly appropriate in the present case not to adjust the data or correlation test procedure for serial dependence. Composite Analysis of Annual Hydrographs Composite analysis (e.g., von Storch and Zwiers, 1999) was performed as a methodologically completely independent check upon the foregoing results, and to examine streamflow seasonality impacts in a robust and straightforward manner. While composite analysis has more lenient data requirements, to ensure full consistency and comparability with the results of the EOF-correlation analysis, the period of common full annual data (1979-1993 with gaps) was again used. A positive or negative event year was defined as one during which the annual DJFM ENSO or A O index was greater or less than half a standard deviation (e.g., Redmond and Koch, 1991) from its 1979-1993 mean. Positive and negative A O event years thus identified are (1989, 1990, 1992, 1993) and (1979, 1981, 1987, 1988), respectively. Positive and negative ENSO years are (1982, 1989) and (1987, 1992), respectively. A O event years apparently have not previously been explicitly published. Our El Nino and La Nina years correspond well to those of, for example, Shabbar et al. (1997) over the same time period once we account for the fact that, in the interest of convenience for streamflow analysis (see previous section), we assign the DJFM SOI average to the JFM year. We thus consider here a hydrologic ENSO event year, which is the year following ENSO event onset; it is the latter to which the title of ENSO event year is more conventionally assigned. In terms of the usual nomenclature for climatological composite analysis (e.g., Shabbar et al, 1997), this is exactly equivalent to performing composite analysis of hydrologic data from the [+1] year, using the [0] year SOI data to 67 identify ENSO events. Note that due to data gaps, we miss the substantial 1982-83 El Nino (following the more conventional labelling of ENSO event years). The composite annual hydrograph of the rf 1 river for positive AO events was then calculated as: n where r = 1,12 is the month, Q is the number of positive A O years over the period of record, and qyxn is the corresponding seasonally-standardized monthly streamflow. The dimensions of the composite are thus standard deviations from the 1979-1993 hydroclimatic normal. Expressions for negative A O and positive and negative ENSO years are analogous. For the four glacial rivers, and separately for each month, r = 1,12, the median value o f + A O C T t t l =1 ,4 is compared to that of 'AOCTi^u using a two-sided nonparametric Mann-Whitney test (see chapter four for a description of this method). The emphasis, therefore, is upon contrasts between opposite climate states rather than deviations from normal conditions. This procedure was also performed for nival streams, and the entire process was then repeated for ENSO events. 5.3 Results +AO 5.3.1 EOF-Correlation Analysis The three resulting LMUPC time series are illustrated in Figure 18, and the leading-mode eigenvectors and corresponding percentages of variance explained are provided in Table 6. In all cases, the leading mode was well-separated from higher modes and significant at a = 0.01, and explained a substantial proportion of the total spatiotemporal variance, providing at least a tentative indication that the EOF results are reliable. Note that individual higher modes were responsible for much lower percentages of total variance. Overall, the leading-mode eigenvector entries were about equal with the exceptions of the Kluane and Wheaton 68 N (a) 0 Jan-75 Jan-81 Jan-87 Jan-93 Jan-75 Jan-81 Jan-87 Jan-93 Jan-75 Jan-81 Jan-87 Jan-93 Figure 18: Leading-mode unrotatedprincipal component time series, (a) regional result (generated by EOF analysis of data from all eight rivers); (b) glacial result (using data only from four glacial rivers) ;(c) nival result (data from four nival rivers). Dimensions are number of standard deviations from 1975-1993 hydroclimatic normal. 69 Rivers. The Kluane River eigenvector entry for both EOF analyses in which it appears is of identical sign, but smaller magnitude, relative to those corresponding to the other watersheds. This indicates that Kluane River discharges are positively but comparatively weakly correlated to the corresponding LMUPC time series, implying that its response to regional-scale climatic driving forces may be attenuated. Such behaviour might be caused by Kluane Lake, which in both area and depth dwarfs the lakes of the other catchments considered, and which may thus possess sufficient storage volume to modify interannual flow variability. The effect of extensive lake cover upon riverine responses to large-scale circulation patterns may itself be a useful direction for future studies (see Arnell, 2002). There is no immediately apparent cause for the somewhat anomalous behaviour of the Wheaton River eigenvector entries. Nevertheless, to a first approximation, the leading mode for each case yielded a roughly uniform spatial pattern, so it is the corresponding LMUPC time series that contain most of the (potentially) interesting behaviour, and the temporal correlation-based approach to assessing ENSO and A O effects discussed above thus seems appropriate. The resulting correlation matrix is presented in Table 7. Our interpretation begins with the glacial and nival LMUPC time series. The correlation results provide two key insights. First, the glacial and nival LMUPC time series do not mutually correlate, indicating a decoupling of interannual flow variability between glacial and nival rivers. This is broadly consistent with previous empirical studies demonstrating that glacial and nival rivers respond to climate variability in moderately to radically different ways, depending on the timescale of investigation and local glacioclimatic conditions {Henshaw, 1933; Meier, 1969; Krimmel and Tangborn, 1974; Zuming, 1982; Fountain and Tangborn, 1985; Ferguson, 1985; Braithwaite and Olesen, 1988; Chen and Ohmura, 1990a; Moore, 1992; Neal et al, 2002; Lafreniere and Sharp, 2003; other chapters, this thesis). However, such apparently strong, systematic independence between the annual flow anomalies of glacial and nival rivers does not appear to have previously been explicitly demonstrated. Second, the results suggest selective teleconnectivity of glacial and nival catchments to Arctic and Pacific circulation patterns, respectively, where we adopt Redmond and Koch's (1991) broader definition of teleconnections as climate linkages over great distances. The 70 glacial LMUPC time series correlates strongly with the A O index (a = 0.02), but not at all with the SOI; in contrast, the nival LMUPC time series is moderately correlated with the SOI (a = 0.10), but not at all with the A O index. Prior work has focussed on the magnitude of streamflow responses to a single, strong, and apparently locally dominant circulation pattern: the PDO in the Alaskan panhandle (Neal et al., 2002), and ENSO in Alberta (Lafreniere and Sharp, 2003). The selective glacial-nival teleconnectivity demonstrated here in a region of somewhat weaker, mixed climate signals has not been previously identified. Table 6: Leading-mode eigenvectors and percentages of variance explained. These are left unitless, with dimensionality ascribed instead to the corresponding LMUPC time series (Figure 18) to maintain consistency with other figures; the converse would be equally technically valid.. run % glacierized basins nival basins var W R A H K R K L ARBR WRA WRC M C R W DRHJ B C M glacial 45 0.56 0.02 0.62 0.54 - - - -nival 48 - - - - -0.12 0.60 0.54 0.58 regional 34 0.37 0.11 0.42 0.40 0.10 0.43 0.46 0.33 Table 7: (Symmetric) correlation matrix between LMUPC time series and climate indices. Annually averaged regional, glacial, and nival LMUPC time series and mean DJFM ENSO and AO indices were compared. Results are expressed as significance levels at which null hypothesis of no linear relationship can be rejected, found using a nonparametric randomization procedure; corresponding correlation coefficients also shown in parentheses. Correlations along the diagonal ("D ") are perfect by definition. regional glacial nival ENSO AO regional D 0.02 (+0.74) 0.01 (+0.82) > 0.25 (-0.29) 0.10 (+0.60) glacial - D > 0.25 (+0.24) > 0.25 (+0.15) 0.02 (+0.76) nival - - D 0.10 (-0.57) > 0.25 (+0.28) ENSO - - - D > 0.25 (+0.03) A O - - - - D 71 The hydrology of snow and glaciers is complex and incompletely understood (see reviews by Fountain and Walder, 1998 and Singh and Singh, 2001), and thus so too are the detailed hydrologic responses of nival and glacial watersheds to climatic forcing (e.g., Moore and Demuth, 2001; Lafreniere and Sharp, 2003). At a coarser level, however, the type of decoupling of glacial and nival streamflows found here may be interpreted in terms of the presence or absence of a large ice reservoir potentially available for melting. Consider for the sake of argument a year in which summer temperatures are anomalously high, with all else held constant. As its hydrologic inputs are fixed in this case, a nival watershed can respond only with a streamflow decrease due to increased evapotranspiration losses; these are balanced or exceeded in a glacierized watershed, however, by increased melt generation from the glacial reservoir. The net response to a climatic signal can thus potentially be diametrically opposite between adjacent and otherwise identical glacial and nival rivers, accounting for both decoupling and differential responses to a given large-scale circulation pattern. For selective teleconnectivity of glacial and nival rivers to different circulation patterns, we have the additional requirement that each pattern locally causes different surface meteorological effects. This constraint does not seem unreasonable for ENSO and the AO, which are physically distinct and statistically linearly uncorrelated (but see Wu and Hsieh, 2004). The local meteorological impacts of ENSO and the A O are addressed in somewhat greater detail in a subsequent section. The regional EOF-correlation results are enigmatic. The leading EOF mode based on all eight rivers exhibits less convincing correlations with the circulation indices than do its glacial or nival counterparts, and explains a substantially lower proportion of total dataset variance than either the glacial or nival leading modes. Moreover, the regional LMUPC time series turns out to be nothing more than a simple linear combination of the glacial and nival LMUPC time series: L M U P C regional = B [LMUPCglacial)* (LMUPC nival) (15) which correlates with the observed regional time series at R - 0.99 for (B, 0.25 (-0.14) 0.25 (+0.41) glacial - D > 0.25 (-0.06) > 0.25 (+0.23) 0.10 (+0.55) nival - - D 0.25 (-0.44) > 0.25 (0.00) ENSO - - - D > 0.25 (+0.03) A O - - - - D 74 5.3.2 Hydrologic Composite Analysis Figure 19 presents the resulting suite of composite annual hydrographs, with a shaded bar denoting months for which median values are different between positive and negative index states at a ~ 0.05. Also illustrated, for purposes of comparison, are composites developed for the stacked glacial and nival time series; these were not incorporated into the Mann-Whitney test as they are obviously not independent samples. Composite annual hydrographs derived from the glacial and nival LMUPC time series (not shown) were found to be virtually identical to the stack-based composites. Streamflow responses to the Arctic Oscillation are clear. During a positive A O event year, glacial streamflows are much greater April through June relative to a negative A O year; the overall effect at the annual level, then, is a streamflow increase. This is fully consistent with the strong positive correlation we found between the A O index and the glacial LMUPC time series. The nival streamflow composites show significant increases and decreases in spring through early summer, and late summer through early fall, respectively, in comparison to a negative A O year. The composite analysis thus demonstrates that nival systems are not, in fact, immune to A O effects as implied by the correlation analyses. Results using the two methodologies are nevertheless fully consistent: the nival monthly streamflow changes associated with A O events largely mutually cancel at the annual level, leading to no significant correlation between annual mean nival discharge and the A O index. From a pragmatic, net water resource perspective, the selective teleconnectivity of glacial and nival rivers to the Arctic Oscillation, identified on the basis of EOF-correlation analysis, is indeed robust. ENSO responses are not as strong or coherent. Nevertheless, relative to negative ENSO events, La Nifias are locally associated with low glacial flows in March and high glacial flows in May. These, however, mutually cancel at the annual level. Similarly, then, glacial rivers are indeed somewhat susceptible to ENSO effects but not in terms of annual net water resources, and the composite results are thus fully consistent with the absence of a significant correlation between the SOI and the annual glacial LMUPC time series. During 75 Figure 19: Composite annual hydrographs of standardized hydrometric station data (following two pages). Black: positive event years; grey: negative event years. Units are number of standard deviations from 1979-1993 monthly normals for each river. Shaded bars denote months for which the median value of the four average annual hydrographs (one for each river) is significantly different between positive and negative event years. Bold lines are composites constructed from glacial and nival stacks; these were not included in the statistical comparison. Different symbols are used to identify watersheds. Glacial rivers (a and c): Wann = triangle, Alsek = square, Kluane = circle, White = diamond; nival rivers (b andd): Big Creek = triangle, Dezadeash = square, M'Clintock = circle, Wheaton = diamond. 76 -1.5 I -i 1 1 1 1 1 1 1 1 1 1 1 2 3 4 5 6 7 8 9 10 11 12 month 1 2 3 4 5 6 7 8 9 10 11 12 month a La Nina year, the significant changes in nival discharges consist of low flows in both July and August, leading to a net annual decrease. Again, this is consistent with the EOF-correlation results. In general, then, the glacial and nival composites show coherent responses to ENSO and the A O consistent with the EOF-correlation and stack-correlation findings. However, considerable variability exists between rivers. This may represent the effects of locally relatively weak circulation pattern effects, and interbasin variability in watershed characteristics such as lake cover and hydrogeology (see chapter two). The composites of the Kluane and Wheaton Rivers deviate most from the other rivers in the glacial and nival groups, respectively, which is consistent with the leading-mode eigenvector entries. Interestingly, Kluane deviations from the glacial river norm occur primarily from about August forward or over January-March, which may be consistent with the previously-proposed lake mechanism for Kluane deviations, insofar as the storage capacity of a large lake may alter both the annual recession curve and winter baseflow. Furthermore, Wheaton deviations from the nival river norm occur primarily during winter baseflow, which in this case may suggest a hydrogeologic control. Overall, however, the responses of all the rivers, including the Kluane and Wheaton, are quite similar to their (glacial or nival) peers during the summer freshet which, as noted previously, is the most important streamflow season in terms of both pragmatic water resource impacts and hydroclimatic information content. Note that the foregoing interpretation focused on those months for which differences in the composites were statistically significant, and that these in turn also largely occur during the freshet. Along similar lines, the contrasts between the circulation pattern responses of all the glacierized rivers, versus all the nonglacierized rivers, are stark. This suggests (together with other considerations; see chapter three) that the presence or absence of glacierization is the dominant watershed-level control, rather than the annual snowpack effects of elevation. The latter, though (generally) greater in glacierized catchments, varies in a more continuous fashion between watersheds and thus seems less likely to yield the observed binary distinction. 79 5.3.3 Composite Analysis of Surface Meteorological Effects Given the physical complexity of the relationships between large-scale circulation patterns, local meteorological effects, and streamflow impacts, it might be unduly optimistic to hope that the hydrologic patterns identified above could be readily explained in terms of a few basic meteorological parameters, such as temperature and precipitation (see, for example, Meier, 1969). Nevertheless, given that ENSO and the A O obviously can only impact streamflow via surface meteorological pathways, it is reasonable to expect that local meteorological anomalies associated with ENSO and A O events are, at a minimum, not inconsistent with the hydrological results. To this end, composite analyses similar to those described above were performed using monthly mean temperature and total precipitation data. Time series were considered for stations of the Historical Canadian Climate Database (HCCD; see Vincent and Gullet, 1999 and Mekis and Hogg, 1999) present within the study area, and possessing adequate temporal coverage over the common period of streamflow record to ensure comparability with the hydrologic investigation. The stations used were Burwash Landing, Haines Junction, and Whitehorse Airport for temperature, and Burwash Landing, Atlin, and Whitehorse Airport for precipitation. Due to the smaller number of stations (three) relative to river discharge (four for each composite analysis), a two-sided Mann-Whitney test at a = 0.05 cannot be performed (e.g., Daniel, 1990); instead, shading on Figure 20 denotes significance using a two-sided test at a = 0.10, which is equivalent to the appropriate (on a month-by-month basis) one-sided test at a = 0.05. Although significant temperature effects are seen for all months but September, positive A O years are most strongly associated with a large positive spring-summer temperature anomaly peaking in June, inducing a substantial increase in glacial melt production and streamflow. On the other hand, evapotranspiration (ET) losses can be substantial due to fairly high typical temperatures and very long day lengths in summer, and the presence of substantial wetlands. In addition, moderate increases in winter snowpack associated with positive A O events, considered alone, could delay the onset of glacial melt production, reducing streamflow due to the far greater albedo of snow (about 0.6 to 0.8) relative to ice (0.3 to 0.4) (see Meier, 1969) (note that glacial and nival systems thus respond oppositely to anomalies not only in summer temperature but also in snowpack, all else held constant). Moreover, 80 Figure 20: Composite annual regimes of mean temperature and total precipitation. Black: positive event years; grey: negative event years. Units are number of standard deviations from 1979-1993 monthly normals for each of the three met stations. Shaded bars denote months for which the median value of the three average annual thermographs and hyetographs (the terms are used loosely here) is significantly different between positive and negative event years. Bold lines are composites constructed from stacked temperature and precipitation time series; these were not included in the statistical comparison. Different symbols are used to identify met stations. Temperature (a and c): Burwash Landing = triangle, Whitehorse - square, Haines Junction = circle; precipitation (b and d): Burwash Landing = triangle, Whitehorse = square, Atlin = circle. 81 -1.5 -i 1 1 1 1 1 1 1 1 1 < 1 1 2 3 4 5 6 7 8 9 10 11 12 month summer precipitation decreases slightly. Judging from the streamflow composites, these three effects are more than offset by the temperature-induced increase in glacial melt production. This seems reasonable, given that mass balances in the region are strongly negative under historical conditions of progressive climatic warming, these glaciers are at least semi-continental, and temperature may in general be the strongest control over glacial melting, which each imply a strongly positive streamflow sensitivity to summer ablation rates (see Hodge et al, 1998; Bitz and Battisti, 1999; Arendt et al, 2002; Alley, 2003; chapter three; chapter six). It is unclear from the meteorological composites, however, why the positive glacial streamflow anomaly terminates by about July. In nival systems, which do not enjoy the benefit of glacial reservoirs, increased ET losses and minor summer rain decreases appear to be balanced by the greater snowpack associated with positive A O events, so the net hydrologic effect consists solely of the forward shift in freshet timing observed in Figure 19(b). ENSO events do not appear to be locally associated with significant snowpack effects, as no precipitation anomalies occur in January or February, and the large observed winter temperature anomalies are largely hydrologically irrelevant from a snow accumulation perspective here due to very cold baseline winter temperatures. However, La Nina years are associated with a positive summer temperature anomaly roughly similar to that experienced in positive A O years, although the summer temperature difference between positive and negative ENSO events is much smaller, and of somewhat different form, than that between positive and negative AO events. As with positive A O events, then, increases in glacial melt production and net glacial streamflow are seen relative to negative index years, but are far less dramatic, more erratic, and statistically significant only for May. Lower observed March glacial river discharge might be due to anomalously low temperatures at this time, delaying by a few weeks the earliest stages of the melt freshet and, more broadly, the onset of spring. Note that while local clirnatological winter (Wahl, 1987) and river ice cover (from Environment Canada, 1999) typically last until April, and the freshet does not generally start in earnest until April or May, the net radiation balance at Whitehorse becomes positive by March and freeze-thaw events may not be uncommon by this time (see Wahl, 1987). The March glaciofluvial regime is therefore no longer completely dominated by aquifer releases as in wintertime, and contemporaneous meteorological anomalies may 84 have some streamflow impact. Note that nival streams also qualitatively appear to exhibit a similar March effect, which is to be expected given the mechanism proposed. In nival rivers, La Nifia-associated increased ET losses do not appear to be balanced by precipitation increases, unlike the case for positive A O events; although spring-through-fall precipitation anomalies occur, these are of opposite sign for different months and seem to mutually cancel. The net hydrological effect of a La Nifia in nival streams is thus the streamflow decrease seen in the hydrograph composites. Clear, robust delineation of linkages between circulation patterns, local surface meteorological effects, and streamflow impacts would obviously be gratifying. However, the ambiguity and nonuniqueness of the foregoing qualitative interpretation of the meteorological middleman's role underscores a broader, and probably more important, point. Streamflow, temperature, and (to a somewhat lesser degree) precipitation are each seen to show distinctive responses to the circulation indices considered. Nonetheless, it would be very difficult to reliably predict the water resource response to a suite of meteorological changes observed to be associated with the phase of a particular circulation index or, by the same token, to infer from an observed streamflow response the corresponding meteorological changes. This may in part arise from data limitations: homogeneous humidity, windspeed, and cloudiness time series, for example, might provide some of the missing pieces of the above puzzle but unfortunately are not included in the HCCD. However, the filtering effects of watershed-specific processes and characteristics -in particular, fundamental differences in the water resource responses of glacial and nival systems - greatly heighten the difficulty. The results thus suggest that, realistically, in variably-glacierized regions the surface hydrological and meteorological responses to variability in large-scale climatic driving forces might best be assessed directly and independently (though preferably in tandem), rather than attempting to infer one from the other. 5.4 Discussion The EOF-correlation, stack-correlation, and composite results all take us in the same direction. Although the record lengths and number of hydrometric stations available are 85 modest, the statistical methods employed account for the limiting effects of sample size upon the reliability of inferences drawn from that data. The statistically significant relationships identified here between circulation patterns and streamflow, and the influences of watershed glacierization upon those relationships, thus appear robust. However, important details vary from place to place. Conceptually, it would seem the differential glacial-nival streamflow impacts of large-scale circulation variability would be location-dependent, because (at a minimum) the local meteorological effects of such variability are also location-dependent. This appears to be borne out in practise. Our records are too short to consider the PDO impacts evaluated by Neal et al. (2002) and no one has previously considered A O effects in the current context, but we can compare our ENSO results to those (found using a very different approach) of Lafreniere and Sharp (2003). In their analysis of two adjacent watersheds in Alberta, they observed a substantial streamflow increase in response to the 1997-98 El Nifio in a glacier-fed river, but a decrease in the nival system; this was attributed to ENSO-related low 1997-98 winter snowfall and high 1998 spring and summer temperatures. However, El Nifio impacts upon temperature and, in particular, precipitation appear much greater in Alberta than in southwest Yukon; moreover, the two regions respond differently to the asymmetries between El Nino and La Nifia events (see Shabbar and Khandehar, 1996; Shabbar et al., 1997). There may also be hydrologically important glaciological differences between the regions (see chapter six), and while both regions possess freshet-dominated fluvial regimes, their baseline climatic conditions are certainly different as well (i.e., subarctic semi-continental versus temperate continental). Such considerations may be responsible for the contrast between our results and those of Lafreniere and Sharp (2003): we found little net water resource impact in glacial rivers from ENSO events, and a modest El Nino-associated streamflow increase in nival rivers relative to La Nifia years. It also bears mentioning, however, that whereas we have sampled the effects of several El Nifio and La Nina events upon a fair number of basins, so that our results may reasonably be considered representative of regional ENSO impacts, Lafreniere and Sharp (2003) captured the effects of a single warm event and no cold-phase episodes upon only two rivers, so their conclusions may thus be less generally representative of local responses. 86 Nevertheless, while the local details no doubt do vary considerably, the work presented here and by Lafreniere and Sharp (2003) does converge to a common fundamental result: glacial and nival rivers respond to ENSO in very different ways. Together with our Arctic Oscillation results and the differential PDO responses identified by Neal et al. (2002) in coastal southeastern Alaska, it can be said with some confidence that such contrasts extend to the fluvial impacts of large-scale atmosphere-ocean circulation patterns in general. 5.5 Conclusions The potential for glacier- and snowmelt-fed rivers to respond differently to two large-scale ocean-atmosphere circulation patterns, El Nifio-Southern Oscillation and the Arctic Oscillation, was empirically investigated. Regionally coherent glacial and nival hydroclimatic signals were first extracted from the station data using both empirical orthogonal function analysis and a stacking procedure. The resulting annually-averaged glacial and nival time series were found to be mutually uncorrelated, suggesting a decoupling of interannual streamflow variability between the two fluvial regimes. In addition, the glacial signal was significantly positively correlated to the A O index and not at all to the SOI, whereas the nival signal was modestly negatively correlated to the SOI and not at all to the A O index, suggesting selective teleconnectivity. Composite analyses of annual hydrographs performed using discharge data from the individual hydrometric stations were fully consistent with the EOF-correlation and stack-correlation results, and provided details of ENSO- and AO-related changes in hydrologic seasonality. Results of parallel composite analyses of the limited homogeneous temperature and precipitation datasets available for the study area were generally consistent with the streamflow interpretations. Both glacial-nival hydroclimatic decoupling and selective teleconnectivity are fresh discoveries; however, both do appear strongly consistent, at least at a conceptual level, with prior studies. 87 CHAPTER 6 IMPACTS OF CLIMATIC CHANGE 6.1 Introduction Perhaps the most widely-suggested, yet least clearly-demonstrated, effect of watershed glacial cover upon fluvial responses to climatic variability consists of streamflow change in glacier-fed rivers due to progressive glacial ice mass loss arising from climate warming (e.g., Chen and Ohmura, 1990b; Konovalov and Shchetinnicov, 1994; Hodgkins, 1997; Pelto and Riedel, 2001). The negative mass balance trends of mountain glaciers under post-mid-nineteenth century warming conditions, although not universal, are nonetheless global in scope and may be accelerating (e.g., Dyurgerov and Meier, 2000). In the worst-case scenario of G C M simulations which predict the demise of a substantial proportion of the world's alpine glaciers by ca. 2100, this might cause streamflow increases during an initial melting phase followed by declines relative to present conditions as glacial reservoirs are progressively exhausted (IPCC, 2001a). Such considerations recently led the United Nations to identify the potential water resource impacts of glacial responses to climate change as a major scientific and policy issue (Showstack, 2002). However, in spite of the gravity of the matter, the amount of overall research it has generated, and the relevance, accessibility, and spatial and temporal coverage of data from existing general-purpose streamgauge networks, little or no effort has been made to specifically assess whether there is empirical evidence for differential streamflow trends between glacierized and glacier-free basins under historical warming conditions. This knowledge gap seems conspicuous, particularly since previous research has a number of limitations. Mass balance studies are innumerable and often presume to have strong water resource implications, insofar as a negative mass balance may indicate an increase in glacial melt production and runoff. However, although mass balance is of obvious hydrological import, it is not a comprehensive and unique descriptor of the potential streamflow effects of glacioclimatic variability; as just one example, melt generation and streamflow can rise due to an increase in the amplitude of the annual mass balance cycle, 88 with no net change in ice mass. In addition, virtually all previous work has focussed tightly upon mass balance-proglacial runoff relationships, often leading to an emphasis upon low-order glacial headwater streams and exclusion from the analysis of nonglacial streamflow influences, such as catchment hydrogeology, direct precipitation-runoff relationships, and vegetative, lake, and wetland cover. In contrast, greatest human needs, biodiversity, and biomass occur far downstream from the glacier terminus (Vannote et al, 1980; Milner and Petts, 1994; Roper and Scarnecchia, 2001), so the wider practical effects of glaciohydrological change might be obscured or mitigated by nonglacial controls and inputs. Results from relatively uncommon quantitative models of glacial runoff responses to 2xCC»2 GCM-predicted climate scenarios (Davidovich and Ananicheva, 1996; Singh and Kumar, 1997) are only as reliable as the climate model output upon which they are based (e.g., Hamlet and Lettenmaier, 1999b; see also Coe, 2000). Actual streamflow data against which theories regarding glaciological-hydrological interactions might be tested are often not incorporated into the analyses. Similarly rare studies based on joint glaciological and hydrological data may provide the most convincing evidence, but these are subject to some of the above restrictions and have also yielded contradictory results: although all their study glaciers exhibit negative mass balance trends, Collins (1989) and Braun and Escher-Vetter (1996) found increasing proglacial discharge trends in continental Europe, whereas Moore and Demuth (2001) and Demuth et al. (2002) found negative glacial river runoff trends in temperate regions of maritime and continental North America. Perhaps most importantly, few if any studies have incorporated a significant control population of nonglacierized catcliments from the same hydroclimatologic region. A discovery of general regional trends in glacial streamflow, for example, could lose much of its significance from a glacioclimatic or, in particular, comparative glacial-nival perspective if it turns out that adjacent nival streams show similar trends. Clearly, no single study suffers all these limitations, but each possesses at least one. This suggests that an alternative approach to the problem might prove instructive. The salient question, particularly but not solely from a water resource perspective, is simply this: do rivers with glaciers in their watersheds respond to climatic warming trends in a manner which is demonstrably different than those without? A sensible approach to addressing this issue, inspired by previous nonparametric statistical trend analyses of 89 historical streamflow data, is to empirically compare observed long-term discharge trends in adjacent rivers, holding factors other than catchment glacierization as constant as practicable between watersheds. We present here such a controlled natural experiment, portions of which are also discussed by Fleming and Clarke (2003) or Fleming (in review), conducted using several decades of total annual flow volume data from our five glacierized and four nival basins, with catchment areas ranging from hundreds to thousands of square kilometres. These are located in a region which is undergoing an observed, approximately regionally uniform, long-term regional warming trend (chapter two). In this statistical framework, nonglacial influences upon streamflow are nonspecifically but directly incorporated into the analysis procedure as noise, and a relationship between streamflow variability and watershed glacierization is concluded only if it is detectable above this noise at a standard confidence level. No assumptions are necessary as to the ultimate origin, nature, and longevity of the driving climatic changes. The results are, of course, subject to their own set of limitations. Nevertheless, this approach offers a unique new perspective on the potential for glacial control of water resource responses to climatic warming. 6.2 Data Processing and Methodology Daily discharge data from the nine study rivers were processed (Jones, 1969) to obtain annual time series of total yearly flow volume, vr: M VT = \0.10 Omed 5.4 1.8 0.05 0.82 0.05 qmean 2.8 -1.5 >0.10 0.75 0.05 Umax 0.0 -2.6 >0.10 0.53 >0.10 qjQR 1.7 -4.6 >0.10 0.69 0.05 tmax -1.1 -0.7 >0.10 -0.17 >0.10 tc -0.9 -1.5 >0.10 0.53 >0.10 interesting to note that both fluvial systems show negative tmax trends, consistent also with the tc trend analysis. The timing metrics thus suggest that the freshet is shifting slowly forward in time, at a rate roughly equal between glacial and nival systems. A few observations regarding the reliability and meaningfulness of the trend estimates listed in Table 9 and used in the comparisons of Table 10 are in order. First, the trend directions determined here appear consistent with those of previous analyses in the region (Zhang et al., 2000; Whitfield and Cannon, 2000a,b; Zhang et al., 2001; Whitfield, 2001). In addition, applications of the two-sided nonparametric Mann-Kendall test for monotonic trend (see appendix two) tend to be suggestive of trend significance for the time series considered here. The tests were performed without deserialization (see again appendix two) and at a 96 significance level of 0.10, which is not uncommon in such applications (e.g., Burn, 1994; Cunderlik and Burn, 2002). It is also vital to note that the slope estimator of Theil (1950) and Sen (1968) returns accurate values for signal-to-noise ratios considerably lower than those for which significance tests corifirm trend at a conventional a, where we define the signal to be trend and noise to be all other time series components (appendix two). This estimator thus reliably quantifies subtle trending processes, even when these were not clearly identifiable in a particular time series using the Mann-Kendall test; recall that it is the trend estimates themselves that form the basis of the central statistical comparison of Table 10 and which are, therefore, our main concern. Finally, and perhaps most importantly, it is highly doubtful that such strong relationships between glacierization and vr trends would be found if the trend estimates were spurious. This notion may be quantified using a binomial distribution-based method as follows. Define a discrete random variate, U, such that for a given river, U - 1 if the observed trend is positive and the basin is glacierized, or if the trend is negative and the basin is glacier-free; and U = 0 if the trend is negative for a glacial river or positive for a nival river. The probability that a trend estimated in the presence of noise is exactly zero to arbitrary precision is infinitesimally small (the values of Table 9 have been truncated to two significant digits for presentation). Consequently, there is a 50/50 chance that the estimated trend in a random time series is positive or negative and, by the same token, that U- 1 or 0. Evaluating the probability that the observed vr result (££ / = 9 successes over n = 9 basins, or trials., with Bernoulli probability, pB = 0.5 for each trial) is a chance occurrence is much like determining the probability of getting nine heads in nine consecutive tosses of a fair coin, and may be evaluated using the binomial distribution: (18) where: ( n 1 n! ZU) {n-IU)!W! (19) 97 Evaluation of equation (18) yields a probability of approximately 0.002 that the aforementioned one-to-one correspondence of positive and negative vr slope estimates to basin glacierization and nonglacierization, respectively, arises from chance trends in random time series. In other words, the null hypothesis of no association between trend direction and watershed glacierization can be rejected at a < 0.01 in favour of the alternative hypothesis that glacial (nival) rivers experience positive (negative) flow volume trends. Strictly speaking, the method requires an absence of spatial correlation at length scales greater than the streamgauge separation distances; comparing the VT trend estimates of Table 9 to the station locations of Figure 1, this assumption is reasonable for total annual flow volume. 6.4 Interpretation 6.4.1 General Historically-observed meteorological and glaciological trends in the study area may be summarized on the basis of previous work as follows. Estimated meteorological trends (Zhang et al, 2000; Whitfield and Cannon, 2000a,b; Whitfield, 2001) vary somewhat depending on the datasets and methodologies used, but climate in the region is clearly growing warmer, particularly but not exclusively in summer, with the possible exception of a narrow window in late winter. Trends in total precipitation are less clear; analyses suggest either increases, or a seasonally-dependent mix of increases and decreases, so that the region is certainly not becoming drier overall, and may be growing somewhat wetter. The sole long-term glacial mass balance record from the region, developed for the Kaskawulsh Glacier using an innovative remote sensing approach, is indicative of strongly negative annual net balance trends (Arendt et al, 2002). Seasonal balance data are not available. Time series from Alaskan glaciers (e.g., data presented by Hodge et al, 1998) also suggest a negative net annual mass balance trend, manifested as positive trends in ablation and balance cycle amplitude with no clear trend in accumulation. The Alaskan glaciers lie too far from the study area to be unambiguously representative of local conditions (see previous chapter), but note that the region over which a roughly uniform long-term glacioclimatic 98 trend is observed seems intuitively likely to be more extensive than that over which a unifoirm meteorological response to a circulation pattern might be expected. This premise is lent some support by the fact that negative mass balance trends are the global norm for alpine glaciers (e.g., Dyurgerov and Meier, 1999, 2000), whereas local responses to a particular ocean-atmosphere circulation phenomenon show spatially complex teleconnectivity patterns. Although available glaciological and meteorological data in the region are inadequate to perform detailed basin-by-basin analysis or modelling, this suite of trends does paint a self-consistent overall picture of linked climatological, glaciological, and hydrological change in southwestern Yukon and northwesternmost British Columbia. Observed trends in temperature, precipitation, net annual mass balance, seasonal balances, and annual balance cycle amplitude are all strongly consistent with increased glacial meltwater production. These may also suggest that shifts in annual mass balance are dominant over balance cycle amplitude effects, with the impacts of warming exceeding those of increased precipitation, if the latter is indeed present. Both of these glaciological effects, however, are associated with a warming climate and both are fully consistent with observed increases in the total annual flow volumes of glacier-fed rivers. Declining vr values in nonglacierized catchments also suggest the dominance of increasing trends in temperature from the perspective of direct streamflow impacts. Lakes and wetlands are common in the study area, summer weather is warm, and summer day lengths are very long (see chapter two), so losses to evapotranspiration can be substantial; these apparently are increasing faster than gains in the annual water input (if positive precipitation trends are indeed present) to strictly nival streams. Our interpretations of increased summer ET losses as the source of negative vr trends in nonglacierized watersheds (which do not enjoy the benefit of increased melt production from large reservoirs of glacial ice), and increased summer glacial melt as the source of positive vr trends in glacier-fed rivers (exceeding the effects of increased ET), are also generally consistent with the qmean, qmed, and qiQR results. In addition, recognizing that freshet timing is dominated by spring and summer temperatures, the forward-moving freshets observed for both glacial and nival systems seem consistent with the warming trend. Note: that potential overall hypsometric differences between glacierized and nonglacierized watersheds, which are obviously static over the period of record, seem unlikely to yield 99 opposite time rates of change in annual river flow volume except insofar as they help determine which catchments host glaciers. 6.4.2 Baseflow Trends in VT are emphasized here because this metric describes the gross behaviour of the rivers under long-term climatic change - i.e., whether the rivers are growing larger or smaller over time. Total annual flow volume is also, of course, the single most relevant hydrometric variable from a traditional water resources perspective. Several of the supplementary metrics considered, such as qiQR, are also closely tied to VT- In contrast, winter baseflow is largely irrelevant to most conventional engineering or geomorphologic issues. However, the complex observed baseflow trend pattern - uniformly positive trends in glacierized catchments, and a mixture of trend directions in nival watersheds - is intriguing from an ecological perspective, as winter flows are of particular significance to fish survival in cold regions (see chapter seven). The likelihood that this pattern arises from chance trends in random time series is assessed using a nonparametric Monte Carlo randomization technique. In a given realization, random baseflow records are generated by scrambling the temporal order of the observations in each of the nine original time series. Trends are then estimated using equation (17). The trend pattern thus generated is taken to match the observed pattern if the following condition is satisfied: {mu >0allu = \,Ng\and {(mv > 0 any v = l,N„)and (mv < 0 any v = l,N„)} (20) where Ng and N„ are the number of glacierized and nival basins, respectively. 105 such realizations are synthesized. The probability that the observed trend pattern occurs by chance is the proportion of the ensemble that satisfies equation (20). This probability was found to be 0.02; i.e., the observed pattern is statistically significant at a < 0.05. 100 Recalling that baseflow is largely derived from aquifer storage, so that baseflow trends reflect hydrogeologic changes, and noting that surface and subsurface hydrologic catchment areas appear to correspond reasonably well within the study area (see Owen, 1967), it is reasonable to posit that the observed trend pattern arises from systematic basin-to-basin variability in the net balance of trends in aquifer recharge and aquifer properties. It was demonstrated in this chapter that total flow volume here is increasing in glacierized watersheds but decreasing in nival watersheds. However, aquifer properties are also believed to be changing due to the permafrost impacts of historical climatic warming. These watersheds lie in areas of sporadic (valley bottoms) to continuous (alpine zones) permafrost, which exerts powerful control over hydrogeologic properties (e.g. Burn, 1992; Owen, 1967). Permafrost may be regionally degrading in various fashions under warming conditions (Burn, 1992), which has been hypothesized to increase the connectivity and storage capacity of aquifers, resulting in greater winter baseflow (Woo, 1990; Michel and van Everdingen, 1994). Warming-induced trends in aquifer recharge and aquifer storage capability are in the same direction for glacierized catchments, resulting in uniformly positive baseflow trends. In contrast, such hydrogeologic trends lie in opposite directions for nival catchments, so that the net direction of groundwater resource impacts may be sensitive to interbasin variability in basic geologic characteristics, leading to a mixture of positive and negative baseflow trends in different nival watersheds. 6.5 Implications and Limitations Glacial responses to climatic variability are demonstrated here to result in opposite discharge trends from one river to the next within an otherwise hydroclimatologically uniform study area. The analysis thus offers clear empirical evidence in support of the common hypothesis that climatic warming can materially affect downstream water resources specifically via glaciological pathways, and also suggests that regional generalizations of interpreted or projected hydrologic trends may not, therefore, be tenable in variably-glacierized regions. While these observed trend patterns are dramatic, particularly for vr, most of the trends themselves are very subtle over the period of record relative to year-to-year variability. One consequence is that first-order nonstationarity is not a substantial issue for the analyses in preceding chapters, particularly given that most of the techniques 101 employed are relatively insensitive to small or moderate departures from stationarity in their current applications (e.g., CV and EOF analyses; see also chapter two). Nevertheless, the resulting cumulative long-term flow volume changes can represent a reasonably substantial shift in baseline conditions in absolute terms. For example, the apparently small values of trend in the normalized time series (Table 9) translate into sizable per annum changes in average total water throughput in the Kluane and Dezadeash Rivers, for example, of +1.2xl07 m3yr"' and -3.9xl06 m 3 yf\ respectively. The final severity of such changes, of course, depends upon the ultimate temporal persistence of the driving climatic trends. Some specific implications of potential concern include exacerbation of existing differences in the ecohydrology (e.g., Ward, 1994) and water chemistry and sediment load (e.g., Singh and Singh, 2001) of glacierized and nival systems, and spatial juxtaposition of opposite trends in riparian community composition (e.g., Law et al, 2000), water supply and hydroelectric power generation capacity (e.g., Gleick, 1987), and contaminant mobilization (e.g., Chang et al, 2001). The analysis is subject to two main limitations. First, our physical interpretations of the local mechanisms whereby given changes in meteorological forcing differentially affect adjacent glacierized and nonglacierized systems (see preceding section), although reasonable, remain largely qualitative and the subject requires closer investigation. Second, the specific details of the results may not be geographically generalizable. Recall the differences in inferred runoff responses between temperate regions of southwestern Canada and continental Europe to which we alluded in the introductory section. These were concluded by Moore and Demuth (2001) to result from differences in firn cover, which may in turn be connected with differences in relative stages of retreat, and in degrees of sensitivity to climatic variability, due to dissimilar basic glacioclimatic environments. It is very unlikely that differential streamflow trends between glacierized and nival watersheds are unique to the southwestern Canadian subarctic. However, the manner in which such differential trends are locally realized could potentially vary considerably due to contrasting climatic, glaciological and hydrological conditions, and the analysis requires replication elsewhere. 102 6.6 Conclusions Nonparametric statistical techniques were applied to historical streamflow data to determine whether rivers with and without catchment glacial cover respond in significantly different ways to a warming climate. It was found that glacier-fed rivers grew larger, whereas nival streams grew progressively smaller, over the historical record under an observed regional warming trend. Although some of these trend effects are subtle, the overall result is statistically significant at restrictive confidence levels. Combined consideration of hydrological, meteorological, and glaciological trends suggests that the streamflow consequences of increasing temperature exceed those from a possible concurrent rise in precipitation in the study area, causing increases in both glacial meltwater production and evapotranspiration; the former appears the dorninant net hydrologic effect in glacierized catchments, and the latter in glacier-free watersheds. By empirically demonstrating that catchment glacial cover can result in opposite trends in total annual flow volume from one river to the next within an otherwise hydroclimatologically uniform study area, the analysis presents strong evidence in support of the hypothesis that climatic warming can materially affect downstream water resources specifically via glaciological pathways, and also implies that regional generalizations of interpreted or projected hydrologic trends may not be tenable in variably-glacierized regions. This is the first study to directly identify glacial-nival hydroclimatological contrasts of this type; however, comparisons with prior work suggest that the local form of such contrasts may differ greatly between glaciohydroclimatic regions. 103 CHAPTER 7 ECOLOGICAL IMPLICATIONS 7.1 Introduction As noted in previous chapters, glacial-nival streamflow differences likely have substantial implications to a range of issues, including availability and forecasting of water supply and hydroelectric power resources, flood prediction and control, and reliable inference of climatic variability from hydrometric data. In an effort to give additional shape to the broader potential consequences of the statistical hydroclimatological results presented in this thesis, we shift our attention here to an additional issue: comparative glacial-nival hydroecology. There has been increasing interest over the last decade in the comparative ecology of glacial and nival rivers. Ward (1994) suggested that the hydrophysical and consequent habitat characteristics of glacier-fed and strictly nival streams were sufficiently different to warrant designation as distinct ecosystems (kryal and rhithral, respectively). Milner and Petts (1994) went on to posit a conceptual model of lotic biomass and biodiversity as a function of glacial influence. The general idea behind this model, which has subsequently seen considerable refinement (see Milner et al, 2001), is that the lower temperatures, higher turbidity, and greater channel instability of glacier-fed rivers tend to result in a relatively impoverished ecosystem. This model largely focuses on water quality and invertebrates, however, whereas more recent work has turned to streamflow characteristics and salmon stocks. Specifically, Dorava and Scott (1998) and Dorava and Milner (2001) noted that the relatively high late-summer to early-fall flows of glacial rivers coincide with major salmon runs, whereas nival streamflows in cold-region environments typically decrease to very low levels at this ecologically crucial time of year. On the basis of these observations and intriguing preliminary analyses of southeast Alaskan salmon escapement data, they suggested that watershed glacial cover engenders larger salmon populations, particularly in combination with lake cover, which is inferred to maintain the habitat benefits of glacierization while mitigating its liabilities. 104 Notwithstanding this exciting progress in the comparative ecology of loyal and rhithral systems, at least three key issues remain unaddressed. First, whereas glacial and nival rivers have been shown here to exhibit a number of substantial differences in their annual streamflow regimes (chapter three), previous consideration of the ecological consequences of seasonal discharge distinctions has been restricted to autumn flows. Second, no attempt has been made to examine the potential ecological consequences of differential glacial-nival responses to interannual climatic variability, the primary subject of this thesis. Finally, previous work has considered glacial-nival contrasts in. the populations and diversity of macroinvertebrates or populations of fish (specifically, salmon), but no assessment has been made of potential distinctions in fish species diversity. Here., I explore the potential ecological consequences of the hydrological results from previous chapters. An assessment is first made of the likely fish habitat implications of systematic differences in the annual streamflow regime, and in interannual streamflow variability at various timescales, of glacial and nival systems. A preliminary comparison of fish taxa richness between the kryal and rhithral systems is then performed using available fish species presence/absence data for the study rivers. A few of the issues considered here have been summarized in Fleming (in press). While this is the least quantitative chapter of the thesis, and strays somewhat into the domain of the life sciences, it is both appropriate and valuable to obtain and present a more concrete understanding of the potential wider implications of the comparative hyeVoclimatology of glacial and nival rivers described in the foregoing analyses. 7.2 Ecological Background 7.2.1 Study Area Abell et al. (2000) place our study region within the Yukon freshwater ecosystem, consisting of arctic river and lake habitat, with roughly uniform aquatic zoogeography. An interesting partial exception is the upper Alsek watershed; the Alsek River upstream of Lowell (Naludi) Glacier, and the Dezadeash River, contain a hybrid Yukon/Pacific species assemblage that 105 apparently arose from locally complex and intertwined glacial, hydrologic, and ecological histories {McPhail and Lindsey, 1970; Lindsey et al, 1981; Lindsey and McPhail, 1986; Pahlke and Etherton, 2001; Fisheries and Oceans Canada, 2002; see also Cruikshank, 2001). Salmonids variably present in the study rivers are arctic grayling (Thymallus arcticus), chinook (spring, king, tyee) salmon (Oncorhynchus tshawytscha), chum (dog) salmon (Oncorhynchus keta), lake trout (Salvelinus namaycush), Dolly Varden (Salvelinus malmd), freshwater sockeye salmon (Kokanee) {Oncorhynchus nerka), inconnu (Stenodus leucichthys), least cisco (Coregonus sardinella), round whitefish (Prosopium cylindraceum), and lake whitefish (Coregonus clupeaformis) (Yukon FISS database: Fisheries and Oceans Canada, 2002). Available fish distribution data for the area do not consistently differentiate between individual coregonine species, so these are treated generically here as general whitefish. Whitefish are considered a subfamily within the salmonidae in this study (e.g., Wooding, 1997); this taxonomic choice has little effect upon the analyses presented. Nonsalmonid species present in one or more of the study rivers are northern pike (Esox lucius), slimy sculpin (Cottus cognatus), longnose sucker (Catostomus catostomus), burbot (Lota lota), and lake chubb (Couesius plumbeus). 7.2.2 Data Species identified in each of the rivers are listed in Table 11. Information was retrieved from Fisheries and Oceans Canada (2002) only for the study rivers, and not from other water bodies in the same watershed, to help ensure appropriateness of comparisons between hydrologic and ecological characteristics for a given stream. In addition, Alsek data in Fisheries and Oceans Canada (2002) pertain to the lower river. Alternative upper Alsek basin biophysical data sources (Lindsey et al, 1981; Lindsey and McPhail, 1986) list species only in particular lakes rather than in the main-stem river. As there are major ecological differences between the upper and lower Alsek (see above), and only the former appears in our study region and is monitored by the ARBR streamgauge (Table 1), Alsek River fish species are not incorporated here. Note also that the Internet-based Yukon FISS database is a compilation of limited available data from disparate sources, for a very remote and 106 incompletely-studied area. Data sources are listed online but are not generally available for scrutiny; fish survey methodologies may have varied between catchments, and it is conceivable that the fish distribution data could be incomplete. Fish abundance datasets sufficiently comprehensive for interbasin comparisons do not appear to be available for the study rivers. Table 11: Fish species presence (X) or absence (—). Data are available for eight of the study rivers. species glacierized basins nival basins Takhini White Kluane Wann Wheaton M'Clintock Dezadeash Big Creek salmonids chinook salmon X X X - - X - X chum salmon X X X - - - - X arctic grayling X X X X X X X X general whitefish X X X - X X X X lake trout X - X - - - X -Dolly Varden - - - - - - X -Kokanee - - - - - - X -subtotal 5 4 5 1 2 3 5 4 nonsalmonids burbot X X X - - X - -slimy sculpin X X X X X X X X lake chub X - - - - X - -longnose sucker X X X - - X X -northern pike X - X - - X X -subtotal 5 3 4 1 1 5 3 1 all fish species total 10 7 9 2 3 8 8 5 107 7.3 Habitat Implications 7.3.1 Freshet Floods One of the starkest streamflow contrasts found here between glacial and nival regimes is the greater annual hydrologic cycle amplitude of the former (chapter three). Considering, for example, Gibson and Myers'' (1988) review of the influences of seasonal river discharge upon salmon survival, and noting that in cold regions the annual flood occurs in summer and the drought in winter, it appears that the ecological implications of discharge during the high-flow season can involve a trade-off between greater habitat availability and excessive habitat disturbance. Although it seems intuitively likely that the systematically larger freshet floods of glacierized watersheds are ecologically significant, the specific fisheries implications thus appear ambiguous. 7.3.2 Spawning Flows The environmental significance of adequate spawning-season flows has been recognized by scientists for at least half a century (Wicket, 1958), and probably much longer by subsistence, recreational, and commercial anglers. In particular, Dorava and Milner (2000) noted that spawning-season discharge controls channel navigability, oxygen availability for eggs, and food and cover abundance. They also suggested that the higher natural flows associated with kryal systems during autumn may thus present a significant advantage to spawning habitat quantity and quality, particularly given that rhithral discharges diminish to very low values in cold regions at this time of year. The spawning significance of the higher fall discharges investigated more thoroughly in chapter three, however, is by no means universal. The salmonids identified in Table 11 variously spawn in late summer through late autumn in the study area, with most being fall spawners (McPhail and Lindsey, 1970). However, arctic grayling is an important exception, spawning immediately after spring break-up. In addition, anadromous and adfluvial species are likely more sensitive to such spawning-season flow differences, although local 108 populations may vary in their migratory characteristics (McPhail and Lindsey, 1970). The ecological significance of glacially-induced high spawning-season flows also clearly differs with regional hydroclimatic characteristics. On the other hand, glacial meltwater inputs at more temperate maritime locations might moderate the effects of interannual variability in the onset timing of heavy autumn rainfall; such kryal inputs tend to be stable from year to year (see below). Finally, recall that a glacier-fed river will produce substantially higher spawning-season flows than a nival river of similar catchment area, but not relative to a nival stream chaining a far larger watershed. 7.3.3 Winter Habitat Availability The winter low-flow regime in cold regions is of particular ecological significance due to the overall physical harshness of such environments (e.g., Gibson and Myers, 1988; Cheng et al, 1993; Smith et al, 2001); the dominance of abiotic ecological controls during cold-regions winters (e.g., Cunjak and Power, 1986; Power et al, 1993; Alfredsen and Tesaker, 2002); and the severity and duration of cold-regions winter low flows (e.g., Gibson and Myers, 1988; Cheng et al, 1993; Clausen and Biggs, 2000), which are particularly acute in the arctic and subarctic (see Reynolds, 1997). The effects of cold-season flow depressions and massive attendant seasonal habitat losses are apparent in the winter biogeography of salmonids within and near the current study area. Juvenile chinook migrate substantial distances through the upper Yukon River basin to reach areas of greater winter aquatic habitat (Bradford and Grout, 2000), and viable overwintering locations for arctic grayling and arctic char in northern Yukon and Alaska are restricted to spring- (i.e., groundwater-) fed pool habitat, which remains unfrozen (Craig and McCart, 1974; Craig and Poulin, 1975). Moreover, significant positive correlations have been established between survival of young salmonids and winter discharge (Chadwick, 1981; Gibson and Myers, 1988; Cunjak et al, 1998) or winter groundwater inputs (Curry et al, 1995) in other cold regions. The larger winter specific discharge of glacier-fed rivers (chapter three), supported by greater aquifer storage and release, thus implies a superior overwintering environment on a unit- catchment-area basis relative to rhithral systems. 109 7.3.4 Attenuation of High-Frequency Interannual Variability There is abundant evidence that, at downstream locations in larger watersheds where the areal proportion of glacial coverage therefore tends to be relatively low, glaciers attenuate high-frequency interannual streamflow variability (see chapter four). The downstream water resources and lotic habitat associated with glacial systems are therefore more stable than those of nival rivers. As alluded to above, this phenomenon may intertwine with the greater autumn flows of glacial systems in an interesting way: not only do glacial rivers offer greater spawning-season habitat, but that spawning-season habitat is also more stable at an interannual level. Such effects may be of broad significance, as salmon run timing is affected in more temperate regions, such as coastal British Columbia, by pluvial inputs that can show considerable year-to-year variability, which in turn might be offset by glacierization. Two aspects of the hydrological results from chapter four add further complexity to the potential ecological implications of glacial attenuation of interannual streamflow variability. First, we discovered that such dampening does not extend to winter baseflow. Glacial systems may therefore offer more, but not more hydroclimatologically stable, winter habitat. Second, we found that the timescale at which this attenuation mechanism operates appears to be relatively short, the consequence being that glacial hydroclimatological dampening does not extend to the lower-frequency forms of variability discussed in the following two sections. The interannual stability of habitat quantity is thus increased only up to a point. 7.3.5 Selective Teleconnectivity Glacial and nival river discharges can respond in very different (indeed, potentially opposite) ways to large-scale ocean-atmosphere circulation patterns, such as ENSO, the PDO, and the A O (see chapter five). This is very likely to have ecological implications via both direct habitat availability impacts and by modifying salient water quality parameters. For example, while fisheries studies which unambiguously reflect solely the freshwater environmental impacts of such climate modes appear to be rare (see Mote et al., 2003), water temperature and discharge anomalies associated with ENSO have indeed been 110 demonstrated (at least in nonglacial environments) to have significant effects upon lotic and lentic fish populations (Gunn, 2002; Smolders et al, 2002). We found that glacial river flows in the study area are coupled to the A O but not ENSO, nival river flows are coupled to ENSO but not the AO, and glacial and nival rivers are decoupled from each other. The implication, then, does not appear to be that the habitat associated with one fluvial regime is more or less susceptible to the vagaries of large-scale ocean-atmosphere circulation patterns than that of the other (see also preceding section). Rather, habitat quality and quantity follow different schedules between glacial and nival rivers. This may have powerful consequences for interbasinal differences in the temporal varialbilities of fish populations, which may in turn have implications, for example, to the regional representativeness of lotic fish survey data. Whether diversity differences could also be expected to arise from this type of environmental control is far less clear. 7.3.6 Differential Long-Term Trends Three basic types of differential long-term trend pattern in response to regional climatic warming were identified in chapter six: (1) opposite trends in total annual flow volume, with river size increasing in glacial rivers but decreasing in nival streams; (2) approximately opposite trends in interquartile range, with freshet magnitudes increasing in glacierized watersheds but decreasing in nival and lightly-glaciered catchments; and (3) a more complex pattern of baseflow trends, being universally positive in glacierized systems but of varying sign and magnitude in nival streams. If the driving warming conditions continue, (1) and (2), and in a less consistent fashion, (3) will reinforce and exacerbate the annual streamflow regime differences identified in chapter three. Until glacial ice reservoirs begin to be exhausted, glacial rivers will thus enjoy further increases in overall water resource productivity, freshet magnitude, and autumn discharges relative to nival rivers, presumably matched by increases in the ecological consequences of these properties as discussed above. Recall, however, that long-term glaciofluvial responses to climatic warming appear to be strongly dependent on local glaciohydrological conditions (see chapter six), and a wide range of studies has also suggested region-to-region spatial variability in nonglacial 111 streamflow trends under historical warming. The trend results and their ecological imputations may therefore be site-specific. 7.4 Fish Species Richness Although some of the ecological implications discussed above, while potentially significant, do not clearly lead to either an amelioration or degradation of overall environmental conditions, a view of lotic habitat as consisting of a balance sheet with both assets and liabilities may nevertheless prove useful. Ecologically positive kryal characteristics are shown here to include substantially higher spawning-season and winter habitat availabilities, which are of potentially great ecological significance in cold regions; a long-term trend toward a deepening of these distinctions between glacial and nival systems, which has probably been in place since at least the end of the neoglacial; and less extreme interannual discharge variability, up to a point. Negative kryal characteristics (section 7.1) also clearly extend to the current study area as well (higher freshet turbidity, Table 1; lower freshet water temperatures, Table 1; and probably greater channel instability, insofar as inspection of to pographic maps shows the glacier-fed study rivers to typically be braided, whereas the nival study rivers are generally not). The ecological effects of negative kryal properties have been identified largely on the basis of macroinvertebrate studies (see Milner and Petts, 1994; Ward, 1994; Milner et al, 2001). Nevertheless, such hydrologic characteristics likely also have substantial negative implications to fish populations and diversity via both the indirect effects of a bottom-up negative trophic cascade, as well as more direct impacts (e.g., Craig and McCart, 1974; Scrimgeour et al, 1994). How do these habitat assets and liabilities compare in terms of net fisheries impacts? As noted previously, available biophysical survey data for the study rivers are restricted to piscine species presence/absence, which is adequate for an exploration of taxa richness, a simple measure (e.g., Dodds, 2002) of biodiversity. The results, consisting simply of the sums listed in Table 11, suggest little evidence for a direct relationship between degree of catcliment glacierization and the number of salmonid or total fish species present, regardless of whether the Dezadeash River and its hybrid species assemblage is included in the comparison. It would thus appear that either piscine species richness is not sensitive to such 112 otherwise ecologically salient hydrophysical consequences of watershed glacial cover, or that positive and negative kryal characteristics effectively mutually cancel from a fish species richness perspective. This negative result is, in fact, interesting, because it suggests that the largely entomologically-derived kryal-rhithral model of decreasing lotic diversity with increasing glacial influence (Milner and Petts, 1994; Milner et al, 2001) might not apply to fish. It is also interesting to note, however, that the two study rivers with the highest total number of fish species, and two of the three rivers tied for the greatest total number of salmonid species (in both cases, the Takhini and Kluane Rivers), possess both glacial influences and substantial lake cover. Like basin glacierization, degree of lake cover considered alone was seen not to correspond to taxa richness. Although they were concerned solely with salmon populations, this overall result appears at least vaguely consistent with the hypothesis of Dorava and Scott (1998) and Dorava and Milner (2000) that the combination of lake and glacial influences is optimal for fish in cold regions, since lakes preserve or accentuate at least one kryal habitat advantage (higher spawning flows) while mitigating at least one kryal habitat liability (higher turbidity) by trapping sediment. 7.5 Conclusions Fish stocks are threatened by a variety of forces including overutilization, habitat destruction, and natural and (it is widely believed) anthropogenic climatic variability. Coming to terms with these threats requires a deeper understanding of the complex hydroecological relationships which in large part shape habitat conditions. Glaciers are known to modify the hydrologic characteristics of the rivers they feed; abiotic controls of diversity and population are typically dominant in physically harsh environments, particularly cold-regions winters; and glacierization is common in several regions of prime (or formerly prime) salmonid habitat, including western North America, Scandinavia, and continental Europe. The hydroecological impacts of watershed glacial cover may assume additional urgency given the observed global-scale and possibly accelerating retreat of alpine glaciers (e.g., Dyurgerov and Meier, 2000) and its effect upon various aspects of river flow (e.g., chapter six) and, therefore, lotic habitat (e.g., Rouse et al, 1997). Linked 113 climafological, glaciological, hydrological, and ecological change may therefore prove to be a central issue in the successful long-term management of fisheries resources. This strongly suggests a need for a far better understanding of the effects of watershed glacierization upon fish habitat than we currently possess. Six types of glacial-nival hydrologic differences, identified or refined in previous chapters, were evaluated here from a piscine habitat perspective. In summary, relative to nival systems, the greater annual hydrologic cycle amplitudes of glacial rivers appear to be of ambiguous ecological significance; the higher autumn flows of glacial rivers are an important habitat asset, but primarily for fall spawners; the higher winter baseflows of glacial rivers may be a crucial environmental advantage to all species in cold regions; the attenuation of high-frequency interannual streamflow variability by moderate levels of watershed glacierization could be a stabilizing force in a variety of hydroclimatological settings; the decoupling and selective teleconnectivity of glacial and nival systems may not clearly favour one environment over the other in an average sense, but do imply potentially independent fish population variability between the two fluvial regimes; and finally, generally opposing long-term hydrologic trends in the glacial and nival systems considered imply that certain of the foregoing contrasts will grow sharper over time, provided that climatic warming continues. Some of the positive environmental conditions arising from watershed glacial cover might compensate for previously-described negative kryal characteristics, such as greater turbidity and channel instability and lower water temperatures. This is consistent with available fish species presence/absence data, which suggest no direct relationship between glacial cover and fish taxa richness; this result also implies that watershed glacierization affects macroinvertebrate and piscine diversity in very different ways. However, a tentative link was identified between species richness and the combination of glacierization and substantial watershed lake cover, which appears consistent with prior studies suggesting that this combination is optimal from a salmon population perspective. 114 CHAPTER 8 SUMMARY AND SYNTHESIS 8.1 Review Prior work has identified two ways in which glacial and nival river flows respond differently to interannual variability in climatic forcing: glacial systems show attenuated responses to year-to-year meteorological fluctuations for moderate levels of watershed glacierization (i.e., at downstream locations); and the two systems respond in contrasting, even opposite, manners to at least two modes of ocean-atmosphere circulation behaviour, ENSO and PDO. A considerable number of studies also imply that glacial and nival streamflows experience contrasting impacts from long-term climatic change, although such effects have never been empirically demonstrated. In addition, the annual cycles of glacial and nival rivers are known to possess distinctions of clear interest to studies of interannual hydroclimatic variability, insofar as interannual variations are, in general, expressed as changes in the annual cycle, particularly in hydroclimatic regions exhibiting strong seasonality. This includes most if not all glacial and nival rivers. The foregoing suite of phenomena correspond to a range of variability timescales, progressively increasing from the seasonal cycle to year-to-year variability modification, circulation pattern responses, and climate change impacts. Our understanding of the comparative hycfroclimatology of glacial and nival rivers is subject to two main, albeit somewhat related, types of limitations. The first consists of a series of knowledge gaps specific to each phenomenon. For example, while it is widely-recognized among cold-regions hydrologists that the characteristic annual flow regimes of glacial and nival rivers show timing distinctions, very little work has been done to further explore the full potential array of glacial-nival contrasts in the annual cycle. Numerous studies have identified the ability of watershed glacial cover to attenuate year-to-year hydrologic fluctuations, but these have been largely restricted to consideration of simple measures of typical yearly flow and are, therefore, rather imprecise as to how such variability modification is expressed. Only two studies have addressed the different manners in which 115 glacier- and snowmelt-fed rivers respond to large-scale ocean-atmosphere circulation patterns, and each considered only one such climate mode. And although differential glacial-nival responses to long-term climatic shifts have been so widely-hypothesized as to perhaps be regarded as a given, such effects have never actually been proven. There are, therefore, a host of unanswered - and apparently, previously unasked - questions regarding each phenomenon. The second limitation is perhaps more subtle, and regards the relationships that each of the four types of streamflow dynamical contrast may bear to both local glaciohydroclimatic conditions, as well as to each other in a given environment. Prior work regarding glacial modification of the seasonal cycle and of the severity of short-term interannual fluctuations are suggestive of a geographic universality, at least to the degree that such glacial-nival distinctions have been identified and explored. The situation for circulation patterns and climate drift, however, is not so straightforward. While the PDO is in some sense a slow-moving analog to ENSO (Zhang et al., 1997), the two modes of large-scale climate behaviour are by no means the same, and the results of Neal et al. (2002; PDO) and Lafreniere and Sharp (2003; ENSO) are therefore difficult to meaningfully compare. Nevertheless, the complex spatiotemporal teleconnections generally associated with most or all circulation patterns strongly imply that the local meteorological effects arising from a given climate mode, and thus the corresponding glaciohydrologic effects and whatever glacM-nival fluvial contrasts might be present, vary between regions. In addition, strong inter-regional differences in the observed responses of glacial river discharges to historical climatic warming (compare Moore and Demuth, 2001 and Demuth et al, 2002 against Collins, 1989 and Braun and Escher-Vetter, 1996) indicate that, if differential long-term glacial-nival trends exist, these patterns are certainly not geographically universal, in spite of the fact that negative mass balance trends (almost) are. One consequent difficulty, then, is that it may be challenging to construct robust generalizations with respect to how glacial and nival rivers respond differently to variability in interannual climatic forcing at longer timescales. A necessity for region-by-region mapping of such glacial-nival distinctions is thus implied, both in the interest of identifying local impacts for pragmatic planning purposes, and as a prerequisite for gaining an understanding of the underlying physics. Such 116 mapping has barely begun for circulation pattern responses, and not at all for climatic change impacts. Given the limitations of previous work, however, a second implication of such complex spatial heterogeneity is that we do not, in fact, know how the glacial and nival rivers in any one region respond to climatic variability. Specifically, with the partial (albeit notable) exception of Fountain and Tangborn (1985), who touched upon annual cycle differences and long-term mass balance trends in their study of interannual variability modification, each prior study has essentially been limited to one of the four phenomena listed above; furthermore, prior studies have generally been conducted in different regions, and certainly using different datasets. As a result, the manner is which a given set of rivers in a particular glaciohydroclirnatic region responds to climatic variability at a full range of timescales has never been ascertained or, indeed, explored - a considerable issue when, as we have seen, the results from one region cannot be readily transferred to another. The purpose of this thesis was to obtain a substantially more thorough understanding of the differences and similarities between glacial and nival responses to climatic variability by addressing the foregoing limitations to our current knowledge. This was accomplished through application of nonparametric statistical and time series analysis methods to historical discharge data from a set of four nival and five glacial rivers in a climatologically homogeneous region of southwest Yukon and northwest British Columbia. Each of the four phenomena discussed above were considered, and the overall approach consisted of deterrnining how streamflows responded to climatic variability at a given timescale, followed by comparison of the results between glacial and nival watersheds. All of the previously-unasked questions listed above, and others, were addressed. In addition, this appears to be the first body of work to trace out the effects of catchment glacial cover upon streamflow responses to climatic variability across a range of timescales in a single region, using a single dataset. To help bring home some of the broader potential implications of the results, ecological ramifications were also considered. The results are summarized and integrated in the following sections. 117 8.2 Summary of Results 8.2.1 Annual Flow Regimes The typical annual flow regimes of glacial and nival rivers were compared and contrasted using three statistical approaches: (i) Kolmogorov-Smirnov analysis of the cumulative distribution functions of annual hydrographs averaged across all glacial and all nival rivers, to assess flow magnitude differences; (ii) Kolmogorov-Smirnov analysis of the cumulative distribution functions, again obtained from averaging across all glacial and all nival streams, of temporal discharge accumulation, to assess flow timing distinctions independent of magnitude characteristics; and (iii) empirical orthogonal function analysis, applied to full non-deseasonalized discharge time series, using a number of data pre-processing techniques, in conjunction with various modal significance tests including a novel nonparametric Monte Carlo randomization approach, and followed by Kendall partial rank correlation, in order to evaluate the role of watershed glacierization in determining overall spatiotemporal patterns in seasonal- and regional-scale discharge variability. The analyses indicated that: (1) relative to nival systems, glacial rivers possess larger, longer freshets that peak later in the melt season; (2) glacial rivers also possess, on a unit-catchment-area basis, substantially higher winter baseflows; and (3) degree of watershed glacial cover is roughly on par with basin area in determining the amplitude of the annual hydrologic cycle, which is a good indicator of overall flow magnitudes in the hydroclimatic environment considered. In addition, comparison of average annual nival and glacial hydrographs to the typical annual temperature regime suggested that their freshet flows are snowpack- and temperature-limited, respectively, which was interpreted to indicate that the aforementioned regime distinctions indeed arise specifically from glacierization impacts, rather than as a result of broader hypsometric contrasts. It was also concluded that differences in the average annual cycles of glacial and nival rivers may be reflective of negative long-term glacial mass balance trends. 118 8.2.2 Interannual Variability Modification In order to robustly evaluate the potential of glacierization to modify year-to-year variability in different components of the annual hydrologic cycle, eight annual hydrometric variables, capturing various aspects of freshet and non-freshet flow magnitude and timing, were defined; corresponding annual time series were constructed from the available discharge data; and quantile-based coefficients of variation were calculated for each metric and river and then compared, on a metric-by-metric basis, to presence and degree of watershed glacial cover using Mann-Whitney median and Spearman rank correlation tests, respectively. It was determined that catchment glacial cover attenuates overall measures of freshet-related flow and timing, but little effect was seen upon extreme events measures, and winter baseflow appeared largely immune to the phenomenon. For those metrics showing an association, the relationship was monotonically decreasing with increasing glacial influence; this result is consistent with prior work, largely restricted to a single measure of typical annual flow, which indicated such behaviour for relatively low levels of glacierization, i.e., locations well downstream from the glacier terminus, which is the case here. The results for freshet-related metrics appear consistent with the widely-suggested albedo-reservoir mechanism for glacial attenuation of interannual hydrologic variability. Failure of winter baseflow to show signs of this phenomenon may result from the additional variability-dampening capability, common to both glacial and nival watersheds, provided by the aquifer storage mechanisms which dominate cold-season flow. The results of preliminary Lomb-Scargle analyses, which must be viewed as tentative due to the combination of modest record lengths and the notorious susceptibility of spectral analysis to time series duration, suggest that the ability of glacierization to stabilize flow is limited to variability timescales of roughly two years or less, and thus to brief, year-to-year fluctuations arising from short-lived meteorological anomalies. 8.2.3 Large-Scale Ocean-Atmosphere Circulation Patterns The impacts of El Nino-Southern Oscillation and the Arctic Oscillation, also known as the Northern Annular Mode, upon glacial and nival river flows were assessed using two 119 independent techniques. First, the regionally-coherent components of streamflow variability were extracted from the available suite of discharge time series via empirical orthogonal function analysis, using methodologies similar to those summarized in section 8.2.1 but implemented in a somewhat different manner, as well as a stacking procedure. The resulting leading-mode unrotated principal component time series, and stacked time series, corresponding to glacial and nival rivers were then linearly correlated against ENSO and A O indices using a nonparametric Monte Carlo randomization procedure. Second, composite analyses of annual hydrographs from each catchment were performed using ±0.5 standard deviation ENSO and A O index values as circulation event triggers, with significance determined on a month-by-month basis between positive and negative climate events using Mann-Whitney tests. In both cases, seasonally standardized discharge data (i.e., with the mean annual cycle removed and standardized to zero mean and unit variance) were used. The results indicate that: (1) the regionally coherent modes of glacial and nival streamflow variability are uncorrelated, suggesting a hydroclimatic decoupling of the two fluvial regimes; (2) glacial and nival annual streamflow anomalies are associated with the A O and ENSO, respectively, but not the reverse, suggesting selective teleconnectivity from a net annual water resource perspective; and (3) both glacial and nival systems are indeed impacted by both climate modes at a monthly level, but seasonal interactions are such that the aforementioned net annual selective teleconnectivity occurs. The results are tentatively but reasonably interpretable in terms of relationships between the local meteorological impacts of ENSO and the A O , which were also investigated using composite analysis, and the presence or absence of a glacial potential melt reservoir within the watershed. 8.2.4 Climate Change Impacts The potential for glacial and nival basins to experience contrasting water resource responses to Idstorically-observed long-term climatic changes, and in particular a regional warming trend, was investigated by determining robust rank-based regression slopes for the total annual flow volume of each study basin, and then comparing these against watershed glacial cover using Mann-Whitney median and Spearman rank correlation tests. It was found that, over the period of record, presence or absence of catchment glacierization controlled whether rivers in the study area grew larger or smaller over time, respectively. An 120 additional binomial distribution-based test provided further confirmation of the observed trend pattern at highly restrictive confidence levels. The results appear readily explainable vis-a-vis increased glacial melt production and higher evapotranspiration losses under the long-term warming trend. Supplementary trend analyses of additional hydrometric variables describing other aspects of the annual flow regime were generally consistent with the yearly flow volume results, but suggested that changes in overall river size are not accommodated by simple, linear hydrograph shifts. In this regard, it is of particular interest that trends in groundwater-supported winter baseflows, investigated more closely using randomization modelling due to their potential ecological import, showed relatively complex patterns which may reflect interactions between climate wanning-induced changes in aquifer recharge and the permafrost regime. 8.2.3 Kryal-Rhithral Hydroecology The potential ecological ramifications of the foregoing glacial-nival streamflow contrasts were first qualitatively explored in terms of fish habitat quality and quantity. Overall glacial rivers offer both advantages and disadvantages relative to nival systems. Impacts upon aquatic diversity were then more quantitatively investigated using available fish species presence/absence data from the study rivers. The results must be viewed as prehminary due to fisheries data limitations, but some intriguing possibilities were raised by the analysis. Little difference in taxa richness was apparent between glacial and nival ecosystems, which contrasts with previous work restricted to entomological diversity, presumably reflecting the somewhat different habitat needs of fish and lotic macroinvertebrates. However, it was found that the combination of watershed glacierization and extensive lake cover may facilitate higher fish species richness relative to both nival systems, regardless of degree of lake cover, and lake-poor glacial watersheds. The result is intuitively consistent with previous studies associating this combination of watershed properties with higher salmon populations, and appears to reflect the ability of lakes to mitigate the fish habitat disadvantages of glacierization while rnaintaining its advantages. 121 8.3 Synthesis The four sets of physical hydrological results summarized in the preceding section, along with results from previous investigators, may be synthesized as follows into a more integrated view of the comparative hycfroclimatology of glacial and nival rivers. 8.3.1 General Synthesis The shortest timescale of variability considered here is the seasonal cycle. The annual flow regimes of glacial and nival rivers exhibit strong differences in essentially every aspect of their flow and timing. Overall, glacial rivers appear considerably more productive from a water resource perspective, reflecting in essence the additional and temporally-extended seasonal contributions of glacial meltwater to the catchment; in particular, they exhibit freshets which are more extensive in both volume and duration, and also support much higher winter baseflows, all else being equal. Every component of the annual hydrograph appears affected by the presence or absence of watershed glacial cover. Most of these contrasts appear geographically universal, although certain newly-identified characteristics, such as baseflow differences, require further confirmation. Stepping up to short timescales of interannual climatic and hydrologic variability, watershed glacierization attenuates the severity of hydrologic effects arising from year-to-year meteorological fluctuations, at least for relatively low levels of watershed glacierization, which correspond to locations well downstream from the glacier terminus and where both human and ecological demands and impacts are greatest. Like contrasts in the typical annual hydrographs of glacial and nival rivers, glacial dampening nominally extends to interannual variability in all flow and timing aspects of the annual flow regime, but is in this case statistically significant only for measures of freshet magnitude and timing, as well as for overall yearly streamflow metrics which are heavily influenced by freshet flows, such as annual mean daily discharge and total annual flow volume. The dampening mechanism which essentially arises from the contrasting albedos of snow and ice and the availability of a glacial potential melt reservoir, seems to be a universal property of alpine glaciers, but appears restricted in its effectiveness to variability timescales of about two years and less. 122 While the precise physical reason for this particular upper limit is unclear, it seems intuitively likely that more protracted climatic forcing eventually simply overwhelms the internal glacial mechanism for variability attenuation. Moreover, although the two-year limit was inferred from graphical comparisons of power spectra estimated from relatively modest-duration hydroclimatic records, and therefore is approximate and tentative, the result is fully consistent with analyses corresponding to longer variability timescales, which show that glacial and nival rivers respond very differently to climate circulation modes and climatic change but do not suggest any glacial attenuation of the severity of such responses. Specifically, it is now reasonably well-established, at least for rivers in northwestern North America, that glacial and nival rivers can respond in radically contrasting ways to large-scale ocean-atmosphere circulation patterns. Modes of large-scale climate variability studied to date are the PDO, ENSO, and the AO. To first-order, such differential responses arise from the combination of a specific suite of meteorologic anomalies associated with particular phases of each climate mode, and the presence or absence of a glacial reservoir within the watershed potentially available for melt production; more intricate internal glacial hydrologic processes may also play some role. Differential responses appear primarily, and perhaps solely, present during the melt freshet. The contrasts can be striking, such as opposite responses to a given circulation pattern between adjacent glacial and nival rivers, or selective net water resource teleconnectivity of various watersheds within a given hycfroclimatic region to different circulation patterns depending on presence or absence of catchment glacierization. The form of such contrasts in glacial and nival responses to a given circulation pattern seem to vary substantially between gkciohydroclimatic regions, primuily as a result of differing meteorological signatures arising locally from that climate mode. At the longest timescales that can be empirically addressed using a finite-length hydroclimatic time series, glacial and nival rivers show different long-term trends in response to historically-observed climate changes. The current study is the first to directly investigate this phenomenon, and diametrically opposite trends in river size over the historical record were observed between glacial and nival rivers under progressive climatic warming, with watershed glacierization determining whether the study rivers grew larger or 123 smaller over time. The contrast in water resource impacts appear to be statistically clearest for total annual flow volume, consistent with the physical interpretation of such differential trends as the product of wairning-induced changes in total glacial melt production and total evapotranspiration losses. The trend direction differences are nevertheless clearly expressed in other, more specifically freshet-related, flow magnitude metrics. This is to be expected as freshet flows largely control total annual flow volume, although measures of flow timing associated with the melt freshet do not appear to be differentially affected between the fluvial regimes. Long-term trends appear to extend to measures of non-freshet discharge, but the more complicated pattern of winter baseflow shifts may reflect a broader array of process controls, including warming-related changes in the permafrost regime which are of less importance to freshet flows and total annual flow volume. Differences in observed glacial streamflow trends between glaciohydroclimatic regions imply that the specific suite of hydrologic trend contrasts observed here are unlikely to appear elsewhere, although more generally, glacial-nival contrasts in long-term trends should likely be expected in other regions. Note that the monotonic progression in the foregoing discussion from shorter to longer timescales of variability in climatic forcing and hydrologic response is, in a sense, also circular. Specifically, the long-term responses of glacial and nival rivers to climatic change have, in turn, been interpreted both here and previously to play an important role in seasonal hydrologic cycle differences between the two fluvial regimes. 8.3.2 Timescale Relationships The generalizations that follow are based on a rather limited number of studies, and are thus tentative and no doubt subject to considerable elaboration, revision, or outright invalidation on the basis of future research. Nevertheless, three broad patterns relating to variability timescales seem at least dimly visible and worthy of note. First, the longer the timescale of the process considered, the less geographically universal the specific local expression of that process is likely to be. While the four phenomena considered here - glacial-nival contrasts in the annual flow regime, responses to high-124 frequency interannual meteorological fluctuations, climate mode effects, and climate change impacts - all appear to be fully general, the manner in which the two shorter-timescale processes are locally manifested appears to be more-or-less independent of location, whereas the two slower processes are realized quite differently in different gkciohydroclimatic regions. It is not clear that this behaviour is entirely a true function of timescale per se. The spatial variability apparent in the direction of glacial streamflow trends under almost-universal warming and negative mass balance conditions likely results solely from the local antecedent and present glacial conditions, such as the size of the glacier and the relative amount of firn exposed. Observed spatial variability in climate mode response likely reflects both the local general glaciohydrologic environment and, in particular, the complex meteorological teleconnectivity patterns associated with ocean-atmosphere circulation modes. The ability of glaciers to dampen short-term interannual hydrologic fluctuations arises from a purely internal and fundamental glacial albedo/reservoir mechanism, seemingly with physical limits that constrain the timescales at which it can operate, and apparently restricted in its effectiveness to catchments wherein degree of glacial cover above the streamgauge is relatively low, but otherwise independent of the local glacial, climatic, or hydrologic environment. Finally, the geographic universality of annual flow regime differences no doubt arises largely from the planetary scales of both the seasons and, perhaps, negative alpine glacier mass balances, as well as a certain commonality to the inner hydrologic workings of mountain glaciers in a very general sense (e.g., all possess some measure and manner of ablation-season internal liquid water storage). Thus, there may be no overarching, fundamental relationship between timescales and spatial scales of variability in the current context. Second, the longer the timescale considered, the more restricted the differential glacial-nival response is to melt-season flows. Differences in the average seasonal cycles of glacial and nival rivers permeate the entire annual hydrograph; glacial attenuation of high-frequency interannual variability seems to extend throughout the entire year, but is significant only for hydrometric variables corresponding to, or heavily influenced by, the melt freshet; and contrasting responses to intermediate-timescale climatic forcing arising from set modes of ocean-atmosphere circulation appear strictly limited to (approximately) March through about September. Long-term trends provide the spanner in the works with respect to this 125 generalization, for winter baseflow trend differences are apparent between glacial and nival systems, in addition to those corresponding to freshet or freshet-dominated streamflow metrics. The observed baseflow trend patterns seem likely to reflect the addition of ground thermal regime influences, which exert proportionately stronger influences upon winter baseflow relative to summer freshet discharges. Overall, however, the majority of the hyfooclimatic information content of cold-regions streamflow time series does appear to reside in measurements taken during, or in hydrometric variables representing or dominated by, the melt freshet. The third generalization hypothesized here with respect to timescale may be the most robust: the longer the timescale of the process, the more subtle the corresponding glacial-nival contrast is likely to be. While elaborate procedures were required to quantitatively and rigorously describe and verify such differences, many of the distinctions between the annual flow regimes of glacial and nival rivers were quite obvious upon visual inspection of the annual hydrograph. The ability of alpine glaciers to moderate the discharge variability of the rivers they feed was not so clear but, with appropriate analysis, was shown to be entirely unambiguous with respect to most freshet-related hydrometric variables. Responses to large-scale ocean-atmosphere circulation patterns, although ultimately found to exhibit a striking selective net teleconnectivity, were far more challenging to identify, describe, and confirm. And while watershed glacierization was seen to impart a similarly dramatic signature upon regional patterns of fluvial responses to historical warming trends, the trends themselves generally represented a rather nuanced aspect of the overall streamflow records. This seeming increase in the subtlety of glacial-nival contrasts with increasing timescale of forcing is undoubtedly a product, in part, of the local gkciohydroclimatic conditions; the relatively weak meteorological (and therefore glacial and hydrological) responses within the current study area to large-scale circulation patterns provide an obvious example. Nevertheless, there is often a very general tendency, perhaps due simply to the artefactual but unavoidable effects of possessing only relatively short historical climate records, for shorter-term fluctuations to overwhelm longer-term variability: a structural trend barely visible in a 30-year time series, for instance, may become the dominant feature of a 300-year record. Consequently, it is reasonable to expect that the apparent strength of glacial-nival 126 contrasts inferred from time series of a given duration should, as a very general rule, decrease as the timescale corresponding to the type of process considered increases. 8.4 Recommendations It is evident from the above summary and synthesis that, although the four types of glacial modification of streamflow variability considered still appear to be distinct processes, the results garnered from each are mutually consistent, and in some ways closely related. Together with the more basic observation that all of these four types of glacial-nival distinction correspond to the response of the same general (glacial or nival) system to the same climatic forcing, the differences fundamentally lying solely in the timescales of forcing and response, this suggests that a comprehensive theory for the comparative hydroclimatology of glacial and nival rivers will ultimately be attainable. To this end, at least two major directions for future research can be defined. Noting that, overall, our knowledge of these atmosphere-cryosphere-hydrosphere relationships remains at a relatively early stage, the first is to further develop our empirical understanding of comparative glacial-nival hydroclimatology through continued application of statistical and time series analysis techniques to a far more comprehensive array of glaciohydroclimatic environments, considering the fullest range of variability timescales observable from historical data. Aside from its obvious inherent merit, such an approach would also help focus future physical modelling efforts by delineating the types and patterns of behaviour they must replicate and explain. The work presented in this thesis can serve as a useful template for empirical studies in other regions. Clearly, certain methodological details may require adjustment on a case-by-case basis. The relatively intricate EOF-correlation and stack-correlation procedures of chapter five, for instance, might be replaced or complemented by much simpler correlation mapping in regions showing stronger responses to large-scale ocean-atmosphere circulation patterns; and additional methods can be incorporated as appropriate, such as split-sample tests for intervention as applied by (for example) Neal et al. (2002) to longer time series in their PDO study. Likewise, some of the new technical developments presented here (e.g., nonparametric testing of EOF modal significance) are general whereas the applicability of 127 others (e.g., baseflow trend assessment by Monte Carlo randomization) are context-dependent; and sophisticated new techniques which have not previously been applied to questions of comparative glacial-nival hydroclimatology (such as nonlinear EOF analysis, inappropriate for the current study due to data availability constraints; see appendix one) might prove highly fruitful in future studies elsewhere. More broadly, however, the general approach and most of the specific methodologies employed here seem fully transportable. It may be tactically useful, to list a few locales to which attention might be particularly directed. Southern hemisphere rivers appear to have been grossly under-represented in such work, so consideration of glacial-nival hydroclimatological contrasts in Andean and New Zealand watersheds may be valuable. The north Cascades of Washington are also an appropriate region for further study, for essentially the opposite reason: a great deal of groundwork has been laid here by previous investigators, including but not limited to Fountain and Tangborn (1985), Pelto and Riedel (2001), and Mote (2003); relatively high spatial densities and long durations of hydroclimatic data are available; a considerable suite of hydrometric stations have already been identified for the purpose of comparing glacial and nival streamflow responses to climatic variability {Fountain and Tangborn, 1985); and the results could have very direct socio-economic consequences for a relatively densely-populated and fast-growing region of North America. In addition, the Torngat Mountains of northern Labrador (see, for example, Rogerson, 1986) may offer a particularly intriguing study site due to the globally anomalous regional climate trends historically observed there, consisting of a strong net annual cooling trend, with little or no long-term shift in ablation-season temperatures (Zhang et al, 2000). Finally, the variably-glacierized rivers of northern India and Pakistan, and of the Indus basin (see Singh and Kumar, 1997) in particular, should be considered due to the extraordinary geopolitical significance of their associated water resources (e.g., de Villiers, 2000). A second general direction for future work consists of examining such statistically-determined phenomena using process-based physical models. Relatively few watershed models directly incorporate glacial runoff effects. Exceptions are the UBC Watershed Model (see Quick, 1995 for a review), and the FIBV model and its derivatives (e.g., Moore, 1993), which could thus provide excellent starting points. A perhaps more philosophically 128 satisfying, potentially more instructive, but certainly more difficult approach could be to couple state-of-the-art watershed and glacial models that simultaneously solve the fundamental mathematical physics of glacial, riparian, and terrestrial hydrology; good candidates might include the models of van der Kwaak (1998) and Flowers (2000). Yet more ambitiously, one might also couple whatever physical glacier-watershed hydrology model is ultimately chosen to a meso- or global-scale climate model. If our purpose is to better understand how and why glacial and nival rivers respond differently to a given set of climatic forcing conditions, however, the vast increase in total model complexity and uncertainty induced by the addition of a climate modelling component might serve only to obfuscate the watershed dynamical processes of interest. Beyond the question of model selection and/or construction, it may also be challenging to obtain reliable external constraints upon the values of the considerable number of physical parameters typically required by such models, particularly in remote, data-sparse regions. In addition, the overall modelling problem would need to involve not only the usual goal of best-fit model calibration to observed hydrographs, but also reproduction of broad and sometimes subtle patterns of system behaviour, as well as provision of detailed physical explanations as to how such behaviour arises. Overall, it would seem sensible to implement the physical modelling approach within a paired research basin framework, a venerable concept for experimentation and analysis in the hydrologic sciences, which in this case would consist of two or more adjacent, previously well-characterized, glacial and nival study catchments. 129 REFERENCES Aalto, R., L. Maurice-Bourgoin, T. Dunne, D.R. Montgomery, C.A. Nittrouer, and J.-L. Guyot, 2003. Episodic sediment accumulation on Amazonian flood plains influenced by El Nifio/Southern Oscillation, Nature, 425,493-497. Abell, R.A., D.M. Olson, E. Dinerstein, P.T. Hurley, J.T. Diggs, W. Eichbaurn, S. Walters, W. WettengeL T. Allnutt, C.J. Loucks, and P. Hedao, 2000. Freshwater Ecoregions of North America: a Conservation Assessment. Island Press, Washington DC, 319 pp. Alfredsen K., and E. Tesaker, 2002. Winter habitat assessment strategies and incorporation of winter habitat in the Norwegian habitat assessment tools, Hydrological Processes, 16, 927-936. Alley, R.B., 2003. Comment on "When earth's freezer door is left ajar," EOS, Transactions of the American Geophysical Union, 84(33), 315. Anderson, J.E., S.-Y. Shiau, and D. Harvey, 1991. Preliminary Investigation of Trend/Patterns in Surface Water Characteristics and Climate Variations, In Kite, G.W., K.D. Harvey (eds.), Using Hydrometric Data to Detect and Monitor Climatic Change, Proceedings of NHRI Workshop No. 8, National Hydrology Research Institute, Saskatoon, 189-201. Anderson, R.L., 1942. Distribution of the serial correlation coefficient, Annals of Mathematical Statistics, 13(1), 1-13. Arendt A.A. , K.A. Echelmeyer, W.D. Harrison, C S . Lingle, V . C . Valentine, 2002. Rapid wastage of Alaska glaciers and their contribution to rising sea level Science, 297, 382-386. Arnell, N. , 2002. Hydrology and Global Environmental Change. Prentice Hall, Harlow, 346 PP-130 Barnett, A.P., 1974. Hydrological studies of the Slims River, Yukon, June-August 1970. In Bushnell, V . C . and M.G. Marcus (eds.), Icefield Ranges Research Project, Scientific Results, Volume 4. American Geographical Society/Arctic Institute of North America, New York/Montreal, 143-150. Bitz, C .M. , and D.S. Battisti, 1999. Interannual to decadal variability in climate and the glacier mass balance in Washington, western Canada, and Alaska, Journal of Climate, 12,3181-3196. Bostock, H.S., 1969. Kluane Lake, Yukon Territory, Its Drainage and Allied Problems, GSC Paper 69-28. Geological Survey of Canada, Ottawa, 19 pp. Box, G.E.P., and D.R. Cox, 1964. An analysis of transformations, Journal of the Royal Statistical Society Series B, 26(2), 211-252. Bradford, M.J., and J.A. Grout, 2000. The importance of small streams as salmon habitat in the upper Yukon River basin. 51st American Association for the Advancement of Science Conference, Whitehorse, Yukon. Available online at http://www.taiga.net/ arctic2000/abstracts/bradford.html. Bralthewaite, R.J., and O.B. Olesen, 1988. Effect of glaciers on annual run-off, Johan Dahl Land, South Greenland, Journal of Glaciology, 34,200-207. Bras, R.L., 1990. Hydrology: An Introduction to Hydrologic Science. Addison-Wesley, Reading, M A , 643 pp. Bras, R.L., and I. Rodriguez-Iturbe, 1993. Random Functions and Hydrology. Dover, New York, 559 pp. Braun, L.N. , and H. Escher-Vetter, 1996. Glacial discharge as affected by climate change, Interpraevent 1996 - Garmisch Partenkirchen, 1, 65-74. 131 Brito-Castillo, L . , A . V . Douglas, A . Leyva-Contreras, and D . Lluch-Belda, 2003. The effect o f large-scale circulation on precipitation and streamflow in the Gul f of California continental watershed, InternationalJournal of Climate, 23, 751-768. Bryan, M . L . , 1972. Variations in quality and quantity of Slims River water, Yukon Territory, Canadian Journal of Earth Sciences, 9, 1469-1478. Bum, C.R. , 1992. Recent ground warming inferred from the temperature in permafrost near Mayo, Yukon Territory. In Dixon, J.C. and A . D . Abrahams, A . D . (eds.), Periglacial Geomorphology, Wiley, Chichester, pp. 327-350. Burn, D . H . , 1994. Hydrologic effects o f climatic change in west-central Canada, Journal of Hydrology, 160, 53-70. Bum, D . H . , and M . A . H . Elnur, 2002. Detection o f hydrologic trends and variability, Journal of Hydrology, 255, 107-122. Bum, D . H . , and E . D . Soulis, 1991. The use o f hydrological variables in detecting climate change: possibilities for single station and regional analysis. In Kite, G .W. and K . D . Harvey (eds.), Using hydrometric Data to Detect and Monitor Climatic Change, Proceedings of NHRI Workshop No.8, National Hydrology Research Institute, Saskatoon, 121-130. Burroughs, W.J . , 1992. Weather Cycles: Real or Imaginary? Cambridge University Press, Cambridge, 201 pp. Cahill, A . T . , 2002. Determination of changes in streamflow variance by means o f a wavelet-based test, Water Resources Research, 38(6), doi: 10.1029/2000WR000192, 1-1-1-14. Cayan, D.R. , and D . H . Peterson, 1989. The influence o f north Pacific atmospheric circulation on streamflow in the west. In Peterson, D . H . (ed.), Aspects of Climate 132 Variability in the Pacific and the Western Americas, Geophysical Monograph 55, American Geophysical Union, Washington DC, 375-397. Chadwick, E.M.P., 1981. Stock-recruitment relationship for Atlantic salmon (Salmo salar) in Newfoundland rivers, Canadian Journal of Fisheries and Aquatic Sciences, 39, 1496-1501. Chang, H. , B .M. Evan, and D.R. Easterling, 2001. The effects of climate change on stream flow and nutrient loading, Journal of the American Water Resources Association, 37(4), 973-985. Chen, H-L, and A.R. Rao, 2002. Testing hydrologic time series for stationarity, Journal of Hydrologic Engineering, 7(2), 129-136. Chen, J., and A. Ohmura, 1990a. On the influence of Alpine glaciers on runoff. In Lang, H. and A. Musy (eds.), Hydrology in Mountainous Regions I, Hydrological Measurements, The Water Cycle, IAHS Publication 193, International Association of Hydrological Sciences, Wallingford, UK, 117-125. Chen, J., and A. Ohmura, 1990b. Estimation of alpine glacier water resources and their change since the 1870s. In Lang, H. and A. Musy (eds.), Hydrology in Mountainous Regions I, Hydrological Measurements, The Water Cycle, IAHS Publication 193, International Association of Hydrological Sciences, Wallingford, UK, 127-135. Cheng, H. , S. Leppinen, and G. Whitley, 1993. Effects of river ice on chemical processes. In Prowse, T.D., and N.C. Gridley (eds.), Environmental Aspects of River Ice, National Hydrology Research Institute, Saskatoon, 75-96. Church, M . , 1974. Hydrology and permafrost with reference to North America. In Demers, J. (ed.), Permafrost Hydrology: Proceedings of Workshop Seminar, Canadian National Committee, International Hydrological Decade, Environment Canada, Ottawa, 7-20. 133 Clague, J.J., and V.N. Rampton, 1982. Neoglacial Lake Alsek, Canadian Journal of Earth Sciences, 19,94-117. Clausen, B., and B.J.F. Biggs, 2000. Flow variables for ecological studies in temperate streams: groupings based on co variance, Journal of Hydrology, 237, 184-197. Coe, M.T., 2000. Modeling terrestrial hydrological systems at the continental scale: testing the accuracy of an atmospheric G C M , Journal of Climate, 13, 686-704. Collins, D.N., 1989. Hydrometeorological conditions, mass balance and runoff from alpine glaciers. In Oerlemans, J. (ed.), Glacier Fluctuations and Climate Change, Kluwer, Dordrecht, 235-260. Conover, W.J., 1971. Practical Nonparametric Statistics. Wiley, New York, 462 pp. Corcoran, C D . , and C R . Mehta, 2002. Exact level and power of permutation, bootstrap, and asymptotic tests of trend, Journal of Modern Applied Statistical Methods, 1(1), 42-51. Coulibaly, P., F. Anctil, P. Rasmussen, and B. Bobee, 2000. A recurrent neural networks approach using indices of low-frequency climatic variability to forecast regional annual runoff, Hydrological Processes, 14, 2755-2777. Craig, P.C., and P.J. McCart, 1974. Fall spawning and overwintering areas of fish populations along routes of proposed pipeline between Prudhoe Bay and the Mackenzie Delta. In McCart, P.J. (ed.), Fisheries Research Associated with Proposed Gas Pipeline Routes in Alaska, Yukon, and Northwest Territories, Biological Report Series, Volume 15, Canadian Arctic Gas Study Limited, Calgary, 1-24. Craig, P . C , and V.A. Poulin, 1975. Movements and growth of arctic grayling (Thymallus arcticus) and juvenile arctic char (Salvelinus alpinus) in a small arctic stream Alaska, Journal of the Fisheries Research Board of Canada, 32, 689-697. 134 Cruikshank, J., 2001. Glaciers and climate change: perspectives from oral traditions, Arctic, 54, 377-393. Cullen, H.M. , A. Kaplan, P.A. Arkin, and P.B. Demenocal, 2002. Impact of the North Atlantic Oscillation on Middle Eastern climate and streamflow, Climatic Change, 55, 315-338. Cunderlik, J.M., and D.H. Burn, 2002. Local and regional trends in monthly maximum flows in southern British Columbia, Canadian Water Resources Journal, 27(2), 191-212. Cunjak, R.A., and G. Power, 1986. Winter habitat utilization by stream resident brook trout (Salvelinus fontinalis) and brown trout (Salmo trutta), Canadian Journal of Fisheries and Aquatic Sciences, 43, 1970-1981. Cunjak, R.A., T.D. Prowse, and D. L. Parrish, 1998. Atlantic salmon (Salmo salar) in winter: "the season of parr discontent"? Canadian Journal of Fisheries and Aquatic Sciences,."55, 161-180. Curry, R.A., D.L.G. Noakes, and G.E. Morgan, 1995. Groundwater and the incubation and emergence of brook trout (Salvelinus fontinalis), Canadian Journal of Fisheries and Aquatic Sciences, 52,1741-1749. Daniel W.W., 1990. Applied Nonparametric Statistics, 2nd Ed. PWS-Kent, Boston, 635 pp. Davidovich, N.V., and M.D. Ananicheva, 1996. Prediction of possible changes in glacio-hydrological characteristics under global warming: southeastern Alaska, USA, Journal ofGlaciology, 42(142), 407-412. Demuth, M.N. , A. Pietroniro, and T.B.M.J. Ouarda, 2002. Streamflow regime shifts resulting from recent glacier fluctuations in the eastern slopes of the Canadian Rocky 135 Mountains. In Towle, J.E. (ed.), Conference Program, CWRA 55th Annual Conference, Canadian Water Resources Association, Winnipeg, 35. Deutsch, C.V., and A . G . Journel, 1998. GSLIB: Geostatistical Software Library and User's Guide, Td Ed. Oxford, New York, 369 pp. de Villiers, M . , 2000. Water. Stoddart, Toronto, 422 pp. Dodds, W.K., 2002. Freshwater Ecology. Academic Press, San Diego, 569 pp. Dorava, J.M., and K . M . Scott, 1998. Role of glaciers and glacial deposits in the Kenai River watershed and the implications for aquatic habitat. In Gray, J.E. and J.R. Riehle (eds.), Geologic Studies in Alaska by the U.S. Geological Survey, 1996, USGS Professional Paper 1595, U.S. Geological Survey, Washington DC, 3-8. Dorava, J.M., and A . M . Milner, 2000. Role of lake regulation on glacier-fed rivers in enhancing salmon productivity: the Cook Inlet watershed, south-central Alaska, USA, Hydrological Processes, 14, 3149-3159. Dracup, J.A., and E. Kahya, 1994. The relationships between U.S. streamflow and La Nifia events, Water Resources Research, 30,2133-2141. Dyurgerov, M.B., and M.F. Meier, 1999. Analysis of winter and summer glacier mass balances, Geografiska Annaler, 81A(4), 541-554. Dyurgerov, M.B., and M.F. Meier, 2000. Twentieth century climate change: evidence from small glaciers, Proceedings of the National Academy of Science, 97, 1406-1411. Eaton, B., M . Church, and D. Ham, 2002. Scaling and regionalization of flood flows in British Columbia, Canada, Hydrological Processes, 16, 3245-3263. Edgington, E.S., 1995. Randomization Tests, 3rdEd. Marcel Dekker, New York, 409 pp. 136 Eltahir, E.A.B., and G. Wang, 1999. Nilometers, El Nino, and climate variability, Geophysical Research Letters, 26, 489-492. Enfield, D.B., A . M . Mestas-Nunez, and P.J. Trimble, 2001. The Atlantic multidecadal oscillation and its relation to rainfall and river flows in the continental U.S., Geophysical Research Letters, 28, 2077-2080. Environment Canada, 1999. HYDAT CD-ROM v.99-2.00. Water Survey of Canada, Ottawa. Ferguson, R.I., 1985. Runoff from glacierized mountains: A model for annual variation and its forecasting, Water Resources Research, 21(5), 702-708. Fisheries and Oceans Canada, 2002. Yukon Fisheries Information Summary System (FISS). Fisheries and Oceans Canada, Pacific Region, Habitat Enhancement Branch, Delta, British Columbia. Available online at http://habitat.rhq.pac.dfo-mpo.gc.ca/cfdocs/fiss /dcfOl.cfrn Fleming, S.W., in review. Impacts of climatic trends upon groundwater resources, aquifer-stream interactions, and freshwater habitat in glacierized watersheds, Yukon Territory, Canada. In Knight, P.G. (ed.), Glaciology and Earth's Changing Environment, Blackwell, Oxford. Fleming, S.W., in press. Comparative analysis of glacial and nival streamflow regimes with implications for lotic habitat quantity and fish species richness, River Research and Applications. Fleming, S.W., and G.K.C. Clarke, 2002. Autoregressive noise, deserialization, and trend detection and quantification in annual river discharge time series, Canadian Water Resources Journal, 27(3), 335-354. 137 Fleming, S.W., and G.K.C. Clarke, 2003. Glacial control of water resource and related environmental responses to climatic warming: empirical analysis using historical streamflow data from northwestern Canada, Canadian Water Resources Journal, 28(1), 69-86. Fleming, S.W., and G.K.C. Clarke, in review. Attenuation of high-frequency interannual streamflow variability by watershed glacial cover, ASCE Journal of Hydraulic Engineering. Fleming, S.W., A . M . Lavenue, A .H. Aly, and A. Adams, 2002. Practical applications of spectral analysis to hydrologic time series, Hydrological Processes, 16, 565-574. Flowers, G.E., 2000. A Multicomponent Coupled Model of Glacier Hydrology. PhD thesis, University of British Columbia, 265 pp. Fountain, A .G. , and W.V. Tangborn, 1985. The effect of glaciers on streamflow variations, Water Resources Research, 21, 579-586. Fountain, A.G. , and J.S. Walder, 1998. Water flow through temperate glaciers, Reviews of Geophysics, 36(3), 299-328. Fowler, H.J., C.G. Kilsby, and P.E. O'Connell, 2003. Modeling the impacts of climatic change and variability on the reliability, resilience, and vulnerability of a water resource system, Water Resources Research, 39(8), doi: 10.1029/2002WR001778, 10-1-10-11. Franks, S.W., and G. Kuczera, 2002. Flood frequency analysis: evidence and implications of secular climate variability, New South Wales, Water Resources Research, 38(5), doi: 10.1029/2001WR000232, 20-1-20-7. Fulton, R.J., 1995. Surflcial Materials of Canada, GSC Map 1880A. Geological Survey of Canada, Ottawa. 138 Gibson, R.J., and R.A. Myers, 1988. Influence of seasonal river discharge on survival of juvenile Atlantic salmon, Salmo salar, Canadian Journal of Fisheries and Aquatic Sciences, 45, 344-348. Gleick, P.H., 1987. Global climate changes and regional hydrology: impacts and responses. In Solomon, S.I., M . Beran, and W. Hogg (eds.), The Influence of Climate Change and Climatic Variability on the Hydrologic Regime and Water Resources, IAHS Publication 168, International Association of Hydrological Sciences, Wallingford, 389-402. Gordey, S.P., and A.J. Makepeace, 2001. Bedrock Geology Map, Yukon Territory, GSC Open File 3754 / DIAND Open File 2001-1. Geological Survey of Canada / Indian and Northern Affairs Canada, Ottawa. Gunn, J.M., 2002. Impact of the 1998 El Nifio event on a lake charr, Salvelinus namaycush, population recovering from acidification, Environmental Biology of Fishes, 64, 343-351. Haggerty, R., S.A. McKenna, and L .C . Meigs, 2000. On the late-time behavior of tracer test breakthrough curves, Water Resources Research, 36(12), 3467-3479. Haggerty R., S.M. Wondzell, and M.A. Johnson, 2002. Power-law residence time distribution in the hyporheic zone of a 2nd-order mountain stream, Geophysical Research Letters, 29(13), doi: 10.1029/2002GL014743, 18-1-18-4. Hamlet, A.F. , and D.P. Lettenmaier, 1999a. Columbia River streamflow forecasting based on ENSO and PDO climate signals, Journal of Water Resource Planning and Management, 125, 333-341. Hamlet, A.F. and D.P. Lettenmaier, 1999b, Effects of climate change on hydrology and water resources in the Columbia River basin, Journal of the American Water Resources Association, 35(6), 1597-1623. 139 Harrison, D.E. , and N.K. Larkin, 1998. El Nino-Southern Oscillation sea surface temperature and wind anomalies, 1946-1993, Reviews of Geophysics, 36, 353-399. Henshaw, F.F., 1933. Notes on variation of runoff on the Pacific slope, Transactions of the American Geophysical Union, 14,431-435 Hirsch, R.M., J.R. Slack, and R.A. Smith, 1982. Techniques of trend analysis for monthly water quality data, Water Resources Research, 18(1), 107-121. Hirsch, R.M. , J.R. Slack, and R. James, 1984. A nonparametric trend test for seasonal data with serial dependence, Water Resources Research, 20(6), 727-732. Hisdal, H. , and O.E. Tveito, 1993. Extension of runoff series using empirical orthogonal functions, Hydrological Sciences Journal, 38(1/2), 33-49. Hodge, S.M., D.C. Trabant, R.M. Krirnmel, T.A. Heinrichs, R.S. March, and E.G. Josberger, 1998. Climate variations and changes in mass of three glaciers in western North America, Journal of Climate, 11, 2161-2179. Hodgkins, R., 1997. Glacier hydrology in Svalbard, Norwegian High Arctic, Quaternary Science Reviews, 16, 957-973. Hoerling, M.P., A. Kumar, M . Zhong, 1997. El Nino, La Nifta, and the nonlinearity of their teleconnections, Journal of Climate, 10,1769-1786. Horel, J.D., and J.M. Wallace, 1981. Planetary-scale atmospheric phenomena associated with the Southern Oscillation, Monthly Weather Review, 109, 813-829. Hsieh, W.W., 2004. Nonlinear multivariate and time series analysis by neural network methods, Reviews of Geophysics, 42, doi:10.1029/2002RG000112, RG1003 1-25. 140 Hsieh, W.W., Yuval, J. Li , A. Shabbar, and S. Smith, 2003. Seasonal prediction with error estimation of Columbia River streamflow in British Columbia, Journal of Water Resource Planning and Management, 129, 146-149. Huang, N.E. , Z. Shen, S.R. Long, M.C. Wu, H.H. Shin, Q. Zheng, N-C Yen, C C . Tung, and H.H. Liu, 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proceedings of the Royal Society of London Series A, 454, 903-995. IPCC, 2001a. Climate Change 2001: Impacts, Adaptation and Vulnerability, Contribution of Working Group II to the Third Assessment Report of the Intergovernmental Panel on Climate Change. United Nations, New York. Available online at http://www.grida.no/ climate/ipcc_tar/wg2/. IPCC, 2001b. Climate Change 2001: The Scientific Basis, Contribution of Working Group I to the Third Assessment Report of the Intergovernmental Panel on Climate Change. United Nations, New York. Available online at http://www.grida.no/climate/ipcc_tar/ wgl/. Jackson, D.A., 1993. Stopping rules in principal components analysis: a comparison of heuristical and statistical approaches, Ecology, 74, 2204-2214. JISAO, 2002. Online Climate Data Archive, Joint Institute for the Study of the Atmosphere and Ocean. University of Washington, Seattle. Available online at http://tau.atmos. washington.edu/data_sets/#time_series. Johnson, P.G., 1986. Holocene paleohydrology of the St. Elias Mountains, British Columbia and Yukon, Geographie Physique et Quaternaire, 40,47-53. Jolliffe, I.T., 1987. Rotation of principal components: some comments, Journal of Climatology, 7, 507-510. 141 JolMe, I.T., 1993. Principal component analysis: a beginner's guide - II, pitfalls, myths, and extensions, Weather, 48,246-253. Jones, R.E., 1969. Approximate integrator of functions tabulated at arbitrarily spaced abscissas, Report SC-M-69-335. Sandia Laboratories, Albuquerque. Code available online at http://www.netlib.org/slatec/src/. Jury., M.R., 1998. Statistical analysis and prediction of KwaZulu-Natal climate, Theoretical and Applied Climatology, 60, 1-10. Jury, W.A., W.R. Gardner, and W.H. Gardner, 1991. Soil Physics, 5th Ed. John Wiley, New York, 328 pp. Kahya, E„ and J.A, Dracup, 1993. U.S. streamflow patterns in relation to the El Nino/Southern Oscillation, Water Resources Research, 29,2491-2503. Karl, T.R., A.J. Koscielny, and H.F. Diaz, 1982. Potential errors in the application of principal component (eigenvector) analysis to geophysical data, Journal of Applied Meteorology, 21,1183-1186. Kendall, M . , 1942. Partial rank correlation, Biometrika, 32,277-283. Kendall M . , and J.D. Gibbons, 1990. Rank Correlation Methods, 5th Ed. Oxford University Press, Oxford, 260 pp. Kirchner, J.W., X. Feng, and C. Neal, 2000. Fractal stream chemistry and its implications for contaminant transport in catchments, Nature, 403, 524-527. Kirchner, J.W., X. Feng, and C. Neal, 2001. Catchment-scale advection and dispersion as a mechanism for fractal scaling in stream tracer concentrations, Journal of Hydrology, 254,82-101. 142 Konovalov, V . G . , and A.S. Shchetinnicov, 1994. Evolution of glaciation in the Parniro-Alai Mountains and its effect on river run-off, Journal of Glaciology, 40(134), 149-157. Krimmel, R.M. , and W.V. Tangborn, 1974. South Cascade Glacier: the moderating effect of glaciers on runoff. In Washichek, J.N (ed.), Proceedings of the Western Snow Conference, 42nd Annual Meeting, Colorado State University, Fort Collins, 9-13. Lafreniere, M . , and M . Sharp, 2003. Wavelet analysis of inter-annual variability in the runoff regimes of glacial and nival stream catchments, Bow Lake, Alberta, Hydrological Processes, 17(6), 1093-1118. Lall, U. , and M . Mann, 1995. The Great Salt Lake: a barometer of low-frequency climatic variability, Water Resources Research, 31, 2503-2515. Latif, M . , D. Anderson, T. Barnett, M . Cane, R. Kleeman, A. Leetmaa, J. O'Brien, A. Rosati, and E. Schneider, 1998. A review of the predictability and prediction of ENSO, Journal of Geophysical Research, 103, 14375-14393. Law, D.J., C.B. Marlow, J.C. Mosley, S. Custer, P. Hook, and B. Leinard, 2000. Water table dynamics and soil texture of three riparian plant communities, Northwest Science, 74(3), 234-241. Leith, R.M., and P.H. Whitfield, 1998. Evidence of climate change effects on the hydrology of streams in south-central British Columbia, Canadian Water Resources Journal, 23, 219-230. Leith, R.M., and P.H. Whitfield, 2000. Some effects of urbanization on streamflow records in a small watershed in the Lower Fraser Valley, B.C., Northwest Science, 74(1), 69-75. Lettenmaier, D.P., 1995. Stochastic modeling of precipitation with applications to climate model downscaling. In von Storch, H. and A. Navarra (eds.), Analysis of Climate Variability: Applications of Statistical Techniques, Springer, Berlin, 197-212. 143 Lettenmaier, D.P., E.F. Wood, and J.R. Wallis, 1994. Hydro-climatological trends in the continental United States, 1948-88, Journal of Climate, 7, 586-607. Lindsey, C C , K. Patalas, R.A. Bodaly, and C P . Archibald, 1981. Glaciation and the physical, chemical, and biological limnology of Yukon lakes, Canadian Technical Report of Fisheries and Aquatic Sciences, 966, i-37. Lindsey, C C , and J.D. McPhail, 1986. Zoogeography of fishes of the Yukon and Mackenzie basins. In Hocutt, C H . and E.O. Wiley (eds.), The Zoogeography of North American Freshwater Fishes, Wiley, New York, 639-674. Lins., H.F., 1985. Streamflow variability in the United States: 1931-78, Journal of Climate and Applied Meteorology, 24,463-471. Lins, H.F., 1997. Regional streamflow regimes and hydroclirnatology of the United States, Water Resources Research, 33(7), 1655-1667. Lins, H.F., and P.J. Michaels, 1994. Increasing U.S. streamflow linked to greenhouse forcing, EOS, Transactions of the American Geophysical Union, 75(25), 281, 284-285. Lorenz, E.N., 1956. Empirical Orthogonal Functions and Statistical Weather Prediction, Technical Report, Statistical Forecast Project Report 1, Department of Meteorology, Massachusetts Institute of Technology. MIT, Boston, 49 pp. Manly, B.F.J., 1997. Randomization, Bootstrap, and Monte Carlo Methods in Biology, 2nd Ed. Chapman and Hall, London, 399 pp. Mantua, N.J., S.R. Hare, Y. Zhang, J.M. Wallace, and R.C. Francis, 1997. A Pacific interdecadal climate oscillation with impacts on salmon production, Bulletin of the American Meteorological Society, 78(6), 1069-1079. 144 McPhail, J.D., and C C . Lindsey, 1970. Freshwater fishes of northwestern Canada and Alaska, Bulletin of the Fisheries Research Board of Canada, 173, 1-381. Meier, M.F., and W.V. Tangborn. 1961. Distinctive Characteristics of Glacier Runoff, USGS Professional Paper 424B, B14-B16. U.S. Geological Survey, Washington DC, 2 pp. Meier, M.F., 1969. Glaciers and water supply, Journal of the American Water Works Association, 61(1), 8-12. Mekis, E . , and W.D. Hogg, 1999. Rehabilitation and analysis of Canadian daily precipitation time series, Atmosphere-Ocean, 37(1), 53-85. Meko, D.M., and C.W. Stockton, 1984. Secular variations in streamflow in the western United States, Journal of Climate and Applied Meteorology, 23, 889-897. Michel, F.A., and R.O. van Everdingen, 1994. Changes in hydrogeologic regimes in permafrost regions due to climatic change, Permafrost and Periglacial Processes, 5, 191-195. Milner, A . M . , and G.E. Petts, 1994. Glacial rivers: physical habitat and ecology, Freshwater Biology, 32,295-307. Milner, A . M . , J.E. Brittain, E . Castella, and G.E. Petts, 2001. Trends of macroinvertebrate community structure in glacier-fed rivers in relation to environmental conditions: a synthesis, Freshwater Biology, 46, 1833-1847. Moore, R.D., 1992. The influence of glacial cover on the variability of annual runoff, Coast Mountains, British Columbia, Canada, Canadian Water Resources Journal, 17, 101-109. 145 Moore, R.D., 1993. Application of a conceptual streamflow model in a glacierized drainage basin, Journal ofHydrology, 150, 151-168. Moore, R.D., and M.N. Demuth, 2001. Mass balance and streamflow variability at Place Glacier, Canada, in response to recent climate fluctuations, Hydrological Processes, 15(18), 3473-3486. Moore, R.D., A.S. Hamilton, and J. Scibek, 2002. Winter streamflow variability, Yukon Territory, Canada, Hydrological Processes, 16(4), 763-778. Mote, P.W., 2003. Trends in temperature and precipitation in the Pacific Northwest during the twentieth century, Northwest Science, 77(4), 271-282. Mote, P.W., E.A. Parson, A.F. Hamlet, W.S. Keeton, D. Lettenmaier, N . Mantua, E.L. Miles, D.W. Peterson, D.L. Peterson, R. Slaughter, and A.K. Snover, 2003. Preparing for climatic change: the water, salmon, and forests of the Pacific Northwest, Climatic Change, 61,45-88. Mysak, L.A. , 1986. El Nifio, interannual variability and fisheries in the northeast Pacific Ocean, Canadian Journal of Fisheries and Aquatic Sciences, 43, 464-497. Neal, E .G. , M.T. Walter, and C. Coffeen, 2002. Linking the Pacific decadal oscillation to seasonal discharge patterns in southeast Alaska, Journal of Hydrology, 263, 188-197. North, G.R., T.L. Bell, R.F. Cahalan, and F.J. Moeng, 1982. Sampling errors in the estimation of empirical orthogonal functions, Monthly Weather Review, 110, 699-706. Owen, E.B., 1967. Northern hydrogeological region. In Brown, I.C. (ed.), Groundwater in Canada, GSC Economic Geology Report No. 24, Geological Survey of Canada, Ottawa, 173-194. 146 Pahlke K.A. , and P. Etherton, 2001. Abundance of the Chinook Salmon Escapement on the Alsek River, 2000, Fishery Data Series No. 01-30. Alaska Department of Fish and Game, Anchorage, 29 pp. Pelletier, P.M., 1990. A review of techniques used by Canada and other northern countries for measurement and computation of streamflow under ice conditions, Nordic Hydrology, 21,317-339. Pelto, M.S., and J. RiedeL 2001. Spatial and temporal variations in annual balance of North Cascade Glaciers, Washington 1984-2001, Hydrological Processes, 15, 3461-3472. Phillips, D., 1990. The Climates of Canada. Environment Canada, Ottawa, 176 pp. Piechota, T.C. , and J.A. Dracup, 1999. Long-range streamflow forecasting using El Nifto-Southern Oscillation Indicators, Journal of Hydrologic Engineering, 4,144-151. Pitman, E.J.G., 1937. Significance tests which may be applied to samples from any populations, II, the correlation coefficient test, Journal of the Royal Statistical Society, B4, 225-232. Power, G., R. Cunjak, J. Flannagan, and C. Katopodis, 1993. Biological effects of river ice. In Prowse, T.D. and N.C. Gridley (eds.), Environmental Aspects of River Ice, National Hydrology Research Institute, Saskatoon, 97-119. Preisendorfer, R.W., 1988. Principal Component Analysis in Meteorology and Oceanography. Elsevier, Amsterdam, 426 pp. Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, 1992. Numerical Recipes in Fortran 77, 2nd Ed. Cambridge University Press, Cambridge, 933 pp. Quick, M.C. , 1995. The UBC watershed model. In Singh, V.P. (ed.), Computer Models of Watershed Hydrology, Water Resources Publications, Littleton, 233-280. 147 Redmond, K.T. , and R.W. Koch, 1991. Surface climate and streamflow variability in the western United States and their relationship to large-scale circulation indices, Water Resources Research, 27, 2381-2399. Reynolds, J.B., 1997. Ecology of overwintering fishes in Alaskan freshwaters. In Milner, A . M . and M.W. Oswood (eds.), Freshwaters of Alaska: Ecological Syntheses, Springer-Verlag, New York, 281-302. Richman, M.B., 1986. Rotation of principal components, Journal of Climatology, 6, 293-335. Roden, G.I., 1989. Analysis and interpretation of long-term climatic variability along the west coast of North America. In Peterson, D.H. (ed.), Aspects of Climate Variability in the Pacific and the Western Americas, Geophysical Monograph 55, American Geophysical Union, Washington DC, 93-111. Rogerson, R.J., 1986. Mass balance of four cirque glaciers in the Torngat Mountains of northern Labrador, Canada, Journal of Glaciology, 32(111), 208-218. Roots, E.F., 1989. Climate change: high-latitude regions, Climatic Change, 15, 223-253. Roper, B.B., and D.L. Scarnecchia, 2001. Patterns of diversity, density, and biomass of ectothermic vertebrates in ten small streams along a North American river continuum Northwest Science, 75(2), 168-175. Rouse, W:R., M.S.V. Douglass, R.E. Hecky, A.E. Hershey, G.W. Kling, L . Lesack, P. Marsh, M . McDonald, B.J. Nicholson, N.T. Roulet, and J.P. Smol, 1997. Effects of climate change on the freshwaters of arctic and subarctic North America, Hydrological Processes, 11, 873-902. 148 Salas, J.D., J.W. Delleur, V. Yevjevich, and W.L. Lane, 1980. Applied Modeling of Hydrologic Time Series. Water Resources Publications, Littleton, 484 pp. Scrimgeour, G.J., T.D. Prowse, J.M. Culp, and P.A. Chambers, 1994. Ecological effects of river ice break-up: a review and perspective, Freshwater Biology, 32,261-275. Sen, P.K., 1968. Estimates of the regression coefficient based on Kendall's tau, Journal of the American Statistics Association, 63, 1379-1389. Shabbar, A., B. Bonsai, and M . Khandekar, 1997. Canadian precipitation patterns associated with the Southern Oscillation, Journal of Climate, 10, 3016-3027. Shabbar, A. , and M . Khandekar, 1996. The impact of El Nino-Southern Oscillation on the temperature field over Canada, Atmosphere-Ocean, 34, 401-416. Showstack, R., 2002. Mountain regions in peril, EOS, Transactions of the American Geophysical Union, 83(6), 54. Simpson, H.J., M.A. Cane, A .L . Herczeg, S.E. Zebiak, and J.H. Simpson, 1993. Annual river discharge in southeastern Australia related to El Nifio-Southern Oscillation forecasts of sea surface temperatures, Water Resources Research, 29, 3671-3680. Singh, P., and N. Kumar, 1997. Impact assessment of climate change on the hydrological response of a snow and glacier melt runoff dominated Himalayan river, Journal of Hydrology, 193, 316-350. Singh, P., and V P . Singh, 2001. Snow and Glacier Hydrology. Kluwer, Dordrecht, 742 pp. Sivapalan, M . , C. Jothityangkoon, and M . Menabde, 2002. Linearity and nonlinearity of basin response as a function of scale: discussion of alternative definitions, Water Resources Research, 38(2), doi: 10.1029/2001WR000482,4-1-4-5. 149 Sk0ien, J.O., G. Bloschl, and A.W. Western, 2003. Characteristic space scales and timescales in hydrology, Water Resources Research, 39(10), doi: 10.1029/2002WR001736, 11-1-11-19. Smith, B.P.G., D .M. Hannah, A . M . Gurnell, and G.E. Petts, 2001. A hydrogeomorphological context for ecological research on alpine glacial rivers, Freshwater Biology, 46, 1579-1596. Smolders, A.J.P., M.A. Guerrero Hiza, G. van der Velde, and J.G.M. Roelofs, 2002. Dynamics of discharge, sediment transport, heavy metal pollution, and Sabalo (Prochilodus lineatus) catches in the lower Pilcomayo River (Bolivia), River Research and Applications, 18,415-427. Spence, C , 2002. Streamflow variability (1965 to 1998) in five Northwest Territories and Nunavut rivers, Canadian Water Resources Journal, 27(2), 135-154. Tang, Y., and W.W. Hsieh, 2002. Hybrid coupled models of the tropical Pacific: II, ENSO prediction, Climate Dynamics, 19, 343-353. Tangborn, W.V., 1984. Prediction of glacier derived runoff for hydroelectric development, Geografiska Annaler, 66A(3), 257-265. Telford, W.M., L.P. Geldart, and R.E. Sheriff, 1990. Applied Geophysics, 2nd Ed. Cambridge University Press, Cambridge, 770 pp. Tempelman-Kluit, D., 1980. Evolution of physiography and drainage in southern Yukon, Canadian Journal of Earth Sciences, 17, 1189-1203. Theil, H. , 1950. A rank-invariant method of linear and polynomial regression analysis, I, II, and III, Proceedings, Koninklijke Nederlandse Akademie van Wetenschappen, 53, 386-92, 521-525, 1397-1412. 150 Thompson, D.W.J., and J.M. Wallace, 1998. The Arctic Oscillation signature in the wintertime geopotential height and temperature fields, Geophysical Research Letters, 25, 1297-1300. Thompson, D.W.J., and J.M Wallace, 2000. Annual modes in the extratropical circulation, part I: month-to-month variability, Journal of Climate, 13, 1000-1016. Thompson, D.W.J., and J.M. Wallace, 2001. Regional climate impacts of the Northern Hemisphere annular mode, Science, 293, 85-89. Trenberth, K .E . , 1977. Relationships between inflow to Clutha lakes, broad-scale atmospheric circulation parameters, and rainfall, New Zealand Journal of Science, 20, 63-71. Trenberth, K.E . , and T.J. Hoar, 1996. The 1990-1995 El Nino-Southern Oscillation event: longest on record, Geophysical Research Letters, 23, 57-60. Tsonis, A.A. , A . G . Hunt, and J.B. Eisner, 2003. On the relation between ENSO and global climate change, Meteorology and Atmospheric Physics, 84, 229-242. Ulrych, T., and S. Kaplan, 2002. Blind Deconvolution via Independent Component Analysis. Consortium for the Development of Specialized Seismic Techniques, Vancouver, 15 pp. van der Kwaak, J.E., 1999. Numerical Simulations of Flow and Chemical Transport in Integrated Surface-Subsurface Hydrologic Systems. PhD thesis, University of Waterloo, 217 pp. Vannote, R.L., G.W. Minshall, K.W. Curnmins, J.R. Sedell, and C.E. Cushing, 1980. The river continuum concept, Canadian Journal of Fisheries and Aquatic Sciences, 37, 130-137. 151 Viles, H.A., and A.S. Goudie, 2003. Interannual, decadal and multidecadal scale climatic variability and geomorphology, Earth-Science Reviews, 61, 105-131. Vincent, L .A, and D.W. Gullett, 1999. Canadian historical and homogeneous temperature datasets for climate change analyses, International Journal of Climatology, 19, 1375-1388. von Storch, H. , 1995. Misuses of statistical analysis in climate research. In von Storch, H. and A. Navarra (eds.), Analysis of Climate Variability: Applications of Statistical Techniques, Springer-Verlag, Berlin, 11-26. von Storch, H. , and F.W. Zwiers, 1999. Statistical Analysis in Climate Research. Cambridge University Press, Cambridge, 484 pp. WaliL H.E. , D.B. Fraser, R.C. Harvey, and J.B. Maxwell, 1987. Climate of Yukon, Climatological Studies Number 40. Atmospheric Environment Service, Ottawa, 323 pp. Walker, J.F., and D. Wang, 1997. Measurement of flow under ice covers in North America, Journal of Hydraulic Engineering, 123( 11), 1037-1040. Wallace, J.M., and D.S. Gutzler, 1981. Teleconnections in the geopotential height field during the Northern Hemisphere winter, Monthly Weather Review, 109, 784-812. Wallace, J.M., E . M . Rasmusson, T.P. Mitchell, V . E . Kousky, E.S. Sarachik, and H. von Storch, 1998. On the structure and evolution of ENSO-related climate variability in the tropical Pacific: lessons from TOGA, Journal of Geophysical Research, 103, 14241-14259. Wallace, J.M., and D.W.J. Thompson, 2002. Annular modes and climate prediction, Physics Today, 55, 28-33. 152 Walters, R.A., and M.F. Meier, 1989. Variability of glacier mass balances in western North America. In Peterson, D.H. (ed.), Aspects of Climate Variability in the Pacific and the Western Americas, Geophysical Monograph 55, American Geophysical Union, Washington DC, 375-397. Wang, X . L . , and V.R. Swail, 2001. Changes of extreme wave heights in northern hemisphere oceans and related atmospheric circulation regimes, Journal of Climate, 14, 2204-2221. Ward, J.V., 1994. Ecology of alpine streams, Freshwater Biology, 32,277-294. Whitfield, P.H., and A.J. Cannon, 2000a. Recent climate moderated shifts in Yukon hydrology. In Kane, D.L. (ed.), AWRA Spring Specialty Conference Proceedings, American Water Resources Association, Anchorage, 257-262. Whitfield, P.H., and A.J. Cannon, 2000b. Recent variations in climate and hydrology in Canada, Canadian Water Resources Journal, 25, 19-65. Whitfield, P.H., 2001. Linked hydrologic and climate variations in British Columbia and Yukon, Environmental Monitoring and Assessment, 67, 217-238. Wicket, W.P., 1958. Review of certain environmental factors affecting the production of pink and chum salmon, Journal of the Fisheries Research Board of Canada, 15, 1103-1126. Woo, M . , 1990. Consequences of climate change for hydrology in permafrost zones, Journal of Cold Regions Engineering, 4, 15-20. Woo, M . , and R. Thorne, 2003. Comment on 'Detection of hydrologic trends and variability' by Burn, D.H. and Hag Elnur, M.A. , 2002, Journal of Hydrology, 111, 150-160. 153 Wooding, F.H., 1997. Lake, River, and Sea-Run Fishes of Canada. Harbour Pubkshing, Madeira Park, 303 pp. Wu, A., and W.W. Hsieh, 2004. The nonlinear Northern Hemisphere winter atmospheric response to ENSO, Geophysical Research Letters, 31, doi: 10.1029/2003GL018885, 1-1-1-4. Yue, S., P. Pilon, and G. Cavadias, 2002. Power of the Mann-Kendall and Spearman's rho tests for detecting monotonic trends in hydrological series, Journal of Hydrology, 259, 254-271. Zaltsberg, E . , 1990. Potential changes in mean annual runoff from a small watershed in Manitoba due to possible climatic changes, Canadian Water Resources Journal, 15(4), 333-344. Zhang, X., K.D. Harvey, W.D. Hogg, and T.R. Yuzyk, 2001. Trends in Canadian streamflow, Water Resources Research, 37(4), 987-998. Zhang X. , L.A. Vincent, W.D. Hogg, and A. Niitsoo, 2000. Temperature and precipitation trends in Canada during the 20th century, Atmosphere-Ocean, 38, 395-429. Zhang, Y. , J.M. Wallace, and D.S. Battisti, 1997. ENSO-like interdecadal variability: 1900-n, Journal of Climate, 10, 1004-1020. Zuming, L . , 1982. A study on the variation coefficient of annual runoff of the rivers in northwest China, In Lang, H. and A. Musy (eds.), Hydrological aspects of alpine and high-mountain areas, IAHS publication 138, International Association of Hydrological Sciences, Birmingham 285-294. 154 APPENDIX 1 NOTES ON EOF ANALYSIS Al.l General Say we have time series of some quantity of interest, y, measured at m = l,M locations at / = IJV times. This quantity could be streamflow hydrographs from different rivers, for example, or the price of a pint of beer at various bars over time. It will be assumed here that at each location, the time series has been centred - that is, it has had its temporal mean removed, so that the y values are departures from the local long-term mean. We will start by viewing these data as a vector, y(tf), them - \,M entries of which correspond to the m = l,M locations. This is in fact a suite of N vectors, each corresponding to a given time, tt. The underscore symbol is used here to denote vector quantities. We wish to expand y_(t^ into a linear combination of maps, and time series, ka(t^: M with: A t K e 7 •y(ti) (22) where the T superscript denotes the transpose and the symbol, •, denotes the vector dot product. Each k refers to a different "mode," and k= 1,M, so that there are as many modes as locations, and each mode has its own vector *g and time series maftj). Like yft$ for a particular th & for a given value of A: is Af-dimensional, and the individual entries of any one of these k = 1,M vectors correspond to the m = 1,M locations. Empirical orthogonal function analysis consists of solving (21) and (22) for ae and va(tj), subject to the following constraints: (i) all the k = \,M k§. must be orthogonal (i.e., Mmutually perpendicular vectors in an M-dimensional hyperspace); (ii) each kg. must be of unit magnitude (usually; other 155 scaling protocols are possible); (iii) all the k = \,Mkaftj) must be mutually uncorrelated (i.e., time-independent); and (iv) the linear combination of *e and ka(tj) for k = 1 must explain as much of the total variance of all the y values as possible, with higher modes (linear combinations of ne and ^ afti) for successively larger values of k) accounting for progressively lower proportions of the total spatiotemporal variance. Together, (i) and (ii) dictate that the k = 1 ,M ice must be orthonormal. The general nature or meaning of this type of expansion can be viewed in a number of (ultimately equivalent) ways. One is that it constitutes the decomposition of a complex dataset into a number of orthogonal basis functions, In comparison to spectral analysis, for example, the suite of vectors, is analogous to a suite of sinusoids of varying frequency; and the eigenspectrum, &i for all k = l,M (to be discussed in due course), is analogous to the Fourier coefficients for each of those sinusoids. It may also be viewed as a transformation of the coordinates of the problem to a more natural geometry: *g over all k = 1,A/ consist of orthogonal unit vectors defining a new coordinate system and for each k, the time: series kd(tj) represents the projection of the original data onto each of those new coordinate axes. The fundamental motivation for EOF analysis, however, is one of reducing dataset dimensionality. There are as many modes as there are locations. However, as noted above, these modes are constructed such that the leading mode (the linear combination corresponding to k = 1) accounts for the greatest amount of the variability in the data possible, and the amount of variability explained by each mode decreases with increasing k. As a result, the series of (21) may be truncated at some k « M, without losing much of the information content of the original data. Given the type of product yielded by the expansion (for a given mode, a simple spatial map dynamically weighted by a single time series), this truncation can enormously but justifiably simplify subsequent data interpretation, analysis, or modelling. Recall from previous chapters that the leading mode alone is often sufficient to capture the dominant characteristics of spatiotemporal variability, even in very large and complex datasets. This is, of course, subject to the caveat that there must be some kind of interesting overall spatiotemporal pattern lurking within the data, waiting to be captured; and that a linear decomposition of the type in (21) and (22) is capable of capturing it. 156 Nevertheless, in practise, the method has proven itself able to extract dominant modes of behaviour in an extremely wide array of problems. Moreover, application of linear EOF analysis to data representing the output of distinctly nonlinear processes has been found to yield sensible and useful results in many cases (e.g., identification of basic ENSO patterns within tropical Pacific SST and SLP data), irrespective of violation of the underlying technical assumptions. The method appears, therefore, to be reasonably robust. The mechanics of implementation to a given dataset, y, are as follows. Note that the vector notation used above is more convenient for illustrating the fundamental premise of EOF analysis, whereas matrix notation is more convenient for describing (and, in practise, implementing) the analysis procedure, but the two are fully consistent. For example, the time, space, and modal indices (/ = \J^,m- l,Af, and k = \Jrf, respectively) are identical; and a column of Y (see below) is y(t^ for a given /, a column of E is & for a given k, and a row of A is kaftj) for a given k. We begin by assembling the observations at i - l^V times from m = \,M locations into a M-row by A -^column data matrix, Y, each row of which contains the observed time series from a particular location: m = I, i = 1 later times -> different locations m = \,i = N Y (23) m :- M,i = 1 m = M,i = N The order in which the locations are listed in the matrix is irrelevant, but must be kept track of to interpret the spatial patterns that will be produced. If, as discussed in the introductory paragraph, the time series have been centred, i.e., the mean value of the row has been subtracted from each of the observations in that row, with the procedure repeated separately for each row of Y, then the data covariance matrix, C, may be obtained as: C = —YYT N (24) 157 C is a square symmetric M by M matrix, the elements of which provide the covariance between each pair of time series in Y; the diagonals consist of the zero-lag autocovariances, i.e., simply the variance of each time series. If in addition to having been centred, each of the time series assembled into Y has also been individually normalized by its own standard deviation, so that all the time series are standardized to zero mean and unit variance, then C becomes the correlation matrix, where the diagonals are by definition unity and the off-diagonal entries provide the linear Pearson product-moment correlation coefficients between each time series and every other time series in Y. The basis functions, or map patterns, consist of the eigenvectors of C. These are contained in a square M by M matrix, denoted E, each column of which contains one eigenvector: m = 1, k = 1 higher EOFs different locations I m = \,k = M (25) m = M, k = 1 m = M, k = M Each column thus corresponds to a different EOF mode. The elements of each eigenvector (i.e., the eigenvector for given mode) correspond to different locations in space and are arranged in precisely the same order as a given column of Y. The corresponding eigenvalues are organized as a column vector: 2=> k = \ higher EOFs (26) 158 Both E and X are arranged such that mode k = 1 explains the most variance (as discussed above) and is associated with the largest eigenvalue, with k increasing and eigenvalue magnitude decreasing downward through the matrices; the k values are consistent between E and X, so that the A* eigenvalue corresponds to the A* eigenvector. The first column of E, therefore, contains the leading-mode eigenvector and the first row of X contains its corresponding eigenvalue. As discussed in previous chapters, eigenanalysis was performed here using subroutines from Press et al. (1992). The precise manner whereby eigenvectors arise in EOF analysis varies between full mathematical developments (compare, for example, the derivations of Lorenz, 1956 versus von Storch and Zwiers, 1999). The various methods by which EOF analysis can be derived (or rediscovered, as the case may be) will not be discussed here. However, it is reassuring to note that, as a general truth of matrix algebra, the eigenvectors of a real symmetric matrix (such as Q are orthogonal, provided that the eigenvalues are not equal or nearly equal (degenerate); and also that in the present case, C possesses a total of M eigenvectors. Thus, it is clear that the eigenvectors of C can serve as orthogonal basis functions in the way prescribed by (21). The descending series of eigenvalues in (26) are known as the eigenspectrum, and while they are not directly used in the expansion, they are a crucial component of the EOF analysis result. Specifically, the proportion of the total spatiotemporal variance of the original data explained by the map pattern and time series corresponding to a particular mode, k = v, is given by: %var = ^ — 1 0 0 % (27) M k=l The: values taken on by *X are key to assessing the physical significance of the various EOF modes, as discussed in detail in previous chapters. Note that while efforts to quantify modal 159 significance must address distributional issues, the EOF analysis technique itself, as described in this appendix, makes no distributional assumptions (e.g., Lins, 1997). The final step in the EOF decomposition is to evaluate the M-row by N-column principal component time series matrix, A: A = E Y (28) which is arranged in the following manner: A k = 1, / = 1 later times higher EOFs I . . m = U = JV~ k = M,i = \ m = M,i - N (29) Thus, each row of A provides the principal component time series for a given mode. In particular, the leading-mode unrotated principal component (LMUPC, in previous chapters) time series is given by the first (top) row of A. A1.2 Gridding Karl et al. (1982) suggested that application of EOF analysis to spatial point measurements of a continuous field, such as ocean temperature or atmospheric pressure, can result in a distortion of the eigenvector patterns such that areas of dense spatial sampling receive disproportionately high loadings; interpolation of station data to a regularly-spaced grid was suggested as a solution. The implicit assumption is that the interpolated data are reasonably reliable, which in turn requires adequate spatial density in the station data being gridded and an interpolation scheme appropriate to the problem. Both conditions may, in certain applications, be difficult to satisfy in practise. Nevertheless, many if not most 160 climatological applications of EOF analysis do indeed employed gridded data, although this may also reflect issues of convenience: certain commonly-used datasets, such as that of the NCEP reanalysis, are made publicly available in gridded form. We are concerned here, however, with the manner in which individual hydrometric stations respond to climatic forcing, and interpolating the hydrometric station data to a grid would, therefore, be obviously inappropriate. More generally, such gridding procedures are likely unsuitable for streamflow data regardless of the context of the study. River discharge records are not point measurements of a continuous field, but spatially integrated measurements of a discontinuous field, and smearing such data over the significant physical discontinuity of the watershed drainage divide through interpolation is likely to cause far more problems than it could conceivably solve (Lins, 1997). Likewise, there are issues of physical reasonableness: spatial interpolation of streamflow records to a regular grid will certainly place river flow values at locations where there is no river. More broadly still, the very concept of spatial interpolation has no meaning in most applications of the parent technique, principal components analysis (see below), and thus certainly cannot be viewed as a standard and requisite data processing procedure. A1.3 Rotation and Nonlinear Approaches Given the fairly considerable attention devoted to rotation in the past, it is summarily addressed in this appendix. The basic notion behind rotation is to manipulate the results of a standard EOF analysis such that data clusters (i.e., correlation patterns) in our M-dimensional data space are more effectively pierced by the rotated axes. The mechanics of the procedure are not, in general, trivial. The seminal work on rotation in a climatological context is that of Richman (1986), who advocated rotation as a means of addressing four problems he claimed to exist with unrotated EOF analysis results: domain shape dependence, which consists of a predictable sequence of eigenvectors for successive modes, determined by the shape of the spatial data domain; subdomain instability, whereby one obtains different overall spatial patterns if the data are split spatially prior to analysis; eigenvector ^determinancy for degenerate eigenvalues, so that if neighbouring modes possess similar eigenvalues, the directions of the corresponding eigenvectors are arbitrary, 161 and the EOF map patterns and PC time series are meaningless; and a lack of physical interpretability. The drawbacks to rotation, however, are even more numerous (e.g., Jolliffe, 1987; von Storch and Zwiers, 1999). At least 19 different rotation procedures exist; the choice of which of these to use is largely arbitrary and subjective, yet often highly pertinent to the result. Such procedures also require the user to select a certain number of modes to rotate, this choice again being somewhat arbitrary, but typically leading to different results. The alleged augmentation of physical interpretability afforded by rotation is another somewhat subjective criterion and relevant only when the unrotated results are uninterpretable; moreover, rotation can render a physically reasonable unrotated solution uninterpretable. In general, the first k rotated modes explain less variance than the first k unrotated modes. Rotation also involves a loss of modal independence, which is one of the most basic and desirable properties of EOF analysis: depending on the rotation technique used, either eigenvector orthogonality or principal component time series independence between modes is sacrificed. The domain shape dependence and subdomain instability issues raised by Richman (1986) were more-or-less dismissed by Jolliffe (1987). More broadly, Jolliffe (1987) suggested that rotation constitutes an attempt to force an EOF analysis to do something other than what it was intended to do: capture as much of the correlation structure of the data as possible in a series of independent modes of progressively decreasing significance. Without wishing malice upon Richman's (1986) work, to my mind it possesses a few other difficulties: Richman (1986) appeared to assume densely-populated (e.g., finely-gridded), rectilinear spatial data domains in his analysis of domain shape dependence and subdomain instability, which does not seem relevant to many applications, such as those considered in this dissertation; virtually no attention was paid to the relative qualities of PC time series corresponding to rotated and unrotated solutions, a considerable limitation given that these can be of equal or greater practical importance to spatial patterns, depending on the specific physical application; and while the results of standard EOF analyses, whatever their limitations might or might not be, are intuitively attractive and accessible, it is considerably less clear what a rotated solution represents, particularly since each rotation technique does something different. 162 Given that nearly two decades have passed since Richman's (1986) work, it is instructive to examine what consensus the hydroclirnatological community may have reached, either explicitly or implicitly, with regard to rotation. A great number of unrotated EOF analyses have yielded results which do not appear plagued by Richman's (1986) modal pattern artefacts or physical uninterpretabihty. The use of unrotated EOF analysis in the determination of widely used climate indices, such as the PDO (Mantua et al., 1997) and A O (Thompson and Wallace, 1998) indices, has raised little or no controversy. Applications of EOF analysis to spatiotemporal streamflow variability are relatively rare. Hisdal and Tveito (1993) and Cullen et al. (2002) did not use rotation, whereas Lins and Michaels (1994) did; more helpfully for our immediate purposes, Lins (1985, 1997) considered both unrotated and rotated (varimax and promax, Lins, 1985; varimax only, Lins, 1997) solutions, and found that rotation did not materially affect, or improve, interpretation of streamflow EOF analysis results. Finally, a recent and fairly authoritative text of statistical climatology by von Storch and Zwiers (1999) indicated that rotation procedures are not to be viewed as a standard post-processing procedure, and suggested that they be applied only if the eigenvalues are degenerate. As discussed in preceding chapters, this was not found to be the case in the analyses presented here for the first and second modes, which are all that are relevant when general attention is restricted to the leading mode. Note also that interpretations of the unrotated EOF analysis results in this thesis were checked against methodologically independent lines of evidence, and were found to be consistent with findings from graphical and distributional analyses (chapter three) or composite analysis (chapter six). An additional potential reason why rotation has apparently fallen largely into disuse is very different and worth a brief discussion in its own right: the recent emergence of nonlinear EOF analysis. A number of nonlinear EOF-like procedures have recently been developed and/or applied in the earth system sciences (e.g., Huang et al, 1998; Ulrych and Kaplan, 2002), but in a climatological context, the dominant approach has been artificial neural network (ANN)-based nonlinear EOF analysis. A review of ANN-based nonlinear analogs to a variety of linear, eigenanalysis-based multivariate statistical techniques is provided by Hsieh (2004). The fundamental advantage of nonlinear EOF analysis over conventional, rotated or unrotated, EOF analysis is that it can properly accommodate a far broader range 163 of process dynamics, and thus may often explain a considerably greater proportion of total dataset variance using fewer modes. Consequently, if one is unhappy with conventional, linear, unrotated EOF analysis, it is generally more sensible now to step up to nonlinear methods rather than to an ad hoc rotation of the linear solution. Unfortunately, nonlinear EOF analysis also had its drawbacks. A N N methods are reknown for both their power and their inscrutability: nonlinear EOF analysis can in some cases be statistically superior, but it may not be entirely clear what the procedure has done or how it did it. In addition, as noted previously, linear methods often prove serviceable when applied to data representing nonlinear phenomena, such as ENSO impacts upon tropical Pacific surface climate. Finally, copious amounts of data (in particular, long records) are generally necessary for adequate ANN training and testing; the data requirements for the method are thus likely well in excess of, for example, streamflow data availability in the current study area. A1.4 EOF Analysis versus PCA A final note on nomenclature is warranted. The terms empirical orthogonal function analysis and principal component analysis (PCA) are typically used interchangeably. This is legitimate, since the two are fundamentally mathematically identical. However, application of the term EOF analysis, which presumably arose from the title of Lorenz's (1956) seminal work in multivariate statistical meteorology, is restricted to the climatological and oceanographic communities, and that term is adopted here for its specificity. PCA results can be used and processed in a number of ways; for example, fisheries biologists often use PCA as an ordination technique, and the end product is often nearly unrecognizable in the context of maps and time series discussed here. In contrast, EOF analysis typically refers to particular space-time problems of the types addressed in this thesis, and implies a certain perspective, drawn from long climatological experience, regarding how data should be processed, how the method should be implemented, and how the results are presented and interpreted. While EOF applications to streamflow data have been relatively limited to date, and meteorological or oceanographic rules for its implementation may not in all cases be fully appropriate to questions of riparian science (see, for example, above and Lins, 1997), overall the particular approaches used by those communities appear to serve as the most relevant guide available for application of the technique to the space-time fluvial 164 hydroclirnatological issues considered here. Adoption of the term EOF analysis thus serves as a useful flag as to the perspectives and protocols that were drawn from during implementation of the method. Note finally that factor analysis, although similar in concept to PCA/EOFA and sometimes used as an additional synonym, is in fact an entirely distinct technique (Jolliffe, 1993). 165 APPENDIX 2 MARKOV CHAINS, DESERIALIZATION, AND TRENDS A2.1 Introduction A2.1.1 Nonparametric Statistical Detection of Monotonic Trend Rigorous statistical identification of long-term trends in hydrological time series is relevant to a multitude of surface water and groundwater applications, and is of particular importance to assessing river discharge records for evidence of climate change signals. Trend identification is usually performed in this context using the nonparametric Mann-Kendall test, which accepts or rejects the null hypothesis of randomness against the alternative of monotonic trend. The advantage of the Mann-Kendall test is that it makes as few assumptions as possible: the potential trend may be linear or nonlinear, and no presumptions are made as to the underlying statistical distribution. This technique is often used in conjunction with the rank-based, outlier-resistant slope estimator of Theil (1950) and Sen (1968), which was introduced and discussed in chapter six. Application of the Mann-Kendall test to a given time series proceeds as follows (e.g., Burn, 1994). The test statistic is calculated as: where JC, for i = 1,JV are the sequential data values of the time series, and the sign function is defined as: (30) 1, 7T>1 sgnfa) =0 , n = 0 -1, n<0 (31) 166 A corrected test statistic, W, is very nearly Gaussian for N > 10, and has a mean of zero under the null hypothesis, so it may therefore be transformed to a z-score to facilitate the hypothesis test: W'=W-sgn(w) (32) and: W (33) z = sw The above techniques are typically applied to an annual time series consisting of some metric derived from daily streamflow data. Such measures include but are not limited to total annual discharge, annual mean daily discharge, calendar day number of peak flow, mean discharge for each month, and freshet timing. Analysis of annual rather than daily time series involves a substantial and undesirable reduction of dataset size. This is nevertheless a generally preferable approach because the fluvial impacts of climate changes, be they natural or anthropogenic, may involve timing shifts or changes in extreme events, or be seasonal in nature; responses of this type, if present, are typically more apparent in the derived annual metrics than in daily flow data. A2.1.2 Markov-Chain Processes and Deserialization It has long been recognized that the Mann-Kendall test is not robust to serial correlation (e.g., Hirsch et al, 1982; Hirsch and Slack, 1984), with serial correlation specifically taken to refer to persistence in the form of AR(1) noise: where JC, is the present (/*) measurement of the hydrological quantity under consideration, is the previous (/-1th) measurement, st is white (serially independent and normalfy-X; - cx i-l (34) 167 distributed) random noise with mean, ju, and standard deviation, a, and c is the lag-1 serial correlation coefficient, given by: N-\ E f o - * X * < + i - * ) where /V is again the number of data. The general presence of such first-order Markov-chain processes in daily discharge data, for instance, is beyond dispute. For example, if it is stormy today and streamflow increases in response, it is quite likely that discharge will also be relatively high tomorrow due to lingering rainfall or interflow. The consequence of such autoregressive noise is that the effective confidence level or probability of a type-I error, cr, of the Mann-Kendall test can be substantially higher than the nominal value, and false positive trend identifications are more likely to result. Another way of looking at this could be to note that Markov-chain processes have a tendency to cause extended upward or downward "runs" in a time series (e.g., von Storch and Zwiers, 1999), which might appear as a trend within a finite-duration record; unlike a structural trend, however, such clear but chance AR(l)-induced trends would eventually disappear if the record duration was extended by a factor of, say, two or three. Although the physical basis for AR(1) noise in annual streamflow-derived time series is far less clear, and Lettenmaier et al. (1994) suggested that the effects of serial correlation upon trend tests in such applications could be assumed negligible, concern about the potential impacts of autoregressive processes upon trend analysis nevertheless appears to be gaining momentum. For example, assessments of the potential effects of urbanization and climate change upon river behaviour undertaken by Leith and Whitfield (2000) and Zhang et al. (2001), respectively, applied the prewhitening or deserialization scheme of von Storch (1995) to deal with serial correlation: p w xt =x, -cx,_, (36) 168 where pwXi is the prewhitened data to be used in the subsequent trend analysis, and c is the lag-1 serial correlation coefficient as determined directly from the data using equation (35) (i.e., by the method of moments). Prewhitening is only utilized if observed c is greater than some critical value, cCru, taken to be 0.1 by von Storch (1995) on the basis of Monte Carlo results. Comparing equations (34) and (36), the prewWtening step essentially filters out AR(1) noise. This method is hereafter referred to in this appendix as single-stage prewhitening. Other recent examples of its application to river flow data include Burn and Elnur (2002), Spence (2002), and Cunderlik and Burn (2002). A more sophisticated deserialization procedure, which assumes that the time series consists of a linear combination of straight-line trend and AR(1) noise and uses an iterative method to simultaneously obtain the most appropriate values of c and slope, m, was introduced by Zhang et al. (2000) and elaborated upon by Wang and Swail (2001). This multi-stage prewhitening method is nevertheless essentially identical to single-stage prewhitening in its basic philosophy. An excellent description of the approach is provided by Wang and Swail (2001). A2.1.3 Problem, Purpose, and Scope The broader ramifications of applying such procedures have never been investigated. A potentially damaging consequence may arise from the assumption that the serial correlation observed in a time series results from an AR(1) process. In truth, large serial correlation coefficients can result not only from other autoregressive stochastic (e.g., ARMA, ARIMA) processes, but also from deterministic processes. These include periodic and quasi-periodic oscillations and, perhaps ironically, the very trends under investigation. If the high serial correlation observed in a particular time series is due to an interannual periodicity or long-term linear or nonlinear trend, application of prewhitening procedures would thus amount to a deliberate corruption of the data. The purpose of this appendix is to evaluate, in support of certain previous portions of the thesis, some of the potential consequences of deserialization and assess the relative merits of single- and multi-stage methods using Monte Carlo experiments. Elements of this work are discussed also in Fleming and Clarke (2002). A number of ancillary issues in the statistical 169 identification and quantification of trend in annual streamflow-derived time series are also addressed. A2.2 Method Time series were generated using four models: (1) white noise; (2) AR(1) noise; (3) a linear trend contaminated with white noise; and (4) a linear trend contaminated with AR(1) noise. Each time series is of length N = 50 and, to provide some physical meaningfulness, is considered here to consist of annual median daily discharge, qmed, over a 50-year period. Over much of Canada, and particularly throughout the north, a 50-year streamflow record is of typical to long duration. For each of the cases specified above, trend analyses were conducted upon the synthetic time; series using the Mann-Kendall test, the Mann-Kendall test preceded by single-stage prewhitening, and the Mann-Kendall test preceded by multi-stage prewhitening. A two-sided test was used, which is typical in such applications. The ccrU value used as the trigger for multi-stage prewWtening was set to 0.1 rather than 0.05 {Wang and Swail, 2001) to maintain consistency with the von Storch (1995) approach. For models (1) and (2), which do not contain a structural trend, the observed confidence level, dobs, was calculated as: «<*» = T ^ ( 3 7 ) MC where NPOS is the number of time series realizations for a particular model which show a positive test for trend at a nominal significance level, anom, of 0.05, and NMC is the total number of time series realizations generated using that model, which is always equal to 104 here. For models (3) and (4), which do contain a structural trend, the observed probability of committing a type-II error (accepting the null hypothesis of no trend when a trend is in fact present), fi0bS, was calculated as: 170 R - \ - N p o s Pobs ~l N MC (38) The statistical power of a test is defined as 1-/?, so that power is synonymous with the ability to detect relationships; note, however, that equations (37) and (38) are applied to different models, and it should not be concluded that a - 1-/?. Both dobs and p\bs would ideally be nil. In practice, we would like to see a0bs equal to the nominal value of 0.05 applied in the Mann-Kendall test, and to minimize fi0bs-In some applications, the estimated value of trend may be as important as the statistical acceptance or rejection of trend in a particular time series. If it is necessary to use a deserialized time series for trend identification, it would also seem necessary, or at least consistent, to use the same pre whitened data for trend estimation. For models (3) and (4), therefore, we also consider the effects of prewhitening upon estimated trend values using the percent mean absolute error in slope, MAE: . . . . 100% ^F. M A E = X\mobs, j ~ mnom (39) Hue >i ' where m0bsj is the slope observed using the estimator of Theil (1950) and Sen (1968) for the 7* time series realization, and mnom is the slope of the linear trend actually applied in the model used to generate that time series. A2.3 Results A2.3.1 Mann-Kendall Test, AR(1) Noise, and Deserialization Figure 22 presents Monte Carlo results for the case of qmed - £i (white noise) with ju, a constant bias, of 200 m3/s and a, year-to-year variability about that mean value, ranging from 5 to 40 m3/s. The particular value of // used in any of the models considered in this paper is irrelevant to the Monte Carlo results and is set only to provide some physicality. 171 Since a simply scales the variability o f the white noise model without altering the structure of the time series, the results o f Figure 22 are thus independent o f both JJ. and cr. In general it is very unlikely that the best-fit trend to any given white noise time series is exactly zero. The Mann-Kendall test (and other tests for trend) are tuned by the user's selection o f a„0m such that only anom ' 100% o f the time wil l the signal-to-noise ratio, where signal is the chance trend and noise refers specifically to white noise, be sufficiently large to result in a positive trend identification. A s noted in chapter two, this is the basic concept behind such significance tests. Perhaps less widely-recognized is the fact that strong serial correlation can occur in a white noise time series, due either to a chance trend or chance autoregressive behaviour. The implications o f this fact are apparent in the results: whereas without prewhitening ctobs = anom - 0.05, with prewhitening a0bs 3 0.04. Results from single- and multi-stage prewhitening are indistinguishable for this model. Such inappropriate but only occasional deserialization o f structurally uncorrelated white-noise time series is largely irrelevant from a practical perspective, but it does provide a preliminary quantitative indication that prewhitening can indeed pose adverse side effects. a obs 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 • r a w d a t a • s i n g l e - s t a g e p r e w h i t e n i n g • m u l t i - s t a g e p r e w h i t e n i n g 10 15 20 25 30 35 40 [m3/s] Figure 22: Case (I) - significance levels found for white noise time series. The same vertical scales are used in all graphs of error frequencies (Figures 23, 24, 25, and 27) to facilitate comparison. 172 Figure 23 presents results for the case o f qmed, t = c qmed, t-i + £i, i.e., an AR(1) noise series. Note that the same vertical axis scales are used in all graphs o f error frequencies in this paper to facilitate comparison. Values o f /j and cr are taken to be a constant 200 and 25 m 3/s, respectively. Again, the value o f cr does not affect the result. The prescribed serial correlation coefficient, c, is varied from 0 to 0.35; this is about the range considered by von Storch (1995). Without prewhitening a0bs increases with c, with a0bs noticeably greater than a„om for c > ~0.1; and with single-stage prewhitening, which assumes that serial correlation is entirely due to an AR(1) process, a0bs = otnom = 0.05. This is consistent with the results o f von Storch (1995) and confirms that single-stage prewhitening works well provided that the structural model it assumes is indeed valid. Multi-stage prewhitening provides an intermediate result, with a0bs slightly higher than a„om for higher values o f c, but substantially lower and thus more appropriate than that obtained without prewWtening. This is not unexpected given that multi-stage deserialization assumes both AR(1) noise and a linear trend to be present. In cases when trend is not present, some portion o f the high observed lag-1 serial correlation coefficient is nonetheless effectively attributed to trend, and the prewhitening procedure is incomplete. • r a w d a t a • s i n g l e - s t a g e p r e w h i t e n i n g • m u l t i - s t a g e p r e w h i t e n i n g 0.6 aobs 0.5 0.4 0.3 0.2 °o W i , T I , h i Lfcn,In,tirtrliH 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 C Figure 23: Case (2) - significance levels for AR(1) time series. 0.9 0.8 173 We next consider the case of qmed,, = mtt+ £i, i.e., a linear trend with white noise. Values of p and a are again taken to be a constant 200 and 25 m3/s, respectively, while slope, m, ranges from 0.125 to 1 m3/s/yr. Since the model contains a trend, we consider not a0bs but fiobs, the frequency with which trend is incorrectly rejected by the Mann-Kendall test. The results (Figure 24) clearly show that multi-stage and, in particular, single-stage deserialization generally increase /? (decrease the power) of the Mann-Kendall test. Overall, the iincreases in associated with the application of single-stage prewhitening to data which violate the assumptions of the method appear at least as substantial and considerably more consistent (Figure 24) than the desirable decreases in dobs obtained when the method is applied to AR(1) time series (Figure 23). For the case of m = 0.75, for example, the probability that the Mann-Kendall test will incorrectly accept the null hypothesis of no association is about 18% without prewhitening, but jumps to 28% if single-stage prewhitening is applied. This arises from the fact that serial correlation is prevalent in these time; series and thus triggers prewhitening, but is in truth due to linear trend, rather than AR(1) noise as assumed by the single-stage method. Multi-stage prewhitening performs substantially better as it effectively attributes at least some portion of the observed serial correlation to linear trend. Nevertheless, some degree of AR(1) noise is still assumed, thus providing results of intermediate quality (/?= 22% for m = 0.75, for instance). An important point with respect to these simulations is that, unlike cases (1) and (2) above, the results are sensitive to a since larger year-to-year variability obscures a constant structural trend. The effects of increasing cr are qualitatively equivalent to decreasing m; this issue is discussed in detail by Yue et al. (2002). Another salient conclusion from case (3) is that application of single-stage prewhitening to time series which contain linear trends but no autoregressive processes can induce large errors in the slope estimates (Figure 25). The effect is minimal for very small trends, presumably since the corresponding observed serial correlation coefficients from equation (35) are low; single-stage prewhitening does not, therefore, greatly corrupt the data (see equation (36)), or perhaps not at all (if c0bs < ccrit). For larger linear slopes, however, serial correlation increases and so does inappropriate deserialization, causing average slope errors of 20% or more. The effect is virtually nonexistent for multi-stage prewhitening, which 174 detrends the time series using the method o f Theil (1950) and Sen (1968) simultaneously with estimation o f AR(l)-induced c and deserialization (see Wang and Swail, 2001 for details). Apparently, incorrect partitioning o f the observed lag-1 serial correlation coefficient between both linear trend and AR(1) noise by the multi-stage method, while decreasing the power o f the Mann-Kendall test, does not materially affect slope estimate accuracy. 1 T 0.9 P obs 0.125 [m7s/yr] Figure 24: Case (3) - type-II error probabilities, linear trends plus white noise. MAE l%] • r a w d a t a • s i n g l e - s t a g e p r e w h i t e n i n g • m u l t i - s t a g e p r e w h i t e n i n g • l l 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1 m [m'/s/yr] Figure 25: Case (3) - slope M A E for time series consisting of linear trends and white noise. 175 Finally, we evaluate the case of qmed, i~mti + c qmed, t-i + £u i-e., a linear trend contaminated with lag-1 autoregressive noise. Values of / / and crare again taken to be a constant 200 and 25 ni3/s, respectively, while m ranges from 0.125 to 1.0 as in case (3) above. In the interest of conciseness and to facilitate comparison to graphical representations of the results found for the other three models considered, the results for this case (Figures 26 and 27) are presented only for c = 0.15; in general, it was found that increasing c is qualitatively equivalent to increasing m. Overall, the p0b$ responses are similar to those found for a linear trend with white noise (Figure 24), except that the differences between the results for raw data, single-stage, and multi-stage prewhitening are slightly exaggerated. Although the structural model assumed by multi-stage deserialization is fully honoured here, the procedure still must effectively apportion serial correlation between trend and AR(1) noise. This, may not be performed with complete accuracy, possibly accounting for the persistence of power losses to case (4). As was found for case (3), average error in slopes estimated using non-prewhitened data and multistage prewhitening are mutually similar, but are considerably higher here. Moreover, single-stage deserialization, which resulted in much larger slope errors in case (3), can actually provide more accurate trend estimates for model (4) for small to moderate slopes. Key to understanding this result is our definition of an accurate slope estimate as one which closely matches the structural value of m specified in the Monte Carlo simulations, rather than the actual trend of the resulting time series over the period of record. As discussed in later sections, these two quantities may be systematically different when autoregressive noise is present. For smaller structural slopes, much of the observed (net) trend results from AR(1) noise, and removal of some or all the AR(l)-generated trend by single-stage prewhitening results in more accurate slope estimates using the above definition. For larger slopes, however, much of the observed (net) trend results from the prescribed structural trend (i.e., m), and single-stage deserialization, which assumes that serial correlation is entirely due to AR(1) noise, results in removal of some substantial portion of that underlying structural trend. 176 0.125 025 0.375 0.5 0.625 0.75 0.875 1 m fV/s/yr] Figure 27: Case (4) - slope errors, linear trends plus AR(1) noise. A2.3.2 Alternative Approaches Increases in ciobs beyond nominal levels in the presence o f AR(1) noise are typically phrased in the literature as a lack o f resilience o f the Mann-Kendall test to serial correlation. This 177 could easily be taken to imply that the problem is specific to the Mann-Kendall test, which in turn begs the question of whether some other test for trend significance might be more appropriate. This issue was investigated by repeating the simulations of case (2) above using both the conventional parametric Pearson test and a novel approach for trend identification first put forward by Leith and Whitfield (1998). This method (referred to here as the Leith-Whitfield test) divides the time series into early and late parts (here we use halves), computes the median value of each population, and compares the medians using the nonparametric Mann-Whitney test; as such, it might be regarded as an interesting extension of split-sample tests for intervention (e.g., Anderson et al, 1991; Mantua et al, 1997). The Monte Carlo simulations revealed that all three tests show virtually indistinguishable responses to AR(1) noise (see again Figure 23). This may be particularly telling for the Leith-Whitfield approach, since the method does not involve the point-by-point comparisons of a correlation test, which would seem to be more vulnerable to lag-1 correlation, but instead compares the measures of central tendency of two populations. Collectively, the simulations suggest that the problem is a general one: it would appear that chance trends are simply more likely to occur in autoregressive noise than in serially independent noise, and the effect increases with the magnitude of the serial correlation coefficient. As is the case with white noise, then, the trends identified in AR(1) time series are real but not reflective of any underlying structural trend; it is only in the frequency and magnitude of spurious trends that serially dependent and independent noise differ. Two other methods of adjusting time series for AR(1) noise common in oceanographic and atmospheric science applications are differencing and determination of effective degrees of freedom. The latter approach involves deterrnining the first zero crossing of the correlogram or, alternatively, the integral time scale of the record. The effective record length or number of degrees of freedom used in the statistical test is then modified accordingly. However, this approach still effectively assumes that serial correlation arises from an autoregressive process, though not necessarily of order one. Differencing involves replacing the time series with the difference between successive data values, which is identical to single-stage prewhitening with c = 1 (equation (36)). However, this procedure is useless for removing 178 AR(1) noise in the context of facilitating accurate statistical trend detection because it also removes trend. Indeed, the method is essentially a rough finite-difference approximation of dX(t), where X(t) is the time series under consideration, or of dX(t)/dt if the data are discretized at unit time (e.g., an annual record with time units of years rather than, say, months). In the latter case, and for an ideal model consisting of linear trend and no noise, the differenced time series simply consists of a constant value equal to the slope of the trend. Neither approach provides any advantages over the single- and multi-stage methods with respect to trend detection reliability. A third, very simple, approach suggested by von Storch (1995) is to "prune" (i.e., decimate) the data. Removing every second data point, for example, might result in a time series which no longer exhibits strong lag-1 autocorrelation. However, since AR(1) processes apparently lead to increased trend detections only because such real but chance trends simply have a greater probability of arising when noise is autocorrelated, it is unclear whether this approach addresses the fundamental issue of inappropriate trend detections or just changes the value of lag-1 c measured from the data. Moreover, in the analysis of annual river discharge records, dataset sizes are very seldom sufficiently long for this data loss to be affordable. A2o3.3 Ancillary Results The above assessment of the susceptibility of the Pearson and Leith-Whitfield trend significance tests to autoregressive noise provided an opportunity to cursorily investigate a few other technical details. First, it may not be immediately clear that acceptance or rejection of median equality at some specific confidence level using the Leith-Whitfield test is exactly equivalent to acceptance or rejection of monotonic trend using a correlation (e.g., Mann-Kendall or Pearson) test at that same confidence level. This equivalency was easily assessed, and confirmed, by repeating the simulations of case (1) above using the Leith-Whitfield approach. Some initial investigation of the relative power of these three trend detection methods was also performed using the case (3) model. The Pearson test provided slightly greater power than the Mann-Kendall test, which is not surprising given that all the models used here assume Gaussian noise and that parametric tests can be somewhat more 179 powerful than nonparametric approaches when the assumption of normality is satisfied (e.g., Daniels, 1990). This may suggest that testing of a given dataset for normally-distributed noise be conducted prior to selection of a method for statistical trend detection, although the condition of normality is typically not well-satisfied for hydroclimatological data (see chapter two). The Leith-Whitfield approach appears somewhat less powerful than the Mann-Kendall test for the case (3) model, but this cannot be taken to indicate inferiority of the method; indeed, we suspect that this test might be more robust to such common difficulties as outliers and nonideal noise distributions. Rather, these results strongly suggest a need to expand upon the work presented here and by Yue et al. (2002) through a comprehensive assessment of the relative merits of the variety of trend significance tests available. Finally, a result of potentially substantial general utility is apparent from comparison of the simulation results of Figures 24 and 25. Defining trend as the signal and all other time series components as noise, the method of Theil (1950) and Sen (1968) returns accurate trend estimates for signal-to-noise ratios considerably lower than those for which the Mann-Kendall test confirms trend at a conventional a. Obviously, caution must be exercised when interpreting and using slope estimates which correspond to trends that have not been statistically verified. Nevertheless, it is clear that failure to identify trend at a typical pre-selected confidence level is far from proof of the absence of trend (see Figure 24 and Yue et al, 2002), and the result bodes well for the use of the Theil (1950) and Sen (1968) estimator in hydrologic applications which are concerned primarily with trend quantification rather than detection per se. A2.4 Discussion A2.4.1 A Practical Example The Monte Carlo results show that, in theory, the application of deserialization procedures substantially reduces the power of the Mann-Kendall test and can cause significant slope estimate errors. A simple illustrative example of the potential practical implications of this problem is provided here using observed streamflow data from the western Canadian 180 subarctic and single-stage deserialization. An annual time series of February median daily discharges of the Yukon River at Whitehorse was constructed from the H Y D A T database (Environment Canada, 1999) and is illustrated in Figure 28. In contrast to our Monte Carlo simulations, we do not know a priori the underlying structural model of this data. Nevertheless, visual inspection of the hydrograph qualitatively but very clearly reveals a substantial positive monotonic trend. Application of the two-sided Mann-Kendall test predictably and unambiguously confirms trend at a = 0.01. However, significant serial correlation (c > 0.1) is present and, using the procedure of von Storch (1995), the adjustment of equation (36) is therefore applied. With single-stage prewhitening, the trend is no longer significant at the 0.01, or even 0.05, level. Although there is indeed a possibility that the observed trend is a chance occurrence, the consequent acceptance of the null hypothesis of no trend seems, in this example, an unreasonable result. Moreover, application of single-stage prewhitening causes the slope estimate to decrease from 1.9 m3/s/yr to only 0.33 m3/s/yr, which again is not consistent with the hydrograph. 200 QFeb [rrfVs] 150 A 100 A 50 TT 1940 1950 1960 1970 1980 1990 2000 Figure 28: Observed February median daily discharge, Yukon River. 181 A2.4.2 Serial Correlation from Other Deterministic Processes Preliminary analyses, not presented here in the interest of conciseness, strongly suggest that a large lag-1 serial correlation coefficient can also arise in annual streamflow records from periodic or quasi-periodic hydrometeorological oscillations, such as ENSO. The implications of such effects for prewhitening have not been previously recognized. Although a detailed analysis is beyond the scope of this appendix, some tentative discussion is warranted as a possible starting point for future research. In practise, AR(1) and monotonic trend time series have indistinguishable correlograms and power spectra, particularly when the latter is contaminated with any significant degree of noise, serially independent or not. Specifically, the correlograms of both show a progressive decay from the lag-0 value of unity, and both periodograms display similar inverse frequency-scaled responses (hence "red noise" as a synonym for autoregressive noise). These methods cannot, therefore, generally be used to separate the two processes, particularly when applied to short, noisy annual hydrological datasets. In contrast, periodic behaviour can easily be distinguished from AR(1) noise on the basis of the autocorrelation function of the data or its Fourier transform, the periodogram (power spectrum). Frequency-domain band-cut filtering could be used to remove the periodicities, and the filtered time-domain signal might then be assessed for lag-1 serial correlation. A possible complication with this approach is that interannual oscillations tend not to be strictly periodic, but typically consist of broad frequency bands; ENSO is a clear example of this (see Burroughs, 1992). Depending on the particular dataset, removal of all significant bands in the power spectrum may thus require zeroing of a substantial proportion of the frequency-domain signal, which seems likely to lead to processing artefacts. At a minimum, this suggests that while the filtered time series might be adequate for estimating a lag-1 serial correlation coefficient theoretically devoid of the effects of interannual oscillations, the original time series should be used for all subsequent analyses. 182 A2.4.3 Deserialization Assumptions and Real Physical Systems The Monte Carlo results presented here clearly demonstrate that the recently-developed multi-stage prewhitening method, although effectively something of a compromise approach, is generally superior to the single-stage method. Combined with deseasonalization, multi-stage prewhitening is an appropriate and effective tool for addressing the effects of AR(1) processes upon trend detection in daily (and perhaps weekly or monthly) streamflow records, which clearly contain autoregressive noise, as well perhaps as the meteorological and oceanographic data for which it was first developed (Zhang et al., 2000; Wang and Swail, 2001). Moreover, the fact that such improvements are possible over the initial work of von Storch (1995) holds promise for substantial further progress in deserialization methods. As discussed in the introductory section of this paper, however, efforts to evaluate long-term trends in river flow data have justifiably focussed not upon daily data but, rather, suites of annual time series. In this specific context, the search for a universal data processing scheme which renders the Mann-Kendall (or any other) test for trend generally reliable and of constant significance level, regardless of the nature of random noise in the time series, may ultimately amount to chasing shadows. When we are faced with the analysis of some observed dataset, we in general do not know its structure: indeed, the basic purpose of the analysis is to identify that structure. In particular, it is not known whether noise is serially independent; if substantial serial correlation is present, whether it is attributable to structural trend or autoregressive noise or both; and if a trend is present, whether it is chance occurrence, an underlying structural feature of the data, or a combination of these as seen in our model (4). This is complicated by the fact that there appears to be no reliable, systematic method for empirically discriminating between the effects of trend and autoregressive noise upon the observed serial correlation, which in turn is used by prewhitening methods as the basis for deserialization. Likewise, we cannot wait for the accumulation of an additional half-century or more of hydrometric data to discriminate between a structural trend, and an AR(1)-induced run which is otherwise indistinguishable from trend over the available record. 183 The fundamental problem with prewhitening is that the underlying assumptions are no less restrictive than those associated with a failure to prewhiten. If we choose not to deserialize the data when observed c > ccru, we assume that: (1) noise is serially independent; and (2) the observed serial correlation in the time series is therefore due to something other than autoregressive noise, such as a trend or cycle. However, an AR(1) process is apparently more likely to generate chance monotonic trends than a white noise process. As a result, when we apply a trend test under the assumption of white noise (which underlies the nominal significance level of the test), we run a greater-than-anticipated risk of falsely identifying trends should AR(1) noise happen to be present. On the other hand, if we choose to deserialize the data when observed c > cCru, we assume that: (1) the noise is serially correlated; and (2) the serial correlation observed in the time series is therefore due in whole (single-stage prewhitening) or in part (multi-stage prewhitening) to autoregressive noise, rather than to a trend or cycle. When we apply a trend test under the assumption of AR(1) noise (which underlies deserialization), we run a substantially greater risk of falsely rejecting trend if the noise is in fact serially independent. Consequently, the application of prewhitening methods by no means renders a trend significance test any more general. Rather, such processing techniques are useful and appropriate when the assumption of AR(1) noise is valid and, as shown by the Monte Carlo results, hairnful when it is not. It would appear, then, that deserialization should only be applied if there is some strong external constraint supporting the assumption of autoregressive noise. A convincing physical basis for strong AR(1) processes in annual streamflow records has not been established. It might be conjectured that various storage mechanisms, if sufficiently large, slow-responding, and well-coupled to the river under consideration, could generate interannual autoregressive noise insofar as the hydrological conditions of one year might thereby influence those of the next. Such mechanisms might include major lakes, aquifers, and wetlands. The detailed physical causes of a true AR(1) process, however, may be quite complex. To the best of our knowledge, no attempt has been made to clearly ascertain whether such natural reservoirs in fact yield formal AR(1) processes on a consistent interannual basis, or simply dampen year-to-year flow variability to varying degrees. It is 184 also conceivable that interannual autoregressive noise might occur in streamflow as a direct result of serial dependence in meteorological driving variables. However, the caveats discussed previously with respect to serial correlation in river discharge also apply to time series of precipitation and temperature: observed autocorrelation must represent AR(1) noise rather than a trend or oscillation if it is to be used as justification for prewhitening. Moreover, many rivers do not always respond in a simple linear fashion to external meteorological forcing, and what true AR(1) noise might be present in precipitation and temperature could conceivably be damped or overshadowed by other processes, interactions between processes, and heterogeneous basinal characteristics. A2.5 Conclusions The evaluation of long-term trends in yearly discharge records, such as annual peak daily flow or total annual runoff, is important to a variety of issues including water resource planning, flood hazard studies, and the assessment of historical data for evidence of anthropogenic climate change effects. Prewhitening or deserialization procedures have recently been developed and applied to adjust statistical tests of monotonic trend, and the nonparametric Mann-Kendall test in particular, for sensitivity to serial dependence. Deserialization attributes much or all of the observed serial correlation in a time series to an autoregressive process; however, deterministic processes can also lead to a large lag-1 serial correlation coefficient, and the physical basis for autoregressive noise may be weaker for annual rather than more finely-discretized (e.g., daily) streamflow records. The potential consequences of using such procedures are investigated here through a suite of Monte Carlo simulations. It is found that prewhitening can substantially and inappropriately reduce the power of trend significance tests and increase slope estimate errors. The choice of whether deserialization is applied is to some degree left to the judgement and conservatism of the individual practitioner, but it is suggested that: (1) such procedures not be applied to a given annual hydrologic time series unless there is a strong site-specific physical basis for the assumption of AR(1) noise; and (2) if deserialization is performed, very recently-developed multi-stage techniques appear preferable. A number of useful ancillary results regarding trend identification in streamflow-derived data are also presented; of particular note is that the rank-based slope estimator of Theil (1950) and Sen (1968) returns accurate values for 185 signal-to-noise ratios substantially lower than those for which trend can be confirmed for a particular time series using currently-available hypothesis-testing procedures. 186