@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix skos: . vivo:departmentOrSchool "Forestry, Faculty of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Szeftel, Pascal"@en ; dcterms:issued "2010-11-10T14:45:10Z"@en, "2010"@en ; vivo:relatedDegree "Doctor of Philosophy - PhD"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description "The concept of hydrologic connectivity provides a temporally dynamic aspect to an otherwise static description of streamflow generation processes. The research presented in this dissertation focuses on the linkages between slopes, riparian zones and streams – the three main constitutive landscape elements of steep montane catchments – and the ways in which hydrologic connectivity affects hydrological processes. The research was conducted in Cotton Creek Experimental Watershed (CCEW), a snow-dominated meso-scale catchment located in the Kootenay Mountains, south-eastern British Columbia, Canada. First, the controls on slope water delivery to the stream were examined across a range of flow regimes. The spatial patterns and their dominant controls shifted from high to low flow. Catchment contributing area and slope length were good predictors of streamflow generation during high flow periods. At low flow, a new topographic index reflecting the spatial organization of flow pathways within a catchment was the dominant control on streamflow generation. Second, the dynamics of coupling between slopes, riparian zone and stream were examined through analysis of streamflow diel fluctuations. Periods representing high, intermediate and low baseflow were selected to investigate the dominant controls on diel fluctuations at nine stream gauges arranged in a nested design. Streamflow diel fluctuations were generated by evapotranspiration, the diurnal pattern of which was represented using vapour pressure deficit measured at a weather station within the catchment. Response times between climatic forcing and the various stream gauges revealed that both in-stream wave advection and dispersion processes, and transient storage processes in the riparian aquifers needed to be accounted for to explain spatio-temporal patterns of streamflow diel fluctuations. Last, the influence of slope water contributions to the stream on modelled solute transient storage and in-stream transport was examined. It was found that calibrated parameters representing transient storage vary systematically with the assumed magnitude and location of water fluxes to and from the stream."@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/29908?expand=metadata"@en ; skos:note " Stream-Catchment Connectivity and Streamflow Dynamics in a Montane Landscape by Pascal Szeftel A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Forestry) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) November 2010 © Pascal Szeftel 2010 ii Abstract The concept of hydrologic connectivity provides a temporally dynamic aspect to an otherwise static description of streamflow generation processes. The research presented in this dissertation focuses on the linkages between slopes, riparian zones and streams – the three main constitutive landscape elements of steep montane catchments – and the ways in which hydrologic connectivity affects hydrological processes. The research was conducted in Cotton Creek Experimental Watershed (CCEW), a snow-dominated meso- scale catchment located in the Kootenay Mountains, south-eastern British Columbia, Canada. First, the controls on slope water delivery to the stream were examined across a range of flow regimes. The spatial patterns and their dominant controls shifted from high to low flow. Catchment contributing area and slope length were good predictors of streamflow generation during high flow periods. At low flow, a new topographic index reflecting the spatial organization of flow pathways within a catchment was the dominant control on streamflow generation. Second, the dynamics of coupling between slopes, riparian zone and stream were examined through analysis of streamflow diel fluctuations. Periods representing high, intermediate and low baseflow were selected to investigate the dominant controls on diel fluctuations at nine stream gauges arranged in a nested design. Streamflow diel fluctuations were generated by evapotranspiration, the diurnal pattern of which was represented using vapour pressure deficit measured at a weather station within the catchment. Response times between climatic forcing and the various stream gauges revealed that both in-stream wave advection and dispersion processes, and transient storage processes in the riparian aquifers needed to be accounted for to explain spatio-temporal patterns of streamflow diel fluctuations. Last, the influence of slope water contributions to the stream on modelled solute transient storage and in-stream transport was examined. It was found that calibrated parameters representing transient storage vary systematically with the assumed magnitude and location of water fluxes to and from the stream. iii Preface Chapters 3 to 5 of this dissertation were prepared as stand-alone manuscripts for submission to peer-reviewed journals. I am the senior author of all these chapters. I was responsible for the conceptualization, design, data collection, data analyses and writing of these chapters. Specific contributions by my supervisors and co-authors are outlined below. Chapter 3 This chapter was co-authored with Professors R. D. (Dan) Moore, Markus Weiler and Younes Alila, who contributed to the design of the data collection, provided suggestions and guidance for the statistical analyses and revised the manuscript prior to submission to a peer-reviewed journal (May 07, 2010). It is currently in review. The manuscript has been modified to conform to the overall dissertation format and avoid redundancies. Chapter 4 This chapter will be co-authored with Professors R. D. (Dan) Moore and Markus Weiler, who contributed to the design of the data collection and provided suggestions for the analysis. The advective and analytical models presented in this chapter were developed by R. D. Moore, but were programmed and run by the author. Chapter 5 This chapter was co-authored with Professors R. D. (Dan) Moore and Markus Weiler, who contributed to the design of the data collection, provided suggestions and guidance for the statistical analyses and revised the manuscript prior to submission. The manuscript has been accepted for publication subject to final revisions in Journal of Hydrology (June 07, 2010). It has been modified to conform to the overall dissertation format and avoid redundancies. iv Table of contents Abstract .............................................................................................................. ii  Preface ............................................................................................................... iii  Table of contents ............................................................................................... iv  List of tables .................................................................................................... viii  List of figures .................................................................................................... ix  Acknowledgments ............................................................................................ xii  Dedication ....................................................................................................... xiii  1  Introduction .............................................................................................. 1  1.1  In search for alternatives to traditional hydrological modelling ......................................................................................................... 1  1.2  Connecting hillslope to watershed scales ............................................. 2  1.3  Hydrologic connectivity – an emergent property ................................. 4  1.4  Hydrologic connectivity and landscape discretization ......................... 6  1.5  Thesis objectives and outline ................................................................ 7  2  Study area and data collection ................................................................. 9  2.1  Land cover .............................................................................................. 9  2.2  Geology ................................................................................................. 10  2.3  Surficial material ................................................................................. 10  2.4  Meteorological and streamflow data ................................................... 11  2.5  Streamflow gauging ............................................................................. 11  v 2.6  Stage-discharge relation ...................................................................... 14  3  Linking landscape and geomorphic organization with water delivery to the stream network ....................................................................... 21  3.1  Introduction ......................................................................................... 21  3.2  Methods ................................................................................................ 25  3.2.1  Mapping stream discharge in space and time ........................................ 25  3.2.2  Terrain analysis ........................................................................................ 31  3.2.3  Statistical analyses ................................................................................... 33  3.3  Results .................................................................................................. 36  3.3.1  Discrete flow increments .......................................................................... 36  3.3.2  Principal component analysis .................................................................. 36  3.3.3  Linear regression ...................................................................................... 37  3.4  Discussion ............................................................................................ 39  3.4.1  Spatial variability and shifting controls on water delivery to the stream ................................................................................................................... 39  3.4.2  Contributing area versus accumulated area of stream feeding cells in relation to DEM quality ........................................................................... 41  3.4.3  Unexplained variance of spatial variability of water delivery .............. 42  3.5  Conclusion ............................................................................................ 44  4  Streamflow diel fluctuations reveal catchment-stream connectivity dynamics .......................................................................................................... 56  4.1  Introduction ......................................................................................... 56  4.2  Methods ................................................................................................ 59  4.2.1  Conceptual model of propagation of diel perturbations ......................... 59  4.2.2  Selection of time periods .......................................................................... 59  4.2.3  Streamflow response time to climate forcing .......................................... 60  4.2.4  Extraction of diel fluctuation characteristics ......................................... 60  4.2.5  Semi-distributed advective model of diel perturbations ........................ 62  4.2.6  Analytical model of linkages between slope, riparian zone and stream ................................................................................................................... 66  4.3  Results .................................................................................................. 68  4.3.1  Whole network assessment of lagged streamflow responses ................. 68  vi 4.3.2  Characteristics of streamflow diel fluctuations ...................................... 69  4.3.3  Simulations from the one-dimensional advective model ....................... 70  4.3.4  Results from the analytical model ........................................................... 72  4.3.5  Area of influence ....................................................................................... 72  4.4  Discussion ............................................................................................ 73  4.4.1  What is the driver of streamflow diel fluctuations in CCEW during the study period? ....................................................................................... 73  4.4.2  Amplitude of streamflow diel fluctuations .............................................. 73  4.4.3  Streamflow response time to climatic forcing ......................................... 75  4.4.4  Origins of streamflow diel fluctuations ................................................... 76  4.5  Conclusion ............................................................................................ 78  5  Influence of distributed flow losses and gains on the estimation of transient storage parameters from stream tracer experiments ..................... 90  5.1  Introduction ......................................................................................... 90  5.2  Study site description .......................................................................... 92  5.3  Methods ................................................................................................ 93  5.3.1  Solute transport modelling ...................................................................... 93  5.3.2  Delineation of study reaches .................................................................... 94  5.3.3  Tracer injection and calculation of reach-scale net inflow ..................... 94  5.3.4  Implementation of lateral fluxes into OTIS ............................................ 96  5.3.5  Estimation of parameters and derived metrics ...................................... 98  5.3.6  Assessment of parameter estimates ........................................................ 99  5.3.7  Statistical analysis ................................................................................. 100  5.4  Results ................................................................................................ 102  5.4.1  Assessment of goodness of fit ................................................................. 102  5.4.2  Comparisons between configurations InBh - InDn - InUp .................. 103  5.4.3  Standard application of OTIS - configuration SdAp ............................ 105  5.4.4  Solute mass balance, timing and location ............................................. 105  5.4.5  Assessment of parameter uncertainty .................................................. 106  5.5  Discussion .......................................................................................... 107  5.5.1  Patterns of parameter estimates among configurations InBh – InDn – InUp ......................................................................................................... 107  5.5.2  Parameter uncertainty ........................................................................... 111  5.5.3  Standard application of OTIS - configuration SdAp ............................ 112  vii 5.5.4  Implications for future stream tracer experiments and modelling with OTIS............................................................................................................. 114  5.6  Conclusion .......................................................................................... 115  6  Conclusion ............................................................................................. 132  6.1  Summary of key findings ................................................................... 132  6.2  Limitations and future directions ..................................................... 134  6.2.1  The need for additional measurements................................................. 134  6.2.2  The merits and challenges of nested designs ....................................... 135  6.2.3  The need to account for hydrologic connectivity ................................... 137  References ...................................................................................................... 138  viii List of tables Table 3-1 : List of descriptors extracted from the terrain analysis. ................ 46  Table 3-2 : Relative weights of the first three principal components. ............ 46  Table 3-3 : Selected variables and their transformation as input to the stepwise linear regression procedure. ................................................... 47  Table 3-4 : Final multiple linear models listed in order of decreasing streamflow. ............................................................................................. 48  Table 4-1 : Periods selected for extraction of diel fluctuation characteristics.79  Table 4-2 : Characteristics of stream reaches as input to the advective model. ................................................................................................................ 79  Table 5-1 : Reach characteristics. ................................................................... 117  Table 5-2 : Lateral fluxes parameterization for each configuration. ............ 117  Table 5-3 : Metrics derived from solute transport parameters and their significance. .......................................................................................... 118  Table 5-4 : RMSE values [mg l-1]. ................................................................... 118  Table 5-5 : Optimized solute transport parameters. ...................................... 119  Table 5-6 : ANOVA table for all five parameters. .......................................... 122  Table 5-7 : p-values arising from pairwise paired t-tests, adjusted according to the Bonferroni procedure. ............................................................... 122  Table 5-8 : Comparison of ܨ݉݁݀200 inter-reach variability from previous studies ................................................................................................... 123  Table 5-9 : Solute mass recovery for the different configurations of spatial organizations of lateral inflow ............................................................. 123  Table 5-10 : Solute mass held in the transient storage zone (MS) at the end of the experiment ..................................................................................... 124  Table 5-11 : Ratios AS/A. ................................................................................. 124  ix List of figures Figure 2.1 : Location of Cotton Creek Experimental Watershed .................... 15  Figure 2.2 : Overview of the catchment from sub-basin ES, looking south- west. ........................................................................................................ 15  Figure 2.3 : Geology in CCEW. ......................................................................... 16  Figure 2.4 : Surficial material in CCEW. ......................................................... 16  Figure 2.5 : Stream reach immediately upstream of the Cotton Creek (CO) gauging station. ...................................................................................... 17  Figure 2.6 : V-notch weir installed at the Elk South gauging station. ........... 17  Figure 2.7 : Typical electric conductivity breakthrough curve. ....................... 18  Figure 2.8 : Salt dilution kit. ............................................................................. 18  Figure 2.9 : Master curves for EC-NaCl concentration relation for five periods. ................................................................................................... 19  Figure 2.10 : Rating curve for gauge a) CM, b) EL, c) EU and d) GL. ............ 20  Figure 3.1 : Location of continuous and discrete sites and associated Incremental Sub-Basin (ISB) and Incremental Contributing Area (ICA). ...................................................................................................... 49  Figure 3.2 : Distribution of ICAs. ..................................................................... 49  Figure 3.3 : Typical response to snowmelt and rainfall inputs at the outlet of the Glory sub-basin (GL00). .................................................................. 50  Figure 3.4 : Mask of stream feeding cells based on flow lines derived from the D8 algorithm. ......................................................................................... 50  Figure 3.5 : Flow accumulation from stream hillslope cells (clear) to stream cells (shaded). ......................................................................................... 51  Figure 3.6 : ICA loadings of the first Principal Component. ........................... 51  Figure 3.7 : Mapping of flow increments. ......................................................... 52  Figure 3.8 : Temporal variability of flow increments. ..................................... 53  Figure 3.9 : Loadings of PC 1 – 3 from the PCA with normalized flow x increments. ............................................................................................. 53  Figure 3.10 : Time series scores for PC 1 against Julian day and against the flow regime. ............................................................................................ 54  Figure 3.11 : Relations between flow increment and contributing area. ........ 55  Figure 4.1 : Conceptual model of generation and travel of diel perturbations. ................................................................................................................ 80  Figure 4.2 : Baseflow periods selected in year 2005, illustrated at gauge ES. ................................................................................................................ 80  Figure 4.3 : Extraction of the diel fluctuations of vapour pressure deficit with a Kolmogorov- Zurbenko low-pass filter. .............................................. 81  Figure 4.4 : Distributions of hourly VPD for the three periods. ...................... 81  Figure 4.5 : Cross-correlation coefficients between VPD and streamflow records .................................................................................................... 82  Figure 4.6 : Response time distribution among the nine gauges .................... 83  Figure 4.7 : Diel fluctuations for the high baseflow period. ............................ 84  Figure 4.8 : Amplitudes of diel fluctuation [l s-1]. ............................................ 85  Figure 4.9 : (a) mean and (b) standard deviation of the amplitude of diel fluctuations as a function of catchment contributing area. ................. 86  Figure 4.10 : Distributions of daily times of troughs streamflow. .................. 86  Figure 4.11 : (a) mean and (b) standard deviation of the transformed time of daily minimum flow as a function of maximum channel length. ........ 87  Figure 4.12 : Observed streamflow at stream gauges EL07 (inflow) and EL01 (outflow), along with modelled streamflow at EL01. ........................... 88  Figure 4.13 : Predicted discharge from the riparian aquifer to the stream modulated by values of the time constant a. ........................................ 88  Figure 4.14 : Interpolated streamflow maxima at stream gauge EL01. ......... 89  Figure 5.1 : Locations of study reaches. ......................................................... 125  Figure 5.2 : Diel discharge fluctuations during the data collection period. . 125  Figure 5.3 : Spatial organization of lateral fluxes ......................................... 126  Figure 5.4 : Simulated stream solute breakthrough curve for reaches EL02 and EL07 .............................................................................................. 127  Figure 5.5 : Parameter estimates for all four configurations ........................ 128  Figure 5.6 : Derived metrics for all four configurations ................................ 128  xi Figure 5.7 : In-stream (solid line) and transient storage zone (dashed line) breakthrough curves ............................................................................ 129  Figure 5.8 : Damköhler indices for all four configurations. .......................... 129  Figure 5.9 : Composite scaled sensitivities .................................................... 130  Figure 5.10 : Parameter correlation coefficients (pcc) ................................... 131  Figure 5.11 : Discharge longitudinal profile at reach EL08 for configurations InBh, InDn and InUp. ......................................................................... 131  xii Acknowledgments This work results from direct and indirect participation of a number of persons. First, I thank my supervisors Dan Moore and Markus Weiler for their support and guidance throughout my doctoral program. I also thank Younes Alila for enlightening philosophical discussions and for granting me the opportunity to start a doctoral program at UBC. Funding for this research was provided by the Forest Investment Initiative of the Forest Science Program of British Columbia through three successive grants allocated to Younes Alila, Dan Moore and Markus Weiler. The work presented in this dissertation was part of a larger research project with multiple stakeholders. I thank our industry partners Kalesnikoff and Tembec for their financial contribution and logistic support, and the British Columbia Ministry of Forest and Range (MoFR) for logistic support. I also appreciate sponsorship by Aquatic Informatics through their provision of the Aquarius software package. Special thanks go to David Gluns (MoFR) for his enthusiasm and for sharing his field experience with me. I extend my thanks to the hydrology crowd at large in the faculty of Forestry at University of British Columbia and especially to Axel Anderson, Dan Bewley, Markus Hrachowitz, Georg Jost and Joel Trubilowicz for their friendship, laughs and often sound advices. I received invaluable help in the office and in the field from Charlotte Argue, Marianne Falda-Buscaiot and Francois Mettra. I am indebted to Mina Kavia and Georg Jost for their warm welcome and hospitality during my first few months in British Columbia. Finally, I acknowledge the (im)patience and unconditional support of my grand-mother, parents and brothers who all helped me see this thesis to its end. xiii Dedication To Hélène and André Szeftel and Janine Dreyfus-Cassin 1 1 Introduction 1.1 In search for alternatives to traditional hydrological modelling The International Association of Hydrological Sciences initiated the 2003-2012 decade on Predictions in Ungauged Basins (PUB). The PUB initiative was motivated by the will to move away from traditional hydrological modelling and search for innovative hydrological theories and reduced prediction uncertainty [Sivapalan et al., 2003]. Diverse approaches to hydrological modelling resulted in the development of a number of predictive tools of diverse natures (e.g., empirical, stochastic, process-based). Process-based models attempt to represent physical processes observed in natural systems. Their structure may vary depending on process representation, and the ways in which they account for landscape spatial heterogeneity. Applications of such models present some limitations that were promptly highlighted in the PUB agenda. First, the model structure is dictated by the modeller’s observations and understanding of a catchment’s internal processes. Moreover, it is difficult to test a posteriori the relevance of a model structure because the dynamics of a hydrological system will never be completely captured [Oreskes et al., 1994]. Second, model calibration remains a necessary step toward hydrological predictions in most applications. Model calibration is the process of adjusting parameter values so that model outputs (simulated values) mimic observed data. Hydrological models rely on calibration not only to compensate for our lack of understanding of actual processes and process interactions [McDonnell et al., 2007], but also to cope with scale-related issues. Process description is inherently scale-related [Klemeš, 1983, Oreskes et al., 1994]; therefore, the traditional, yet incorrect, implementation of micro-scale physics in hydrological models for catchment-scale predictions will inevitably lead to a 2 reliance on model calibration and the use of effective parameters [Kirchner, 2006; McDonnell et al., 2007]. To address these issues, field hydrologists have worked toward a better understanding of hydrological processes so that model structure could be improved. The field of hillslope hydrology has especially provided a wealth of detailed process descriptions and deepened our understanding of small-scale hydrological processes. However, the diversity of landscape properties and climate forcing under which hydrological processes have been investigated has inevitably led to a myriad of idiosyncratic behaviours that are at odds with the rather rigid model structure formulation of current catchment models [McDonnell et al., 2007]. It is questionable whether additional overly detailed descriptions of small-scale processes will produce the long awaited breakthrough in hydrological understanding. Rather, hydrologists should strive to develop new theories and concepts to describe processes directly at the appropriate scale – that of the whole catchment – or to allow for a transfer of knowledge gained from hillslope studies to larger scales. 1.2 Connecting hillslope to watershed scales Considering the relatively poor understanding of hydrological processes at the watershed scale compared to that at the plot and hillslope scales, attempting to upscale these processes seems a logical way forward. However, simplification of hillslope process descriptions is required, in order to “whittle down unnecessary details” [Sivapalan, 2003] prior to upscaling. Klemeš [1983], describing the work from Eagleson [1972], first coined the terms upward and downward routes as complementary approaches to develop theories that transcend spatial scales. Klemeš’ [1983] vision has since resonated in the hydrological community through repeated calls for a search for organizing principles governing landscape properties and ultimately hydrological response characteristics [Dooge, 1986; Kirchner, 2003; Sivapalan, 2005; McDonnell et al., 2007, Troch et al., 2009]. More recently, hydrologists have sought emergent properties, by studying spatial variability 3 through inter-hillslope comparisons [Uchida et al., 2006] and variability through time at a set hillslope under a broad range of climate forcing [Tromp- van Meerveld and McDonnell, 2006a; Lehmann et al., 2007]. The term emergent refers to a characteristic large scale hydrological behaviour arising from concurrent realizations of small scale, often stochastic, processes. Non- linear dynamics, and especially threshold behaviour, has often been identified as an emergent property of a catchment hydrological response (e.g., [Spence and Woo, 2003; Tromp‐van  Meerveld  and  McDonnell,  2006b; Hooke, 2007; James and Roulet, 2007; Harman et al., 2009]). Spence and Woo [2003] proposed the “fill-and-spill” mechanism to explain the temporal variability of the hydrological response of a 5 ha catchment within the Canadian Shield landscape (comprising bedrock uplands and soil-filled valleys). They demonstrated that the water storage capacity of valley segments needed to be progressively filled with antecedent inputs in order to connect uplands to valleys, and generate a catchment-scale response. A similar threshold mechanism was also observed in the Panola Mountain Research Watershed [Tromp-van Meerveld and McDonnell, 2006b]. Harman et al. [2009], focusing on recession flow, demonstrated that a non-linear catchment-scale recession behaviour, here a power law relation, may emerge from the aggregation of heterogeneous, yet linear responses of multiple hillslopes (sensu linear reservoir theory [Sivapalan et al., 2002]). Emergent behaviours are not intuitively predictable from the direct observation of small-scale processes. They are only the large-scale expressions of these processes. Therefore, the concept of emergence provides a convenient avenue for the exploration of hydrological processes: a seemingly complex non-linear response may be explained from simpler processes occurring at smaller scale. More importantly, emergent behaviours highlight the critical need to move away from observing small-scale processes in isolation and account explicitly for their interactions when examining large-scale responses. 4 1.3 Hydrologic connectivity – an emergent property Because flow pathways can be defined from the pore scale to the watershed scale they are an essential component of the notion of emergent properties in hydrology. Flow pathways have the potential for transmitting matter and energy across the landscape. When this potential is realized, an elevated state of connectivity is reached and individual responses of landscape elements aggregate into a macroscale response. Therefore, understanding the spatiotemporal variability of connectivity between landscape elements is critical in defining a continuum of hydrological response across spatial scales [Bracken and Croke, 2007; Ali and Roy, 2008]. The role of connectivity in streamflow generation has been documented for decades (e.g., [Hewlett and Nutter, 1970; Taylor, 1982; Roulet, 1990; Taylor et al., 1993; Devito et al., 1996]). However, despite an upsurge in the use of the term “hydrologic connectivity” in recent years, there is still a lack of consensus for a unified definition [Ali and Roy, 2008]. Pringle [2001] defined hydrologic connectivity as the “water-mediated transfer of matter, energy, and/or organisms within or between elements of the hydrologic cycle”. This definition makes no assumptions with regards to spatial and temporal scales, which results in hydrologic connectivity being assessed in diversified contexts. Numerous hillslope–scale studies identified subsurface flow as the dominant mechanism of runoff generation in shallow forest soils [Dunne and Black, 1970; Mosley, 1979; McGlynn et al., 2002]. However, only recently has hillslope hydrological response been viewed as the emergent expression of a multitude of sub-scale individual feature properties. Hillslope excavation revealed the critical role of hydrologic connectivity among preferential flow features in delivering water to the lower face of the hillslope [Sidle et al., 2000; Anderson et al., 2009]. Stochastic representation of macropore density and direction [Weiler et al., 2003] or small-scale percolation processes [Lehmann et al., 2007] into numerical models captured the range of connected versus disconnected states and allowed for a more faithful account of the hillslope hydrological behaviour. 5 In contrast, large-scale studies usually examine the influence of anthropic activity on hydrologic connectivity, or lack thereof, in extended stream networks. Dams typically play a major role in the inhibition of downstream sediment transport, ultimately leading to coastal erosion in deltas and estuaries [Pringle, 2001]. More diffuse and seemingly inoffensive anthropic riverine developments may also have, when accumulated, severe repercussions on downstream ecosystem viability due to hydrologic connectivity. For instance, the cumulative channelization and diversion through pipes of headwater streams within large basins reduced on-site nutrient processing, which resulted in increased nutrient loadings and subsequent eutrophication in higher order streams and rivers [Freeman et al., 2007]. Headwater streams have been the subject of close attention not only because they account for two-thirds or more of the total stream length in a hydrographic network [Leopold et al., 1995], but also because they are direct recipients of water and solute fluxes originating on slopes [Gomi et al., 2002]. Traditionally, headwater catchments have been discretized into a mosaic of constitutive elements that may or may not generate runoff or contribute to the catchment hydrological response. In essence, this is the distinction between active versus contributing area [Ambroise, 2004]. At a set time, an active area generates runoff but will become contributing only if the area is hydrologically connected to the catchment outlet. The transition of various areas within the catchment from being active to being contributing is usually assessed indirectly via traditional hydrograph separation methods. In contrast with hillslope-scale assessments of hydrologic connectivity, the question is no longer what mechanism causes hydrologic connectivity but more what the implications are, in terms of water yield and solute mobilization, when hydrologic connectivity changes throughout the catchment [Stieglitz et al., 2003]. Temporal patterns of stream water dissolved organic carbon [McGlynn and McDonnell, 2003; Pacific et al., 2010], alkalinity [Tetlzaff et al., 2007] or nitrate [Ocampo et al., 2006] have been key to understanding the dynamics of hydrologic connectivity in the slope- riparian transition of meso-scale catchments. 6 1.4 Hydrologic connectivity and landscape discretization Hydrologists have employed a range of approaches to discretizing a landscape for the purposes of describing hydrologic connectivity. Terminology varies depending on the morphological settings but all approaches differentiate terrestrial from aquatic entities with a more or less well defined transition zone in between. “Hillslope, floodplain and active channel” were proposed as generic terms that could apply to any water-land interface from an ecosystem perspective [Gregory et al., 1991]. The representation of hydrologic connectivity in the landscape discretization process is not mundane; it is often associated with the development of a conceptual model to explain large-scale hydrological behaviour. For instance, Hydrological Response Units (HRU) are at the heart of several hydrological models including Precipitation-Runoff Modelling System [Leavesley et al. 1983; 2002] and the Cold Regions Hydrological Model [Pomeroy et al., 2007]. Within a model, HRU are defined as spatial units within which processes and states can be adequately represented by single sets of parameters, state variables, and fluxes “but having a place in a landscape sequence or water/snow cascade” [Pomeroy et al., 2007]. Landscape discretization can also be found as a critical component of conceptual models. Spence and Woo [2003] investigated hydrologic connectivity between “bedrock upland and soil- filled valley” in the Canadian Shield landscape. The function of bedrock uplands in their conceptual model was to collect water inputs and transmit them to downslope valleys. Soil-filled valleys acted as a switch (the “fill-and- spill” concept discussed in section 1.2) and delivered water to the stream only when cumulative upslope contributions exceeded a certain threshold. In a different environment, Sidle et al. [2000] proposed the “hydrogeomorphic paradigm” to explain the hydrological response of a small headwater catchment discretized into “zero-order hollow, riparian zone and stream”. The activation of a network of preferential flow pathways was vital to linking all responses of individual landscape elements and deliver water to the stream. More recently, Buttle et al. [2004] and Jencso et al. [2009] focused on the linkages between “hillslope, riparian zone and stream” in two contrasting 7 environments. Landscape organization was key to sustaining slope water delivery to the stream. The schema for landscape discretization used by Jencso et al. [2009] in a steep montane catchment will be retained in this dissertation, as this is the one that relates the most closely to the site in which data were collected for the present study. 1.5 Thesis objectives and outline The literature review of the previous sections has highlighted the following points: i. the future of hydrological science does not lie solely in the study of processes in isolation ii. traditional reductionist approaches to modelling a catchment’s hydrological response often fail to account for emergent behaviors iii. hydrologic connectivity is an appropriate avenue for a holistic approach to studying catchment-scale hydrological response, because it allows for the expression of emergent behaviors The present work explores the temporal dynamics of linkages between upslope areas, riparian zones and streams in a meso-scale catchment. It focuses on three aspects of watershed functions from a stream perspective – delivery of water from upslope areas, diel streamflow fluctuations resulting from evapotranspiration, and downstream solute transport – and illustrates the emergence and influence of connectivity at the reach- or catchment-scale. This study draws upon nested streamflow measurements to capture the spatial variability of stream-catchment connectivity. A nested study design is powerful, but has not often been employed in previous studies due to the substantial resources required to install stream gauges or to conduct spatially distributed streamflow measurement campaigns. A detailed review and evaluation of the use of nested streamflow measurements in hydrologic studies will be provided in Chapter 3. The remainder of this thesis comprises five chapters. Chapter 2 8 includes a description of the study site and the data collection installations and procedures. Chapter 3 addresses the variability in space and time of slope water contributions to the stream network. This chapter tests the relevance of contributing area and landscape organization as first order controls on streamflow generation and its spatio-temporal variability. The linkages and interactions among lateral hillslope discharge, hyporheic exchange, evapotranspiration and streamflow are explored in Chapter 4. In this chapter, catchment-wide hydrologic linkages are explored through a multiple-site analysis of streamflow diel fluctuations during baseflow periods. Observations from Chapter 3 set the basis for a new solute transport modelling framework presented in Chapter 5. This chapter explores the influence of hydrologic coupling between slopes, the riparian zone and the stream on transient storage and downstream transport of solutes. More specifically, it investigates the sensitivity of solute transport parameters to the specification of spatially variable boundary conditions. Together, these three complementary chapters provide a comprehensive assessment of hydrological connectivity between the dominant elements of a montane landscape with a focus on its implications for streamflow dynamics and downstream solute transport. The thesis closes with a summary of key findings and recommendations for future research (Chapter 6). Figures and tables are inserted at the end of each chapter. 9 2 Study area and data collection The work presented in this dissertation is part of a larger research project investigating the response of snow accumulation and ablation, runoff generation and sediment transport processes to forest disturbance such as logging practices and ongoing mountain pine beetle epidemic. Research was conducted in the Cotton Creek Experimental Watershed (CCEW) (49º21’25”N, 115º48’00”W), a 17.4 km2 catchment located in the Kootenay Mountains, south-eastern British Columbia, Canada (Figure 2.1). A picture of the catchment is provided in Figure 2.2. Climate normals computed at a nearby long-term meteorological station (Cranbrook airport, 49º36’36”N, 115º47’00”W, 939 m a. s. l.) are characteristic of a snow-dominated environment. Monthly maximum temperature varied between -4 ºC in winter and 26 ºC in summer during the period 1961-1990. During the same period, total precipitation averaged 386.5 mm, falling mainly as snow (32%) from November to March and as rain (68%) during the rest of the year. 2.1 Land cover Harvesting activities in the CCEW date back 70 years and currently 34% of the watershed has been harvested, creating a mosaic of forested, forest regeneration and clear-cut areas (Figure 2.1). The areal coverage of species is 51% Lodgepole pine (Pinus contorta), 6% Engelmann spruce (Picea engelmannii), and 5% subalpine fir (Abies lasiocarpa). Along most of the channel network, the riparian vegetation has been preserved in its original state. It is mainly composed of Engelmann spruce (Picea engelmannii), with balsam fir (Abies lasiocarpa) dominating at higher elevations where the foothills are steep and the riparian zone poorly developed, and western red cedar (Thuja plicata) in the wetter and flatter lower part of the watershed. 10 Deciduous trees such as alder (Alnus tenuifolia), trembling aspen (Populus tremuloides) and willow (Salix sp.) are often the dominant understory riparian vegetation, especially in disturbed areas such as old logging roadbeds and skid trails. 2.2 Geology Geological maps were provided by the BC Geological Survey (www.em.gov.bc.ca/geology) at the 1:250000 scale. In CCEW, the geology is dominated by two Pre-Cambrian formations: Creston and Kirchener (Figure 2.3). Both formations are composed of homogeneous sediments consisting of fine-grained quarzites, argillites and limestones. 2.3 Surficial material A distributed assessment of surficial material was conducted by Halleran [2004] at 168 soil survey sites located throughout the whole catchment. This resulted in a simplified map of locally dominant surficial material (Figure 2.4). The vast majority of CCEW was glaciated, leaving a thick compact basal till layer. In some locations, a relatively looser ablation till mixed with glacio-fluvial sediments cover linear segments of the watershed. These deposits have been identified as potential evidence of paleo- fluvial systems. Most of these systems make up the current stream network. Areas covered by pure ablation till are usually found at mid-slope. These areas are characterized by a large water storage capacity relative to areas covered in basal till. Bedrock outcrops are most often located on ridge tops and along the watershed boundaries, but also along the main divides within the watershed. Colluvial zones of moderate size are often found downslope of nearby bedrock outcrops. Alluvial deposits are found along the lower reaches of CCEW stream network. Except for the deeply incised stream channels and bedrock outcrops, surficial materials are overlaid by a soil layer about 70 cm deep. 11 2.4 Meteorological and streamflow data Meteorological variables were continuously recorded at two climate stations located in pine regeneration stands at 1390 m and 1780 m a. s. l., denoted as the lower and upper stations, respectively (Figure 2.1). The stations recorded hourly air temperature, relative humidity, wind speed, shortwave radiation, rainfall and snow depth. Distributed snowmelt inputs were calculated from stratified snow surveys over the snow accumulation and ablation periods [Jost et al., 2007; 2009]. CCEW encompasses two second order streams, Elk Creek (EL) and Cotton Creek, (CO). Elk and Cotton creeks each encompass two first-order headwater streams, Elk North (EN) and Elk South (ES), and Cotton North (CN) and Cotton South (CS). A fifth headwater (Glory Creek, GL) is part of the CO sub-catchment. A nested network of nine stream gauges was installed in autumn 2004 in order to record streamflow for catchments with contributing areas ranging from about 1 km2 to 17.4 km2 (Figure 2.1). For reference, Figure 2.5 illustrates the stream reach immediately upstream of the CO stream gauge (natural cross-section) while the V-notch weir installed at the ES stream gauge is shown in Figure 2.6. Stage was continuously recorded with Odyssey capacitance water level probes on a 15 min time step at the nine stream gauges (two V-notch weirs and seven naturally constricted cross-sections). According to the manufacturer, the resolution of Odyssey probes is about 0.8 mm. The loggers were encased in PVC standpipes (32 mm inner diameter) screened at their base with geotextile to avoid the accumulation of fine sediments. Discrete discharge measurements were conducted at a range of flow regimes and rating curves were subsequently established for every stream gauge to convert stage records to streamflow. 2.5 Streamflow gauging Distributed assessments of streamflow are an essential aspect of the present work. Some field campaigns required over 25 streamflow measurements to be conducted in a single day. This unusually large number 12 of measurements introduced time constraints to the stream gauging exercise. Moreover, motorized access was limited so that the complete “salt dilution kit” had to fit in a large backpack. The relatively large set of data that resulted from multiple field campaigns was, as much as possible, processed automatically with custom-built macros. Field protocol Discharge typically ranged from about 1 l s-1 in the headwater stream at low flow to about 1 m3 s-1 at the outlet of CCEW during the peak of the freshet season. Salt dilution gauging is a suitable method for this range of discharge values and was therefore adopted for all streamflow measurements conducted within CCEW. A number of measurement campaigns were conducted during years 2005-2008 to gauge streamflow at a range of flow regimes and at multiple locations within CCEW. Discharge was measured using a slug injection of sodium chloride (NaCl) following standard procedures described by Day [1976] and Moore [2005]. The technique involved the release of a known mass of salt at the injection location and the recording of an electrical conductivity (EC) breakthrough curve (BTC) at the sampling location. Electrical conductivity was used as a surrogate for the concentration of NaCl in stream water, which necessitated further calibration to convert EC readings into salt concentration and derive the salt concentration BTC. Electric conductivity readings are sensitive to stream water temperature. Therefore, they were corrected for temperature with the loggers’ built-in compensation curve. Two EC meters were used: the WTW Multi 350i and the Topac Consort C931. Both meters logged on a time step of five seconds and allowed for the storage of 1200 to 1800 data points. A typical EC breakthrough curve is presented in Figure 2.7. If complete tracer mass recovery and a steady flow regime over the duration of the stream gauging measurement are assumed, discharge is computed as follows, using a simple mass balance approach: 13 ܳ ൌ ெ ׬ ஼ሺఛሻ·ௗఛ ೟ బ    2‐1 where M is the mass of NaCl [g] injected into the stream, C [g l-1] is its concentration, Q [l s-1] is the discharge estimate and t [s] the duration of the BTC. EC-salt concentration calibration Ideally, the slope of the locally linear relation between electric conductivity and salt concentration (calibration slope) is determined for every stream gauging trial. This is, however, time consuming and was deemed incompatible with the time constraints of the field campaigns. Instead, a master calibration curve that relates background EC and individual calibration slopes was derived for every streamflow measurement campaign from a subset of all the sites at which streamflow was gauged within a set campaign. An illustration of the master curves for five successive field campaigns is given in Figure 2.9. For instance, during the third field campaign, streamflow was gauged at 69 locations within CCEW. Among these 69 sites (some measurements were replicated), calibration slopes were measured 30 times. Calibration slopes at the 39 un-calibrated sites were derived from the appropriate master calibration curve from their value of background EC. Uncertainty of discharge estimates For a reach with efficient lateral mixing, uncertainty in an individual measurement is typically 5% [Day, 1976; Moore, 2005]. Two measurement campaigns at low flow (September 5, 2006, n=50) and high flow (June 3, 2007, n=53) were conducted using two EC probes placed on either side of the sampling cross- section to test the assumption of complete lateral mixing of the injected tracer. In 97% of all measurements, the paired discharge 14 estimates were within 7% of each other, consistent with an uncertainty of 5% for each measurement at 95% confidence. In addition, the precision with which background EC (EC of the stream water prior to the injection of salt) was determined is equal to the EC meter round-off error (1 µs cm-1). Depending on the shape of the EC breakthrough curve (especially, its flatness), an uncertainty of 1 µs cm-1 on the determination of the background EC value may result in a large uncertainty in the derivation of the discharge estimate (up to 20%). For this reason, the mass of salt injected in the stream was estimated in order to record an EC breakthrough curve that peaked between 60 and 100 µs cm-1 above background EC. 2.6 Stage-discharge relation A stage-discharge relation or rating curve is required to convert a time series of stream water level into a hydrograph. An empirical relation was fit to discrete streamflow measurements with a non-linear least squares regression. The relation between water level and discharge is as follow: lnሺܳሻ ൌ ܽ ൅ ܾ · ݈݊ ሺ݄ ൅ ܿሻ 2‐2 where Q [l s-1] and h [cm] are the discharge and stream water level, respectively, a, b and c are the model parameters and ln(·) is the natural logarithm operator. Equation 2-2 actually describes a log-transformed power law relation between water level and discharge. Examples of rating curves are provided in Figure 2.10. Extrapolation of rating curves beyond the range of data used to fit them can involve substantial uncertainty. Therefore, streams were gauged extensively around the peak of the freshet and during low base flow in order to cover the range of observed streamflow and minimize extrapolation. 15 Figure 2.1 : Location of Cotton Creek Experimental Watershed and main instrumentation. Figure 2.2 : Overview of the catchment from sub-basin ES, looking south- west. 16 Figure 2.3 : Geology in CCEW. Figure 2.4 : Surficial material in CCEW. 17 Figure 2.5 : Stream reach immediately upstream of the Cotton Creek (CO) gauging station. The channel is about 3 m wide. Figure 2.6 : V-notch weir installed at the Elk South gauging station. The Odyssey water level logger encased in a screened PVC standpipe and the staff gauge (about 1.5 m tall) are visible on the left bank of the controlled cross-section. 18 Figure 2.7 : Typical electric conductivity breakthrough curve. Background EC is equal to 187.0 µs cm-1. Figure 2.8 : Salt dilution kit. From left to right and top to bottom, 1 l plastic volumetric flask, 1 l plastic mixing bottle, concentrated solution of NaCl (4.86 g l-1),WTW 350i EC meter, injection plastic bottle, 5 ml glass pipette and sealed bag of NaCl (15.004 g). 19 Figure 2.9 : Master curves for EC-NaCl concentration relation for five periods. Points are individual calibration slopes. Regression lines are master curves for distinct periods. 20 Figure 2.10 : Rating curve for gauge a) CM, b) EL, c) EU and d) GL. Dots are individual discharge measurements. Root Mean Square Error is equal to 26.1, 13.2, 7.80 and 1.73 l s-1 for CM, EU, EL and GL, respectively. a b c d 21 3 Linking landscape and geomorphic organization with water delivery to the stream network 3.1 Introduction An ongoing theme in hydrology is the search for organizing principles that could explain the linkage between landscape heterogeneity and hydrologic response [Klemeš, 1983; Dooge, 1986; Sivapalan, 2005; McDonnell et al., 2007]. In particular, there is a long tradition of research linking streamflow dynamics to the geomorphic structure of a catchment. For example, Horton formulated geomorphic laws that were further developed by Shreve [1966] and ultimately led to the Geomorphologic Instantaneous Unit Hydrograph framework (GIUH) [Rodríguez-Iturbe and Valdés, 1979]. However, the GIUH framework relies on the assumption that travel time on slopes is negligible compared to that within the catchment’s stream network and therefore cannot be applied to small or meso-scale basins. At the meso-scale, numerous studies focused on the relation between a catchment’s contributing area and its hydrologic response [Pilgrim et al., 1982]. Contributing area has been found to be a poor predictor of the catchment’s hydrological response in semi-arid and arid catchments because of the high spatial variability of convective storm rainfall and the influence of channel transmission losses on the integrated hydrological response [Renard, 1977]. However, the relation between a catchment’s contributing area and its hydrological response is stronger in more humid, temperate environments [Allis, 1962; Sopper and Lull, 1965; Kincaid et al., 1966; Pilgrim et al., 1982]. In such environments, a topographic index introduced by Beven and Kirkby [1979], which combines upslope contributing area and local slope gradient, has been used to analyze and model runoff generation dynamics (e.g., [Moore and Thompson, 1996; Hjerdt et al., 2004; Monteith et al., 2006]). However, 22 considerable scatter exists in empirical relations between streamflow and contributing area, evidence that other factors must be considered. At smaller scales, those of the small catchment and the hillslope, detailed characterization of runoff generation and water delivery processes has been sought after, in order to gain a better understanding of precipitation-runoff relations. In temperate forested catchments, subsurface flow has been identified as the dominant process for water delivery to the stream [Hoover and Hursh, 1944; Hursh, 1944 cited in Hewlett and Hibbert, 1967; Kirkby, 1978, cited in Tani, 1997; Mosley, 1979; Peters et al., 1995]. For over four decades, hillslope scale studies have focused on two key questions: what are the dominant mechanisms that drive subsurface flow; and where does subsurface flow occur both laterally and vertically? These studies documented the amount, timing and chemistry of subsurface flow, be it matrix flow or preferential flow in macropores and soil pipes [Mosley, 1979; Tsuboyama et al., 1994; Sidle et al., 2000; Weiler and Naef, 2003; Weiler and Fluhler, 2004; Anderson et al., 2009]. However, perhaps due to the small size of the experimental site or the nature of the research questions, no general conclusions with regards to universal controls on subsurface flow mechanisms have been drawn. Rather, a myriad of idiosyncratic subsurface flow drivers have been identified. At sites with relatively simple geometry or well defined boundary conditions, surface topography has been identified as the dominant control on runoff generation. In their benchmark study, Dunne and Black [1970] found that variable saturated areas, sustained by upslope subsurface flow, were the principal sources of water delivery to the stream. Such saturated areas were only observed on a concave hillslope (as opposed to planar and convex hillslopes), which indicates that surface topography and especially curvature is a first order control on runoff generation and delivery to the stream. Anderson and Burt [1978] also identified convergent topography as a major control on subsurface flow generation. In shallow soils with a low permeability substratum, hillslope hydrology is predominantly controlled by surface topography, because in this specific case, surface topography is a reasonable proxy for the bedrock (or confining layer) topography upon which 23 saturated flow pathways typically develop [McDonnell, 1990; McGlynn and Seibert, 2003], and also for the hydraulic gradients driving flow. However, working on a glaciated slope in south coastal British Columbia, Hutchinson and Moore [2000] showed that the surface topography is a good proxy for water table gradients and resulting subsurface flow paths at higher flows; at lower flows, the basal confining layer became a better proxy for the water table surface and thus hydraulic gradients. The role of confining layer topography has been well studied at the Panola Mountain Research Watershed, where the spatial variability of subsurface storm flow at the bottom of a trenched hillslope was related to bedrock topography indices [Freer et al., 1997; 2002]. In contrast to the results of Hutchinson and Moore [2000], this relation became stronger as hillslope wetness increased due to a higher hydrological connectivity between flow pathways. At the same site, Tromp-van Meerveld and McDonnell [2006b] showed that under any magnitude of rainfall event, saturated areas developed in various regions of the hillslope, characterized by shallow soils and lower soil moisture deficit, but this did not necessarily generate subsurface storm flow at the trench face. However, increased connectivity of the flow pathways and significant water delivery to the trench face, at the bottom of the hillslope was achieved for cumulative precipitation exceeding a threshold of 55 mm. This led to the confirmation of the “fill-and-spill” theory (initially proposed by Spence and Woo [2003]) as being the dominant temporal control on water delivery, in which the filling of bedrock depressions works as a connectivity switch and controls delivery of subsurface storm flow to the trench face. These recent studies illustrate a shift from the focus on runoff generation processes alone to their study in conjunction with the hydrologic connectivity of flow pathways to deliver water to the stream. Intuitively, the potential for disconnected flow pathways or hydraulic barriers increases with increasing spatial scales and as we move from the hillslope to the catchment scale, the importance of hydrologic connectivity cannot be ignored. Although direct measurement of hydrologic connectivity has been studied at the hillslope scale, where it was mainly restricted to discrete preferential flow 24 features [Sidle et al., 2001; Anderson et al., 2009], connectivity is usually investigated at larger scales to sample enough spatial heterogeneity. For instance, snowmelt and rainwater at various locations in a catchment may not contribute significantly to the hydrograph due to poor connectivity with the downslope stream [McNamara et al., 2005; Ocampo et al., 2006]. Hydrologic connectivity may be particularly affected during spring snowmelt, when it is common for lower portions of slopes to lose their snow cover first. As the lower slopes dry out, the connectivity of flow paths linking snowmelt from upslope areas to the stream channel may become less effective. The persistence of hydrologic connectivity between hillslope and stream throughout the catchment was recently related to distributed upslope accumulated area along the stream [Jencso et al., 2009]. In other studies, patterns of soil moisture deficit [Devito et al., 1996] and bedrock storage capacity [Buttle et al., 2004] were found to control catchment scale runoff delivery. A number of studies have used incremental discharge measurements to investigate the controls on runoff generation. Anderson and Burt [1978] related discharge increments to surface topographic convergence. Huff et al. [1982] related lateral contributions to the stream to bedrock outcrops dipping into the stream. They hypothesized that water collected along bedding planes, then moved laterally along strikes to finally discharge into the stream. The use of incremental discharge measurements has also been used to highlight the influence of a permeable karstic bedrock [Genereux et al., 1993] and the effect of spatially heterogeneous inputs [Kuraś et al., 2008] on the spatial variability of contributions to the stream network. Findings from Kuraś et al. [2008] confirmed that in snow-dominated environments, local temperature and incoming radiation, topography and land cover characteristics all interact to control spatial and temporal dynamics of snowmelt processes, which in turn, affect the dynamics of runoff generation and delivery to the stream [Cazorzi and DallaFontana, 1996; Ferguson, 1999; Williams and Tarboton, 1999; McCartney et al., 2006]. The present study examines the spatio-temporal variability of runoff delivery to the stream in a partly forested, snow-dominated catchment. Based 25 on the analysis of a spatially intensive set of streamflow data, the following questions are addressed: (1) Can the spatial variability of water delivery be described with simple, GIS-derived predictors? (2) Does the nature of these predictors shift along with dominant runoff processes throughout the year? (3) Does the dominance of contributing area, a standard control on water delivery under relatively homogeneous temperate forcing, still hold in climatically and physiographically complex environments? A new index that represents the topographic complexity of the landscape is proposed as dominant control on water delivery to the stream when hydrologic connectivity between upslope areas and stream is limited. 3.2 Methods 3.2.1 Mapping stream discharge in space and time Differential gauging Water delivery to the stream was assessed using the differential gauging method. Discharge was measured at various locations along the stream network of CCEW over two consecutive days, during a number of field campaigns. The discharge increment along the stream reach bounded by two neighboring sites is computed as the difference in discharge between the sites. The incremental contributing area is defined as the area drained by the reach segment. Water is delivered to the reach segment in the form of gross lateral inflows that encompass discrete overland flow contributions such as tributaries or neighboring springs, diffuse subsurface contributions seeping from the channel banks and discharge from hydrologically connected aquifers. The same reach may lose water in the form of gross lateral outflow that encompasses discharge to hydrologically connected aquifers, evaporation or via hyporheic flow pathways that extend beyond the study reach. A discharge increment is therefore the net integrated measurement of water fluxes from and to the stream along the reach segment. This quantity may be negative when gross losses exceed gross contributions. 26 As for any integrated measurement, the flow increment summarizes the sub-reach scale spatial variability and, due to its ease of measurement, has the advantage of being applicable to basin scales up to a few tens of square kilometers. Previous studies demonstrated the relevance of the differential gauging method to draw inferences on subsurface hydrological processes, or as a sampling strategy for water quality (e.g., [Huff et al., 1982; Genereux et al., 1993; Grayson et al., 1997; Ruehl et al., 2006; Kuraś et al., 2008]). The magnitude and spatial pattern of gross fluxes from and to the stream is of utmost importance when investigating water budget, solute mobilization and transport processes, and biogeochemical cycling. However, measuring gross fluxes integrated over the length of the reach involves the analysis of a reach- scale tracer breakthrough curve [Payn et al., 2009], which is time consuming and was not feasible given the broader scopes of the current study. Flow increment uncertainty as a criterion for the determination of the locations of discharge measurement sites Discrete streamflow measurements were conducted at fifty four sites distributed throughout the perennial stream network of CCEW (Figure 3.1). The distance between sites, roughly 250 m, was set a priori to ensure sufficient flow gain between two neighboring sites. Fine tuning of the final location of every site was based on consideration of sampling channel cross- section and mixing reach morphology. At every site, the length of the mixing reach (between injection and sampling locations) was set a priori to 25 times the channel width following recommendations by Day [1976] and then adjusted depending on flow conditions (spring freshet versus autumn baseflow). Mixing reaches ideally showed no or small pools and a least one constriction to facilitate lateral mixing of the tracer. When necessary (most often during less turbulent baseflow periods), small artificial constrictions were built between injection and sampling locations to enhance lateral mixing of the tracer. Typically, a discharge measurement campaign lasted two consecutive days, gauging both Cotton and Elk creeks (Figure 2.1). Either Elk or Cotton 27 creek was gauged entirely over one day, starting at its outlet and moving upstream to avoid tracer contamination. Early analysis of the magnitude of flow increments raised concerns about their precision. Because a flow increment is essentially a difference between two, sometimes three (in case of a confluence) discharge measurements, the flow increment magnitude is often small relative to streamflow measurement magnitude, and the resulting relative error in flow increment is large. For a reach with efficient lateral mixing, uncertainty in an individual measurement is typically 5% [Day, 1976; Moore, 2005]. For two independent discharge measurements, Q1 and Q2 [l s-1], the variance (var(·)) of their difference is: ݒܽݎሺܳଶ െ ܳଵሻ ൌ ݒܽݎሺܳଶሻ ൅ ݒܽݎሺܳଵሻ 3‐1 Application of Equation 3-1 revealed that the uncertainty in flow increments was often greater than their magnitude, especially for lower elevation reaches where streamflow measurements were greater and flow increments smaller relative to their headwater counterparts. While streamflow generally increases downstream as contributing area is being accumulated, it may decrease when local or diffuse gross stream losses caused by evaporation and discharge to the hyporheic zone and aquifers exceed gross stream gains (e.g., [Hamilton and Moore, 1996; Story et al., 2003; Ruehl et al., 2006]). Nonetheless, in order to reduce the uncertainty in flow increments and reduce the influence of local processes compared to their catchment-scale integrated counterparts, every other site was discarded (Figure 3.1). Doing so lengthened stream reaches, increased flow increments and decreased their error bound. In the following, continuous site refers to the nine streamflow gauges equipped with a water level sensor and for which continuous records were available. Discrete site refers to every site that was not permanently equipped with a water level sensor but at which streamflow was gauged 28 during measurement campaigns. An Incremental Sub-Basin (ISB) is defined as the area drained between two (sometimes three, in the case of a confluence) continuous sites, while an Incremental Contributing Area (ICA) is defined as the area drained between two (sometimes three, in the case of a confluence) neighboring discrete sites. The whole stream network was divided into 23 reaches ranging from 226 to 911 m in length, and each draining an incremental contributing area ranging from 0.05 to 1.1 km2. Three reaches drained substantially larger areas, with incremental contributing areas ranging from 1.7 to 2.7 km2 (Figure 3.1 and Figure 3.2). Temporal mapping of discharge measurements The streamflow regime is characterized by a snowmelt-induced freshet with peak flows typically occurring in May, followed by a prolonged recession limb ending between September and October. Responses to rainfall events are superimposed on the main recession limb, with time scales ranging between a day and a week (Figure 3.3). The months of September and October are normally dominated by a baseflow period with no major rain events and fairly constant streamflow. Water stored as soil moisture and groundwater generally decreases following the end of snowmelt, resulting in a decline in discharge due to storage depletion. Therefore, discharge at the outlet of the watershed was used as an index of overall catchment wetness. Discharge measurement campaigns were conducted during years 2005- 2007. Dates were set in order to capture the seasonal variability of streamflow. However, during the freshet season (high flow) the upper part of the drainage network was still covered in snow, thus impeding stream gauging efforts. Therefore, all measurement campaigns occurred after the annual peak flow. This resulted in nine measurement campaigns, with high flow being sampled once, intermediate flow three times and baseflow five times (Figure 3.3). 29 Synchronization of discharge measurements Each main sub-basin (Elk or Cotton sub-basin) was gauged over an entire day during a typical salt dilution campaign. Diel fluctuations were observed at all continuous streamflow gauges. These fluctuations were caused by radiation-induced snowmelt during the freshet and evapo- transpiration during the baseflow period and will be further examined in Chapter 4. Since every streamflow measurement campaign extended over two entire days, the data not only captured the spatial variability of streamflow among discrete sites but also its temporal variability. While Kuraś et al. [2008] considered that: “adjusting for this desynchronization with continuously measured outlet flow, however, was not a possibility, as the assumption of a direct relation between these measures contradicts the premise of our research,” the temporal adjustment of discharge measurements was deemed necessary to remove the temporal source of variability from the data and better extract the streamflow spatial variability throughout the catchment. The amplitude of diel fluctuations often reached 10% and up to 23% of the daily streamflow during low flow periods. Therefore, failing to correct for diel fluctuations would have resulted in substantial errors for streamflow increment estimates. Model of diel streamflow fluctuations Observations of diel streamflow fluctuations led to the selection of a sinusoidal model with time-varying amplitude that could account for asymmetry between the rising and falling limbs of the fluctuations. The following model: 30 ܳ ൌ ܳ଴ ൅ ൫ܣ଴ ൅ ܣଵሺݐ െ ݐ଴ሻ൯ . cos ቀ ଶగ ் ሺݐ െ ݐ଴ሻ ൅ |߮|ቁ 3‐2 was fit to observed hydrographs over a period of time ranging from 8 to 15 h and centered on the time of actual discharge measurements at discrete sites. t and t0 are the time of measurement and initial time, respectively, T is the period of the signal (equal to one day), Q0 is the initial discharge value, A0 and A1 are two amplitude terms and φ is the phase. The time-dependent amplitude A1 was included in the model to account for dissymmetry in the diel fluctuations. The four parameters (Q0, A0, A1 and φ) were optimized using a nonlinear optimization routine implemented in the programming language R (http://www.r-project.org/index.html). The absolute value of the phase was used in the model in order to “lock” its polarity and get consistent parameter estimates. Since the sampling time was 15 minutes, between 32 and 60 observations were available to fit each model. Information transfer from continuous to discrete sites Streamflow measurements at all sites (discrete and continuous) were corrected to 09:00 of the first day of the measurement campaign, a reference time chosen arbitrarily. At continuous sites, the corrected discharge value was that recorded at 09:00. At discrete sites, the nearest continuous site was used to correct its value as defined by Equation 3-3: ܳ௧೘೐ೌೞ ௗ௜௦ ൌ ொ೟೘೐ೌೞ ೎೚೙೟ ொ೟ೝ೐೑ ೎೚೙೟ · ܳ௧ೝ೐೑ ௗ௜௦    3‐3  where tmeas and tref are the actual measurement time and the reference time (09:00), respectively, and the superscripts cont and dis refer to discharge measurements at continuous and discrete sites, respectively. For instance, the streamflow measurement was equal to 17.5 l s-1 at 10:05 at the discrete site EL02. The nearest continuously gauged site, EL01 recorded 24 l s-1 at 31 10:05 and 23.5 l s-1 at the reference time 09:00. The ratio 23.5/24=0.98 was used as a multiplicative factor to correct the discharge measurement at EL02 (0.98.17.5=17 l s-1). The use of this method implies that streamflow time series at two distinct locations along the same reach exhibit similar temporal patterns as long as they are reasonably close to each other. Here, no assumption was made with regards to the dependence of streamflow on contributing area to transfer information among locations. In contrast to Kuraś et al. [2008], who worked in a 4.74 km2 experimental watershed that was only continuously gauged at its outlet, the density of stream gauges in CCEW was 2.5 times larger and the nearest continuous site was never farther than 1 km along the stream network. Considering the generally close proximity of continuous and discrete sites, we are confident in our approach to synchronize streamflow measurements. 3.2.2 Terrain analysis Derivation of DEM-based variables and indices A 10 m Digital Elevation Model (DEM) was generated from ortho- rectified stereo photographs. Mean vertical error and precision under the forest canopy were 2.4 m and 6.4 m, respectively. The System for Automated Geoscientific Analyses (SAGA, http://www.saga-gis.org), a Geographic Information System (GIS) platform, coupled with the programming language R, was used to extract topographic information from the DEM. The D8 algorithm [O'Callaghan and Mark, 1984] as implemented in SAGA was used to derive the flow accumulation grid and related stream network and flow pathways because it provides an unambiguous catchment delineation solution. First, flow direction and stream network were computed using the flow accumulation algorithm and an accumulated area threshold for channel initiation set equal to 13.5 ha, based on field observations. Then, ICAs were delineated given the locations of discrete and continuous sites. This delineation was based on surface topography as represented in the DEM and led to the derivation of a mask that designates every grid cell as belonging to 32 a unique ICA. The mask of ICAs was used to derive zonal statistics for every descriptor and every ICA. Some descriptors were computed from the population of grid cell values within every ICA (e.g., the mean distance to channel); others were derived directly at the scale of the ICA (e.g., percentage of clear-cut area). When a variable was derived from the population of cells, the following characteristics were computed: mean, standard deviation, median, minimum and maximum, and lower and upper quartiles. These will subsequently be referred to as m, sd, q50, q0, q100, q25 and q75, respectively. Physiographic variables that were extracted from the DEM analysis are listed in Table 3-1. The DEM was also resampled to 25 m and 50 m to capture the larger scale topographic structures (e.g., contour curvature). The nearest neighbor algorithm was used as a resampling method. The same descriptors derived from the 10 m DEM were also derived from the 25 m and 50 m DEMs. Stream feeding cells Hydrologic connectivity between upslope areas and stream is requisite for the delivery of water to the stream. When connectivity is established, water is delivered through the cells that drain directly into the stream network. The expression stream feeding cells will be used to designate these cells. They have to be distinguished from both the remaining hillslope cells that are not directly draining into the stream, and the stream network cells (Figure 3.4). The steps involved in identifying the cells that drain directly into the stream network were: i. Flag all cells that belong to the stream network, ii. Create a buffer of one cell width around the stream network, iii. Flag all cells within this buffer for which the flow direction points at a stream cell. The distribution of accumulated areas derived from all stream feeding cells of an ICA represents the population of contiguous sub-catchments discharging along the stream reach of this ICA. It is of particular interest 33 because it reflects the topographic complexity of the landscape, and therefore of the flow concentration of upslope areas in connection with the stream. The mean accumulated area is sensitive to the number of stream feeding cells in a given ICA but not to the value of accumulated area they bear (Figure 3.5a versus Figure 3.5c). Alternatively, the standard deviation of accumulated area of stream feeding cells was used in order to reflect the spatial heterogeneity of landscape elements within the ICA (Figure 3.5a versus Figure 3.5b). It is hypothesized that the standard deviation of the accumulated area of stream feeding cells increases with increasing flow concentration and thus reflects a propensity forhydrologic connectivity between upslope areas and stream and ultimately, water delivery to the stream. 3.2.3 Statistical analyses Principal Component Analysis A principal component analysis (PCA) was first run on the flow increment data to extract the dominant sources of variability contained in the dataset. In this analysis, there are 23 variables (discrete sites) and nine observations per variable (streamflow measurement campaigns).In this case, because the number of variables is greater than that of observations, the calculation was based on a singular value decomposition (SVD) of the data matrix using the R function “prcomp”. As applied here, PCA generates a set of vectors in which the components for each site (loadings) reflect a spatial pattern; each vector is associated with a time series of scores that indicate the temporal variability of the strength and polarity of the corresponding spatial pattern. Principal components (PCs) with eigenvalues substantially greater than unity represent dominant patterns of spatio-temporal variability within the data set, and thus potentially represent an underlying structure. This analysis is based on the hypothesis that the spatial pattern of connectivity, and thus of water delivery to the stream channel, should vary with 34 catchment wetness. If this hypothesis is valid, then the time series scores associated with at least one of the dominant PCs should be correlated with catchment wetness. Discrete flow increments measured at high flow were larger than those measured at baseflow. Therefore, most of the variability contained in the dataset was concentrated in the measurements collected during the freshet, rather than in those collected during baseflow periods. We first ran a PCA on the non-normalized flow increment dataset. Most of the variability (94%) was contained in the first principal component (PC 1). Moreover, the loadings of PC 1 did not show any contrast (all values were negative) (Figure 3.6). The contrast in hydrological regime between freshet and baseflow periods was clearly of sufficient magnitude that it hid any other source of variability, spatial or temporal. In order to filter out this source of variability, equal weights needed to be given to all measurement campaigns in the PCA. To serve this purpose, all flow increments belonging to the same measurement campaign, at a given date, were divided by their sum. Consequently, relatively small increments measured at baseflow, for instance, were normalized by a relatively small sum and the resulting normalized flow increments were of similar magnitude compared to their normalized counterparts measured during other flow regimes. Note that while the main effect of this normalization is to remove the overall seasonal cycle, the expression of this cycle remains in the spatial pattern of flow increments. It is common to standardize variables (mean equal to zero and standard deviation equal to one) prior to running a PCA because they might be of different magnitudes or have different units. In the present case, all variables were flow increments and it was desirable to maximize the contrasts between them. Therefore, flow increments were centered (by subtracting the mean) but their standard deviations were not standardized prior to running the PCA. 35 Linear regression models After the dominant sources of variability were identified via the PCA, linear regression models were implemented to explore the controls on the variability contained in each flow increment measurement campaign. Since it is commonly acknowledged that contributing area is a dominant control on runoff delivery, the linear relation between flow increment and ICA contributing area alone was first examined. Then, multiple linear regression models were built to explain the spatial variability of runoff delivery to the stream. The development of such models involved two steps. First, a restricted number of predictors (transformed when necessary, to increase linearity) was selected based on visual assessment of univariate plots. These predictors were meant to contain enough information to explain the spatial variability of flow increments under any flow regime. Then, an automated stepwise regression algorithm was implemented to fit a specific linear regression model to individual flow increment campaigns. Only the most significant predictors were retained in the final statistical models to make them as robust and parsimonious as possible. Therefore, depending on the flow increment campaign, both the nature and number of significant predictors may vary. In the aforementioned models, ICA contributing area was systematically retained in the final statistical models even when not found to be statistically significant in the presence of other predictors. This approach was preferred to using specific flow increments as response variable, as the very notion of spatial uniformity underlying the concept of specific discharge contradicts our research hypotheses in relation to spatial heterogeneity of runoff generation and delivery processes. Both response variable (flow increments) and predictors were standardized (mean equal to zero and standard deviation equal to one). Doing so removes all scale effects among predictors and the respective standardized regression coefficients are a reflection of their dominance in a given model. The standardized regression coefficient measures the change of the standardized response variable for a unit change of the standardized predictor, all other standardized predictors of the model being held constant. 36 3.3 Results 3.3.1 Discrete flow increments The highest flows and largest flow increments occurred on June 3-4, 2007. Streamflow ranged from 15.4 to 134.6 l s-1 in the Elk sub-basin and from 61.7 to 385 l s-1 in the Cotton sub-basin. All reaches were net gaining with the exception of reaches EL09 (-0.19 l s-1) and EN13 (-0.22 l s-1), which were marginally losing water. The lowest flows and smallest flow increments occurred 0n October 7-8, 2007. Streamflow ranged from 0.5 to 11.8 l s-1 in the Elk sub-basin and from 3.8 to 34.2 l s-1 in the Cotton sub-basin. All reaches were net gaining with the exception of reach ES11 (-0.2 l s-1). The spatial variability of flow increments throughout the whole catchment is presented in Figure 3.7. Light (blue) and dark (red) solid circles indicate net gaining and losing reaches, respectively. The temporal variability of flow increments among flow regimes can be visualized from the map legends in Figure 3.7 but is best highlighted in Figure 3.8. 3.3.2 Principal component analysis Kaiser’s rule [1961] was used to determine the number of significant PCs. According to that rule, the cut-off value of proportion of variance was set equal to 1/23≈0.04. The first three PCs were kept in the analysis, the fourth one being of marginal power (proportion of variance equal to 0.05). They accounted for 85% of the variance with more than 50% contained in the first PC (Table 3-2). PC loadings are mapped in Figure 3.9. For all three PCs, loadings showed a strong variability among ICAs. PC 1 was positively expressed in the ICA containing the Glory headwater sub-basin (CO5) while ICAs belonging to Elk South and Cotton North sub-basins usually exhibited negative loadings (Figure 3.1). Two other ICA loadings containing major tributaries to the Cotton Creek main channel (CO3 and CO7) showed negative loadings. The time series score for PC 1 plotted in a near linear fashion with the 37 date and in a downward curvilinear relation with the sum of flow increments (Figure 3.10). PC 1 was negative during high flow and gradually increased to reach positive values at baseflow. This designates PC1 as the expression of the seasonal cycle and separates ICAs that are highly responsive to snowmelt input from those that show a more constant response throughout the year. Loadings for PC 2 and PC3 did not appear to follow a simple pattern, and regression analysis did not reveal any clear relations with topographic variables. 3.3.3 Linear regression Contributing area as single predictor The relation between flow increment and incremental contributing area was generally weak but improved at higher flow regimes, compared to flow increment campaigns during baseflow periods (Figure 3.11). Patterns of flow increments during baseflow were similar across years (Figure 3.11g, h and i). This suggests that there is a stationary baseline response, which is likely insensitive to antecedent soil moisture conditions, as snowmelt input in 2006 was nearly three times as large as that in 2005 (289.3 mm versus 107.5 mm [Jost et al., 2007]). Multiple linear regression The process of selecting explanatory variables for input to linear regression models was not only based on visual assessment of univariate plots between a predictor and flow increments but also on a conceptual representation of the dynamics of runoff generation and routing to the stream network. Based on the results from previous smaller scale studies, water delivery to the stream was expected to be influenced by 1) the amount of water infiltrating into the ground and 2) the likelihood that infiltrated water would reach the stream. Contributing area and vegetation were expected to affect the amount of infiltrated water, while slope length and curvature were expected to provide the flow concentration required to sustain water delivery 38 to the stream. Predictors that were retained as input to the stepwise procedure regression are summarized in Table 3-3. In some cases, the relation between the response variable and some of the predictors was slightly curvilinear. These predictors were transformed accordingly and both the original variable and its transformation were input to the stepwise regression procedure. Since the primary objective of this study was to determine the controls on water delivery to the stream network, we did not look for regression models that would give maximal predictive power. Rather, simple models based on plausible causal relations between dependent and independent variables were sought after. For that reason, every multiple linear regression model was screened for multicollinearity by computing the Variance Inflation Factor (VIF). It is usually acknowledged that a VIF value of five or higher is evidence of high correlation between two or more predictors [Kutner et al., 2004]. The number of predictors in every statistical model was further reduced in order to ensure that VIF values were lower than this threshold, thus keeping multicollinearity as low as possible in the final regression models. In practice, predictors with the lowest explanatory power (highest p- value) were removed from the list of explanatory variables in an iterative process until the VIF of a given regression model fell below a value of five. In general, regression fits were reasonable among flow increment campaigns with multiple R2 ranging from 0.48 to 0.80 (adjusted R2 ranged from 0.40 to 0.76). However, only 33% (18% adjusted) of the variance of the first campaign of fall 2006 (05-09-2006) was explained by the linear regression model (Table 3-4). There was no clear pattern between flow regime and goodness of fit of the regression as expressed by the adjusted R2 (p = 0.44). In other words, the selection of predictors did not appear to introduce any bias in the resulting model goodness of fit with regards to flow regime. Contributing area was found to be not significant at the 95% confidence level in all but three flow increment campaigns, two of them at high flow. The median overland flow distance to the stream was found to be a significant predictor in three out of nine occurrences, all three at high flow regime. The standard deviation of accumulated area of stream feeding cells 39 (sfc.aa.sd) was a significant predictor in six out of the nine regression models, most often at low flow. Median contour curvature computed from the 10 m DEM was never a significant predictor but its counterpart computed from the 50 m DEM was marginally found to be a significant predictor, usually at intermediate flow regimes. Forest cover, as a fraction of the ICA contributing area, was only once found to be a significant predictor. Percentages of selected surficial materials were never found to be significant predictors (Table 3-3). 3.4 Discussion 3.4.1 Spatial variability and shifting controls on water delivery to the stream The spatial variability of water delivery to the stream was substantial and dynamic in time. The PCA confirmed the hypothesis that the spatial pattern of water delivery varies systematically with catchment wetness. The multiple regression analysis further revealed that, during high flows, incremental contributing area and median overland flow distance to channel controls the spatial variability of lateral inflow. During intermediate to low flows, the topographic complexity of the landscape controlled hydrologic connectivity between upslope areas and stream, and therefore water delivery to the stream. These statistical results are discussed below in light of our general understanding of snowmelt runoff processes and field observations at CCEW. During the freshet season, the spatial variability of runoff delivery to the stream is, to a large extent, driven by the spatial variability of snow accumulation and melt. Snowmelt is controlled by the local energy budget which is a function of incoming radiation, aspect, vegetation and atmospheric conditions. The complex interactions between all these variables yield patterns of runoff generation that are spatially heterogeneous and highly variable in time. The catchment-wide hydrological response is the result of 40 snowmelt occurring over an area for which size and location are dynamic in time. Snowmelt is then redistributed downslope, decreasing the water deficit over much larger areas than the area over which snowmelt occurs. Field observations in CCEW revealed that during the snowmelt period, overland flow generated at the soil-snowpack interface emerges immediately downslope of the snowline or isolated snow patches, then infiltrates further downslope and reaches the stream network via subsurface flow pathways. In addition, the concurrent activation of ephemeral tributaries to route runoff as concentrated overland flow to the perennial stream network (e.g., [Reggiani et al., 1998]) was also observed. Sustained concentrated overland flow over several hundreds of meters was observed, especially in the large regenerating clearcut in the Elk South headwater (Figure 2.1), where the soil structure has been disturbed by relatively recent logging activity. It is proposed that the dominance of incremental contributing area and median overland flow distance to channel as controls on water delivery to the stream is the direct reflection of the redistribution of water through subsurface and concentrated surface flow paths. As soil water deficit increases following the freshet, the relative contributions of concentrated surface flow to streamflow decrease while those of subsurface flow and groundwater discharge increase. Conjointly, the heterogeneity of flow convergence on slopes, as defined by the standard deviation of accumulated area of stream feeding cells, becomes a driver of water delivery to the stream. The persistence of hydrologic connectivity between stream and upslope areas in relation with upslope accumulated area has recently been demonstrated by Jencso et al. [2009]. However, in their study, upslope water contributions to the stream were not quantified. Results of the present study reveal that rather than ICA contributing area as a whole, its internal structure and decomposition into concave, convex and planar areal elements affects slope water redistribution and ultimately water delivery to the stream, especially at intermediate to low flows. It is further proposed that a quantitative relation exists between water delivery to the stream and upslope water redistribution within the ICA. 41 3.4.2 Contributing area versus accumulated area of stream feeding cells in relation to DEM quality Extracting properties from the stream feeding cells alone appears to be a promising alternative to extracting these properties from the entire population of cells making up the ICA, as it implicitly includes information about the spatial organization of flow pathways within the ICA. This is especially relevant for the property of accumulated area, since the contributing area of an ICA is simply the sum of accumulated area of all stream feeding cells plus the areal coverage of all stream cells. The contributing area of an ICA, as an integrated value of accumulated areas, only depends on the delineation process and not on the distribution and spatial organization of individual flow pathways that eventually merge at its outlet. The delineation of an ICA depends on the DEM precision and its ability to capture topographic features. It also depends on the type of flow accumulation algorithm used to derive flow pathways. If the DEM fails to capture adequately the true topographic features far enough from the ICA boundaries, this will not affect the resulting value of the ICA contributing area. Conversely, the accumulated area of stream feeding cells is sensitive to the trajectory of individual flow pathways within the ICA. Their location depends on the surface topography as defined in the DEM in the vicinity of the stream network. It is therefore critical to the determination of location and accumulated area of stream feeding cells that the DEM faithfully captures true topographic features throughout the ICA and especially near its stream network. In the present study, accumulated area was derived from a 10 m DEM with vertical error up to 6.4 m under forest canopy. In CCEW, most of the riparian corridor is densely vegetated, even in the vicinity of large clearcuts, because of the preservation of a riparian buffer strip. This most certainly affected the accuracy of the inferred trajectories of upslope flow pathways and especially the derivation of locations and accumulated area values of stream feeding cells, thus introducing some noise in the linear regression 42 models. Although still expensive and not affordable in this study, LiDAR technology has become increasingly popular over the past few years. DEMs derived from such techniques are characterized with a greater precision than their counterparts derived from stereoscopic photography. This is especially true under forest cover [Adams and Chandler, 2002; Su and Bork, 2006]. LiDAR data would provide a refined definition of flow pathways and locations of stream feeding cells. Alternatively, ground survey of the riparian area would at least provide a better representation of surface topography in the vicinity of the stream network. This approach would refine the location of stream feeding cells but is unlikely to significantly affect accumulated area values as they mainly depend on upslope converging patterns. 3.4.3 Unexplained variance of spatial variability of water delivery Between 24 and 58% of the variance could not be explained by the regression models. Although some of this unexplained variance is undoubtedly due to errors inherent in the calculation of flow increments, and errors in topographic indices due to errors in the DEM, three potential process-related sources of unexplained variance are discussed below, with attention to the unexpected non-significance of geology and land cover as predictors. Spatial variability of geologic characteristics Geologic indices and those derived from the map of surficial material were not found to be significant predictors of flow increments in the present study. The present results contrast with a number of previous studies that reported relations between discharge to the stream and geologic rock formations [Armbruster, 1976 cited in Smakhtin, 2001; Genereux et al., 1993; Mulholland, 1993; Tague and Grant, 2004]. The current study design allowed for a refined spatial coverage of streamflow measurements to the detriment of a high temporal sampling frequency. Bedrock nature and structure are usually associated with groundwater flow systems, deeper flow pathways and 43 longer time scales than shallow sub-surface flow systems. Therefore, persistent characteristics of the hydrological response driven by geology should have been captured despite the low temporal frequency characterizing our dataset. The delineation of ICAs was based on surface topography. It is usually accepted that in steep catchments dominated by shallow subsurface flow, flowpath trajectories can be inferred from surface topography [Freeze, 1972]. However, when the groundwater contribution to streamflow cannot be neglected, flow boundaries as defined by surface topography become irrelevant. Due to the anisotropic nature of the hydraulic conductivity of the underlying bedrock, subsurface water can either be redistributed within the catchment or come from neighboring catchments via regional aquifers [Genereux et al., 1993]. Observations of sustained streamflow during baseflow, even after virtually no input to the watershed for several months, support the hypothesis that groundwater delivered substantial amounts of water to the stream (Figure 3.3). However, the spatial resolution of geologic information available for this study was too coarse and did not allow for the analysis of groundwater contributions so that the significance of groundwater fluxes across CCEW divides and ICA boundaries could not be tested. Land cover Soil structure is usually more complex in forests than in disturbed areas. The hydrology of forest soils is strongly influenced by preferential flow pathways which, depending on their hydrologic connectivity, may substantially alter the hydrological response from the plot to the catchment scale [Mosley, 1979; Weiler and Naef, 2003; Anderson et al., 2009]. Moreover, forest cover affects snow accumulation and snowmelt timing [Pomeroy et al., 1998; Baldocchi et al., 2000] and also drives transpiration patterns [Kelliher et al., 1993]. Therefore, vegetation can be a first order control on the spatial distribution of input of water to the soil and its movement, especially in snow-dominated environments [Verry et al., 1983; Jones and Post, 2004; Schnorbus and Alila, 2004]. Jost et al. [2007] demonstrated that in CCEW, under similar topographic settings (elevation, aspect, surface gradient), 44 snowmelt input to the soil was lower for forested areas and was more spatially variable compared to clearcut areas. They also showed that snowmelt was delayed under forest cover compared to clearcut areas due to lower melt rates. Nonetheless, forest cover was not found to be a significant predictor of water delivery to the stream in the present study. In part, this lack of significance might have resulted from the fact that ICAs were delineated based on the locations of streamflow measurement sites along the stream network alone without trying to delineate homogeneous hydrological response units based on topography, geology and forest cover. As a result, nearly every ICA was comprised of both forested and clearcut areas, weakening any signal present in the discharge increments. Local processes The influence of local sub-reach scale processes has been identified as a possible source of noise in the flow increments dataset. Streamflow typically integrates all runoff generation and transport processes occurring at the catchment scale. It also includes local processes such as hyporheic exchange or discharge from the stream to groundwater that will locally modify this basin scale integrated discharge measurement [Hamilton, 2008]. However, despite the fact that there was evidence of such processes as some reaches were consistently losing water (e.g., CN11, Figure 3.1 and Figure 3.7), the influence of these local processes on the calculation of discharge increments cannot be quantified using the differential gauging approach. In the present study, net flow increments were equated with gross lateral inflows, neglecting gross losses to first order. However, under certain circumstances, gross and net inflow may show substantial differences and this approximation is no longer appropriate [Payn et al., 2009]. 3.5 Conclusion In CCEW, a snow-dominated, partly forested catchment, the dominant controls on water delivery to the stream network shift throughout the annual 45 cycle, depending on the overall wetness state of the catchment as measured by the stream discharge at its outlet. During the snowmelt period, when large portions of the ICAs contribute to water delivery via distinct yet concurrent runoff processes, hydrologic response scales, in part, with contributing area. However, as soil moisture deficit increases following freshet, the standard deviation of accumulated area of stream feeding cells, reflecting the spatial organization of flow pathways within the ICA, replaces contributing area as a dominant control on water delivery. Predictors related to land cover and geology did not explain the spatial patterns of water delivery to the stream. The portion of unexplained variance is attributed to the influence of local processes such as hyporheic exchange or stream discharge to connected aquifers, as well as measurement uncertainty. Additional measurements in relation to gross lateral fluxes would help better quantify water delivery to the stream. The present study sought to identify topographic drivers of a catchment’s hydrological response as expressed in terms of statistical characteristics of the distribution of topographic variables computed at the scale of the ICA. This approach contrasts with the more traditional approach in distributed modelling, in which an attempt is made to delineate flow paths explicitly, and is consistent with calls by Koutsoyiannis et al. [2009] and others to move away from purely deterministic representations of physical processes in hydrologic modelling. 46 Table 3-1 : List of descriptors extracted from the terrain analysis. Cell-scale variable (derived from a population characteristic) ICA-scale variable (derived at the scale of the ICA) Elevation [m] Contributing area [m2] Accumulated area [m2] Reach length [m] Topographic gradient [-] Clear-cut area [m2] Aspect [-] Areal percentage of clear-cut [-] Total, contour and vertical curvatures [m-1] Areal coverage of a given surficial material class [m2] Total, vertical and horizontal distance to the stream following flow lines [m] Areal percentage of a given surficial material class [-] Gradient to the stream following flow lines [-] Table 3-2 : Relative weights of the first three principal components. PC 1 PC 2 PC 3 Standard deviation [l s-1] 81.75 57.3 34.22 Proportion of variance [-] 0.51 0.25 0.09 Cumulative proportion [-] 0.51 0.76 0.85 47 Table 3-3 : Selected variables and their transformation as input to the stepwise linear regression procedure. Variables that are derived at the scale of the ICA are highlighted in bold. Other variables are characteristics of the population of cells in a given ICA. All predictors were derived from the 10 m DEM with the exception of cc50.q50. Symbol Predictor description mean min max ca ICA contributing area [km2] 0.691 0.154 2.66 fc.rel Fraction of forested land [%] 28 0 99.7 L Reach length [m] 500 226 912 brk.rel Areal percentage of bedrock outcrop [%] 14.2 0 69.4 bas.rel Areal fraction of pure basal till layer [%] 66.4 21.2 93.7 abl.rel Areal fraction of pure ablation till layer [%] 5.9 0 27.7 ofd.q50 Median overland flow distance to the nearest stream cell, following flow lines [m] 304 95 545 log.ofd.q50 Log (ofd.q50) [-] 2.45 1.98 2.74 cc.q50 Median contour curvature from the 10 m DEM [m-1] 3.6 10-4 -1.63 10-4 1.02 10-3 fg50.q50 Median gradient to the stream following flow lines, from the 50 m DEM [-] 0.27 0.17 0.35 ti.mean Mean topographic index [-] 8.33 7.76 8.70 cc50.q50 Median contour curvature, from the 50 m DEM [m-1] 2.42 10-4 -2.20 10-4 1.11 10-3 sfc.aa.sd Standard deviation of accumulated area from stream feeding cells [km2] 63.59 0.16 254.21 48 Table 3-4 : Final multiple linear models listed in order of decreasing streamflow. Standardized regression coefficients and their associated p-value (bold). A missing value indicates that the predictor was not retained in the final model. Variables that were not found to be significant in any of the flow increment campaigns were omitted. Symbols are identical to that presented in Table 3-3. m, sd and q50 refer to mean, standard deviation and median, respectively. date flow [l s-1] R2 ca fc.rel sfc.aa sd ofd q50 ti m cc50 q50 fg50 q50 log(ofd) q50 03/06/2007 545 0.58 0.689 - - 0.664 -0.555 - - - <0.001 - - 0.002 0.01 - - - 01/07/2005 272 0.76 -0.056 -0.245 0.885 0.326 - - - - - 0.034 <0.001 0.011 - - - - 26/06/2007 203 0.4 0.571 - - 0.643 -0.604 - - - 0.004 - - 0.008 0.017 - - - 24/08/2005 81 0.69 -0.388 - 1.104 - - -0.427 - - 0.097 - <0.001 - - 0.003 - - 05/09/2006 66 0.57 -0.501 - 1.152 - - -0.293 - - 0.073 - <0.001 - - 0.069 - - 10/09/2006 63 0.18 -0.087 - 0.536 - - -0.619 - -0.434 0.822 - 0.163 - - 0.07 - 0.212 11/10/2005 56 0.51 -0.219 - 0.918 - - -0.242 - - 0.447 - 0.004 - - 0.152 - - 26/10/2006 48 0.55 -0.451 - 1.085 - 0.325 - 0.355 - 0.146 - <0.001 - 0.099 - 0.051 - 15/10/2007 45 0.74 -0.618 - 1.326 - - - - 0.268 0.01 - <0.001 - - - - 0.043 49 Figure 3.1 : Location of continuous and discrete sites and associated Incremental Sub-Basin (ISB) and Incremental Contributing Area (ICA). Figure 3.2 : Distribution of ICAs. 50 Figure 3.3 : Typical response to snowmelt and rainfall inputs at the outlet of the Glory sub-basin (GL00). Open circles designate times at which streamflow measurement campaigns were achieved. Figure 3.4 : Mask of stream feeding cells based on flow lines derived from the D8 algorithm. 51 Figure 3.5 : Flow accumulation from stream hillslope cells (clear) to stream cells (shaded). m and sd refer to mean and standard deviation of accumulated area for the population of stream feeding cells. Figure 3.6 : ICA loadings of the first Principal Component. Dot area is scaled with loading value. 52 Figure 3.7 : Mapping of flow increments. Blue (light) and red (dark) solid circles designate net flow gains and losses, respectively. Dot area is scaled with flow increment magnitude. Legends are expressed in l s-1. Maps are ordered by decreasing streamflow magnitude. Measurement date is indicated above every map with format year-month-day. 53 Figure 3.8 : Temporal variability of flow increments. Measurement date is indicated on the x axis with format year-month-day. Figure 3.9 : Loadings of PC 1 – 3 from the PCA with normalized flow increments. Blue (light) and red (dark) solid circles represent positive and negative loadings, respectively. Dot area is scaled with loading magnitude. 54 Figure 3.10 : Time series scores for PC 1 against Julian day and against the flow regime. Julian day is expressed in months. Flow regime is expressed as the sum of flow increments at a given date. Dashed lines are local polynomial regression fits to the coordinates of PC 1. 55 Figure 3.11 : Relations between flow increment and contributing area. Plots are ordered by decreasing flow magnitude. Measurement date is indicated above every map with format year-month-day. 56 4 Streamflow diel fluctuations reveal catchment- stream connectivity dynamics 4.1 Introduction Diel fluctuations of aquifer water table elevations and streamflow have been used to estimate evapotranspiration. This method provides an integrated assessment and, therefore, has often been preferred to point measurement alternatives such as sapflow measurements (e.g., [Granier et al., 2000]) and eddy correlation (e.g., [Baldocchi et al., 1997]), which are difficult to extrapolate to larger areas or apply in complex terrain. White [1932] introduced a method to relate diel fluctuations of the aquifer water table to evapotranspiration. It was later modified to account for differences in bed sediments and soil texture [Loheide et al., 2005; Butler et al., 2007] but the essence of the method remained unchanged. These methods provide a quantification of vegetation water uptake integrated over the zone of influence of piezometers or wells under study. Traditionally, water table records originating from the same study sites are compared to one another to evaluate the spatial variability of vegetation water uptake [Butler et al., 2007]. In a similar fashion, evapotranspiration has been quantified through the magnitude of diel streamflow fluctuations. Evapotranspiration is equated to the “missing streamflow” computed either as the difference between the time series connecting all the daily streamflow maxima and the recorded streamflow time series [Bond et al., 2002], or as the difference between every local daily streamflow maximum and the recorded streamflow time series [Boronina et al., 2005]. Bren [1997] first hypothesized that streamflow diel fluctuations reflect “deep slope processes” or, in other words, the connection between streams 57 and their catchments via subsurface flow pathways. His study resulted in the rejection of this hypothesis as he demonstrated that upslope vegetation water uptake had little influence on streamflow diel fluctuations and the zone of influence was limited to the riparian zone alone. A number of studies followed that did not aim at quantifying evapotranspiration from the study of streamflow or aquifer water table records. Rather, they sought insights into catchment-scale or regional behaviour from comparisons of diel fluctuation characteristics through time and space [Lundquist and Cayan, 2002; Czikowsky and Fitzjarrald, 2004]. The main characteristics of diel fluctuations — namely, amplitude and the timing of troughs and peaks — were used as indicators reflecting the dominant drivers of diel streamflow fluctuations. However, because no previous study used nested streamflow gauges within a single catchment, generation and propagation of evapotranspiration perturbations at multiple scales within a catchment could not be addressed. A number of conceptual models have been developed to explain streamflow diel fluctuations. In one model, lateral water contributions originate on slopes and propagate downslope to reach the channel. Upon reaching the riparian zone, the upslope signal is perturbed by vegetation water uptake, which introduces a daily fluctuating component (e.g., [Bren, 1997; Bond et al., 2002]). As water deficit increases throughout the catchment, the hydrologic connectivity between upslope area and the stream is maintained only through deeper, longer hillslope flow pathways, causing an increased time lag between the lateral slope contribution perturbed by riparian vegetation water uptake (excitation signal) and the streamflow record at the downstream gauge (response signal) [Bond et al., 2002]. In the model proposed by Bren [1997], water fluxes are always directed toward the stream; only their magnitude fluctuates during the day and shows a general seasonal decrease as a catchment’s water deficit increases. Bond et al. [2002] also allowed for evapotranspiration perturbations to reach the stream through hyporheic flow pathways. The hyporheic zone is the saturated zone surrounding the channel where stream water mixes with water from shallow riparian aquifers along hyporheic flow pathways and ultimately discharges 58 back into the stream. Because field observations were not consistent with the conceptual model of Bond et al. [2002], Wondzell et al. [2010] proposed an alternative model, in which riparian vegetation taps directly into the hyporheic zone, thus increasing movement of stream water into the hyporheic zone (downwelling). During periods of low or negligible evapotranspiration, typically at night, reversed head gradients drive hyporheic water back to the stream (upwelling) or at least decrease the magnitude of hyporheic recharge via stream water infiltration. However, their data did not support either the original or their alternative model, and they concluded that local processes occurring within the riparian zone were insufficient, on their own, to explain the spatial pattern of diel fluctuations within the well and piezometer network. Moreover, they suggested that processes occurring at larger scales, such as that of the whole stream network or adjacent hillslopes, should be considered in relation to water table and streamflow diel fluctuations. Szilagyi et al. [2008] also suggested that both “local” and “overall” head gradients drive water fluxes within the shallow aquifer and to the stream. Their implementation of a mechanistic model reproduced observations from Gribovszki et al. [2008], in which diel minima of water levels in riparian wells lagged “anomalously” behind those of streamflow. In this study, it is hypothesized that streamflow diel fluctuations reflect the linkages between slope water delivery, riparian vegetation uptake, hyporheic exchange and in-channel streamflow, and should be investigated to deepen our understanding of hydrological processes controlling hydrologic connectivity within catchments. Traditionally, streamflow diel fluctuations have been observed at a single location with a focus their seasonal variations. The present study is based on the premise that comparison of diel fluctuations at multiple sites along a stream network can shed new light on these processes. A conceptual model of propagation of diel perturbations from upslope areas to the stream is presented. Within this framework, the characteristics of diel fluctuations at nine stream gauges nested within CCEW are explored. Two models – a simple, semi-distributed one- dimensional advective model and a lumped analytical model – were formulated and applied to assist in separating the effects of in-channel 59 processes (advection and dispersion) from those related to hyporheic exchange and lateral inflow. More specifically, the following questions are addressed: (1) What are the propagation mechanisms of diel fluctuations generated by evapotranspiration processes? (2) What can be learnt on catchment-stream connectivity from the analysis of these diel fluctuations? 4.2 Methods 4.2.1 Conceptual model of propagation of diel perturbations Water was assumed to be delivered from the upslope areas, free of perturbations, transits through the riparian aquifer to reach the stream and travels downstream. Streamflow diel fluctuations were assumed to be generated by evapotranspiration in the riparian aquifer. Evidence to support this assumption will be discussed in section 4.3.1. This model, represented in Figure 4.1, was designed to focus on 1) the generation and propagation of perturbations within the riparian aquifer and to the stream and 2) the propagation of perturbation down the stream. 4.2.2 Selection of time periods Three periods, roughly one month apart and representing distinct flow regimes along the overall recession period during summer 2005, were selected based on meteorological records (Table 4-1). The start of each period was chosen to be at least two days following the end of any prior rain event. Periods were also chosen on the basis that they had relatively constant diel fluctuations of air temperature. These periods will be referred to as high, intermediate and low baseflow periods. The probability of exceedance of discharge at the start of each period was equal to 51%, 82% and 87% for high, intermediate and low baseflow periods, respectively (Figure 4.2). 60 4.2.3 Streamflow response time to climate forcing Bond et al. [2002] used sapflow records to quantify evapotranspiration demand. However, no such measurements were available in CCEW and vapour pressure deficit (VPD) was used instead as suggested by Wondzell et al. [2010]. VPD was computed from air temperature and relative humidity records measured at the lower meteorological station (Figure 2.1). The climate station was no farther than 2000 m from any portion of the channel and at an elevation of 1390 m a.s.l., while the entire stream network lies between 1230 and 1800 m a.s.l. Cross-correlation functions were computed for the three selected time periods between the VPD time series and the streamflow records at each of the nine stream gauges. Lag increments were set equal to one hour and the overall lag ranged between zero and twenty-three hours. Since an increase in VPD is accompanied by an increase in vegetation water demand and a subsequent lagged streamflow decrease, the streamflow response time to evapotranspiration was defined to be the time lag associated with minimum (most negative) correlation. 4.2.4 Extraction of diel fluctuation characteristics Streamflow time series were processed to extract three characteristics: mean daily discharge; daily fluctuation amplitude, equal to the difference between daily peak (maximum) and trough (minimum); and time of daily trough. A Kolmogorov-Zurbenko (KZ) filter was applied to VPD and streamflow time series to isolate diel fluctuations The KZ filter involves iterative applications of a moving average window (e.g., [Zurbenko, 1986; Rao and Zurbenko, 1994; Fleming et al., 2009]). It is based on two parameters: the half-window size q, set to be equal to 10 h; and the number of iterations, K, set to be equal to 2. A KZ filter is low-pass and extracts signal components with periods greater than the cutoff value (Pc): 61 ௖ܲ ൌ 2 · ݍ · √ܭ   4‐1 The chosen values for q and K resulted in a cutoff period equal to about 28 h (half period equal to 14.14 h). Residuals with periods smaller than Pc were extracted from the difference between the initial time series and the filtered one (Figure 4.3). Mean daily values for VPD and streamflow were extracted from the unfiltered records, while amplitudes and times of trough were extracted from the time series of residuals. Time series of residuals were further smoothed with a binomial filter in order to remove spurious artefacts associated with measurement noise in the stage records. A binomial filter is a discrete approximation of a Gaussian filter parameterized with one variable, its order. The order of the binomial filter determines its range and smoothing effect. In the present study, it was set equal to 40. This represents a temporal window of 10 h since streamflow was recorded every 15 min. A time wrap transformation was applied to the time series of times of daily minimum discharge so that the transformed values account for the fact that the time of a trough occurring at 23:57 should not be very different from that of a trough occurring at 00:03 the next day. Preliminary analysis of the data indicated that daily times of minimum discharge occurred between 15:00 and 7:00. Therefore, times of trough occurring between 15:00 and 24:00 were mirrored to 24:00 in order to provide a temporal continuum (Equation 4-2). This transformation was preferred to a sinusoidal alternative due to its linear nature. It can be formally represented as follows: ௖ܶሺݐ௖ሻ ൌ ቐ ݐ௖                ݐ௖ א ሾ24: 00, 7: 00ሿ  ݐ௖ െ 24     ݐ௖ א ሾ15: 00, 24: 00ሾ          4‐2  where tc, and Tc are the characteristic time and its transformation, respectively. 62 4.2.5 Semi-distributed advective model of diel perturbations Data required to compute wave celerities were only available for the period between August 31, 2006 and September 09, 2006. Therefore, model simulations were only conducted during this fourth period. Conceptual and mathematical formulation Similar to the approach used by Wondzell et al. [2007], a Lagrangian one-dimensional advective model that accounts for lateral inflows distributed along the whole channel network and downstream attenuation of signal amplitude was developed. Distributed wave celerities and lateral inflows along the whole network were required as inputs to the model. In addition, wave attenuation due to in-channel energy losses had to be accounted for as they may play an increasingly significant role with increasing stream network length: the phase and amplitude of the water wave subject to advection down the stream network can interfere with the phase and amplitude of lateral inflows distributed along the stream network. These interferences may enhance or dampen the wave depending on the time, magnitude and location down the channel and, hence, affect the coherence of the streamflow signal at the downstream gauge [Wondzell et al., 2010]. The water wave resulting from diel variations in evapotranspiration is represented in terms of the spatial and temporal distribution of streamflow, denoted Q(t, i), where t is time and i represents nodes located along the channel network. Time is discretized into increments of magnitude ∆ݐ. As described below, the node spacing was set at ∆ݐ · ܿ, where c is the wave celerity for the reach, so that the wave was advected one node downstream in each time step. The upstream boundary condition, Q(t, 1), was set equal to streamflow recorded at a weir. Calculations progressed downstream from node to node. At each node, the first step was to lag the hydrograph from the upstream node by an interval Δt : 63 ܳ௟௔௚ሺݐ, ݅ሻ ൌ ܳሺݐ െ ∆ݐ, ݅ െ 1ሻ 4‐3 This lagged signal was then smoothed with a Gaussian filter to simulate wave attenuation by in-channel energy losses. The Gaussian filter is parameterized by s, the standard deviation and np, the spatial window over which the filter is truncated: ܩܨሺ߬, ݏሻ ൌ ଵ √ଶగ·௦ ݁ି ഓమ మೞ   4‐4 where τ is the number of nodes up and down the channel to the wave location of interest and was allowed to vary between 1 and (np-1)/2. Formally, the filtering is a convolution in time: where ܳ௟௔௚ሺݐ, ݅ ሻ is the input (lagged) function, which is convoluted with GF to produce ܳ௦ሺݐ, ݅ሻ, the lagged and smoothed discharge that represents the effects of advection and dispersion. The final step in the calculation at each node was to add the local lateral inflow into the reach: where ݍሺݐ, ݅ሻis the lateral inflow into reach i at time t. ܳ௦ሺݐ, ݅ሻ ൌ ෍ ܩܨሺ߬, ݏሻ ఛୀ ௡೛ିଵ ଶ   ఛୀି  ௡೛ିଵ ଶ   ·  ܳ௟௔௚ሺݐ െ ߬, ݅ሻ   4‐5  ܳሺݐ, ݅ሻ ൌ ܳ௦ሺݐ, ݅ሻ ൅ ݍሺݐ, ݅ሻ 4‐6  64 Discretization of the stream network The stream network was divided into seven reaches of length Li. Values of Li were derived from a GIS analysis. Standard salt dilution gauging protocol was followed to measure discharge [Day, 1976; Moore, 2005] at both ends of every reach on September 05, 2006. Net lateral contributions integrated over the length of every reach were computed by taking the difference between outflow and inflow of every reach. On September 1-2, 2006, reach-scale tracer slug injections were carried out for every reach. Further details of the data collection protocol are presented in section 5.3.2. The time to peak of the breakthrough curve was assumed to reflect a representative water velocity over the length of the reach. Ultimately, every reach was characterized with a net lateral inflow and a representative velocity (Table 4-2). In the context of propagation of perturbations down the channel, wave celerity (c) is meaningful; water velocity (v) relates to the travel time of water molecules, not a wave. Under the assumptions of a channel with rectangular cross-section and shallow depth, wave celerity is approximately 1.5 times the velocity [Henderson, 1966]. The extent to which these assumptions have been met at CCEW is unknown. Nonetheless, this relation was used to derive a representative wave celerity value for every reach. Reach lengths were chosen so that every adjusted Li was a multiple of ܿ௜ · ∆ݐ. Doing so ensured that reach boundaries were aligned with space increments ܿ௜ · ∆ݐ and artificially increased the stream network length by only 1.2%. The whole stream network is therefore represented as a one- dimensional sequence of nodes that are ܿ௜ · ∆ݐ meters apart. Distributed lateral inflows Every lateral inflow contribution, q(i, t) was defined as the sum of a fixed term qf(i, t), reflecting water slope contribution before being altered by evapotranspiration, and a variable term qv(i, t), reflecting evapotranspiration uptake. The term qf(i, t) was constrained by measured reach-scale streamflow increments (Table 4-2). The time-varying term qv(i, t) was set to be a function 65 of VPD. For a number of tree species, a quasi-linear relation can be assumed between transpiration rate and VPD. For example, a quasi-linear relation holds for Western Red Cedar for VPD values ranging from 0 to 4 kPa [Kavanagh et al., 2007]. For trembling aspen, the dominant riparian understory species in CCEW, a linear relation between transpiration rate and VPD can be assumed for VPD values lower than 1 kPa [Hogg and Hurdle, 1997]. During the simulation period, VPD was lower than 1 kPa, 81% of the time, and never greater than 1.49 kPa. Therefore, qv(i, t) was parameterized as being proportional to VPD. The coefficient of proportionality, k, is a parameter of the advective model. The resulting expression of the lateral contributions is as follows: The total travel time of any perturbation from catchment locations that are hydrologically connected to the stream is the sum of a travel time along the flow pathway that connects the source of the perturbation to the channel (typically, along slopes and through riparian areas) and a travel time along the channel, from the perturbation channel entry location to the downstream gauge (Figure 4.1). Perturbation travel time from its source to the point of channel entry was not implemented in the advective model; instead, perturbations were lumped along the slope or riparian flow pathways and applied directly along the channel. Lateral perturbations were delivered to the stream at every node of the stream network. Sensitivity analysis Several simulation runs were conducted to assess the sensitivity of the model output to a change in parameter values. For instance, the amplitude of diel streamflow fluctuations of the modelled outflow depends on the coefficient of proportionality, k, between VPD and ET. It also depends on the smoothing effect of the Gaussian filter. A model parameterized with a large k ݍሺݐ, ݅ሻ ൌ ݍ௙ሺݐ, ݅ሻ ൅ ݇ · ܸܲܦሺݐ, ݅ሻ 4‐7  66 simulates diel fluctuations with large amplitudes; a model parameterized with a large Gaussian filter span, np, yields simulated diel fluctuations with dampened amplitudes. 4.2.6 Analytical model of linkages between slope, riparian zone and stream The differences in response time to climatic forcing among stream gauges were examined more precisely during the fourth period (August 31, 2006 – September 09, 2006), using a lumped analytical model, which is based on the assumption that diel streamflow fluctuations arise purely due to transpirative extraction of water from the riparian aquifer and the resulting decrease in discharge from the riparian aquifer to the stream. The model was formulated as follows: i. Water is delivered from the hillslopes to the riparian aquifer at a rate ܪሺݐሻ ൌ ܾ െ ܿ · ݐ, where b [l s-1] and c [l s-1 s-1] are parameters, and t is the time [s]. A linear decline should be a reasonable approximation of streamflow recession for periods of a few days. ii. Evapotranspiration draws water from the riparian aquifer with a diurnal fluctuation given by ܧܶሺݐሻ ൌ ݀ · sin ሺ߱ · ݐ  ൅  ߠሻ, where ߱ ൌ 2ߨ/ܶ [s-1] is the angular frequency , with T, the period, set equal to one day, θ [-] is a phase lag, and d [l s-1], is a parameter. iii. The riparian aquifer, with storage denoted as S(t) [l], discharges water to the stream at a rate equal to aS(t), where a [s-1] is a parameter. That is, the riparian aquifer behaves like a linear reservoir with characteristic time, 1/a. Note that t is, in fact, the time elapsed since a time origin t0. This time origin t0 was chosen so that θ=0, thus simplifying the expression of ET(t) and the integration of Equation 4-8. Given the three assumptions provided above, the water balance for the riparian aquifer can be expressed as: 67 Equation 4-8 can be solved to yield the following expression for the riparian aquifer: where G [l] is a constant of integration. After transformations, the expression of the riparian aquifer’s contribution to streamflow is: where φ is a phase lag, given by ߮ ൌ ܽݎܿݐܽ݊ሺെఠ ௔ ሻ. The analytical model does not account for a lag time between the source of the perturbation and its point of entry in the channel. This approximation is based on the assumption that perturbation travel time along flow pathways in slopes and riparian aquifers is negligible compared to the characteristic time, 1/a, of riparian aquifers. Although this assumption may not be strictly correct, the interpretation of the role of riparian aquifer dynamics in lagging ET perturbations to the downstream gauge remains unchanged, as presented in section 4.4.4. ݀ܵሺݐሻ ݀ݐ ൅ ܽܵሺݐሻ ൌ ܾ െ ܿ · ݐ ൅ ݀ · sinሺ߱ · ݐሻ  4‐8  ܵሺݐሻ ൌ 1 ܽ ቂቀܾ ൅ ܿ ܽ ቁ െ ܿ · ݐቃ ൅ ݀ ܽଶ ൅ ߱ଶ ሾܽ · sinሺ߱ · ݐሻ െ ߱ · cosሺ߱ · ݐሻሿ ൅ Geି௔௧  4‐9 ܽܵሺݐሻ ൌ ቂቀܾ ൅ ܿ ܽ ቁ െ ܿ · ݐቃ ൅ ܽ · ݀ ܽଶ ൅ ߱ଶ sinሺ߱ · ݐ ൅ ߮ሻ ൅ Geି௔௧  4‐10 68 4.3 Results 4.3.1 Whole network assessment of lagged streamflow responses VPD generally decreased from high to intermediate to low baseflow periods. Mean values computed for each period were equal to 0.75 kPa, 0.46 kPa and 0.28 kPa for high, intermediate and low baseflow periods, respectively (Figure 4.4). The temporal variability of VPD records also exhibited a decrease from high to low baseflow period, indicating that the amplitude of diel fluctuations decreased from high to low baseflow periods. Time series of residuals following KZ filtering were used instead of the unfiltered time series to isolate diel fluctuations. Cross-correlation functions between VPD time series and all nine stream gauge records for all three selected periods are shown in Figure 4.5. In general, streamflow was least cross-correlated with VPD during the high baseflow period and most cross- correlated during the low baseflow period. Response time ranged from five to sixteen hours. The median response time increased with decreasing flow regime, ranging from 11, 12 and 14 hours for high, intermediate and low base flow periods, respectively. The variability of response time among stream gauges decreased with decreasing flow regime: the inter-quartile range (IQR) was equal to 4, 3 and 1 hour for high, intermediate and low base flow, respectively (Figure 4.6). Response time did not correlate significantly with catchment contributing area for any of the selected periods: correlation coefficients were equal to 0.015, 0.4 and 0.57, and p-values were equal to 0.97, 0.28 and 0.11, for high, intermediate and low base flow, respectively (data not shown). A corollary of this lack of correlation is that there was no consistent increase in response time going downstream. For instance, during the high baseflow period, streamflow diel fluctuations recorded at gauge EL07 lagged five hours behind those recorded at gauge EL01, located 1540 m downstream along Elk Creek (Figure 3.1). This characteristic, although reduced, persisted throughout the intermediate baseflow period. For any of the three periods, there was no significant discrete source of lateral flow between these two streamflow gauges which signal could have drowned that 69 recorded at the upstream gauge EL07. A similar behaviour was observed along Cotton Creek. During the intermediate baseflow period, the streamflow at gauge CN10 (outlet of the headwater Cotton North) lagged behind that at the downstream gauge CO01 (outlet of the Cotton Creek sub-basin) (Figure 4.5b). Topographic variables describing the characteristics of each reach and its contributing slopes were available from a terrain analysis (section 3.2.2). Correlation analysis revealed no statistically significant relations between response time and maximum channel length, maximum slope length and mean slope gradient, each derived at the scale of the appropriate incremental sub-basin. 4.3.2 Characteristics of streamflow diel fluctuations General characteristics A sample of streamflow diel fluctuations is presented in Figure 4.7 for the high baseflow period (July 18 – August 05, 2005). The amplitude and shape of diel streamflow fluctuations varied through time, even within periods. At most sites, diel fluctuations were characterized by steep falling limbs and progressive rising limbs. Amplitude variations through time and space are highlighted in Figure 4.8. Patterns of amplitude for a set time period were generally variable. At high baseflow, the relation between amplitude and catchment scale was strongly positive but weakened at intermediate and low base flow periods. Similar trends were observed in the relation between the standard deviation of diel amplitude and catchment scale (Figure 4.9). Streamflow minima occurred between 15:00 and 7:00, mainly between 20:00 and 4:00 (Figure 4.10). They generally occurred earlier during the day during the high baseflow period and later during the low baseflow period. Analysis of variance supported the fact that differences in the timing of daily minima among periods were statistically significant (p-value < 0.001). There was no clear relation between the timing of daily streamflow minima and 70 maximum channel length (Figure 4.11). Singularities of observed streamflow time series The reach outflow (gauge EL01) observed during the period over which simulations were run showed double peaks in September 4-6, 2006 (Figure 4.12). A more comprehensive visual exploration of the dataset revealed that this transient phenomenon occurred at other stream gauges and would never last more than few days. A similar behaviour was that of the observed reach inflow (gauge EL07) which presented a “bend” in the daily rising limb during various days throughout the simulation period. This bend may be the result of the addition of two peaks, close enough in time that the overlap is significant and the peaks are no longer discernible. 4.3.3 Simulations from the one-dimensional advective model Distributed water wave celerities were only available for the lower section of Elk Creek (reaches EL01-EL09, Figure 3.1). Attempts were made to derive celerity values for the whole CCEW stream network in order to run simulations over an extended stream network. However, from the available sample of celerity values, celerity correlated only mildly with discharge (r2=0.56) and poorly with channel gradient (r2=0.07). Wave celerity depends on complex interactions between channel morphology, channel gradient, streambed roughness and water level. The use of predicted celerity values to run simulations throughout the whole stream network was judged too uncertain and simulations were run only for the stream reach that stretched between EL07 and EL01 (Table 4-2). Sensitivity analysis Simulation results are presented in Figure 4.12. The smoothing effect of the Gaussian filter depends on the standard deviation, s. The length of the span of the filter, np, was adjusted so that it extended three times the value of s, both upstream and downstream the water wave. Doing so ensured that the 71 filter was not truncated (the value of the weights at both extremities of the filter was close to zero). The model parameters, k and s, were calibrated manually. The parameter k [l s-1 m-1 kPa-1] is the coefficient of proportionality between VPD forcing and lateral inflow per unit reach length. The modelled reach outflow in Figure 4.12a is characterised by a light Gaussian smoothing (s = 20 min, np = 2 h) and a small effect of ET forcing on streamflow (k = 0.0046 l s-1 m-1 kPa- 1). The amplitude of the diel fluctuations of the modelled reach outflow originated from the reach inflow EL07, with little downstream attenuation. The modelled reach outflow in Figure 4.12b is characterised by a strong Gaussian smoothing (s = 3 h, np = 18 h) and no effect of ET forcing on streamflow (k = 0 l s-1 m-1 kPa-1). The amplitude of the diel fluctuations of the modelled reach outflow originated from the reach inflow EL07, with a strong downstream attenuation. Lateral contributions along the reach were restricted to their non-varying term, qf(t,i) and did not generate diel fluctuations in the modelled reach outflow. The modelled reach outflow in Figure 4.12c is characterised by a strong Gaussian smoothing (s = 3 h, np = 18 h) and a strong effect of ET forcing on streamflow (k = 0.24 l s-1 m-1 kPa-1). The amplitude of the diel fluctuations of the modelled reach outflow originated mainly from the lateral contributions to the stream and, to a lesser extent, from the reach inflow of stream gauge EL07, with a strong downstream attenuation. Model fit The overall fit between observed and simulated outflow time series was poor regardless of the simulation run. Observed and simulated outflow time series were clearly out of phase, with a near 180° phase shift for all simulation runs. In the present case, the inflow time series lagged behind the observed outflow time series, a peculiar behaviour described more extensively in section 4.3.1. Since the advective model formulation could only account for a time lag due to downstream wave travel time, it was impossible to simulate the observed outflow time series successfully. Further, the model did not 72 reproduce the double peaks and \"bends\" in the diurnal patterns that were observed on some days. 4.3.4 Results from the analytical model Manual calibration of the analytical model was difficult due to the relatively large number of parameters and its use was eventually limited to a sensitivity analysis to parameter 1/a, the characteristic time of the discharge rate from the riparian aquifer to the stream. Equation 4-9 reveals that the phase of the sinusoidal component of aS(t) is modulated by a. The influence of a on the phase of the riparian aquifer discharge to the stream is illustrated in Figure 4.13. For values of 1/a ranging from 15 min to a day, the phase of the lateral riparian contribution to the stream shifted in time by about 6 h. In other words, a poorly responsive riparian aquifer, characterised by a large characteristic time 1/a, will delay ET perturbations, prior to transmitting them to the stream. Inversely, a responsive riparian aquifer will transmit ET perturbations to the stream with nearly no phase shift. The lateral riparian contribution with phase φ is then added to the water wave transported down the channel with phase Φ. A difference of phase |Φ െ ߮| between 0 and ߨ will result in an increased amplitude of the channel wave, while a difference of phase |Φ െ ߮| between ߨ and 2ߨ will result in a decreased amplitude of the channel wave. 4.3.5 Area of influence The area of influence is the area encompassing all sources of perturbations that are hydrologically connected to the stream. It is computed as the ratio between the “missing streamflow” and the transpiration rate per unit ground area. Missing streamflow was computed as the cumulative difference between the interpolated daily streamflow maxima and the actual streamflow on a time step of 15 min. following Bond et al. [2002] (Figure 4.14). The resulting daily averaged missing streamflow was found equal to 41.5 m3 day-1, equivalent to an averaged transpiration outflow of 0.48 l s-1. In 73 the absence of data describing tree species and stem density, the slope of the linear relation between VPD and transpiration rate averaged at the stand level was bracketed between 0.3 and 3 mm day-1 kPa-1 based on the literature [Hogg et al., 1997; Moore et al., 2004]. This resulted in daily averaged transpiration rate per unit ground area ranging between 0.25 and 2.5 mm m-2 day-1. Ultimately, the resulting area of influence ranged between 16600 m2 and 166000 m2. Considering a stream buffer which width has been set equal to 80 m, this represents a zone of influence that ranges from 208 to 2075 m in length. 4.4 Discussion 4.4.1 What is the driver of streamflow diel fluctuations in CCEW during the study period? Both snowmelt and evapotranspiration have been reported to generate diel streamflow fluctuations. At CCEW, the snowpack completely disappeared at the upper climate station (1390 m a.s.l.) around April 27, 2005 [Jost et al., 2007], and we are confident that all snow had disappeared prior to the first (high baseflow) period. Moreover, the steep falling limbs and progressive rising limbs as observed in Figure 4.7 are the signature of diel streamflow fluctuations caused by evapotranspiration [Lundquist and Cayan, 2002; Czikowsky and Fitzjarrald, 2004]. These observations support the fact that evapotranspiration was the only driver of diel streamflow fluctuations during the periods of interest. 4.4.2 Amplitude of streamflow diel fluctuations Decrease with baseflow recession The decrease in the amplitude of streamflow diel fluctuations with baseflow recession (from high to intermediate to low baseflow periods) may reflect a progressive decrease in soil water availability for trees on slopes and in the riparian zone during the recession flow. Wondzell et al. [2009] found 74 that the rate of average daily drawdown of water table elevations in riparian wells during the baseflow recession was likely too slow to cause the water table to fall below the root zone over that period, and significantly affect vegetation water uptake. This was supported by the fact that the amplitude of diel fluctuations of water table elevations throughout their network of riparian wells increased with decreasing baseflow, suggesting that rates of evapotranspiration did not decrease. In the same study, the increase of the amplitude of diel fluctuations of water table elevations with baseflow recession was accompanied by a decrease of the amplitude of streamflow diel fluctuations. A similar decrease in the amplitude of streamflow diel fluctuations was observed at CCEW (Figure 4.8 and Figure 4.9a), which indicates that the mechanisms that operate in CCEW may be similar to those inferred by Wondzell et al. [2009]. However, no well or piezometer data were available to confirm this hypothesis. An alternative mechanism may be that the marked decline in mean VPD and in the amplitude of diel fluctuations of VPD with baseflow recession (Figure 4.4) reflects a decrease in meteorological forcing through time, which in turn may induce a decrease of ET rates and a subsequent decrease of the amplitude of streamflow diel fluctuations through time. Increase with catchment scale A surprising result is the increase of the amplitude of streamflow diel fluctuations with catchment scale. One would expect that as catchment scale increases, more sources of ET perturbations are sampled, advected in the stream, and aggregated at the downstream gauge. First, this aggregation increases the likelihood of obtaining desynchronised perturbations, which would reduce the amplitude of diel fluctuations at the downstream gauge. Second, larger spatial scales provide more opportunity for hydrologic connectivity through flow pathways to shift at the daily and seasonal time scales, depending on local soil water availability and ET demand. This, in turn, may increase the temporal variability of the amplitude of streamflow diel fluctuations, as illustrated in Figure 4.9b. 75 Alternatively, the amplitude of diel fluctuations normalized by mean daily streamflow could be examined. Normalized amplitudes decreased with catchment scale (data not shown), confirming results from Lundquist and Cayan [2002] and Czikowsky and Fitzjarrald [2004]. We believe, however, that this decrease is due primarily to the increase of the normalizing term – the average daily streamflow – with catchment scale, rather than to the decrease of non-normalized amplitudes of streamflow diel fluctuations. Therefore, the relation between normalized amplitude of streamflow diel fluctuations and catchment scale may tell little about the linkages between meteorological forcing, vegetation and stream. 4.4.3 Streamflow response time to climatic forcing Streamflow response times to climatic forcing varied substantially among streamflow gauges and through time but are comparable to the response times varying from four to nine hours with baseflow recession that was observed at the stream gauge of a 100 ha catchment [Bond et al., 2002]. After collecting further data in the same catchment, Wondzell et al. [2007] demonstrated that response time increased with decreasing flow velocity. The results of the present study confirm their findings as median response time among all nine streamflow gauges increased from high to intermediate to low baseflow. Moreover, the time of daily minimum streamflow occurred significantly earlier at high baseflow than at low baseflow. In contrast to the increase in mean response times with decreasing flow, the response time variability among stream gauges decreased from high to intermediate to low baseflow regimes, possibly due to shifts in flow pathways that are hydrologically connected to the stream (Figure 4.6). At high baseflow, travel times along shallow flow pathways that depend on local surface topography and near-surface preferred pathways control the streamflow response time to climate forcing. At low baseflow, slope water contributions are only being delivered through deeper flow pathways that may be more homogeneously distributed throughout the whole catchment. Similar inferences were drawn from the analysis of baseflow records among 76 22 catchments ranging in size from 1 to 533 km2 by Brutsaert and Lopez [1998], who identified two successive recession behaviours: a relatively wet short-time response followed by a relatively dry long-time response, each characterised by a recession coefficient. The recession coefficient for the short-time response was much more variable among gauges (coefficient of variation = 288%) than its long-time response counterpart (coefficient of variation = 53%). This may indicate that shallow flow pathways that only contribute to a catchment’s hydrological response during the short-time response are more spatially variable than deeper flow pathways through which slope water reaches the channel during the long-time response. It is hypothesized that a similar behaviour occurs in CCEW. 4.4.4 Origins of streamflow diel fluctuations Simulation results demonstrated that the advective model was capable of reproducing the dampening of diel streamflow fluctuations by downstream energy loss, leading to weak diel streamflow fluctuations at the downstream gauge (Figure 4.12b). They further showed that diel fluctuations at the downstream gauge may originate 1) mainly from the upstream gauge in reaches with minimal in-stream energy loss, or 2) mainly from nearby lateral inflows in cases where the upstream signal (also including upstream lateral inflows) was effectively dampened (Figure 4.12a and Figure 4.12c). There are also alternatives to these two extreme cases, in which in-stream perturbations advected downstream compete with perturbations originating from lateral contributions. In case 1), streamflow diel fluctuations would be generated at the channel heads and would be advected downstream with minimal energy loss. This would induce a systematic increase in response time to evapotranspiration (as indexed by VPD) with catchment contributing area or maximum stream length. However, for each of the periods, no clear relation was found between observed streamflow response time to evapotranspiration and catchment size. Moreover, daily times of minimum streamflow did not vary systematically with maximum channel length. Therefore, differences in 77 streamflow response time to evapotranspiration throughout the whole catchment cannot be explained by in-channel wave travel time alone; the hypothesis that in-channel wave travel time controls the streamflow response time to evapotranspiration is therefore rejected. In case 2), perturbations recorded at the downstream gauge would only be generated by lateral contributions in the close vicinity of the stream gauge. However, the length of the stream buffer of influence derived from mass balance calculations was estimated to range from about 200 m to 2 km. Although this result should be treated with caution, considering the uncertainty that characterizes the type of vegetation and its density in CCEW, it suggests nonetheless that perturbations from sources that are relatively far from the stream gauge may be advected downstream or downslope and generate diel fluctuations at that stream gauge. Therefore, it is unlikely that diel fluctuations recorded at various stream gauges originate solely from nearby sources. The complex and temporally variable shape of diel fluctuations (Figure 4.7 and Figure 4.12) suggests that streamflow diel fluctuations result from the concurrent expression of a number of processes. Preliminary results from the analytical model showed that the characteristics of a single riparian aquifer could greatly influence the phase of streamflow diel fluctuations. This coupling between riparian aquifer and stream could explain why there is no systematic increase of streamflow response time to climatic forcing with catchment contributing area. Depending on the topography of the slope- riparian zone transition and the extent of the riparian zone, several independent riparian aquifers may exist along the channel, each with its own characteristic time. Moreover, and in contrast with the results presented in Figure 4.13, the storage of riparian aquifers varies through time, as a function of upslope inflow, evapotranspiration and outflow to the stream, which further alters the amplitude and phase of perturbations that are transmitted to the stream. The occurrence of dual peaks or “bends” in the rising limbs of streamflow diel fluctuations may simply be the result of two or more competing processes with distinct phases that originate either from in-stream 78 advection of perturbations or from riparian aquifers. Water table elevations in the riparian aquifers and in the stream control the direction of head gradients and associated fluxes between stream and riparian zone. The daily to seasonal variations of these head gradients likely explain the transient nature of some characteristic patterns in the streamflow diel fluctuations (such as the double peak, for instance). 4.5 Conclusion From an analysis of patterns of diel fluctuation characteristics both in time and space, it was found that variations in streamflow timing along the channel network cannot be explained solely by in-stream advection and dispersion. Conversely, the zone of influence of transpiration losses was too large to support the hypothesis that streamflow diel fluctuations originate solely from neighbouring riparian aquifers. Both the numerical and analytical models provided insight into the possible origins of perturbations that are advected to the downstream gauge. While in-stream advection processes can only explain an increase in response time of streamflow fluctuations to climatic forcing with increasing stream length, the dynamics of riparian aquifers, as controlled by the storage-discharge relation, can explain a larger response time at an upstream gauge. This study demonstrated that the information contained in the diel fluctuations of a nested network of stream gauges may reveal novel aspects of the hydrologic connectivity between stream and riparian aquifers. Future studies would benefit from continuous monitoring of riparian water table elevations at multiple locations along the channel. This would shed some light on the synchroneity between the stream and the riparian aquifers, and help delineate independent riparian aquifers. The analytical model presented in this study could be extended with the implementation of multiple riparian aquifers distributed along the channel. Riparian well data could also help constrain the parameterization of the fluxes between stream and riparian aquifers. 79 Table 4-1 : Periods selected for extraction of diel fluctuation characteristics. Period Start End High baseflow 18-07-2005 05-08-2005 Intermediate baseflow 24-08-2005 08-09-2005 Low baseflow 20-09-2005 29-09-2005 Table 4-2 : Characteristics of stream reaches as input to the advective model. Reach Length [m] Celerity [m h-1] Mean gradient [%] QOUT [l s-1] Net budget [l s-1] EL01 248 619 9.5 18.5 1.5 EL02 294 742 10.8 17 -0.2 EL03 263 673 6.5 17.2 0.9 EL04 355 605 10.0 16.3 0.9 EL05 251 508 6.2 15.4 0.5 EL06 256 490 5.7 15 0.7 EL07 82 587 15.9 14.2 1.7 F F D igure 4.1 : C igure 4.2 : ischarge va onceptual Baseflow pe lues are plo model of ge riods selec tted as a lo neration an ted in year garithmic s d travel of 2005, illus cale. diel perturb trated at g 80 ations. auge ES. 81 Figure 4.3 : Extraction of the diel fluctuations of vapour pressure deficit with a Kolmogorov- Zurbenko low-pass filter. Figure 4.4 : Distributions of hourly VPD for the three periods. Cross symbols represent the mean values. 82 Figure 4.5 : Cross-correlation coefficients between VPD and streamflow records for the three periods. The nine stream gauges are separated between three graphs for clarity. (a) Stream gauges of the whole basin and the two main sub-basins. (b) Stream gauges within the Cotton Creek sub-basin. (c) Stream gauges within the Elk Creek sub-basin. 83 Figure 4.6 : Response time distribution among the nine gauges for the three periods. Date format for the x axis is day-month-year. 84 Figure 4.7 : Diel fluctuations for the high baseflow period. Streamflow records were filtered (KZ filter) and smoothed (binomial filter) to isolate diel fluctuation and eliminate spurious artefacts. Units are in l s-1. 85 Figure 4.8 : Amplitudes of diel fluctuation [l s-1]. 86 Figure 4.9 : (a) mean and (b) standard deviation of the amplitude of diel fluctuations as a function of catchment contributing area. Mean and standard deviation were computed for a set time period and a set streamflow gauge. Number of observations depended on the duration of the time period (19, 16 and 10 for the 1st, 2nd and 3rd periods, respectively). Contributing area values are plotted as a logarithmic scale. Smoothers are a lowess fit with span = 1. Figure 4.10 : Distributions of daily times of troughs streamflow. Each boxplot represents daily values across all streamflow gauges. 87 Figure 4.11 : (a) mean and (b) standard deviation of the transformed time of daily minimum flow as a function of maximum channel length. Mean and standard deviation were computed for a set time period and a set streamflow gauge. Number of observations depended on the duration of the time period (19, 16 and 10 days for the 1st, 2nd and 3rd periods, respectively). Channel length values are plotted in the logarithmic scale. Smoothers are a lowess fit with span = 1. 88 Figure 4.12 : Observed streamflow at stream gauges EL07 (inflow) and EL01 (outflow), along with modelled streamflow at EL01. VPD is plotted against the secondary axis. Each graph represents a simulation run with a distinct set of parameters. Parameter values are specified in text. Figure 4.13 : Predicted discharge from the riparian aquifer to the stream modulated by values of the time constant a.Vertical dashed lines are plotted every 6 h. 89 Figure 4.14 : Interpolated streamflow maxima at stream gauge EL01. The “missing streamflow” is computed as the cumulative difference between the time series of interpolated streamflow maxima and the time series of actual streamflow. 90 5 Influence of distributed flow losses and gains on the estimation of transient storage parameters from stream tracer experiments 5.1 Introduction Transient storage plays a key role in the transport and fate of solutes in streams. It encompasses hyporheic exchange and in-stream storage in pools or other regions of slow-moving water. Solute transport dynamics have often been characterized using tracer experiments at the reach scale. This technique involves the estimation of parameters in a transient storage model (TSM) to provide an optimal fit between the simulated tracer breakthrough curve and an experimental tracer breakthrough curve resulting from a controlled injection of dissolved tracer. Fitted model parameters reflect reach- integrated solute transport dynamics and are typically interpreted in relation to the magnitude of transient storage and its variation among streams and with stream discharge [Worman, 1998; Wondzell, 2006; Lautz and Siegel, 2007]. Exchange between the stream water and the different transient storage zones has been modelled with a first order mass transfer [Bencala and Walters, 1983], a one-dimensional diffusive flux [Jackman et al., 1984; Worman, 2000] or various probability density functions for residence time distribution such as power law or gamma distributions [Haggerty et al., 2000]. All models show similar results for shorter time scales, but power law or gamma residence time distributions perform better for longer time periods [Haggerty et al., 2000]. In recent years, many studies have modelled transient storage effects via the One-dimensional Transport with Inflow and Storage model structure (OTIS) [Runkel, 1998]. Due to their spatially integrated nature, reach-scale tracer experiments do not provide information on the sub-reach scale variability of solute transport dynamics. For instance, it is difficult to separate the effects of the 91 hyporheic zone from that of in-stream pools on solute retention. Recent studies have attempted to differentiate the magnitude and timing of exchange with in-stream pools versus the hyporheic zone using one or more of the following approaches: tracer experiments on bedrock reaches, where hyporheic exchange can be assumed to be negligible; tracer experiments at the scale of individual channel units; and by using hydrometric measurements to provide an independent estimate hyporheic exchange via Darcy’s law [Gooseff et al., 2005; Gooseff et al., 2008; Hall et al., 2002; Scordo and Moore, 2009; Stofleth et al., 2008]. In addition to the processes of advection, dispersion and transient storage, TSMs have to account for lateral inflow and outflow, especially for simulating nutrient transport. Water flow and solute mobilization usually originate on hillslopes with downslope flow and solute transport by surface and/or subsurface flow paths. Water and solutes then enter the stream directly or via the hyporheic zone, ultimately to be transported in the advective part of the stream. Lateral inflow is defined as the discharge from the hillslope (or deeper groundwater) to the stream and lateral outflow as the flow from the stream to deep infiltration into the groundwater and substratum or the flow through hyporheic flow paths that extend beyond the (spatial or temporal) domain of study. Mechanisms driving water delivery from or to the stream are highly variable in time and space across a broad range of landscapes and land uses [Genereux et al., 1993; Huff et al., 1982; Kuraś et al., 2008; Shaman et al., 2004]. The extent to which lateral inflow and outflow affect TSM parameterization has been widely overlooked. OTIS is most suited to the present study because its structure accounts for both lateral inflow and lateral outflow. A standard parameterization of lateral fluxes in OTIS first involves the computation of the reach water balance from the difference in measured discharge between the downstream and upstream boundaries of the experimental reach. Either lateral inflow or later outflow is then applied to the reach depending on whether it is net gaining or losing water, respectively. However, it is possible for a stream reach to alternate between gaining and losing segments over distances of tens to hundreds of meters 92 (e.g., [Story et al., 2003, Ruehl et al., 2006]). In this situation, the standard approach would not provide an accurate representation of the solute mass balance, potentially distorting the parameters representing downstream transport and transient storage. This chapter assesses the sensitivity of a TSM parameterization to the specified spatial pattern of lateral fluxes from and to the stream using tracer experiments at the reach scale. Solute transport dynamics were simulated under three contrasting spatial configurations of lateral fluxes: i. The reach is concurrently gaining and losing water over its whole length ii. The upstream half of the reach is gaining water while its downstream half is losing water iii. The upstream half of the reach is losing water while its downstream half is gaining water Simulation results were compared to those from the standard application of the TSM. From the comparison of transient storage model simulations specific to each spatial configuration of lateral fluxes, the following questions are addressed: (1) Do fitted transient storage model parameters depend on the specification of the spatial configuration of lateral fluxes? (2) How large is the uncertainty caused by lack of knowledge of the spatial organization of lateral fluxes relative to the variability in fitted model parameters among reaches? 5.2 Study site description Reach-scale tracer experiments were conducted in Elk Creek (EL), a 5.7 km2 sub-basin within CCEW (Figure 2.1). It encompasses two first-order streams draining two headwaters (Elk North and Elk South) which merge into a second-order stream (Elk Creek). A detailed map of the study area is provided in Figure 5.1. The stream channel features a step-pool-riffle morphology in the two headwaters and pool-riffle morphology in the lower main channel. Despite intensive logging activity in CCEW, most of the 93 riparian vegetation has been preserved in its original state. It is mainly composed of Engelmann spruce (Picea engelmannii), and balsam fir (Abies lasiocarpa) at higher elevations where the foothills are steep and the riparian zone poorly developed, and western red cedar (Thuja plicata) in the wetter and flatter lower part of the sub-basin. Deciduous trees such as alder (Alnus tenuifolia), trembling aspen (Populus tremuloide) and willow (Salix sp.) are typically the dominant understory riparian vegetation, especially in disturbed areas (old logging roadbeds and skid trails). In-stream wood pieces exert a strong influence on the stream morphology. 5.3 Methods 5.3.1 Solute transport modelling We applied OTIS, a one dimensional finite-difference solute transport model that solves for solute concentrations in the stream and transient storage zone [Runkel, 1998]. OTIS accounts for advection, dispersion, lateral fluxes and first order mass transfer between stream and transient storage domains, as represented by the following equations: డ஼ డ௧ ൌ െொడ஼ ஺డ௫ ൅ ଵ ஺ డ డ௫ ቀܣܦ డ஼ డ௫ ቁ ൅ ௤ಽೌ೟಺೙ ஺ ሺܥ௅ െ ܥሻ൅ן ሺܥ௦ െ ܥሻ ௗ஼ೄ ௗ௧ ൌן ஺ ஺ೄ ሺܥ െ ܥௌሻ   5‐1  5‐2  where C, CS and CL are the solute concentrations in the stream, transient storage zone and lateral inflow [g m-3], respectively; A and AS are the stream and transient storage zone cross-sections [m2], respectively; D is the channel dispersion coefficient [m2 s-1]; qLatIn is the lateral inflow per reach unit length [m3 s-1 m-1] to the stream; and qLatOut, the losses to deeper groundwater recharge per reach unit length [m3 s-1 m-1]. The latter does not alter the solute concentration in the stream and thus does not appear in Equation 5-1. 94 However, both lateral fluxes are accounted for when computing the water mass balance. A positive value for qLatIn indicates lateral inflow to the stream and subsequent in-stream solute dilution (assuming CL < C). A positive value for qLatOut indicates losses of stream water and solute mass to deeper groundwater recharge. qLatIn and qLatOut are not interchangeable and in the vast majority of applications will both take positive values. qLatIn may take negative values only to model the effects of stream evaporation on solute concentrations (e.g., [Keefe et al., 2004]). 5.3.2 Delineation of study reaches The active stream network at the time of the freshet in 2005 was divided into 24 contiguous reaches, each roughly 250 m long (Figure 5.1). Tracer tests were conducted over the first two consecutive days of September 2006 to measure solute transport parameters. The study period fell within the baseflow period characterized by minimum fluctuations in streamflow and therefore more reliable estimates of discharge. The upper part of the stream network dried out over the summer and the analysis focuses on the nine lowest adjacent reaches that composed the active stream network at that time of the year (Figure 5.1). We experienced technical issues with the electric conductivity meter while recording data related to the first reach EL01. These data were judged unreliable and thus discarded. For the eight remaining reaches, EL02-EL09, mean longitudinal gradient ranged from 5.9 to 15.9%. The bed sediments were composed of 10-20 cm of alluvial deposit, overlying a compact yet permeable glacial till layer. Characteristics of channel morphology and streamflow for each reach are provided in Table 5-1. 5.3.3 Tracer injection and calculation of reach-scale net inflow During each field campaign, one breakthrough curve (BTC) was recorded for each reach following these steps: i. A tracer pulse was created at the upper end of the lowest reach with a slug injection of sodium chloride (NaCl) lasting between 10 and 20 95 seconds (time needed to rinse the salt off the injection bucket). ii. The BTC was recorded continuously with an electric conductivity meter every 15 seconds until readings dropped back to background. iii. Steps i. and ii. were repeated over the next reach upstream until the eight reaches were covered. Electrical conductivity was used as a surrogate for salt concentration. All measurements were corrected for temperature with the logger’s built-in compensation curve. A global rating curve between electrical conductivity and salt concentration was developed to compute the salt BTC for every reach. Sodium chloride was chosen over organic dye alternatives because of its insensitivity to light and its low adsorption onto bed sediments and organic matter. Slug injection of salt was chosen over constant rate injection to minimize the overall duration of the experimental campaign and the magnitude of flow changes over the duration of the experiment. We acknowledge that the information available for parameter determination in the BTC resulting from a pulse injection is lower than that contained in its constant rate counterpart [Wagner and Harvey, 1997]. Nonetheless, pulse tracer experiments provide sufficient information to derive solute residence time distribution or infer the magnitude and timing of solute transport dynamics from transient storage model parameterization [Gooseff et al., 2003; Haggerty et al., 2002; Zarnetske et al., 2007]. To compute net inflow for each reach, discharge at both ends of every reach was measured on September 5, 2006, using slug injections of table salt [Day, 1976; Moore, 2005] (Table 5-1). There was no substantial change in streamflow patterns between September 1 and September 5, which supported the use of flow data from September 5 for the September 1-2 measurements (Figure 5.2). No rain was recorded during the experiment and streamflow diel fluctuations were mainly due to evapotranspiration (Chapter 4). Moreover, the time between discharge measurements at each end of a given reach never exceeded 25 min. In that time span, streamflow varied less than 1.5%, which is within the discharge measurement uncertainty. 96 5.3.4 Implementation of lateral fluxes into OTIS Spatial patterns of lateral fluxes Due to the integrated nature of flow increment calculations, the gross lateral fluxes, inflow and outflow, and their spatial organization were unknown at sub-reach scales. To understand the sensitivity of fitted TSM parameters to uncertainty in the spatial organization of inflow and outflow, we bracketed the possible range using three hypothetical configurations. The first configuration, called InBh (Inflow Both sub-reaches), is based on the assumption that both lateral inflow and outflow occur concurrently distributed along the entire reach (Figure 5.3a). However, fluxes entering or exiting the channel are driven by local head gradients and studies of hyporheic exchange flow at the channel unit scale have shown that both downwelling and upwelling rarely occur at the same time and location [Anderson et al., 2005]. This consideration led to two alternative configurations for the reach unit. One called InUp (Inflow Upstream), in which the lateral inflow is concentrated in the upper half reach while the lower half experiences only lateral outflow. Conversely, the third configuration, called InDn (Inflow Downstream), characterized by lateral inflow concentrated in the lower half reach and lateral outflow in its upper half (Figure 5.3b and Figure 5.3c). Configurations InDn and InUp were designed so as to maximize the influence of lateral fluxes on solute concentration at the sampling location. Simulation results for the three configurations (InBh, InUp and InDn) were compared to that of a fourth configuration of lateral fluxes that represents a standard application of OTIS (SdAp). In configuration SdAp, parameterization of lateral fluxes is identical for the two sub-reaches. However, and unlike configuration InBh, if the reach to be modelled is a net gaining reach, only qLatIn will be implemented (qLatOut = 0) and if it is a net losing reach, only qLatOut will be implemented (qLatIn = 0) (Figure 5.3d). 97 Representation of the reach unit in OTIS In OTIS, solute transport dynamics are parameterized with the vector θ = (A, D, AS, α) and the two lateral fluxes qLatIn and qLatOut. To account for the spatial variability of lateral fluxes in OTIS and represent the three configurations of the spatial organization of lateral fluxes, the physical reach of length L has been divided into two numerical sub-reaches of length L/2. Subscripts 1 and 2 refer respectively to the upstream and downstream sub- reaches. We can then define a parameter vector for each sub-reach, θ1 and θ2. Since the two sub-reaches belong to the same physical reach, they have same geometric, hydraulic and transient storage properties. Therefore θ1 = θ2. The only difference in parameterization between the two sub-reaches depends on the chosen configuration of spatial organization of lateral fluxes (Figure 5.3, Table 5-2). Reach scale mass balance and relationship between qLatIn and qLatOut As noted in section 5.3.3, the net water gain or loss was known for every reach from discharge increment calculations. If QIN and QOUT are respectively the reach inflow and the reach outflow (Figure 5.3), then the reach scale water balance is: ܳை௎்  െ ܳூே  ൌ  ܮሺݍ௅௔௧ூ௡  െ ݍ௅௔௧ை௨௧ሻ 5‐3 Equation 5-3 implies that a variation in qLatOut will be accompanied by a coordinated variation in qLatIn since the difference QOUT - QIN is fixed and derived from flow increment calculations. There is an infinite number of pairs (qLatIn, qLatOut) that satisfy Equation 5-3 for a given discharge increment QOUT - QIN. However, since qLatOut reflects a loss of solute mass from the stream, consideration of the solute mass balance helps constrain the problem and promotes convergence towards a solution. Further explanations will be provided in the next section, relative to the determination of solute transport parameters. 98 5.3.5 Estimation of parameters and derived metrics Model parameters were determined by optimizing the fit between the simulated and experimental solute breakthrough curves. Best parameter estimates were computed using UCODE_2005, a non-linear optimization package developed by Poeter et al. [2005]. It was first aimed at solving groundwater inverse modelling problems but has been further developed to accommodate a wide range of process models. Scott et al. [2003] used UCODE_2005 in its former version (UCODE), coupled with OTIS, to estimate solute transport parameters from a tracer experiment. They demonstrated the power of the optimization package and its compatibility with OTIS. UCODE_2005 has several advantages compared to alternative calibration methods. It is much less computationally intensive than Monte Carlo simulations and more objective and faster than manual calibration. Moreover, UCODE_2005 yields information on the optimization process such as parameter confidence intervals, sensitivity of the model parameters to observations, correlation between parameters and analysis of residuals, which are powerful tools to assess the robustness of parameters. In the present study, solute transport is characterized by four parameters that are common to both numerical sub-reaches (A, D, AS and α), plus separate specification of qLatOut and qLatIn for each of the two sub-reaches. Because qLatOut and qLatIn are linked by the reach-scale mass balance (Equation 5-3), there are five independent parameters to be estimated: Θ ൌ  ൫ܣ௜, ܦ௜, ܣௌ೔, ߙ௜, ݍ௅௔௧ை௨௧೔൯ 5‐4 with subscript i (that of the sub-reach) taking values of either 1 or 2 for configuration InBh, 2 for configuration InUp and 1 for configuration InDn. Configuration SdAp is parameterized with fixed estimates for lateral fluxes and the model is calibrated over parameters A, D, AS and α only. The goodness of fit between experimental and modelled BTCs was assessed with the root mean squared error (RMSE) objective function: 99 ܴܯܵܧ ൌ  ටଵ ௡ ∑ ቀܿ௜ െ ܿ௜ ᇱሺΘሻቁ ଶ ௡ ଵ    5‐5  where ci is the ith concentration reading of the experimental BTC; c’i is the simulated counterpart of ci and n the number of observations. Best parameter estimates resulting from the optimization process could be used directly to assess the sensitivity of model parameters to the spatial organization of lateral fluxes. However, metrics based on solute transport parameters have been preferred over the parameters themselves because they provide interpretation of exchange processes. [Harvey et al., 1996; Runkel, 2002] (Table 5-3). The analysis was based on both the OTIS parameters and derived metrics. All four configurations were implemented for each reach, yielding four sets of optimized parameters for each BTC: ΘInBh, ΘInDn, ΘInUp, ΘSdAp. 5.3.6 Assessment of parameter estimates To assess the reliability of estimates of solute transport parameters and more specifically exchange parameters between the stream and the transient storage zone, Damköhler indices (DaI) were computed, as suggested by Wagner and Harvey [1997]. The DaI is the ratio between in-stream travel time and the time scale of solute exchange with the transient storage zone (Table 5-3). A DaI much smaller than unity indicates that little solute is exchanged with the transient storage zone, leading to unreliable estimates of exchange parameters. A value of DaI much larger than unity indicates that the solute completely enters transient storage and the effects of in-stream dispersion become similar to that of exchange processes with the transient storage zones. Either way, a departure from unity of more than one order of magnitude is associated with a dramatic increase in the uncertainty of exchange parameter estimates [Wagner and Harvey, 1997]. In addition to the Damköhler index, reliability and robustness of 100 parameter estimates were assessed using outputs from UCODE_2005, including composite scale sensitivity (css), and parameter correlation. The css is a measure of the total information provided by the regression (here, the BTC) about a parameter value [Hill and Tiedeman, 2007]. Here, the css value for each parameter is expressed as a percentage of the css of the parameter to which the model is the most sensitive. Insensitive or problematic parameters, in an optimization sense, show css values of about or less than 1% relative to the most sensitive parameter [Hill and Tiedeman, 2007]. While scaled sensitivity is an essential tool to investigate the model structure and quality of the fit, it does not account for interactions between parameters or the fact that coordinated changes between two or more of them may produce the same fit to experimental BTC. Depending on the complexity of the experimental BTC, the number of physical processes to be represented in the theoretical model and the number of parameters, models can be over- parameterized. In that case, high correlation between parameters arises and there is no longer a unique optimum value for every parameter. A standard output from UCODE_2005 is a set of parameter correlation coefficients (pcc) which can be used not only to assess the uniqueness of the optimized parameter vector but also to identify any ill-posed problems such as incorrect model structure or model over-parameterization. Parameters are likely to be estimated uniquely if pcc absolute values are less than about 0.95 [Hill and Tiedman, 2007]. The quality of parameter estimates was assessed in the light of DaI, css and pcc values. 5.3.7 Statistical analysis The statistical analysis addressed whether the fitted parameters for the different inflow patterns belong to the same population. Formally, the analysis tested the following null and alternative hypotheses: H0: There is no influence of spatial pattern of lateral fluxes on estimated transient storage parameters. Ha: The estimated transient storage parameters vary systematically 101 with the assumed configuration of lateral fluxes. The hypotheses were tested using an analysis of variance (ANOVA) based on a Randomized Complete Block design (RCB) with one fixed factor: the spatial organization of lateral fluxes. Stream reach was considered a random blocking factor. A special application of RCB called Repeated Measures Design was used, as different configurations of the spatial organization of lateral fluxes were applied to the same individuals (the reaches) [Kutner et al., 2004]. Where the ANOVA indicated a significant effect of lateral flow pattern on the parameter estimates, we followed up with pairwise paired t-tests to identify which configurations produced estimates that differed from the others. Because there was no a priori reason to hypothesize which configuration yielded the largest or the smallest parameter estimates, all pairwise comparisons were based on two-tailed t-tests. Moreover, we did not test for any block effect (i.e., differences among reaches) as the purpose of blocking was to decrease the noise and highlight potential differences across configurations. The general form of the statistical model is: ݕ௜௝ ൌ ߚ௜ ൅ ௝ܾ ൅ ߝ௜௝   5‐6 where yij is the parameter estimate from configuration i and reach j; βi is the fixed effect, [InBh, InDn, InUp, SdAp]; bj is the block (stream reach), [EL02- EL09]; and ߝ௜௝ ൎ ࣨሺ0, ߪଶሻ represents the random errors, which are assumed to be independent normal random deviates with a mean of 0 and a variance of σ2. The RCB design treats the experimental reaches as a sample from a broader population about which we want to make inferences. When the ANOVA test leads to rejection of the null hypothesis for a given parameter, the inferences regarding the relation between the parameter value and the 102 configuration of inflow and outflow can be extrapolated to the broader population. All statistical analysis and plotting were conducted with the programming language R. Whiskers in all boxplots presented in this chapter extend to the most extreme data point that is no more than 1.5 times the interquartile range. The statistical significance of differences in estimated parameter values was assessed using a significance level of αC = 0.05. For parameters for which the ANOVA revealed a significant difference among lateral flow patterns, pairwise t-tests were used to assess which configurations differed from each other. Because this procedure requires multiple t-tests for each parameter, the p-values resulting from the t-tests were adjusted using the Bonferroni procedure and then compared to αC/2 [Kutner et al., 2004]. 5.4 Results 5.4.1 Assessment of goodness of fit For each of the 32 simulations (four configurations of lateral fluxes for each of eight reaches), the non-linear optimization algorithm implemented in UCODE_2005 successfully converged toward a solution (best estimate for Θ). Two reaches have been selected in Figure 5.4 to illustrate the goodness of fit of the simulated BTCs from all four configurations to the experimental BTC. Semi-log plots of the same BTCs highlight the departure of simulated BTCs from the experimental ones at the tail of the BTC. Goodness of fit was assessed in a more objective manner using RMSE (Table 5-4). BTCs resulting from InBh, InDn and InUp simulations consistently showed minimal differences among each other. They all departed from the experimental BTCs in the tail of the BTC. BTCs resulting from configuration SdAp consistently departed from experimental BTCs mainly in the tail area but also at the beginning of the rising limb. If none of the simulations reproduced observed behaviour in the tail of the BTC, those from 103 configurations InBh, InDn and InUp consistently underestimated tracer concentration while that from configuration SdAp consistently overestimated observed tracer concentrations (Figure 5.4). For all reaches, configurations InBh, InDn and InUp improved the goodness-of-fit compared to the standard OTIS application SdAp. Accounting explicitly for both lateral inflow and outflow, regardless of their spatial arrangement, reduced RMSE two to three fold (Table 5-4). Note that five parameters were optimized for configurations InBh, InDn and InUp while only four were optimized for configuration SdAp. Optimized parameters are presented in Table 5-5. For configurations InUp and InDn, either qLatOut or qLatIn is each concentrated over a single sub-reach (either the upstream or downstream sub-reach). Conversely, lateral fluxes are applied to the whole length of the reach for configurations InBh and SdAp. For ease of comparison with estimates of qLatIn and qLatOut from configurations InBh and SdAp, those of configurations InDn and InUp were halved. 5.4.2 Comparisons between configurations InBh - InDn - InUp In a first step, patterns of parameter values were compared among configurations InBh, InDn and InUp because these three configurations share the same model structure, based on five parameters. In the next section, we focus on parameter estimates from configuration SdAp, based on a four-parameter model structure and compare them to those obtained from the five-parameter model structure. Visual inspection of the whisker plots suggests that the distributions of A, AS and qLatOut estimates were indeed influenced by the choice of the spatial organization of lateral fluxes along the reach (Figure 5.5). A, AS and qLatOut were systematically larger for configuration InUp and smaller for configuration InDn. Parameters D and α were seemingly not affected by the choice of the spatial arrangement of lateral fluxes. However, ANOVA tests for all five parameters led to the rejection of the null hypothesis and acceptance of the alternative hypothesis, i.e., that the assumed spatial pattern of lateral fluxes has an influence on estimated transient storage parameters at a 104 significance level of 0.05, marginally so for parameter α (p-value = 0.03) (Table 5-6). For parameters A, AS and qLatOut, there were statistically significant differences among all pairs of configurations (all adjusted p-values were lower than 0.0023, Table 5-7). Differences among configurations were small for parameters D and α compared to the variability among reaches (Figure 5.5). The pairwise differences in α were not found to be statistically significant (adjusted p-values greater than 0.44). Estimates of parameter D from configuration InBh was significantly different from those of the two other configurations (both adjusted p-values equal to 0.0027). However, estimates of parameter D between configurations InUp and InDn were not found to be statistically different (adjusted p-value=0.57). Among the inflow- outflow configurations, parameter estimates were most variable across reaches for configuration InUp. Marked differences were found among configurations for the fraction of median travel time due to transient storage, ܨ௠௘ௗଶ଴଴ (Figure 5.6). Estimates were usually larger for configuration InUp (3.0 to 8.0%), intermediate for configuration InBh (3.4 to 7.9%) and systematically smaller for configuration InDn (2.5 to 6.9%). Differences were found to be statistically significant between pairs of configurations InBh-InDn and InUp-InDn but not between configurations InBh and InUp (Table 5-7). Means and coefficients of variation (CV) of ܨ௠௘ௗଶ଴଴ across all three configurations and for every reach are summarized in Table 5-8. Average values ranged from 3% to 7.2% and CVs ranged from 4.2 to 20%, depending on the reach. Visual differences among configurations were small for the transient storage zone mean residence time (TS) and its variability among reaches was similar from one configuration to another (Figure 5.6). However, TS estimates from configuration InBh were statistically different from those from configuration InDn (p-value = 0.0145, Table 5-7). The solute exchange length, LS, was systematically larger for configuration InDn and smaller for configuration InUp. Although characterized by a large variability between reaches, significant differences among configurations were found for LS 105 (adjusted p-values < 0.01). 5.4.3 Standard application of OTIS - configuration SdAp Parameter estimates from configuration SdAp were most often markedly different from those from the other three configurations. Differences were especially striking for parameters AS, D and α (Figure 5.5). Estimates of AS for configuration SdAp were nearly ten times larger than those estimated for configurations InBh, InDn and InUp. Estimates of α for configuration SdAp were usually smaller than those estimated for configurations InBh, InDn and InUp and especially much less variable among reaches. This resulted in relatively larger TS and LS (Figure 5.6). The metric TS was about ten times larger for configuration SdAp than for the other three configurations, while LS was consistently greater than the physical length of the reach (Table 5-1, Figure 5.6). ܨ݉݁݀200 estimates from configuration SdAp were more variable and three to four times larger than estimates from the three other configurations. Consistent and strong differences in parameter estimates from configurations InBh, InDn and InUp as a whole and configuration SdAp were especially striking and motivated further efforts to elucidate the interactions between physical processes and numerical optimization. In order to so, separate mass balance analyses were conducted in both the stream and the transient storage zone. 5.4.4 Solute mass balance, timing and location The stream solute concentration was integrated in time in order to compare solute mass recovery among configurations (Table 5-9). Configurations InBh, InDn and InUp performed equally well in term of solute mass recovery but more solute mass was consistently recovered for configuration SdAp even in the case of reach EL02, a net losing reach, for which qLatOut was implemented, thus allowing solute mass to exit to stream to groundwater. 106 Regardless of the reach or the configuration of spatial arrangement of lateral fluxes, mass recovery reached at best 95.5%, which is an indication that solute mass was either lost to groundwater or remained in the transient storage zone for the duration of the experiment. In order to shed some light on the location where solute mass is being held, the transient storage zone BTC output from OTIS was analysed for a net losing reach, EL02 and a net gaining reach, EL07. For both reaches, solute concentration in the transient storage zone at the end of the experiment was higher with configuration SdAp than with the three other configurations (Figure 5.7). Moreover, the amount of solute held in the transient storage zone at the end of the experiment was systematically greater for simulations with configuration SdAp and up to 25.3% of the injected solute mass for reach EL02 (Table 5-10) 5.4.5 Assessment of parameter uncertainty Damköhler indices from all simulations are presented in Figure 5.8. Among configurations InBh, InDn and InUp, DaI were consistently larger for configuration InUp and smaller for configuration InDn. All 24 values ranged from 4.3 to 31.1 and the three configurations showed marked variability. Median values were 9.25, 11.05 and 13.10, respectively, for configurations InDn, InBh and InUp. DaI values from configuration SdAp were substantially smaller and less variable across reaches than those from the three other configurations, ranging from 1.1 to 1.8. Based on the composite scaled sensitivities, OTIS simulations were most sensitive to cross-sectional area A, followed by qLatOut, AS, D and α (Figure 5.9). For configurations InBh, InDn and InUp, relative sensitivity to qLatOut ranged from 3 to 35% with a mean value of 14% and relative sensitivity to AS ranged from 3 to 8%. Although all relative css values were greater than 1%, relative css values for α remained low (3.3-4.5%). A different pattern was observed for configuration SdAp, for which OTIS simulations were more sensitive to α (4.2-8.5%) than to AS (1.4-2.6%). Absolute values of parameter correlation coefficients (pcc) were binned into two intervals of unequal width: [0, 0.95] and [0.95, 1] (Figure 5.10). For 107 configurations InBh, InDn and InUp and out of the 240 pcc values, (10 pairwise correlations, 24 experimental runs), 31 indicated a dependency between two parameters for a given run. Most often (19 out of 31), α was highly correlated with other parameters. For configuration SdAp, absolute values of correlation coefficients between parameters were usually low and never exceeded 0.95. 5.5 Discussion 5.5.1 Patterns of parameter estimates among configurations InBh – InDn – InUp It remains unclear why estimates of D and α did not show strong contrasts among configurations, although we suspect the lack of contrast resulted from tradeoffs between parameter values during the optimization process to accommodate the shape of the observed BTC. In the following, we focus on interpreting the patterns of A, qLatOut and AS, for which strong contrasts exist among configurations. The systematic variations in the parameter estimates among configurations can be interpreted in relation to the need for all simulations to reproduce the same BTC in terms of timing and stream solute mass balance for the reach. We also cast the results in terms of ܨ௠௘ௗଶ଴଴ to describe the overall effect of transient storage on downstream solute transport. Stream cross section, A Solute advection time as presented in Equation 5-1 depends on the ratio L·A/Q(x). As noted above, this advection time remains invariant among configurations. Consequently, the differences in patterns of parameter A will be discussed in relation to discharge profiles Q(x) among configurations, considering the case of reach EL08 as an example. Discharge profiles along reach EL08 for the configurations InBh, InDn and InUp are shown in Figure 5.11. From this figure, it appears that at any location along the reach, to the 108 exception of its upper and lower boundaries, the simulated discharge will always be greatest for configuration InUp and lowest for configuration InDn. In OTIS, water particle velocity U(x) averaged over the stream cross-section is defined as the ratio Q(x)/A. Water particle velocity controls the advection term in Equation 5-1. Its harmonic mean (MHU) along the reach is computed as follows: ܯܪܷ ൌ ׬ ௗ௫ ௎ሺ௫ሻ ௅ ଴    5‐7 and is directly related to mean solute travel time and, therefore, to time to peak of the BTC. Although the longitudinal profiles of U(x) varied, depending on the configuration of lateral fluxes, their harmonic means must be invariant across configurations because they were constrained by the time to peak of the same experimental BTC. Since:  ܣ ׬ ௗ௫ ொሺ௫ሻ ௅ ଴ ൌ ׬ ௗ௫ ௎ሺ௫ሻ ൎ ܿ݋݊ݏݐܽ݊ݐ ௅ ଴ 5‐8 configuration InUp, characterized by greater values for Q(x) relative to both configurations InBh and InDn, would yield a coordinated increase in A (Figure 5.5). Therefore, estimates of A were greatest for configuration InUp, followed by that of configurations InBh and then InDn. Lateral outflow, qLatOut The parameter qLatOut controls the loss of solute along the reach, the mass of which is given by: 109  ܯ௟௢௦௧ ൌ ׬ ׬ ܥሺݔ, ݐሻݍ௅௔௧ை௨௧ሺݔ, ݐሻ݀ݔ ݀ݐ ் ଴ ௅ ଴ 5‐9 where T is the duration of the experiment. The mass balance for the reach requires that Mlost = Minj-Mrec (solute mass lost = injected solute mass – solute mass recovered), where the solute mass recovered at the sampling site is: ܯ௥௘௖ ൌ ܳை௎் ׬ ܥሺݐሻ݀ݐ ் ଴ 5‐10 The mass recovered only depends on QOUT and C(t), and therefore is independent of the configuration of lateral fluxes. We showed above that at any location along the reach, discharge was always greatest (smallest) for configuration InUp (InDn). The in-stream solute concentration will therefore always be smallest (more diluted) for configuration InUp and largest for configuration InDn. Since the mass of solute lost to lateral outflow is directly related to C(x,t), qLatOut has to be greatest (smallest) for configuration InUp (InDn), which is in agreement with results presented in Figure 5.5. Transient storage zone cross section, AS To understand the sensitivity of AS to the assumed configuration of lateral fluxes, it is useful to consider how it influences the volumetric solute residence time: ௩ܶ௢௟ ൌ ܮሺܣ ൅ ܣௌሻ/ܳ   5‐11 Tvol is the first moment of the BTC and is therefore invariant among configurations since the same experimental BTC is being simulated. Moreover, we showed earlier that the harmonic mean water velocity was also 110 invariant among configurations. If we re-write Equation 5-11 as: ௩ܶ௢௟ ൌ ሺ1 ൅ ܣௌ/ܣሻ · ܮ/ܷ 5‐12 the ratio AS/A must also remain invariant across configurations to verify the equation. The coefficient of variation (CV) across configurations of the ratio AS/A never exceeded 4% depending on the reach (Table 5-11). As a result, since estimates of A were greatest (smallest) for configuration InUp (InDn), estimates of AS were also greatest (smallest) for configuration InUp (InDn) (Figure 5.5). Fraction of median transport time due to transient storage, ࡲ࢓ࢋࢊ૛૙૙ The assumed spatial pattern of lateral fluxes had a significant systematic influence on the estimates of ܨ௠௘ௗଶ଴଴ . However, even though statistically significant, the difference in ܨ௠௘ௗଶ଴଴ among configurations remained small. The coefficient of variation (CV), a measure of the difference between metric estimates across the three configurations, ranged from 4.2% to 20% depending on the reach. A number of previous studies were selected in order to compare the inter-configuration variability of ܨ௠௘ௗଶ଴଴ estimates (in terms of CV) to those found from inter-reach comparisons (Table 5-8). The CVs ranged from 56% to 149%, depending on the study. These values are far beyond the range of inter-configuration variability found in the present study, which might indicate that, even though statistically significant, the variability in ܨ௠௘ௗଶ଴଴ estimates related to the choice of lateral flux configuration is small compared to the inherent inter-reach variability of ܨ௠௘ௗଶ଴଴ estimates. It is important to note that, in the previous studies, reach units were usually chosen so as to maximize the difference among them, in terms of hydraulic or morphological properties. Therefore, CVs of ܨ௠௘ௗଶ଴଴ estimates were probably artificially increased by the search for maximum contrasts among reaches. Moreover, the eight contiguous reaches selected in the CCEW were likely to be under a geographic and climatic “blocking effect”. Individual reaches 111 showed distinct morphological characteristics, but all reaches have a similar geological setting and all experiments were conducted under similar hydrological conditions. The low CVs of ܨ௠௘ௗଶ଴଴ estimates found across configurations could also be attributed to the generally low ܨ௠௘ௗଶ଴଴  values characterizing this study. Estimates ranged from 2.5 to 7.9% across all configurations and reaches, indicating an overall weak influence of transient storage on downstream solute transport. Values of ܨ௠௘ௗଶ଴଴ found in the above- mentioned studies were generally larger, indicating a stronger influence of transient storage on solute transport than in the present study. It is possible that the small differences in ܨ௠௘ௗଶ଴଴  estimates across configurations were due to a generally weak influence of transient storage on downstream solute transport. For reaches showing a much stronger coupling with their neighbouring transient storage zones, the difference in ܨ௠௘ௗଶ଴଴ estimates across configurations may increase. In other words, the variability in ܨ௠௘ௗଶ଴଴  estimates due to various figurations of lateral fluxes may increase with an increase of ܨ௠௘ௗ ଶ଴଴ . 5.5.2 Parameter uncertainty In this study, all estimated parameters and metrics, with the exception of α, were significantly influenced by the configuration of lateral fluxes. It is believed that the lack of a significant difference among configurations for α results from the difficulty in generating reliable estimates from analysis of the BTCs rather than a true lack of sensitivity to the specified pattern of inflows and outflows. Of all of the parameters, OTIS simulations were least sensitive to α, even though relative composite scale sensitivity values were above the critical threshold of 1% (Figure 5.9) [Hill and Tiedeman, 2007]. This lack of sensitivity, coupled with the fact that α was often highly correlated with other parameters during the optimization process (pcc > 0.95, Figure 5.10), support the hypothesis that estimates of α might not have been determined uniquely, even though UCODE_2005 was successful in converging toward a solution. Moreover, Damköhler indices were far greater 112 than 1 and hence were outside of the range suggested by Wagner and Harvey [1997]. This implies that transient storage effects on solute transport dynamics cannot be separated from those of in-stream dispersion processes. As a result, and even though large values of DaI do not necessarily indicate a bias in parameter estimates, they reflect a substantial increase in parameter uncertainty, especially for parameters related to solute exchange with the transient storage zone, α and AS. Given these apparent problems in estimating α, it is likely that, even if the estimates of α were truly related to the choice of the configuration of lateral fluxes, this relation would have been hidden by the uncertainty characterizing estimates of α. 5.5.3 Standard application of OTIS - configuration SdAp In an inverse modelling context, physical processes are often assessed and compared via the values of optimized model parameters. Here, we demonstrated that different simulations based solely on variations of the same model structure, in which both qLatIn and qLatOut were implemented, performed equally well yet yielded significantly different sets of parameters. Further, parameter estimates from configuration SdAp were markedly different from those of the three other configurations of lateral fluxes, with the exception of α. Despite the fact that there was enough information in the signal to determine optimal model parameters, and despite the fact that tight correlations between model parameters were rare, our results demonstrate the potential for equifinality, which undermines our ability to infer physical processes from parameter values for inter-reach comparisons and other purposes. We cannot reject any of the specific model structures because goodness-of-fit was equally satisfactory for all simulations, and optimized parameter values were reasonable. Yet, can we come up with a novel standard for model structure that would be most suitable when lateral fluxes are poorly defined? OTIS models solute fate in both stream and transient storage zone but most often is calibrated only against the stream BTC (versus the transient storage zone BTC). Solute mass recovery is a sensible measure of model quality that may provide decisive information when selecting the 113 most adequate model structure. Consistently more solute mass was recovered for configuration SdAp. This is encouraging as a single model structure comes forward, which would allow further interpretation of model parameters and related solute transport dynamics. However, stream solute concentration was systematically over-estimated in the tail of the BTC with configuration SdAp while all three other configurations led to underestimating observed solute concentration in the tail of the BTC (Figure 5.7). Seemingly higher mass recovery is, in fact, an artefact of the lack of fit in the tail of the BTC. Simulated mass recovery ranged between 59 and 96%, which raised the question of solute fate. Solute decay rates or sorption-desorption dynamics should certainly be accounted for when discussing solute mass recovery. These processes are, however, beyond the scope of the present study, especially given that we used sodium chloride over relatively short time and spatial scales on purpose, so as to mitigate the magnitude of these processes and neglect them in a first order approximation. The stream solute BTC reflects the amount of solute mass to be transported downstream at a given time. The remaining mass of solute can either be stored in the transient storage zone, to later be released in the stream and contribute to the tail of the BTC, or exit the stream via lateral outflow, never to contribute to the BTC. Note that from an experimental point of view, both options are viable. In practice, the stream solute BTC stops being recorded when solute concentration returns to background readings. In the first case, solute would be released from the transient storage zone to the stream at sub-detection concentration levels, while in the second case, solute would effectively exit the stream reach of interest. For a model structure in which both qLatIn and qLatOut are implemented (configurations InBh, InDn and InUp), the dominant mechanism to remove solute mass from the stream is to release it via lateral flow to groundwater reservoirs or hyporheic flowpaths that extend downstream of the sampling location. As a result, the model requires limited amounts of solute to enter the transient storage zone which are more rapidly mobilized back to the stream in the tail of the BTC (Figure 5.7). This also implies that transient 114 storage will have a weak influence on median solute transfer time, ܨ௠௘ௗଶ଴଴ (Figure 5.6). Conversely, for a model structure in which only qLatIn or qLatOut is implemented (configuration SdAp), the preferred mechanism to remove solute mass from the stream is to store it in the transient storage zone. This means that solute is being held in flow paths that will eventually discharge in the reach of interest, but at time scales that extend beyond that of the experiment. From a model parameterization point of view, the transient storage zone is increased accordingly (and sometimes disproportionately) to accommodate the solute mass. Solute will be released at extremely slow rates (low α) resulting in an extended BTC tail (Figure 5.5 and Figure 5.7). Interestingly, this mechanism is observed even for reach EL02, a net losing reach, in which qLatOut was implemented and potentially would have allowed solute mass to exit the stream with lateral outflow (Figure 5.7). This is evidence that model parameterization and subsequent interpretation of solute transport dynamics cannot be dissociated from the model structure. 5.5.4 Implications for future stream tracer experiments and modelling with OTIS Information contained in the stream BTC is insufficient to choose the model structure that would simulate physical processes most faithfully. However, while solute concentration in the transient storage zone fell below detection limit for simulations with configurations InBh, InDn and InUp, it remained above detection limit for simulations with configuration SdAp (Figure 8). Adequate sampling of the transient storage zones, by providing information about the solute concentration dynamics both in time and space, could help constrain the model parameterization and choose the proper model structure. Alternatively, using another tracer/sensor combination to substantially lower the solute concentration detection limit would provide a better definition of the tail of the experimental BTC and related transient storage processes. However, doing so would only assist in selecting the most appropriate model structure; parameter estimation would still be limited by lack of information on the spatial configuration of lateral fluxes, as these 115 estimates were, in most cases, significantly different among configurations InBh, InDn and InUp. This observation highlights the need for a better characterisation of lateral fluxes both in space and time, prior to modelling solute transport dynamics. Downwelling and upwelling can vary along the same reach over relatively short spatial scales, but direct measurements of such fluxes are difficult and can only be bracketed using discharge and solute concentrations at the upper and lower reach boundaries (e.g., [Payn et al., 2009]). Deployment of piezometers to quantify hydraulic gradients in the streambed sediments, wells in the riparian zone, and longitudinal surveys of water temperature and electrical conductivity or water chemistry can provide at least qualitative information on the location of flow losses and gains (e.g., [Story et al. 2003]), which could help guide the specification of flow configuration. This information would also be useful when applying 3-D groundwater models to simulate fluxes between the stream and the transient storage domains, in order to better define the boundary conditions. Recognizing the difficulty in quantifying spatial patterns and fluxes of lateral inflows and outflows in field studies (e.g., [Scordo and Moore, 2009]), we recommend further research employing controlled tracer experiments on small-scale laboratory models, for which boundary conditions can be varied in a deliberate manner. 5.6 Conclusion This study has shown that comparisons of transient storage model parameters among reaches or as a function of time at a single reach may be confounded not only by the currently recognized sources of uncertainty in analysing BTCs, but also by incorrect specification of the spatial patterns of inflow and outflows. Model parameterization is thus highly dependent on the model structure. Moreover, the stream BTC alone does not provide enough information to select the appropriate model structure and then allow for a reliable interpretation of solute transport dynamics from model parameters. We recommend that future studies of transient storage processes include 116 detailed measurements that will provide information on the spatial patterns of flow losses and gains, including piezometers inserted into the stream bed and wells in the riparian zone, and longitudinal profiles of electrical conductivity, water chemistry and temperature. In addition, we call for additional small-scale laboratory tracer experiments, in which boundary conditions are sufficiently well defined (e.g., via piezometry and water chemistry, as well as by the physical design of the laboratory model) to constrain the numerical model structure and boundary conditions. 117 Table 5-1 : Reach characteristics. QOUT refers to the discharge at the lower end of the reach as measured on September 5, 2006, while net budget refers to the integrated value of both gross lateral inflow and gross lateral outflow along the reach. Gaining reaches are characterized by positive net budget values; losing reaches are identified by negative values. Reach Length [m] Sinuosity [-] Mean gradient [%] QOUT [l s-1] Net budget [l s-1] EL01 248 1.12 9.5 18.5 1.5 EL02 294 1.09 10.8 17 -0.2 EL03 263 1.35 6.5 17.2 0.9 EL04 355 1.09 10.0 16.3 0.9 EL05 251 1.09 6.2 15.4 0.5 EL06 256 1.09 5.7 15 0.7 EL07 82 1.05 15.9 14.2 1.7 EL08 249 1.12 9.5 12.6 3.6 EL09 277 1.13 8.3 8.9 1.9 Table 5-2 : Lateral fluxes parameterization for each configuration. Configuration Sub-reach 1 Sub-reach 2 InBh ½ qLatOut, ½ qLatIn ½ qLatOut, ½ qLatIn InDn qLatOut, qLatIn =0 qLatOut =0, qLatIn InUp qLatOut =0, qLatIn qLatOut, qLatIn =0 SdAp ½ qLatOut, qLatIn =0 ½ qLatOut, qLatIn =0 or ½ qLatIn, qLatOut =0 ½ qLatIn, qLatOut =0 118 Table 5-3 : Metrics derived from solute transport parameters and their significance. ܨ௠௘ௗଶ଴଴ has been normalized to a reach length of 200m for inter- reach comparison purposes. Metrics Formulation Signification TS ܣௌ ߙܣ [T] Average solute residence time in the transient storage zone LS ܳ ߙܣ [L] In stream solute travel length before passing into the transient storage zone DaI ߙ ቀ1 ൅ ܣ ܣௌൗ ቁ ܮ ܳ ܣൗ [-] Ratio of in stream travel time to exchange processes time scales ܨ௠௘ௗ ଶ଴଴ ܣௌ ܣ ൅ ܣௌ ቆ1 െ ݁ ିଶ଴଴ቀఈ஺ொ ቁቇ [%] Fraction of median transport time due to transient storage normalized to a reach length of 200 m Table 5-4 : RMSE values [mg l-1]. For reference, the precision of NaCl concentration was 0.54 mg l-1. InBh InDn InUp SdAp EL02 0.198 0.196 0.196 0.429 EL03 0.195 0.193 0.194 0.452 EL04 0.140 0.138 0.139 0.355 EL05 0.163 0.161 0.161 0.434 EL06 0.078 0.078 0.078 0.155 EL07 0.370 0.368 0.369 0.956 EL08 0.225 0.220 0.218 0.483 EL09 0.236 0.232 0.231 0.553 119 Table 5-5 : Optimized solute transport parameters. qLatIn is included in the table for reference despite the fact that it was not an optimized parameter sensu stricto but rather varied co-ordinately with qLatOut or was fixed. Reach Configuration QIN [m3 s-1] A [m2] D [m2 s-1] AS [m2] α [s-1] Net reach budget [m3 s-1] qLatOut [m3 s-1 m-1] qLatIn [m3 s-1 m-1] E02 InBh 1.72E-02 0.125 0.141 0.0096 3.76E-04 -2.27E-04 2.27E-05 2.19E-05 InDn 0.104 0.138 0.0079 3.65E-04 1.89E-05 1.81E-05 InUp 0.153 0.138 0.0116 3.63E-04 2.79E-05 2.71E-05 SdAp 0.130 0.203 0.1076 2.36E-04 7.72E-07 - E03 InBh 1.63E-02 0.132 0.106 0.0116 7.71E-04 8.72E-04 2.67E-05 3.14E-05 InDn 0.114 0.105 0.0098 7.49E-04 2.2E-05 2.67E-05 InUp 0.156 0.105 0.0134 7.49E-04 3.22E-05 3.69E-05 SdAp 0.140 0.153 0.0734 2.78E-04 - 4.74E-06 E04 InBh 1.54E-02 0.140 0.117 0.0124 3.46E-04 8.88E-04 1.92E-05 2.17E-05 InDn 0.113 0.113 0.0098 3.36E-04 1.51E-05 1.76E-05 InUp 0.181 0.114 0.0157 3.37E-04 2.26E-05 2.5E-05 SdAp 0.146 0.174 0.0868 2.03E-04 - 2.49E-06 120 Reach Configuration QIN [m3 s-1] A [m2] D [m2 s-1] AS [m2] α [s-1] Net reach budget [m3 s-1] qLatOut [m3 s-1 m-1] qLatIn [m3 s-1 m-1] E05 InBh 1.50E-02 0.162 0.082 0.0127 2.98E-04 4.61E-04 1.86E-05 2.05E-05 InDn 0.139 0.080 0.0109 2.93E-04 1.57E-05 1.76E-05 InUp 0.191 0.080 0.0149 2.93E-04 2.23E-05 2.41E-05 SdAp 0.167 0.109 0.0705 1.91E-04 - 1.86E-06 E06 InBh 1.42E-02 0.151 0.094 0.0157 8.97E-04 7.35E-04 3.12E-05 3.43E-05 InDn 0.118 0.091 0.0114 8.25E-04 2.37E-05 2.67E-05 InUp 0.202 0.091 0.0196 8.24E-04 4.26E-05 4.56E-05 SdAp 0.164 0.139 0.1861 2.30E-04 - 3.06E-06 E07 InBh 1.26E-02 0.135 0.067 0.0130 8.61E-04 1.65E-03 2.11E-05 4.07E-05 InDn 0.123 0.067 0.0118 8.50E-04 1.79E-05 3.75E-05 InUp 0.149 0.067 0.0143 8.52E-04 2.46E-05 4.42E-05 SdAp 0.142 0.093 0.0397 2.59E-04 - 1.96E-05 121 Reach Configuration QIN [m3 s-1] A [m2] D [m2 s-1] AS [m2] α [s-1] Net reach budget [m3 s-1] qLatOut [m3 s-1 m-1] qLatIn [m3 s-1 m-1] E08 InBh 8.94E-03 0.127 0.078 0.0154 2.25E-04 3.65E-03 1.61E-05 3.18E-05 InDn 0.099 0.074 0.0118 2.22E-04 1.12E-05 2.69E-05 InUp 0.168 0.074 0.0201 2.20E-04 2.35E-05 3.92E-05 SdAp 0.131 0.101 0.0969 2.29E-04 - 1.57E-05 E09 InBh 7.07E-03 0.099 0.082 0.0123 2.10E-04 1.87E-03 8.84E-06 1.73E-05 InDn 0.083 0.080 0.0103 2.08E-04 6.85E-06 1.53E-05 InUp 0.119 0.080 0.0148 2.07E-04 1.15E-05 1.99E-05 SdAp 0.101 0.099 0.0477 2.02E-04 - 8.42E-06 122 Table 5-6 : ANOVA table for all five parameters. The degrees of freedom of the error term for all models are equal to 14. Parameter F-Value p (>F) A 59.993 <0.001 D 27.238 <0.001 AS 43.4815 <0.001 α 4.3909 0.0331 qLatOut 51.074 <0.001 Table 5-7 : p-values arising from pairwise paired t-tests, adjusted according to the Bonferroni procedure. Since we used two-tailed t-tests, they have to be compared to αC/2=0.025. The degree of freedom of the error term for all models is equal to 14. p-values leading to failure to reject H0 are highlighted in bold. Comparison A D AS α qLatOut TS LS ܨ௠௘ௗଶ଴଴ InDn-InBh <0.001 0.0072 0.0023 0.46 <0.001 0.0145 0.0072 0.0026 InUp-InBh 0.0011 0.0072 0.0018 0.44 <0.001 0.0781 0.0099 0.0395 InUp-InDn <0.001 0.5699 0.0018 1 0.00167 1 0.008 <0.001 SdAp-InBh 0.0075 0.002 0.012 0.16 - 0.0035 0.0245 0.0068 SdAp-InDn <0.001 0.0017 0.0112 0.16 - 0.0035 0.4606 0.0059 SdAp-InUp 0.0025 0.0018 0.0136 0.17 - 0.0035 0.0022 0.0083 123 Table 5-8 : Comparison of ܨ௠௘ௗଶ଴଴ inter-reach variability from previous studies and variability across configurations for every reach of the present study. ܨ௠௘ௗ ଶ଴଴ n Mean [%] CV [%] [Runkel, 2002] 53 10.2 149 [Ensign and Doyle, 2005] 8 15.2 63.2 [Lautz et al., 2007] 3 14.8 91.2 [Wondzell, 2006] 8 46.9 56.1 EL02 3 3 16.7 EL03 3 5.6 8.9 EL04 3 3.7 18.9 EL05 3 3.4 11.8 EL06 3 7.5 10.7 EL07 3 7.2 4.2 EL08 3 4.5 20 EL09 3 4.5 13.3 Table 5-9 : Solute mass recovery for the different configurations of spatial organizations of lateral inflow as a percentage of the initial solute mass injected in the stream (Minj). Mass recovery [%] Minj [g] Obs InBh InDn InUp SdAp EL02 298.2 67.6 67.4 67.4 67.4 68.8 EL03 299.8 74.0 74.1 74.1 74.1 75.0 EL04 397.7 63.3 63.4 63.4 63.4 64.2 EL05 400.5 72.4 72.6 72.6 72.6 73.5 EL06 300.2 58.7 58.8 58.8 58.8 59.1 EL07 300.3 89.7 87.8 87.8 87.8 95.5 EL08 306.2 69.0 69.1 69.1 69.1 70.7 EL09 292.8 76.4 76.4 76.4 76.4 78.1 124 Table 5-10 : Solute mass held in the transient storage zone (MS) at the end of the experiment as a percentage of the initial solute mass injected in the stream (Minj). InBh InDn InUp SdAp EL02 CS [mg l-1] 0.576 0.591 0.603 2.38 MS [g] 1.6 1.4 2.1 75.3 MS/Minj [%] 0.55 0.46 0.69 25.3 EL07 CS [mg l-1] 7.95E-04 8.31E-04 8.25E-04 3.89 MS [g] 0.00 0.00 0.00 0.08 MS/Minj [%] 0.00 0.00 0.00 0.03 Table 5-11 : Ratios AS/A. Coefficients of variation (CVs) are calculated within reaches and across configurations. EL02 EL03 EL04 EL05 EL06 EL07 EL08 EL09 InBh 7.68E-02 8.79E-02 8.86E-02 7.84E-02 1.04E-01 9.63E-02 1.21E-01 1.25E-01 InDn 7.57E-02 8.60E-02 8.66E-02 7.84E-02 9.66E-02 9.59E-02 1.20E-01 1.24E-01 InUp 7.58E-02 8.59E-02 8.67E-02 7.80E-02 9.70E-02 9.60E-02 1.20E-01 1.25E-01 CV (%) 0.8 1.3 1.2 0.3 4.2 0.2 0.6 0.1 125 Figure 5.1 : Locations of study reaches. The shaded area encompasses the two headwaters and the ephemeral fraction of the stream network. Figure 5.2 : Diel discharge fluctuations during the data collection period. 126 Figure 5.3 : Spatial organization of lateral fluxes as implemented in OTIS with configurations (a) Inflow Both (InBh), (b) Inflow Upstream (InUp), (c) Inflow Downstream (InDn) and (d) Standard Application (SdAp). QIN, QOUT, qLatIn and qLatOut refer respectively to reach inflow, outflow, gross lateral inflow and gross lateral outflow. qLatIn a QIN QOUT qLatOut InBh c InDn b InUp d or SdAp 127 Figure 5.4 : Simulated stream solute breakthrough curve for reaches EL02 and EL07 in the arithmetic (a, b) and semi- logarithmic (c, d) spaces from configuration SdAp (black) and configurations InBh, InDn and InUp (grey). The experimental BTC is plotted with open circles. Note that the simulated BTCs for InBh, InDn and InUp are visually indistinguishable. 128 Figure 5.5 : Parameter estimates for all four configurations (qLatOut is not defined for configuration SdAp). n=8 for all boxes. AS values from configuration SdAp are plotted against the secondary axis because they are an order of magnitude greater than for the other configurations. Figure 5.6 : Derived metrics for all four configurations (n = 8 reaches for each parameter). 129 Figure 5.7 : In-stream (solid line) and transient storage zone (dashed line) breakthrough curves from InBh-InDn-InUp simulations (grey) and SdAp simulation (black). Experimental BTC is plotted with clear circles. Figure 5.8 : Damköhler indices for all four configurations. 130 Figure 5.9 : Composite scaled sensitivities across eight reaches and configuration SdAp (small solid dots) and the three other configurations (large clear circles). Values are expressed relative to the sensitivity of the parameter to which the model was the most sensitive (A). 131 Figure 5.10 : Parameter correlation coefficients (pcc) across eight reaches for configuration SdAp (small solid dots) and the three other configurations (large clear circles). Figure 5.11 : Discharge longitudinal profile at reach EL08 for configurations InBh, InDn and InUp. 132 6 Conclusion This research examined the dynamics of the linkages between slope, riparian zone and stream in a meso-scale montane catchment responding dominantly to snowmelt. It focused in particular on the spatial distribution of slope water delivery to the stream at a range of flow regimes, and then examined how these lateral contributions are affected by vegetation water demand and affect transient storage and in-stream solute transport processes during the baseflow period. 6.1 Summary of key findings The analysis of discrete flow increments presented in Chapter 3 provided an original perspective on the exploration of a catchment’s dynamics. This chapter highlighted the spatial variability of slope water delivery to the stream, with spatial patterns changing from high to low flows. As soil water deficit increased, and hydrologic connectivity between stream and upslope areas was sustained by fewer flow pathways, dominant controls of slope water delivery to the stream shifted. During the wet state, macro- scale descriptors such as accumulated area and median distance to the stream controlled slope water delivery to the stream. Hydrologic connectivity between stream and upland areas was then sustained relatively homogeneously throughout the whole catchment. During baseflow, an index representing the spatial arrangement of accumulated area in connection with the stream was the dominant driver, thereby suggesting that hydrologic connectivity was sustained by fewer, more specific flow pathways. The coupling between slopes, riparian aquifer and stream was examined in Chapter 4 through the analysis of streamflow diel fluctuations during baseflow recession. Evapotranspiration (as indexed by vapour pressure deficit, VPD) was the sole driver of streamflow diel fluctuations in CCEW during this period. The combined use of a one-dimensional model accounting for in-stream advective and dispersive processes and distributed 133 lateral inflows, and an analytical model accounting for fluxes between upslope area, riparian aquifer and stream, shed some light on the origins of diel perturbations that were recorded at the nine nested stream gauges within CCEW. Streamflow response time to VPD increased with baseflow recession, suggesting that perturbations transitioned from relatively short and shallow flow pathways to relatively longer and deeper as flow receded. However, streamflow response time to VPD forcing did not show a consistent increase with catchment scale or maximum stream length, indicating that in- stream advection was not the only process responsible for delaying evapotranspiration perturbations to a downstream gauge. Simulations were capable of effectively reproducing a dampening of perturbations due to in- stream energy loss. Therefore, streamflow diel fluctuations were likely the aggregated response to local perturbations originating in the vicinity of the stream gauge, and to perturbations originating from more distant sources that were advected downstream with no significant dampening. Further simulations indicated that the recession characteristics of riparian aquifers can modulate the response time of perturbations from their source to their point of entry in the channel. This mechanism provided a plausible explanation to observations of diel fluctuations at a stream gauge lagging behind those of another stream gauge located downstream. Research conducted in Chapter 5 focused on the parameterization of solute transient storage and in-stream transport problems. It was found that incorrect specification of the spatial patterns of lateral inflow and outflows may affect the model structure and confound comparisons of transient storage model parameters among reaches or as a function of time at a single reach. Moreover, the stream BTC alone does not provide enough information to select the appropriate model structure and then allow for a reliable interpretation of solute transport dynamics from model parameters. 134 6.2 Limitations and future directions 6.2.1 The need for additional measurements While the analysis of flow increments allowed for the identification of dominant drivers of slope water delivery to the stream, between 24% and 58% of the variance remained unaccounted for. One main source of uncertainty stemmed from the assumption that gross losses were negligible, and net gains were equal to gross gains (except for net losing reaches). Recent findings from Payn et al. [2009] demonstrated that gross fluxes at the reach scale could be up to one order of magnitude larger than their resulting net budget. It is, however, rather impractical to apply their method for measuring gross lateral fluxes in meso-scale catchments due to time constraints. Our approach could be greatly improved by selecting a stream network in which losses to groundwater aquifers or to the hyporheic zone are negligible (e.g., [Gooseff et al., 2005]). Alternatively, water chemistry from the stream and along upslope and riparian flow pathways should constrain reach-scale ionic mass balances and assist in the derivation of gross flux estimates. Moreover, water chemistry could be used to further explain the patterns of lateral contributions to the stream by drawing inferences on the sources of water and trajectories of flow pathways in relation to bedrock composition. In Chapter 3, detailed knowledge of surficial materials throughout CCEW did not explain the spatial patterns of slope water delivery to the stream. However, surficial material may not always reflect the underlying bedrock due to subsequent alteration processes. Therefore, accounting for bedrock characteristics may prove to be a useful addition to the set of predictors for explaining slope water delivery to the stream, especially during baseflow periods. Simulation results from the analysis of streamflow diel fluctuations highlighted the need for additional well and piezometer data in the riparian zone and lower slopes in order to better understand the dynamics of diel fluctuations from their source to the downstream gauge. Water table elevation records along transversal transects should prove useful to assess 135 the hydrologic connectivity between slope, riparian aquifer and stream, and the extent of the zone of influence away from the stream (e.g., [Jencso et al., 2009]). Water table elevation records along longitudinal transects would shed some light on the synchroneity between the stream and the riparian aquifers. The analytical model presented in Chapter 4 could benefit from these additional data. First, multiple independent riparian aquifers distributed along the channel could be implemented, each with its own characteristics. Second, water table elevation records would greatly help parameterize these riparian aquifers. Research presented in Chapter 5 has demonstrated that parameterization and structure of a Transient Storage Model were dependent on the specification of lateral fluxes from and to the stream. These fluxes are often overlooked in standard applications of the model, which results in ill- defined problems and biased parameter estimates. Here again, water table elevation records in the riparian aquifers may assist in the calculation of the field of hydraulic head in the vicinity of the stream network and in the derivation of fluxes from and to the stream. Alternatively, further small-scale tracer experiments, in which boundary conditions are sufficiently well defined, are necessary to assess the influence of varying boundary conditions on resulting mode parameters. Without such constrained tracer experiments, our ability to investigate further the effect of lateral fluxes on solute transport dynamics will be limited. This dissertation examined catchment dynamics from a stream perspective. All three research chapters focused on the analysis of streamflow data. However, it appears that all of them would have benefited from additional piezometer or well data on slopes and in riparian aquifers. 6.2.2 The merits and challenges of nested designs The preponderance of spatial and temporal scales in observing and interpreting hydrological processes has been acknowledged earlier in this dissertation, and widely by the hydrologic community. A number of studies have focused on the dependence of specific characteristics of the hydrological 136 response on catchment scale (e.g., [Pilgrim et al., 1982; Goodrich et al., 1997; Brutsaert and Lopez, 1998; Lundquist and Cayan, 2002]). The catchments of interest may be nested, or independent (with no spatial overlap). While a sample made up of independent catchments meets the assumption of independence of observations for further statistical analysis, nested designs present some powerful advantages. Nested designs offer a continuum in spatial scales over which hydrological processes can be examined. By considering successively the next basin scale up and its hydrological response in a bottom-up aggregative approach, small-scale description of processes remains (in the inner smallest basins), while the large-scale response of the outer basins is also available. Therefore, it is the ideal design to investigate scale-related issues and especially emergent behaviours, one of the focuses of the “future of hydrology”, as advocated by Wagener et al. [2010]. Moreover, in-stream transport processes are conveniently examined in nested designs. The advection of diel streamflow fluctuations was examined in Chapter 4, and that of solute in Chapter 5. Additional experiments conducted in CCEW involved a long term (12 days) constant rate injection of fluorescent dye within the Elk Creek sub-basin. The dye breakthrough curve was continuously recorded at the various stream gauges down the injection points and discrete synoptic sampling were conducted during the steady phase. While the data need further work and were not presented in this dissertation, they are expected to reveal transient storage and solute transport dynamics at a range of nested spatial scales. As mentioned earlier, nested designs pose the problem of the lack of independence of observations. This is surmountable by simply extracting independent sub-catchments such as Incremental Contributing Areas (ICA) or Incremental Sub-Basins (ISB), as illustrated in Chapter 3. This, however, presents some limitations, mainly with regards to uncertainty. Differences in streamflow measurements had to be computed to derive streamflow contribution from individual ICAs which, in turn, increased the uncertainty of streamflow increments. Mixed-effect models that can handle hierarchical data structures provide an alternative method in which sub-basins are not 137 required to be independent. 6.2.3 The need to account for hydrologic connectivity Hydrologic connectivity was a fundamental element of the present work. While its mechanisms were not examined in detail, it acted nonetheless as an underlying thread upon which various aspects of catchment dynamics could be explored. The relevance of studying both its mechanisms and the ways in which it affects other hydrological processes has been reinforced by recent publications. Wagener et al. [2007] proposed that catchment functions – partition, storage and release of water – should always be examined in conjunction with hydrologic connectivity in their proposed theory on catchment classification. Broadening the picture, Wagener et al. [2010] reiterated their call for a more holistic approach in hydrology with “enhanced integration across disciplines”, inviting hydrologists, ecologists, biologists and chemists to collaborate toward a coherent understanding of catchment dynamics. 138 References Adams, J., and J. Chandler (2002), Evaluation of lidar and medium scale photogrammetry for detecting soft-cliff coastal change, Photogrammetric Record, 17(99), 405-418. Ali, G. A., and A. G. Roy (2008), Revisiting hydrologic sampling strategies for an accurate assessment of hydrologic connectivity in humid temperate systems, Geography Compass, 3(1), 350-374, doi: 10.1111/j.1749- 8198.2008.00180.x. Allis, J. A. (1962), Comparison of storm runoff volumes from small, single- crop watersheds, Agricultural Engineering, (43), 220-223. Ambroise, B. (2004), Variable, 'active' versus 'contributing' areas or periods: a necessary distinction, Hydrological Processes, 18(6), 1149-1155, doi: 10.1002/hyp.5536 ER. Anderson, M., and T. Burt (1978), The role of topography in controlling throughflow generation, Earth Surface Processes, 3(4), 331-344. Anderson JK, Wondzell SM, Gooseff MN, Haggerty R. Patterns in stream longitudinal profiles and implications for hyporheic exchange flow at the H.J. Andrews Experimental Forest, Oregon, USA, Hydrological Processes, 2005;19(15):2931-49. Anderson, A. E., M. Weiler, Y. Alila, and R. O. Hudson (2009), Dye staining and excavation of a lateral preferential flow network, Hydrology and Earth System Sciences, 13(6), 935-944. Armbruster, J. T. (1976), Infiltration index useful in estimating low-flow characteristics of drainage basins, Journal of Research of the US Geological Survey, 4(5), 533-538. Baldocchi, D. D., C. A. Vogel, and B. Hall (1997), Seasonal variation of energy and water vapor exchange rates above and below a boreal jack pine forest canopy, Journal of Geophysical Research-Atmospheres, 102(D24), 28939- 139 28951. Baldocchi, D., F. M. Kelliher, T. A. Black, and P. Jarvis (2000), Climate and vegetation controls on boreal zone energy exchange, Global Change Biology, 6, 69-83. Bencala, K. E., and R. A. Walters (1983), Simulation of solute transport in a mountain pool-and-riffle stream - a Transient Storage Model, Water Resources Research, 19(3), 718-724. Beven, K. J., and M. J. Kirkby (1979), A physically based, variable contributing area model of basin hydrology / Un modèle à base physique de zone d'appel variable de l'hydrologie du bassin versant, Hydrological Sciences Journal, 24(1), 43-69. Bond, B. J., J. A. Jones, G. Moore, N. Phillips, D. Post, and J. J. McDonnell (2002), The zone of vegetation influence on baseflow revealed by diet patterns of streamflow and vegetation water use in a headwater basin, Hydrological Processes, 16(8), 1671-1677, doi: 10.1002/hyp.5022. Boronina, A., S. Golubev, and W. Balderer (2005), Estimation of actual evapotranspiration from an alluvial aquifer of the Kouris catchment (Cyprus) using continuous streamflow records, Hydrological Processes, 19(20), 4055-4068, doi: 10.1002/hyp.5871. Bracken, L. J., and J. Croke (2007), The concept of hydrological connectivity and its contribution to understanding runoff-dominated geomorphic systems, Hydrological Processes, 21(13), 1749-1763, doi: 10.1002/hyp.6313. Bren, L. J. (1997), Effects of slope vegetation removal on the diurnal variations of a small mountain stream, Water Resources Research, 33(2), 321-331. Brutsaert, W., and J. P. Lopez (1998), Basin-scale geohydrologic drought flow features of riparian aquifers in the southern Great Plains, Water Resources Research, 34(2), 233-240. Butler, J. J., Jr., G. J. Kluitenberg, D. O. Whittemore, Loheide,Steven P.,,II, W. Jin, M. A. Billinger, and X. Zhan (2007), A field investigation of phreatophyte-induced fluctuations in the water table, Water Resources Research, 43(2), W02404, doi: 10.1029/2005WR004627. 140 Buttle, J. M., P. J. Dillon, and G. R. Eerkes (2004), Hydrologic coupling of slopes, riparian zones and streams: an example from the Canadian Shield, Journal of Hydrology, 287(1-4), 161-177, doi: 10.1016/j.jhydrol.2003.09.022. Cazorzi, F., and G. DallaFontana (1996), Snowmelt modelling by combining air temperature and a distributed radiation index, Journal of Hydrology, 181(1-4), 169-187. Czikowsky, M. J., and D. R. Fitzjarrald (2004), Evidence of seasonal changes in evapotranspiration in eastern US hydrological records, Journal of Hydrometeorology., 5(5), 974-988. Day, T. J. (1976), On the precision of salt dilution gauging, Journal of Hydrology, 31(3-4), 293-306. Devito, K. J., A. R. Hill, and N. Roulet (1996), Groundwater-surface water interactions in headwater forested wetlands of the Canadian shield, Journal of Hydrology, 181(1-4), 127-147, doi: 10.1016/0022- 1694(95)02912-5. Dooge, J. C. I. (1986), Looking for hydrologic laws, Water Resources Research, 22(9), 465-585, doi: 10.1029/WR022i09Sp0046S. Dunne, T., and R. D. Black (1970), An experimental investigation of runoff production in permeable soils, Water Resources Research, 6(2), 478-490, doi: 10.1029/WR006i002p00478. Eagleson, P. S. (1972), Dynamics of flood frequency, Water Resources Research, 8(4), 878. Ensign, S. H., and M. W. Doyle (2005), In-channel transient storage and associated nutrient retention: Evidence from experimental manipulations, Limnology and Oceanography, 50(6), 1740-1751. Ferguson, R. I. (1999), Snowmelt runoff models, Progress in Physical Geography, 23(2), 205-227. Fleming, S. W., P. Hudson, and E. J. Quilty (2009), Interpreting nonstationary environmental cycles as amplitude-modulated (AM) signals, Canadian Journal of Civil Engineering, 36(4), 720-731, doi: 10.1139/S08-051. Freeman, M. C., C. M. Pringle, and C. R. Jackson (2007), Hydrologic 141 connectivity and the contribution of stream headwaters to ecological integrity at regional scales, Journal of the American Water Resources Association., 43(1), 5-14, doi: 10.1111/j.1752-1688.2007.00002.x ER. Freer, J., J. McDonnell, K. J. Beven, D. Brammer, D. Burns, R. P. Hooper, and C. Kendall (1997), Topographic controls on subsurface storm flow at the hillslope scale for two hydrologically distinct small catchments, Hydrological Processes, 11(9), 1347-1352. Freer, J., J. J. McDonnell, K. J. Beven, N. E. Peters, D. A. Burns, R. P. Hooper, B. Aulenbach, and C. Kendall (2002), The role of bedrock topography on subsurface storm flow, Water Resources Research, 38(12). Freeze, R. A. (1972), Role of subsurface flow in generating surface runoff 2. Upstream source areas, Water Resources Research, 8(5), 1272-1283. Genereux, D. P., H. F. Hemond, and P. J. Mulholland (1993), Spatial and temporal variability in streamflow generation on the west fork of walker branch watershed, Journal of Hydrology, 142(1-4), 137-166, doi: 10.1016/0022-1694(93)90009-X.. Gomi, T., R. C. Sidle, and J. S. Richardson (2002), Understanding processes and downstream linkages of headwater systems, Bioscience, 52(10), 905- 916. Goodrich, D. C., L. J. Lane, R. M. Shillito, S. N. Miller, K. H. Syed, and D. A. Woolhiser (1997), Linearity of basin response as a function of scale in a semiarid watershed, Water Resources Research, 33(12), 2951-2966, doi: doi:10.1029/97WR01422. Gooseff, M. N., J. LaNier, R. Haggerty, and K. Kokkeler (2005), Determining in-channel (dead zone) transient storage by comparing solute transport in a bedrock channel-alluvial channel sequence, Oregon, Water Resources Research, 41(6), W06014, doi: 10.1029/2004WR003513. Gooseff, M. N., R. A. Payn, J. P. Zarnetske, W. B. Bowden, J. P. McNamara, and J. H. Bradford (2008), Comparison of in-channel mobile-immobile zone exchange during instantaneous and constant rate stream tracer additions: Implications for design and interpretation of non-conservative tracer experiments, Journal of Hydrology, 357(1-2), 112-124, doi: 10.1016/j.jhydrol.2008.05.006. 142 Gooseff, M. N., S. M. Wondzell, R. Haggerty, and J. Anderson (2003), Comparing transient storage modeling and residence time distribution (RTD) analysis in geomorphically varied reaches in the Lookout Creek basin, Oregon, USA, Advances in Water Resources, 26(9), 925-937. Granier, A., P. Biron, and D. Lemoine (2000), Water balance, transpiration and canopy conductance in two beech stands, Agricultural and Forest Meteorology., 100(4), 291-308. Grayson, R. B., C. J. Gippel, B. L. Finlayson, and B. T. Hart (1997), Catchment-wide impacts on water quality: the use of `snapshot' sampling during stable flow, Journal of Hydrology, 199(1-2), 121-134. Gregory, S. V., F. J. Swanson, W. A. Mckee, and K. W. Cummins (1991), An ecosystem perspective of riparian zones, Bioscience, 41(8), 540-551. Gribovszki, Z., P. Kalicz, J. Szilagyi, and M. Kucsara (2008), Riparian zone evapotranspiration estimation from diurnal groundwater level fluctuations, Journal of Hydrology, 349(1-2), 6-17, doi: 10.1016/j.jhydrol.2007.10.049. 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-3480, doi: doi:10.1029/2000WR900214. 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). Hall, R. O., E. S. Bernhardt, and G. E. Likens (2002), Relating nutrient uptake with transient storage in forested mountain streams, Limnology and Oceanography, 47(1), 255-265. Halleran, W. H., (2004), Terrain and surficial geology map of Cotton Creek Experimental Watershed. Unpublished report prepared for Tembec Inc. B.C. Division, Cranbrook, BC. Hamilton, A. S., and R. D. Moore (1996), Winter streamflow variability in two groundwater-fed sub-Arctic rivers, Yukon Territory, Canada, Canadian Journal of Civil Engineering, 23(6), 1249-1259. Hamilton, A. S. (2008), Sources of uncertainty in Canadian low flow hydrometric data, Canadian Water Resources Journal, 33(2), 125-136. 143 Harman, C. J., M. Sivapalan, and P. Kumar (2009), Power law catchment- scale recessions arising from heterogeneous linear small-scale dynamics, Water Resources Research, 45, doi: 10.1029/2008WR007392 ER. Harvey, J. W., B. J. Wagner, and K. E. Bencala (1996), Evaluating the reliability of the stream tracer approach to characterize stream- subsurface water exchange, Water Resources Research, 32(8), 2441-2451. Henderson, F. M. (1966), Open Channel Flow, 522 pp., Macmillan, New York. Hewlett, J. D., and A. R. Hibbert (1967), Factors affecting the response of small watersheds to precipitation in humid areas. edited by W. E. Sopper and H. W. Lull, pp. 275-290, Pergamon Press, New York. Hewlett, J. D., and W. L. Nutter (1970), The varying source area of streamflow from upland basins, Interdisciplinary Aspects of Watershed Management. Hill, M. C., and C. R. Tiedeman (Eds.) (2007), Effective Groundwater Model Calibration: With Analysis of Data, Sensitivities, Predictions, and Uncertainty, 464 pp., Wiley and Sons, Inc., Hoboken, New Jersey. Hjerdt, K. N., J. J. McDonnell, J. Seibert, and A. Rodhe (2004), A new topographic index to quantify downslope controls on local drainage, Water Resources Research, 40(5), W05602. Hogg, E. H., and P. A. Hurdle (1997), Sap flow in trembling aspen: implications for stomatal responses to vapor pressure deficit, Tree Physiology, 17(8-9), 501-509. Hooke, J. M. (2007), Complexity, self-organisation and variation in behaviour in meandering rivers, Geomorphology, 91(3-4), 236-258, doi: 10.1016/j.geomorph.2007.04.021. Hoover, M. D., and C. R. Hursh (1943), Influence of topography and soil- depth on runoff from forest land, Transactions-American Geophysical Union, 24, 693-698. Huff, D. D., R. V. O'Neill, W. R. Emanuel, J. W. Elwood, and J. D. Newbold (1982), Flow variability and hillslope hydrology, Earth Surface Processes and Landforms, 7(1), 91-94. Hursh, C. R. (1944), Report of the subcommittee on subsurface flow. Eos, Transaction of American Geophysical Union, 25, 743-746. 144 Hutchinson, D. G., and R. D. Moore (2000), Throughflow variability on a forested hillslope underlain by compacted glacial till, Hydrological Processes, 14(10), 1751-1766. James, A. L., and N. T. Roulet (2007), Investigating hydrologic connectivity and its association with threshold change in runoff response in a temperate forested watershed, Hydrological Processes, 21(25), 3391-3408, doi: 10.1002/hyp.6554 ER. Jackman, A. P., R. A. Walters, and V. C. Kennedy (1984), Transport and concentration controls for chloride, strontium, potassium and lead in Uvas Creek, a small cobble-bed stream in Santa Clara County, California, U.S.A. : 2. Mathematical modeling, Journal of Hydrology, 75(1-4), 111-141. Jencso, K. G., B. L. McGlynn, M. N. Gooseff, S. M. Wondzell, K. E. Bencala, and L. A. Marshall (2009), Hydrologic connectivity between landscapes and streams: Transferring reach-and plot-scale understanding to the catchment scale, Water Resources Research, 45, doi: 10.1029/2008WR007225 ER. Jost, G., R. D. Moore, M. Weiler, D. R. Gluns, and Y. Alila (2009), Use of distributed snow measurements to test and improve a snowmelt model for predicting the effect of forest clear-cutting, Journal of Hydrology, 376(1-2), 94-106. Jost, G., M. Weiler, D. R. Gluns, and Y. Alila (2007), The influence of forest and topography on snow accumulation and melt at the watershed-scale, Journal of Hydrology, 347(1-2), 101-115. Kaiser, H. F. (1961), A note on Guttman lower bound for the number of common factors, British Journal of Statistical Psychology, 14(1), 1-2. Kavanagh, K. L., R. Pangle, and A. D. Schotzko (2007), Nocturnal transpiration causing disequilibrium between soil and stem predawn water potential in mixed conifer forests of Idaho, Tree Physiol., 27(4), 621-629. Keefe, S. H., L. B. Barber, R. L. Runkel, J. N. Ryan, D. M. McKnight, and R. D. Wass (2004), Conservative and reactive solute transport in constructed wetlands, Water Resources Research, 40(1), doi: 10.1029/2003WR002130 145 ER. Kelliher, F. M., R. Leuning, and E. D. Schulze (1993), Evaporation and canopy characteristics of coniferous forests and grasslands, Oecologia, 95(2), 153-163. Kincaid, D. R., H. B. Osborn, and J. R. Gardner (1966), Use of unit-source watershed for hydrologic investigations in the semiarid southwest, Water Resources Research, (2), 381-392. Kirchner, J. W. (2006), Getting the right answers for the right reasons: Linking measurements, analyses, and models to advance the science of hydrology, Water Resources Research., 42(3). Kirchner, J. W. (2003), A double paradox in catchment hydrology and geochemistry, Hydrological Processes, 17(4), 871-874. Kirkby, M. J. (1978), Hillslope Hydrology, 389 pp., John Wiley, Chichester, UK. Klemeš, V. (1983), Conceptualization and scale in hydrology, Journal of Hydrology, 65(1-3), 1-23, doi: 10.1016/0022-1694(83)90208-1. Koutsoyiannis, D., C. Makropoulos, A. Langousis, S. Baki, A. Efstratiadis, A. Christofides, G. Karavokiros, and N. Mamassis (2009), HESS Opinions: 'Climate, hydrology, energy, water: recognizing uncertainty and seeking sustainability', Hydrology and Earth System Sciences, 13(2), 247-257. Kuraś, P. K., M. Weiler, and Y. Alila (2008), The spatiotemporal variability of runoff generation and groundwater dynamics in a snow-dominated catchment, Journal of Hydrology, 352(1-2), 50-66. Kutner, M., C. Nachtsheim, J. Neter, and W. Li (2004), Applied linear statistical models, 5th ed., 1396 pp., McGraw-Hill, New York, NY. Lautz, L. K., and D. I. Siegel (2007), The effect of transient storage on nitrate uptake lengths in streams: an inter-site comparison, Hydrological Processes, 21(26), 3533-3548, doi: 10.1002/hyp.6569 ER. Leavesley, G. H., R. W. Litchy, B. M. Troutman, and L. G. Saindon (1983), Precipitation-Runoff Modeling System: User's Manual. , 83 pp., Denver, Colorado. Leavesley, G. H., S. L. Markstrom, P. J. Restrepo, and R. J. Viger (2002), A modular approach to addressing model design, scale, and parameter 146 estimation issues in distributed hydrological modelling, Hydrological Processes, 16(2), 173-187. Lehmann, P., C. Hinz, G. McGrath, H. J. Tromp-van Meerveld, and J. J. McDonnell (2007), Rainfall threshold for hillslope outflow: an emergent property of flow pathway connectivity, Hydrology and Earth System Sciences, 11(2), 1047-1063. Leopold, L. B., M. G. Wolman, and J. P. Miller (1995), Fluvial processes in geomorphology, Dover Pubns. Loheide, S. P., J. J. Butler, and S. M. Gorelick (2005), Estimation of groundwater consumption by phreatophytes using diurnal water table fluctuations: A saturated-unsaturated flow assessment, Water Resources Research, 41(7), W07030, doi: 10.1029/2005WR003942. Lundquist, J. D., and D. R. Cayan (2002), Seasonal and spatial patterns in diurnal cycles in streamflow in the western United States, Journal of Hydrometeorology, 3(5), 591-603. McDonnell, J. J. (1990), A rationale for old water discharge through macropores in a steep, humid catchment, Water Resources Research, 26(11), 2821-2832. McDonnell, J. J., M. Sivapalan, K. Vache´, S. Dunn, G. Grant, R. Haggerty, C. Hinz, R. Hooper, J. Kirchner, M. L. Roderick, J. Selker, and M. Weiler (2007), Moving beyond heterogeneity and process complexity: A new vision for watershed hydrology, Water Resources Research, 43, doi: 10.1029/2006WR005467 ER. McGlynn, B. L., J. J. McDonnell, and D. D. Brammer (2002), A review of the evolving perceptual model of hillslope flowpaths at the Maimai catchments, New Zealand, Journal of Hydrology (Amsterdam), 257(1/4), 1-26. McGlynn, B. L., and J. J. McDonnell (2003), Role of discrete landscape units in controlling catchment dissolved organic carbon dynamics, Water Resources Research, 39(4), doi: 10.1029/2002WR001525 ER. McGlynn, B. L., and J. Seibert (2003), Distributed assessment of contributing area and riparian buffering along stream networks, Water Resources Research, 39(4), 1082-1088. 147 McNamara, J. P., D. Chandler, M. Seyfried, and S. Achet (2005), Soil moisture states, lateral flow, and streamflow generation in a semi-arid, snowmelt-driven catchment, Hydrological Processes, 19(20), 4023-4038. Monteith, S. S., J. M. Buttle, P. W. Hazlett, F. D. Beall, R. G. Semkin, and D. S. Jeffries (2006), Paired-basin comparison of hydrological response in harvested and undisturbed hardwood forests during snowmelt in central Ontario: I. Streamflow, groundwater and flowpath behaviour, Hydrological Processes, 20(5), 1095-1116. Moore, G. W., B. J. Bond, J. A. Jones, N. Phillips, and F. C. Meinzer (2004), Structural and compositional controls on transpiration in 40-and 450- year-old riparian forests in western Oregon, USA, Tree Physiology, 24(5), 481-491. Moore, R. D., and J. C. Thompson (1996), Are water table variations in a shallow forest soil consistent with the TOPMODEL concept? Water Resources Research, 32(3), 663-669. Moore, R. D. (2005), Slug injection using salt in solution, Streamline Watershed Management Bulletin, 8(2), 1-6. Mosley, M. P. (1979), Streamflow generation in a forested watershed, New- Zealand, Water Resources Research, 15(4), 795-806. Mulholland, P. J. (1993), Hydrometric and stream chemistry evidence of three storm flowpaths in walker branch watershed, Journal of Hydrology, 151(2-4), 291-316. O'Callaghan, J. F., and D. M. Mark (1984), The extraction of drainage networks from digital elevation data, Computer Vision Graphics and Image Processing, 28(3), 323-344. Ocampo, C. J., M. Sivapalan, and C. Oldham (2006), Hydrological connectivity of upland-riparian zones in agricultural catchments: Implications for runoff generation and nitrate transport, Journal of Hydrology, 331(3-4), 643-658, doi: 10.1016/j.jhydrol.2006.06.010 ER. Oreskes, N., K. Shraderfrechette, and K. Belitz (1994), Verification, validation, and confirmation of numerical models in the earth-sciences, Science, 263(5147), 641-646. Pacific, V. J., K. G. Jencso, and B. L. McGlynn (2010), Variable flushing 148 mechanisms and landscape structure control stream DOC export during snowmelt in a set of nested catchments, Biogeochemistry, , 1-19. Payn, R. A., M. N. Gooseff, B. L. McGlynn, K. E. Bencala, and S. M. Wondzell (2009), Channel water balance and exchange with subsurface flow along a mountain headwater stream in Montana, United States, Water Resources Research, 45. Peters, D. L., J. M. Buttle, C. H. Taylor, and B. D. LaZerte (1995), Runoff production in a forested, shallow soil, Canadian Shield basin, Water Resources Research, 31(5), 1291-1304. Pilgrim, D. H., I. Cordery, and B. C. Baron (1982), Effects of catchment size on runoff relationships, Journal of Hydrology, 58(3-4), 205-221. Poeter, E. P., M. C. Hill, E. R. Banta, S. Mehl, and S. Christensen (2005), UCODE_2005 and Six Other Computer Codes for Universal Sensitivity Analysis, Calibration, and Uncertainty Evaluation, U. S. Geological Survey, Techniques and Methods, 6-A11, 299p. Pomeroy, J. W., D. M. Gray, K. R. Shook, B. Toth, R. L. H. Essery, A. Pietroniro, and N. Hedstrom (1998), An evaluation of snow accumulation and ablation processes for land surface modelling, Hydrological Processes, 12(15), 2339-2367. Pomeroy, J. W., D. M. Gray, T. Brown, N. R. Hedstrom, W. L. Quinton, R. J. Granger, and S. K. Carey (2007), The cold regions hydrological process representation and model: a platform for basing model structure on physical evidence, Hydrological Processes, 21(19), 2650-2667, doi: 10.1002/hyp.6787. Pringle, C. M. (2001), Hydrologic connectivity and the management of biological reserves: A global perspective, Ecological Applications, 11(4), 981-998. Rao, S. T., and I. G. Zurbenko (1994), Detecting and tracking changes in ozone air-quality, Journal of the Air and Waste Management Association, 44(9), 1089-1092. Reggiani, P., M. Sivapalan, and S. M. Hassanizadeh (1998), A unifying framework for watershed thermodynamics: balance equations for mass, momentum, energy and entropy, and the second law of thermodynamics, 149 Advances in Water Resources, 22(4), 367-398. Renard, K. G. (1977), Past, present and future water resources research in arid and semiarid areas of the southwestern United States, The Institution of Engineers, Australia, National Conference Publication, 77(5), 1-29. Rodríguez-Iturbe, I., and J. B. Valdés (1979), The geomorphologic structure of hydrologic response, Water Resources Research, 15(6), 1409-1420. Roulet, N. T. (1990), The hydrological role of peat-covered wetlands, Canadian Geographer/Le Géographe canadien, 34(1), 82-83. Ruehl, C., A. T. Fisher, C. Hatch, M. Los Huertos, G. Stemler, and C. Shennan (2006), Differential gauging and tracer tests resolve seepage fluxes in a strongly-losing stream, Journal of Hydrology, 330(1-2), 235- 248. Runkel, R. L. (1998), One dimensional transport with inflow and storage (OTIS): A solute transport model for streams and rivers, U.S. Geological Survey Water-Resources Investigation, Report 98-4018. Runkel, R. L. (2002), A new metric for determining the importance of transient storage, Journal of the North American Benthological Society, 21(4), 529-543. Schnorbus, M., and Y. Alila (2004), Forest harvesting impacts on the peak flow regime in the Columbia Mountains of southeastern British Columbia: An investigation using long-term numerical modeling, Water Resources Research, 40(5). Scordo, E. B., and R. D. Moore (2009), Transient storage processes in a steep headwater stream, Hydrological Processes, 23(18), 2671-2685, doi: 10.1002/hyp.7345 ER. Scott, D. T., M. N. Gooseff, K. E. Bencala, and R. L. Runkel (2003), Automated calibration of a stream solute transport model: implications for interpretation of biogeochemical parameters, Journal of the North American Benthological Society, 22(4), 492-510. Shaman, J., M. Stieglitz, and D. Burns (2004), Are big basins just the sum of small catchments? Hydrological Processes, 18(16), 3195-3206. Shreve, R. L. (1966), Statistical law of stream numbers, Journal of Geology, 150 74(1), 17-37. Sidle, R. C., Y. Tsuboyama, S. Noguchi, I. Hosoda, M. Fujieda, and T. Shimizu (2000), Stormflow generation in steep forested headwaters: a linked hydrogeomorphic paradigm, Hydrological Processes, 14(3), 369- 385. Sivapalan, M. (2005), Pattern, process and function: elements of a unified theory of hydrology at the catchment scale, in edited by M. G. Anderson, pp. 193-220, John Wiley & Sons, Chichester, UK. 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). Sivapalan, M., K. Takeuchi, S. W. Franks, V. K. Gupta, H. Karambiri, V. Lakshmi, X. Liang, J. J. Mcdonnell, E. M. Mendiondo, P. E. O’connell, T. Oki, J. W. Pomeroy, D. Schertzer, S. Uhlenbrook and E. Zehe, (2003), IAHS decade on Predictions in Ungauged Basins (PUB), 2003-2012: Shaping an exciting future for the hydrological sciences, Hydrological Sciences Journal-Journal Des Sciences Hydrologiques, 48(6), 857-880. Sivapalan, M. (2003), Process complexity at hillslope scale, process simplicity at the watershed scale: is there a connection? Hydrological Processes, 17(5), 1037-1041. Smakhtin, V. U. (2001), Low flow hydrology: a review, Journal of Hydrology, 240(3-4), 147-186. Sopper, W. E., and H. W. Lull (1965), The representativeness of small forested experimental watersheds in northeastern United States, I.A.H.S., (66), 441-456. Spence, C., and M. K. Woo (2003), Hydrology of subarctic Canadian shield: soil-filled valleys, Journal of Hydrology, 279(1-4), 151-166, doi: 10.1016/S0022-1694(03)00175-6 ER. Stieglitz, M., J. Shaman, J. McNamara, V. Engel, J. Shanley, and G. W. Kling (2003), An approach to understanding hydrologic connectivity on the hillslope and the implications for nutrient transport, Global Biogeochemical Cycles, 17(4), 16-1 - 16-15, doi:10.1029/2003GB002041. Stofleth, J. M., F. D. Shields, and G. A. Fox (2008), Hyporheic and total 151 transient storage in small, sand-bed streams, Hydrological Processes, 22(12), 1885-1894, doi: 10.1002/hyp.6773 ER. Story, A., R. D. Moore, and J. S. Macdonald (2003), Stream temperatures in two shaded reaches below cutblocks and logging roads: downstream cooling linked to subsurface hydrology, Canadian Journal of Forest Research, 33(8), 1383-1396. Su, J., and E. Bork (2006), Influence of vegetation, slope, and lidar sampling angle on DEM accuracy, Photogrammetric Engineering & Remote Sensing, 72(11), 1265-1274. Szilágyi, J., Z. Gribovszki, P. Kalicz, and M. Kucsara (2008), On diurnal riparian zone groundwater-level and streamflow fluctuations, Journal of Hydrology, 349(1-2), 1-5, doi: 10.1016/j.jhydrol.2007.09.014 ER. Tague, C., and G. E. Grant (2004), A geological framework for interpreting the low-flow regimes of Cascade streams, Willamette River Basin, Oregon, Water Resources Research, 40(4), W04303. Tani, M. (1997), Runoff generation processes estimated from hydrological observations on a steep forested hillslope with a thin soil layer, Journal of Hydrology, 200(1-4), 84-109. Taylor, C. H. (1982), The effect on storm runoff response of seasonal variations in contributing zones in small watersheds, Nordic Hydrology, 13(3), 165-182. Taylor, P. D., L. Fahrig, K. Henein, and G. Merriam (1993), Connectivity is a vital element of landscape structure, Oikos, 68(3), 571-573. Tetzlaff, D., C. Soulsby, S. Waldron, I. A. Malcolm, P. J. Bacon, S. M. Dunn, A. Lilly, and A. F. Youngson (2007), Conceptualization of runoff processes using a geographical information system and tracers in a nested mesoscale catchment, Hydrological Processes, 21(10), 1289-1307, doi: 10.1002/hyp.6309 ER. Troch, P. A., G. A. Carrillo, I. Heidbüchel, S. Rajagopal, M. Switanek, T. H. M. Volkmann, and M. Yaeger (2009), Dealing with landscape heterogeneity in watershed hydrology: a review of recent progress toward new hydrological theory, Geography Compass, 3(1), 375-392. Tromp-van Meerveld, H. J., and J. J. McDonnell (2006a), Threshold relations 152 in subsurface stormflow: 1. A 147-storm analysis of the Panola hillslope, Water Resources Research, 42(2), doi: 10.1029/2004WR003778 ER. Tromp-van Meerveld, H. J., and J. J. McDonnell (2006b), Threshold relations in subsurface stormflow: 2. The fill and spill hypothesis, Water Resources Research, 42(2), doi: 10.1029/2004WR003800 ER. Tsuboyama, Y., R. C. Sidle, S. Noguchi, and I. Hosoda (1994), Flow and solute transport through the soil matrix and macropores of a hillslope segment, Water Resources Research, 30(4), 879-890. Uchida, T., J. J. McDonnell, and Y. Asano (2006), Functional intercomparison of hillslopes and small catchments by examining water source, flowpath and mean residence time, Journal of Hydrology (Amsterdam), 327(3/4), 627-642. Wagener, T., M. Sivapalan, P. Troch, and R. Woods (2007), Catchment classification and hydrologic similarity, Geography Compass, 1(4), 901- 931, doi: 10.1111/j.1749-8198.2007.00039.x. Wagener, T., M. Sivapalan, P. A. Troch, B. L. McGlynn, C. J. Harman, H. V. Gupta, P. Kumar, P. S. C. Rao, N. B. Basu, and J. S. Wilson (2010), The future of hydrology: An evolving science for a changing world, Water Resources Research, 46, doi: 10.1029/2009WR008906 ER. Wagner, B. J., and J. W. Harvey (1997), Experimental design for estimating parameters of rate-limited mass transfer: Analysis of stream tracer studies, Water Resources Research, 33(7), 1731-1741. Weiler, M., and H. Fluhler (2004), Inferring flow types from dye patterns in macroporous soils, Geoderma, 120(1-2), 137-153. Weiler, M., and F. Naef (2003), An experimental tracer study of the role of macropores in infiltration in grassland soils, Hydrological Processes, 17(2), 477-493. Weiler, M., T. Uchida, and J. McDonnell (2003), Connectivity due to preferential flow controls water flow and solute transport at the hillslope scale, Proceedings of MODSIM, Townsville, Australia. White, W. N. (1932), A method of estimating ground-water supplies based on discharge by plants and evaporation from soil: results of investigations in Escalante Valley, Utah, USGPO. 153 Williams, K. S., and D. G. Tarboton (1999), The ABC's of snowmelt: a topographically factorized energy component snowmelt model, Hydrological Processes 13(12-13), 1905-1920. Wondzell, S. M. (2006), Effect of morphology and discharge on hyporheic exchange flows in two small streams in the Cascade Mountains of Oregon, USA, Hydrological Processes, 20(2), 267-287. Wondzell, S. M., M. N. Gooseff, and B. L. McGlynn (2007), Flow velocity and the hydrologic behavior of streams during baseflow, Geophysical Research Letter, 34(24), L24404, doi: 10.1029/2007GL031256. Wondzell, S. M., M. N. Gooseff, and B. L. McGlynn (2010), An analysis of alternative conceptual models relating hyporheic exchange flow to diel fluctuations in discharge during baseflow recession, Hydrological Processes, 24(6), 686-694, doi: 10.1002/hyp.7507. Worman, A. (1998), Analytical solution and timescale for transport of reacting solutes in rivers and streams, Water Resources Research, 34(10), 2703-2716. Worman, A. (2000), Comparison of models for transient storage of solutes in small streams, Water Resources Research, 36(2), 455-468. Zarnetske, J. P., M. N. Gooseff, T. R. Brosten, J. H. Bradford, J. P. McNamara, and W. B. Bowden (2007), Transient storage as a function of geomorphology, discharge, and permafrost active layer conditions in Arctic tundra streams, Water Resources Research, 43(7), W07410, doi: 10.1029/2005WR004816. Zurbenko, I. G. (1986), The spectral analysis of time series, 247 pp. Holland."@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "2011-05"@en ; edm:isShownAt "10.14288/1.0071446"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Forestry"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "Attribution-NonCommercial 3.0 Unported"@en ; ns0:rightsURI "http://creativecommons.org/licenses/by-nc/3.0/"@en ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Stream-catchment connectivity and streamflow dynamics in a montane landscape"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/29908"@en .