SIMULATION OF GROUNDWATER FLOW IN MOUNTAIN WATERSHEDS DEVIN W. CAIRNS M.Sc. Thesis, University of Lethbridge, 2014 A Thesis Submitted to the School of Graduate Studies of the University of Lethbridge in Partial Fulfilment of the Requirements for the Degree [Master of Science] Department of Geography University of Lethbridge LETHBRIDGE, ALBERTA, CANADA © Devin W. Cairns, 2014 SIMULATION OF GROUNDWATER FLOW IN MOUNTAIN WATERSHEDS DEVIN W. CAIRNS Date of Defence: December 9, 2014 Dr. James Byrne Professor Ph.D. Co-supervisor Dr. Dan Johnson Professor Ph.D. Co-supervisor Dr. Hester Jiskoot Associate Professor Ph.D. Thesis Examination Committee Member Dr. Jeff McKenzie Associate Professor Ph.D. Thesis Examination Committee Member Dr. Ivan Townsend Professor Ph.D. Chair, Thesis Examination Committee i Acknowledgements “Thanks for this…Good Work!” - Dr. Jim Byrne “I think that it looks like a really nice contribution to the fields of geomorphometric classification, surface geology classification, and groundwater modelling.” - Dr. Hester Jiskoot “The final product, the permeability map, in that document is really cool.” - Dr. Jeff McKenzie “Interesting - and reads well.” - Dr. Dan Johnson “Yay, you’re done!” - Tanya Byrne, M.Sc. “It is looking pretty good.” - Dr. Ryan MacDonald Words of encouragement are worth the world while tasked with challenging work and should never be forgotten. Each individual above contributed significantly to the completion of this thesis, and I am truly grateful for every second of time committed towards my success. Other academic support was generously provided by Tom Rientjes, Maciek Lubczynski and the University of Lethbridge Department of Geography. In addition, I likely would not have been able to maintain the required work ethic and possessed the necessary resources without the help of Lorene Cairns, Martin Cairns, Cindy Byrne, Paul Byrne, Isaac Cairns, Mackenzie Cairns, Sierra Cairns, Shayanne Cairns, Kelcy Byrne, Justin Byrne, Dr. William Ebert, Gail Ebert, Phyllis Cairns and Albert Cairns. ii Abstract Many mountain watersheds provide a reliable source of freshwater for habitat and human use downstream. Hydrodynamics of these basins can be particularly sensitive to change, which may arise from climate change, natural/anthropogenic alteration and other forcing. The timing and intensity of runoff and the recharge of groundwater must be at least partially understood to simulate the potential effects of change. The effects of groundwater on basin drainage are often neglected in simulations due to a lack of empirical data. This study focused on the integration of remote sensing data, geomorphic principles, theoretical distributions of heterogeneity, basin discretization methods, and saturated flow computation to apply a novel technique to understand groundwater behaviour in mountain watersheds. Methods included a geomorphometric analysis to computationally simulate the distribution of geomorphic landforms, which were used to estimate heterogeneity in the shallow subsurface and provide opportunity to evaluate groundwater flow. Morphometric attributes of various landforms were studied and compared to their genetic origin to identify potential landforms. The resulting landforms were subsequently divided into equivalent porous media units (EMUs) based on the theoretical distribution of heterogeneity within landform types. EMUs were evaluated as irregular units used to discretize a saturated groundwater flow model. Groundwater flow was calculated using recharge simulated by any hydrometeorologic model and was routed using Darcian flow from EMU to EMU. Methods of simulating groundwater flow in this study were found to be well suited for the basin type of the study area used (St. Mary Watershed, Montana, USA), albeit with limitations. Results of the geomorphometric analysis compared well with published surficial geology data. The basin discretization method presented in this research would benefit from implicit groundwater flow solving, and application in a basin where abundant data exist. An implicit scheme would allow faster computation and provide the means for a quantitative comparison of basin outflow and water table elevations, which would be useful to further evaluate the suitability of these techniques. iii Table of Contents Acknowledgements .............................................................................................................. ii Table of Contents .................................................................................................................iv List of Figures ..................................................................................................................... vii List of Tables ........................................................................................................................ ix Chapter 1: Introduction ........................................................................................................ 1 1.1. Research Objectives ........................................................................................................... 1 1.2. Background ......................................................................................................................... 1 Chapter 2: Previous Research and Opportunity for Advancement ............................. 4 GENESYS Hydrometeorological Model .......................................................................... 4 Surface and Groundwater Modeling in Alpine Watersheds ........................................ 4 Chapter 3 Abstract ................................................................................................................ 6 Chapter 3: Novel Geomorphometric Classification of Heterogeneity of Unconsolidated Sediments for Hydrogeologic Characterization in Mountain Watersheds ............................................................................................................................. 8 Introduction ............................................................................................................................. 8 Study Area and Dataset ...................................................................................................... 10 Methodology ......................................................................................................................... 12 3.3.1. Geomorphometric Analysis .................................................................................. 12 3.3.2. Model Indices and Quantitative Verification ................................................... 21 iv 3.3.3. Hydrogeological Characterization ...................................................................... 24 Results .................................................................................................................................... 26 3.4.1. Geomorphometric Model Results....................................................................... 26 3.4.2. Hydrogeological Characterization Results ....................................................... 29 Discussion .............................................................................................................................. 32 3.5.1. Geomorphometric Modeling ............................................................................... 32 3.5.2. Hydrogeologic Characterization ......................................................................... 33 3.5.3. Assumptions and Limitations .............................................................................. 35 Conclusion and Recommendations .................................................................................. 37 Chapter 4: Groundwater Model Discretization and Computation Using Geomorphic-Derived Data for Estimating Unconfined Topographically Driven Groundwater Flow in Mountain Watersheds ............................................................... 40 Introduction .......................................................................................................................... 40 4.1 Study Area .................................................................................................................... 42 Methodology ......................................................................................................................... 42 4.2.1. Spatial Discretization of Domain ........................................................................ 43 4.2.2. Groundwater Flow Simulation and Calibration .............................................. 45 Results .................................................................................................................................... 50 Assumptions and Limitations ........................................................................................... 57 Discussion .............................................................................................................................. 59 v Conclusions ........................................................................................................................... 60 Chapter 5: Thesis Conclusions and Recommendations .............................................. 62 References ............................................................................................................................ 65 vi List of Figures Figure 1 Study Area: Headwaters of the St. Mary River, Montana, USA. Hillshaded and elevation coloured DEM (USGS, 2009) ..................................................................................... 11 Figure 2 Geomorphic landform classification relationships .................................................. 14 Figure 3 Identification of alluvial fan/outwash regions and their delineation ................. 15 Figure 4 Process memberships used in model to identify fluvial features ......................... 16 Figure 5 Process memberships used in the model to identify remaining (sub) surface features (see also Figure 4 and Table 1) ..................................................................................... 18 Figure 6 Map of Geomorphometric Model Results: the overall representation of major geomorphometric features complies with the logical progression of the sediment cascade in glaciated terrain .......................................................................................................................... 27 Figure 7 Simulated distribution of hydrogeological facies. Higher hydraulic conductivity is generally noted at higher elevations, with an attenuation toward the bottom of valleys and the lower reaches due to lacustrine and glacial sediment dominance. ........................ 31 Figure 8 Conceptualization of methods used to assign geometric attributes to each EMU object (separated by colour in figure) for computation. Calculated widths between EMUs are denoted by red lines. ................................................................................................... 45 Figure 9 Modelled depth to bedrock using the approximation of the catenary valley shape and a maximum depth of 20 m at the valley bottom. During calibration, the maximum depth may be optimized and the depth re-calculated. ......................................... 51 Figure 10 Comparison of numeric instability prior to compliance with the von Neumann Criterion using applied fully saturated conditions throughout the watershed ................. 52 Figure 11 Attenuation of fully saturated conditions over 2000 days .................................. 54 vii Figure 12 At later stages in attenuation the presence of freshet is apparent in the hydrograph annually. ...................................................................................................................... 54 Figure 13 Distribution of cumulative groundwater discharge to the surface (not including basin outlet) after 5 years of iterative basin volume attenuation. A satellite image is included for comparison (ESRI, 2014) ....................................................................... 55 Figure 14 Comparison of groundwater discharge methods during iterative attenuation of fully saturated conditions after 14235 days with discharged groundwater being recharged downslope ...................................................................................................................... 56 Figure 15 Comparison of groundwater discharge methods during iterative attenuation of fully saturated conditions after 14235 days with no additional recharge ...................... 56 Figure 16 Water table elevation before (left) and after (right) the attenuation of fully saturated pressure conditions throughout valleys ................................................................... 57 viii List of Tables Table 1 Boundary indices chosen for geomorphometric model ........................................... 21 Table 2 Landform reclassification to create comparable datasets for model verification 23 Table 3 Results of the comparison between simulated and published geomorphology datasets ............................................................................................................................................... 28 Table 4 Results of the comparison between simulated and published geomorphology datasets- total area .......................................................................................................................... 29 Table 5 Estimates of effective porosity and hydraulic conductivity used in model runs (Batu, 1998) ....................................................................................................................................... 53 ix Chapter 1: Introduction 1.1. Research Objectives Techniques for coupling surface-groundwater models could benefit from alternative means of model discretization. Opportunity exists particularly in basin-scale hydrologic modelling, which is an important contributor to the evaluation of water security and climate change related challenges. This research addresses these issues with the introduction of novel techniques to approach basin-scale groundwater simulation. The approach presented in this study includes the conceptualization and implementation of two sequential studies, which includes:  Novel geomorphometric classification of heterogeneity of unconsolidated sediments for hydrogeologic characterization in mountain watersheds; and  Groundwater model discretization and computation using geomorphic-derived data for estimating unconfined topographically driven groundwater flow in mountain watersheds. This thesis is structured around the main thesis research outcomes: two separate papers addressing the above topics, which are intended to be published. The proposed papers are presented herein as Chapters 3 and 4. 1.2. Background Watershed management has become increasingly important in the wake of rapid global population growth which will affect food security and ecosystem health (UNEP, 2007). Continued global industrial growth and urbanization will also place increasing stresses on water resource management (Bell et al., 2012). Other factors such as climate change will only exacerbate these effects (van Leeuwen et al., 2012). Sustainable watershed management relies on robust knowledge of water at various stages in the hydrological cycle, which should never omit groundwater. Global groundwater abstraction rates are increasing with time, which will 1 place additional stresses on aquifers, while also contributing significantly to sea level rise (Wada et al., 2010). Surface water quality and quantity is also greatly dependent on groundwater at local and regional scales (Winter, 1999). Numerical simulations are a proven method of characterizing water resources at the watershed scale, especially when considering the importance of mountains and their contribution to watershed runoff (Marquis et al., 2012, Beniston et al., 1997). Groundwater recharge is typically considered in conceptual surface runoff models; however methods of approximating this phenomenon are often omitted and could benefit from development. Research focused on groundwater recharge at the watershed mass balance scale has exhibited a high degree of ingenuity and variety of approaches due to its challenges (Hayashi and van der Kamp, 2009, Kao et al., 2012). A recurring challenge for researchers has been the availability of data to utilize in the development of watershed groundwater models, as these data are crucial for their success (McCloskey and Finnemore, 1996, Scanlon et al., 2002). Hydrogeological data, such as hydraulic conductivity and the distribution of heterogeneity, often do not exist or are spatially interspersed within a watershed. The mountainous headwaters of watersheds are examples of areas that are large, complex and commonly lack data. As such, these regions are readily available for advancement in groundwater recharge and modeling research (Wilson and Guan, 2004). Geological characterizations are often achieved through the use of geophysical and borehole data. Borehole logs are used to interpolate a distribution of heterogeneity in the subsurface using data retrieved from a preinstalled network. The accuracy of geological characterizations is dependent on many factors, especially the density and distribution of available borehole data. Alternative approaches to subsurface characterization are necessary when borehole data are sparse, especially when study areas are not easily accessible (McCloskey and Finnemore, 1996). Successful intrusive studies are often difficult to complete due to resource limitations. Geophysical data collection, surficial geology maps and geochemical means are useful alternatives; however these methods are predicated on 2 pre-existing data or those that are costly to obtain. Given that hydrogeological characterizations are approximations of the subsurface; this research explores other methods of groundwater recharge approximation in large alpine watersheds. The merit of this research is confounded by the notion that simulations are a viable option in watershed management despite a lack of physical data. The ability to apply a spatial characterization of groundwater recharge with limited data will be a useful tool in linking surface hydrological models and groundwater models (Hevesi et al., 2002, Stone et al., 2001). 3 Chapter 2: Previous Research and Opportunity for Advancement GENESYS Hydrometeorological Model Generate Earth Systems Science (GENESYS) is a physical-based hydrometeorological model created under the supervision of Dr. James Byrne at the University of Lethbridge and has been continuously improved upon, such as the work by MacDonald and Byrne (2009). GENESYS simulates hydrological processes at a high spatial and temporal resolution in mountainous terrain using limited data from local weather stations. Various physical processes are integrated with ecohydrological parameters to simulate temporal water mass balance changes at the watershed scale (MacDonald et al., 2009, 2011, Lapp et al., 2005). GENESYS uses terrain categories, or hydrologic response units, to account for spatially variable ecohydrologic parameters in the watershed while performing the surface water mass balance simulations. Given that the snowpack is of paramount importance in alpine watersheds, GENESYS is highly proficient at simulating snowpack processes. Validation modeling was completed on the Saint Mary River watershed, which receives roughly 70% of its precipitation in the form of snow (MacDonald et al., 2009). GENESYS lacks groundwater processes in its routines, and also lacks water routing capabilities beyond the water budget calculations (MacDonald et al., 2009, 2011). Based on previous application, GENESYS has proven to be suitable for solving the spatially distributed water balance in diverse terrain with high computational efficiency. Surface and Groundwater Modeling in Alpine Watersheds Several recent studies have provided methods of approximating groundwater processes in alpine watersheds. Research regarding the physical principles associated with infiltration and groundwater derived from mountain watersheds is plentiful (Hayashi and van der Kamp, 2009, Gleeson and Manning, 2008). Wilson and Guan (2004) provided a literature review regarding mountain block hydrology and advancements made in quantifying groundwater recharge at the mountain block scale. Their study reviewed methods (water balance, precipitation runoff regression, chloride mass balance, modelling, tracer approaches and others) of quantifying 4 groundwater recharge in alpine watersheds; mostly in the arid and semi-arid regions of the southern United States. The climatic differences between these regions and the alpine region used in the proposed research will likely limit their relevance. This also holds true for more recent studies in other climatic conditions, such as that of Kao et al. (2012) in sub-tropical conditions. Despite the geographic differences in locations of these other studies, the innovative methods used may be beneficial for estimating baseflow. Research conducted by Magruder et al. (2009) is similar to the proposed research with respect to surface process modeling. Magruder et al. (2009) used an ecohydrologic model to calculate water mass balance in a mountain watershed in Montana to estimate groundwater recharge in the mountain block, which they subsequently linked to a groundwater model in MODFLOW. The concept of soil water output is used by Magruder et al. (2009) to describe groundwater recharge thresholds, which are used within a finite-difference solving scheme. This thesis research employs similar methods by using a hydrological model to estimate recharge to the groundwater system; however its primary purposed is to introduce alternative methods of groundwater flow model discretization. 5 Chapter 3 Abstract Novel geomorphometric classification of heterogeneity of unconsolidated sediments for hydrogeologic characterization in mountain watersheds Groundwater controls many aspects of water quantity and quality in mountain watersheds. Yet, groundwater recharge and flow originating in mountain watersheds are difficult to quantify due to challenges in the characterization of local geology, as subsurface data are sparse and difficult to collect. Remote sensing data are more readily available and are beneficial for the characterization of watershed hydrodynamics. We present an automated geomorphometric model to identify the approximate spatial distribution of geomorphic features, considered to represent (sub)surface geology, and to segment each of these features based on relative hydrostratigraphic differences. A digital elevation model (DEM) dataset and predefined indices were used as inputs in a mountain watershed. The model uses periglacial, glacial, fluvial, slope evolution and lacustrine processes to identify regions that are subsequently delineated using morphometric principles. A 10-m cell size DEM from the headwaters of the St. Mary River watershed in Glacier National Park, Montana, was used to extract five morphometric parameters to calibrate the model: elevation, slope, flow direction, flow accumulation, and surface roughness. Algorithms were developed to utilize these parameters to delineate the distributions of bedrock outcrops, periglacial landforms, alluvial channels, fans and outwash plains, glacial depositional features, talus slopes, and other mass wasted material. Theoretical differences in substrate material and hydrofacies associated with each of the geomorphic features were used to segment the watershed into units reflecting similar hydrogeologic properties such as hydraulic conductivity and thickness. The results of the model were verified by comparing the distribution of geomorphic features with published surface geology maps. Although agreement in semantics between datasets caused difficulties, a consensus yielded a comparison Dice Coefficient of 0.65. Our novel method and results can be used to delineate geomorphology 6 and surface geology in remote regions, estimate spatial differences in near-surface groundwater behaviour, and assist in groundwater model calibration. 7 Chapter 3: Novel Geomorphometric Classification of Heterogeneity of Unconsolidated Sediments for Hydrogeologic Characterization in Mountain Watersheds Introduction Advancement in hydrologic modeling is vital for predicting the effects of natural and anthropogenic effects on water resources (UNEP, 2007, Bell et al., 2012, van Leeuwen et al., 2012, Wada et al., 2010). Effective water resource modeling is especially important for mountain watersheds as these provide large quantities of fresh water to users and ecosystems on every continent, either directly from rain runoff, or delayed from melting of snow and ice (Wohl, 2006, Beniston et al., 1997, Marquis et al., 2012). Surface water quality and quantity are greatly dependent on groundwater at local and regional scales, but groundwater behavior is difficult to predict in integrated surface-groundwater hydrologic models (Smerdon et al., 2009), as hydrogeologic data are often absent in mountainous catchments, or are spatially interspersed due to the topographic complexity and remoteness of these areas. As such, new research on groundwater recharge and modeling in mountain catchments is needed, especially for exploring alternatives to data collection and integration (Wilson and Guan, 2004, Kao et al., 2012, Scanlon et al., 2002, Smerdon et al., 2009, Magruder et al., 2009, Manning and Solomon, 2005). Research focused on groundwater behavior at the watershed scale requires ingenuity (Hayashi and van der Kamp, 2009), which includes alternative approaches and data integration procedures in the context of groundwater flow modeling (Motaghian and Mohammadi, 2011), especially when hydrogeologic observations are sparse. Hydrogeologic subsurface interpretation is typically achieved through the use of boreholes (Koltermann and Gorelick, 1996), hydrogeophysics (Hubbard and Linde, 2011), or estimated hydrochemically (Baraer et al., 2009). Limitations in hydrostratigraphic characterization often arise from factors such as the sampling density and distribution, or the availability of hydrogeophysical data, as they may not 8 reflect the scale-dependent needs of many groundwater modeling techniques (McCloskey and Finnemore, 1996, Eaton, 2006). This study examines an alternative approach to hydrostratigraphic characterization using multi-source data and integration methods, by combining remote sensing data with Digital Elevation Models (DEM). Remote sensing data are becoming increasingly available globally with the aid of aircraft and Earth Observation Satellites and the addition of lower cost UAV platforms. These data are valuable for watershed-scale descriptive subsurface characterizations of regions with diverse hydrodynamic properties (Clow et al., 2003). DEMs are now commonly used for the extraction of geomorphic information in the form of geomorphometry (Schmidt and Hewitt, 2004, Napieralski et al., 2007, Pike, 2002, Dymond et al., 1995). Hydrostratigraphic information can be obtained from geomorphic landforms based on principles of sedimentation (de Marsily et al., 1998) and facies-hydrofacies correlations (McCloskey and Finnemore, 1996, Fisher and Barnhill, 1998). Use of morphometric and geomorphic process-based principles, in combination with numerical modeling, has become a widely used method (Bishop et al., 2012) to complete geospatial analyses of DEM data at small scales. Studies similar to that of Saadat et al. (2008) have successfully performed analyses of DEM data to generate geomorphic classifications using factors such as slope, elevations and stream network patterns, which have proven applicable to mountain catchments (Trevisani et al., 2012). In this study, we use surface segmentation to systematically delineate features by observing and classifying data containing mutual attributes, such as flow grids, slope thresholds and trends, surface roughness, etc. (Minár and Evans, 2008), and thus integrate a landform-based approach to coarsely estimate hydrostratigraphy. This is an alternative to complex simulations of micro- scale sedimentary processes or grain size distributions which require larger amounts of data and computation (Koltermann and Gorelick, 1992, Coulthard and Van de Wiel, 2012) and are sometimes only accurate within an order of magnitude (Fisher and Barnhill, 1998, Vienken and Dietrich, 2011). The novel geomorphometric classification we propose is developed for a small 9 mountain watershed bordering the Rocky Mountains, but is designed to be easily applied in complex mountain watersheds using the manipulation of geomorphic feature boundary associations (indices) to calibrate an automated geomorphic characterization. The results of this study will assist in efficient delineation of heterogeneity for groundwater model calibrations. Study Area and Dataset The St. Mary River watershed in Montana, United States, and Alberta, Canada (Figure 1) was chosen as a study area due to its lack of subsurface data. However, coarse surficial geology datasets and borehole logs were available for verification purposes. Previous studies in this area also suggest the need for groundwater modeling to more accurately simulate the hydrodynamics of the watershed (MacDonald et al., 2011, MacDonald et al., 2009). 10 Figure 1 Study Area: Headwaters of the St. Mary River, Montana, USA. Hillshaded and elevation coloured DEM (USGS, 2009) The headwaters of the St. Mary River watershed lie on the eastern flank of the Rocky Mountains in Glacier National Park, Montana. Runoff from the St. Mary River flows to the northeast to its confluence with the Oldman River in southern Alberta. The delineated watershed chosen for this study consists of the headwaters above a stream-gauging station at the International Boundary site (Water Survey of Canada station 05AE027). The total area of the watershed is 1,195 km2, with an elevation range of approximately 3,031 to 1,249 m above sea level (a.s.l.) (MacDonald et al., 2009). 11 The near-surface geology is comprised of consolidated and unconsolidated units ranging in age from Precambrian to Holocene. The thickest of the units is the Precambrian Belt series, of which the Ravalli and Missoula groups occur near the surface in the southwest (upper) portion of the watershed. The remainder of the watershed to the northeast is underlain by undifferentiated Cretaceous shales and sandstones (Ross, 1954, Zientek et al., 2005). Unconsolidated deposits within the watershed are comprised of Quaternary sediments, representative of recent geomorphic processes. Till sheets and moraines present near the current glacier and snowfields date from 10,000 BP and later, and there is an abundance of recent colluvial and alluvial deposits consisting of unvegetated sediments (Carrara, 1987). Methodology 3.3.1. Geomorphometric Analysis Geographic Information Science (GIS) software and its geospatial tools were used to process gridded data and to examine the relationships between geomorphic principles and terrain attributes for model calibration. Gridded data consisted of a no-cost 10 metre (m) cell size DEM (USGS, 2009). These data were considered adequate, as Napieralski and Nalepa (2010) found that DEM spatial resolutions finer than 15 m are significantly more effective in extracting landforms through DEM surface segmentation than at coarser resolutions. Terrain attributes associated with the genetic principles of geomorphic landforms were examined and utilized to calibrate a set of algorithms to simulate the distribution of the landforms. The genetic principles utilized were chosen predominantly based on logic related to morphometry, surface processes, sediment transfer regimes, and hydrology. Chosen terrain attributes (herein referred to as indices) were predefined, optimised, and used as variables within the computationally-derived classification. The basis of the selection of indices are described in further detail in the remainder of this section. Region-growing algorithms (Bishop et al., 2012) utilized to delineate and geomorphic features were written in Python as part of the automated landform delineation. 12 Geographic Information Science (GIS) software and its geospatial tools were used to process gridded data and to examine the relationships between geomorphic principles and terrain attributes for model calibration. Gridded data utilized consisted of a no-cost 10 metre (m) cell size DEM (USGS, 2009) which was downloaded for use. These data were considered adequate, as Napieralski and Nalepa (2010) found that DEM spatial resolutions finer than 15 m are significantly more effective in extracting landforms through DEM surface segmentation than at coarser resolutions. Terrain attributes associated with the genetic principles of geomorphic landforms were examined and utilized to calibrate a set of algorithms to simulate the distribution of the landforms. The genetic principles utilized were chosen predominantly based on logic related to morphometry, surface processes, sediment transfer regimes, and hydrology. Chosen terrain attributes (herein referred to as indices) were predefined, optimised, and used as variables within the computationally-derived classification. The basis of the selection of indices are described in further detail in the remainder of this section. Region-growing algorithms (Bishop et al., 2012) utilized to delineate and geomorphic features were written in Python as part of the automated landform delineation. An evaluation of principles of geomorphologic semantics (Dehn et al., 2001) aided in the selection of criteria to identify relevant landforms. The depositional and erosional regimes in the watershed were characterized based on the gradual changes of terrain attributes derived from the DEM, which is an axiom of surface segmentation described by Minár and Evans (2008). For example, the upper reaches of glaciated valleys are steep enough to be void of significant unconsolidated material as the force of gravity exceeds the friction required to retain unconsolidated sediment on the surface (Maeda and Yuasa, 2009), and coarse-grained material is unlikely to exist on slopes exceeding 40 to 45o (Wichmann et al., 2009). These regions provide detectable boundaries between bedrock outcrop surfaces and colluvium, especially at high elevations where sediment cascades (Becht et al., 2005) dominate over fluvial and biological 13 weathering (Wichmann et al., 2009). These members of the sediment cascade regime have characteristic morphometric attributes that can be used in classification (Wichmann et al., 2009, Marquinez et al., 2003), which were utilized in this study. More specifically, geomorphic regimes were associated using genetic principles that consist of glacial processes, periglacial processes, slope evolution, fluvial processes, and lacustrine depressions (Figure 2). Glacial Processes Alpine Glacial Valleys Bedrock Lacustrine Depressions PeakSlope Outcrops Evolution Tal s Slopes Channelized Regions Fluvial Colluvium Alluvial Fans Processes Outwash Plains Periglacial Outside of Valley Processes Figure 2 Geomorphic landform classification relationships The location and extent of surface landforms were predicted using genetic (process-based) and morphometric analyses, respectively. This sequential approach avoids scale dependence in locating geomorphic features using explicit geometric calculations. Following the evaluation of genetic regimes in the watershed, a morphometric analysis was used to delineate the extent of geomorphic landforms. The final extent of each landform is a result of the optimization of indices used to control rule sets for region-growing algorithms. The geomorphometric characterization of each genetic regime is described in the following subsections. Fluvial Processes The DEM was initially processed to segregate channels and locate regions of fluvial erosion and deposition by computing a surface slope dataset, and applying flow direction and flow accumulation tools in ESRI ArcGIS version 10.1 (Tarboton et al., 1991). Changes in the slopes 14 of stream networks were used as a primary indicator of erosion and deposition intensity- thresholds, as deposition is minor on steep regions of valleys (Benda et al., 2005). A stream slope dataset was generated using elevation differences and slopes of grid cells comprising the channelized regions. The maximum change in stream slope over the four adjacent stream grid cells was subsequently calculated to create a stream slope derivative dataset. This dataset was used to obtain the indices required to link fluvial erosion and alluvial fan/outwash plain deposition thresholds. Visual inspections of publicly accessible satellite imagery (LIST Google Earth/Landsat/etc.) and hillshade representations of the DEM were used to relate the onsets of alluvial fans with the optimal stream slope change indices. Additional indices were optimised for delineating the extent of alluvial fans/outwash plains by examining the surface slope dataset. The first set of indices was used as seed points to initiate the region-growing algorithm, in which the second set of indices were used as variables to calibrate the rule sets that grow the area of the feature (Figure 3). a) Seed Points b) Region Growing c) Final Extent Figure 3 Identification of alluvial fan/outwash regions and their delineation The region-growing algorithm used morphometric factors that include changes in the slope of the alluvial deposit surface (Saito and Oguchi, 2005, Hashimoto et al., 2008), which were mapped in a topographically downward direction over an eight cell neighborhood. The locations of channelized regions were also used to identify remaining fluvial features using a regime of geomorphic process memberships (Figure 4). 15 Channelized Region Channel Slope True False Decreases Onset of Alluvial Fan If Low Perpendicular Slope If Flow Direction within 180o of Topographic Gradient And Low Slope Alluvial Fan Alluvial Outwash Plain Figure 4 Process memberships used in model to identify fluvial features (see also Table 1) Grid cells of channelized regions were also used as seed points for an additional morphometric region-growing algorithm to delineate stepped fluvial profiles and relict channelized regions. This algorithm also used indices derived from the topographic slope dataset to determine the lateral extent of channelized regions. Flat-lying regions surrounding streams were delineated to reflect the areas widened by fluvial processes, and the stepped terraces of mountain headwater streams that are caused by episodic blockages by boulder-sized sediment and organic material (Benda et al., 2005). Slope Evolution and Periglacial Processes The boundaries of landforms within glaciated valleys were determined through characteristic transition zones that reflect variability in genetic mechanisms (Webster et al., 2006, Wichmann 16 et al., 2009). Slope values from within the valley profiles were chosen based on an examination of satellite imagery, the hillshade of the DEM, and the surface slope dataset, to determine appropriate indices to reflect transition zones of slope evolution genetic processes. The indices were optimized by qualitatively comparing the modeled features with satellite imagery, and were then used as inputs in the geomorphic process membership workflow that was used to complete the remainder of the geomorphometric analysis (Figure 5). 17 If Not Channel And High Slope And Planar And Low Slope And High-Moderate Slope And Above Bedrock Outcrop Lacustrine Alpine Elevation Minor Till If Elevation Below Bottom of Bedrock Outcrop Thick Till 2 And Topographic And Flow Direction within 180o of Topographic Gradient 1 And Moderate Slope Talus Slope 1 2 If Elevation Below Bottom of Talus Slope Moderate And Low-Moderate Thickness Till Slope And Flow Direction within 180o of And High Surface PeriglacialTill Topographic Gradient Roughness Colluvium Figure 5 Process memberships used in the model to identify remaining (sub) surface features (see also Figure 4 and Table 1) To identify bedrock outcrop regions, a threshold slope value representing the absence of significant unconsolidated surface material was chosen. Any regions in the DEM exceeding this slope were segmented into bedrock outcrops. The lowest elevations of these bedrock outcrop 18 regions were considered the upper boundary of talus, and were subsequently used as seed points for a talus slope morphometric region-growing algorithm. To reflect the concave nature of talus slope surfaces, downhill decreasing slopes until terminal angle values were used as indices to be applied to the rule sets of the talus region-growing algorithm. Talus features were delineated in a downward topographic direction and computed over an eight grid cell neighborhood. A surface roughness dataset was created to assist in the morphometric analysis of colluvium. Three-dimensional vectors perpendicular to the terrain surface were generated, as described by McKean and Roering (2004). The degree of scattering of the vectors was computed using local standard deviation over a nine cell neighborhood. Surface roughness values indicative of colluvium were identified through an inspection of largely unvegetated alpine regions through satellite imagery (World Imagery from the ESRI online MapServer). The generated distributions of bedrock, talus, and fluvial features were used as indicators of other genetic processes on valley slopes, which were delineated progressively downslope. The final extents of the delineated talus slope regions were used as seed points for the colluvium region-growing algorithm. This algorithm first computed the densities of threshold surface roughness values over an 11-by-11 cell neighborhood and then grew the colluvium regions in a downslope direction. The growing of these regions was terminated at regions of lower surface roughness and/or previously delineated fluvial deposits. Periglacial processes delineated in the classification were limited to scales in the range of freeze- thaw weathering in regions not glacierized during the Holocene. Process and product means of modeling these features was not considered feasible due to the inability to isolate a host of other geomorphic processes (Hall and Thorn, 2011). Regions potentially existing outside of glacial valleys were used to delineate extraglacial regions, given that glaciated valleys have been shown to be the primary control on alpine peak topography (Mitchell and Montgomery, 2006), and that these regions have been subjected to climatic mechanisms related to periglacial processes 19 for extended durations (Whipple et al., 1999). Peak topographic regions were identified by the selection of elevation indices associated with alpine regions. The approximate lowest elevation of alpine regions was identified using satellite imagery to construct estimates of unvegetated elevations, and was used as a threshold for an algorithm to propagate through cells upslope until peak regions were encountered. Topographic peaks were used as seed points for a region- growing algorithm that used changes in slope over a nine cell neighborhood to delineate the flat-lying regions above formerly glacierized terrain, termed ‘glacial valleys’ in this paper. Glacial Landforms and Lacustrine Features Lakes are distributed throughout the study area in alpine regions and lower reaches of valleys. Planar surfaces with an area greater than nine contiguous grid cells were segmented into lacustrine regions, which reflected the areas subjected to lacustrine processes at the point in time at which the DEM imagery was collected. When examining glacial processes, classification methods based on genetic principles associated with the sediment cascade regime were more elusive. The complexity of glacial landforms in mountain environments is exacerbated by peri- and paraglacial geomorphic processes (Slaymaker, 2009; Evans, 2012). Glacierization was once the most powerful genetic regime in mountain watersheds (Mitchell and Montgomery, 2006) but much of the glacial surficial sediment has been covered and redistributed. The retreat of the glaciers left an environment highly susceptible to weathering, slopes failure and fluvial processes. Intact glacial landforms such as moraines exhibit slope and curvature-related attributes that can be evaluated morphometrically (Napieralski et al., 2007), however, their complexity is often overlooked (Dunlop and Clark, 2006) and is subject to scale-related error. This study is primarily relevant to the potential distribution of hydro-facies, and the complexity of discrete glacial landform classification is considered outside its scope due to the lack of fully intact glacial features. A process of elimination was used to overcome the complexity of glacial landforms in the surface segmentation. The glacial surface segmentation was the last to be modeled, and was 20 completed in areas of the DEM not previously segmented into fluvial, mass wasting, periglacial or lacustrine processes. The shapes of the glacial valley profiles and protruding glacial landforms were used as indicators of glacial sediment distributions. Regions with terrain slope values similar to those of talus slopes and colluvium were considered to be draped in a thin layer of till, as these slopes were more likely attributed to relatively thin glacial cover of the lower side walls of U-shaped glacial valleys (Harbor, 1995). Regions with slopes steeper than alluvial sediments, but less steep than other features, were considered bounding regions of glacial landforms, which were classified as medium-thickness till. Regions with shallow slopes characteristic of valley bottom fill were considered thick draped till. 3.3.2. Model Indices and Quantitative Verification Geomorphic interpretations and iterative evaluation of the simulated results yielded a range of optimized indices used in the model, which are presented in Table 1Table 1 Boundary indices chosen for geomorphometric model. Indices Value Alpine minimum elevation. 2100 mASL Periglacial slope range. 0 - 18° Bedrock outcrop slope threshold. 35° Talus slope termination angle. 18° Colluvium surface roughness threshold. 1.85 Alluvial fan onset stream slope change threshold. 2.6 Alluvial fan slope range. 0.2 - 9° Minor glacial sediment slope range. 0.001 - 9° Medium glacial sediment slope range. 9 - 18° Thick glacial sediment slope range. 18 - 35° Table 1 Boundary indices chosen for geomorphometric model The qualitative nature of the geomorphologic interpretation to obtain the indices warranted the comparison of the modeled distribution of geomorphic features with a published dataset. A surficial geology map of Glacier National Park (Carrara, 1990) and a geology map of the 21 Blackfeet Indian Reservation (Cannon, 1996) were the most detailed surface geology maps of the study area available, and had been previously georeferenced (Larsen, 2004). The maps were combined to cover the study region watershed, and were manually digitized and gridded using a 10 m cell size to be used as a comparison dataset. Minor rescaling and reclassification of features in both the modeled dataset, and the surficial geology map was required to create a dataset that reflected differences in scaling and geomorphological semantics (classification terminology and meaning: c.f. Dehn et al., 2001). The modeled dataset was completed using a classification at a 10 m scale, which yielded a greater number of discrete landforms than the published maps. The spatial precision associated with the simulated geomorphic landforms was not inherent in the broad classification, which warranted a spatial generalization of the simulated map. To upscale the modeled dataset for verification, the grid was smoothed using a majority filter over a seven-by-seven cell neighborhood. Uncertainty between names of geomorphic features was also present between the two comparison datasets. One key difference, which included alpine till and colluvium, was likely attributed to the concurrent deposition and highly heterogeneous nature of these features. Several landforms were combined into single categories to create a compatible classification between the datasets (Table 2). 22 Modeled Bedrock Till- Till- Lacustrine Till- Periglacial Talus/ Channel Fan and Mass Classification: Outcrop Thick, Valley coarse Outwash Waste/fine Slope Minor draped Alluvium Colluvium Alluvium colluvium Deposits Published Compatible Classifications: Maps Lacustrine Lacustrine Bedrock Bedrock Diamicton Till- Till- Contine Continental ntal Solifluction Periglacial Frost Rubble Talus Talus/ Colluvium Colluvium Alluvium (Holocene) Alluvium Alluvial Fan Landslide/Till Till (Holocene) Colluvium /Alpine Till- Piedmont Till Rock Glacier Deposit Table 2 Landform reclassification to create comparable datasets for model verification Overlap between datasets was evaluated for the entire watershed and for each feature separately. To evaluate the overlap throughout the watershed, the Dice similarity coefficient (DSC) was used. The DSC is a metric that is used to statistically quantify the spatial similarity of two datasets. The DSC measures the presence or absence of spatially overlapping values between two datasets and ranges from 0 (no similarity) to 1 (complete similarity). The DSC is computed by 2𝐴∩𝐵 DSC(A, B) = , (1) 𝐴+𝐵 23 which is the ratio of the sum of the spatial units where dataset A (simulated) and dataset B (comparison) contain the same classification, to the sum of the total number of spatial units in both datasets (A and B) (Zou et al., 2004). 3.3.3. Hydrogeological Characterization Hydraulic conductivity values were obtained from literature for all classified units to be applicable in models of shallow unconfined groundwater flow (Batu, 1998). Several simplifications were made to accommodate the extensive uncertainty in such hydrogeological characterization. Bedrock was considered to have an overall lower hydraulic conductivity than the overlying unconsolidated sediments. The assignment of secondary permeability and conductivity values (e.g., due to fractures) falls outside the scope of this study, but could be accounted for by other means when conceptualizing a groundwater model. Remaining geomorphic landforms delineated in the geomorphometric analysis were further segmented into sub-units that reflect hydrodynamic differences. These segmentations were characterised by theoretical distributions of unconsolidated sedimentation thickness and hydrostratigraphy. Alluvium Using their total size and slope, alluvial fans and outwash plains were divided into three segments that delineate relatively proximal, medial and distal sections. These section sequences commonly exhibit decreases in hydraulic conductivity and fan thickness (McCloskey and Finnemore, 1996). The approximate triangular area downslope of each grid cell of the fan was first computed to segment the alluvial fans. Secondly, each fan was divided into proximal, medial and distal sections based on decreasing triangular area. The proximal section of larger alluvial fans in the study area (areas found to exceed 0.5 km2) was assigned more area to reflect the higher theoretical volumes of sediment transport at higher angles (Roering et al., 2001). 24 Talus Slopes and Colluvium Hydrogeologic differences were assigned to talus slopes (coarse colluvium) and other mass- wasted material (finer colluvium) using the slope dataset. Talus slopes were divided into two sections based on slope values that reflect positions on the talus slope. Upper steep sections of talus slopes typically contain smaller grain sizes that are more poorly sorted, while the lower regions are coarser and likely transmit higher quantities of groundwater (Clow et al., 2003). As such, steeper regions near the apices of talus slopes were assigned lower relative hydraulic conductivity, while the lower toe regions were assigned higher hydraulic conductivity. Regions of modeled colluvium were assigned hydrogeologic properties in a similar fashion as talus slopes. Less steep regions likely exhibit more intermixing between slope evolution and glacial or fluvial processes, and were assigned lower relative hydraulic conductivity with respect to sediments originating from slope failure on steeper slopes. Till Regions characterized as till were not segmented further to adhere to the respective limitations on hydrostratigraphic scaling. Holocene till found on the steeper slopes of high relief glaciated valleys (alpine till) was assigned a higher relative hydraulic conductivity that is similar to colluvium. Draped till on flat lying regions was assigned an overall lower hydraulic conductivity to reflect sub-glacial till deposits or continental glacier deposits outside of mountainous terrain. The heterogeneous properties of till warrant a larger range of hydraulic conductivity values be assigned to these features. Lacustrine Regions beneath existing lakes were assigned a relative hydraulic conductivity that is dependent on whether they are present above or below the alpine elevation to reflect annual hydrodynamic variations. At higher elevations in glacial valleys where zones of geomorphic processes produce predominantly sediment cascades (Wichmann et al., 2009), the lacustrine regions were assigned 25 higher hydraulic conductivities relative to lower lacustrine regions. Regions predominantly exhibiting lacustrine characteristics in the bottoms of valleys were assigned a lower (the lowest) hydraulic conductivity relative to all other geomorphic landforms. Results 3.4.1. Geomorphometric Model Results The geomorphological segmentation classification of the DEM resulted in the division of the watershed into 10 geomorphic classes (Figure 6). The quantification of overlapping (coincident) landform classifications and percentages of areas contained within are included in Table 3 and Table 4. 26 Figure 6 Map of Geomorphometric Model Results: the overall representation of major geomorphometric features complies with the logical progression of the sediment cascade in glaciated terrain 27 Percentage of area of simulated Most commonly Number of Number of landform misclassified grid cells for coincident within landform within each unit in simulated grid published map published map Surface Geology published map cells landform landform Lacustrine 404,920 360,904 55 Alluvium Talus/Colluvium, Bedrock 2,557,757 1,270,963 38 undifferentiated Talus/Colluvium, Colluvium/Piedmont 1,363,798 794,677 60 undifferentiated Till Colluvium/Piedmont 5,109,518 2,382,200 40 Till- Continental Till Colluvium/Piedmont Till- Continental 354,630 183,285 25 Till Colluvium/Piedmont Alluvium 359,521 277,046 56 Till Talus/Colluvium, Periglacial 218,277 78,009 35 undifferentiated Table 3 Results of the comparison between simulated and published geomorphology datasets 28 Total area Total area percentage of percentage of simulated published Surface Geology landform landform Lacustrine 4 4 Bedrock 15 25 Talus/Colluvium, 22 13 undifferentiated Colluvium/Piedmont 29 49 Till Till- Continental 13 3 Alluvium 14 3 Periglacial 2 2 Table 4 Results of the comparison between simulated and published geomorphology datasets- total area A spatial overlap comparison between the modelled and reference datasets yielded a DSC of 0.63 over the entire watershed. Distributions of lacustrine, bedrock, and colluvium features in the modeled dataset were in relative agreement, with DSC values of 0.85, 0.62 and 0.58, respectively, while a DSC of 0.43 indicated less spatial overlap between talus/colluvium. Till, alluvium, and periglacial features did not overlap well, with DSC values of 0.21, 0.30 and 0.37, respectively. 3.4.2. Hydrogeological Characterization Results The spatial distribution of relative heterogeneity is dependent on the chosen relative facies distributions, which were deductively chosen in this study based on theoretical concepts. Figure 6 is an example of units likely to exhibit bulk hydrostratigraphic differences based on several boundary conditions, as described in the methodology. These conditions are easily modifiable within the model, and may be assigned indices that progressively improve the performance of groundwater models. In effect, the results of this dataset are not an end-member of the 29 hydrogeological characterization, rather an example of an applicable option. Theoretical hydraulic conductivity ranges from approximately 0.1 x 10-8 to 7.1 x 10-4 cm/s for the lower values to 6.6 x 10-1 to 3.1 cm/s for the higher range (Batu, 1998). 30 Figure 7 Simulated distribution of hydrogeological facies. Higher hydraulic conductivity is generally noted at higher elevations, with an attenuation toward the bottom of valleys and the lower reaches due to lacustrine and glacial sediment dominance. 31 Discussion 3.5.1. Geomorphometric Modeling The calibration indices and automatical delineation of geomorphic landforms proposed in this study originated through hydrologically significant genetic and morphometric principles. The geomorphometric analysis is completed through the perception of hydrological model calibration, which includes geomorphic semantics that commonly relate to hydrodynamics (Dehn et al., 2001). The indices established in this study are selected through large-scale fluid erosion and deposition. This scale does not represent every explicit geomorphologic process, it is an a priori initial characterization for large, complex watersheds(Schmidt and Hewitt, 2004). The most salient example of geomorphic processes described in the analysis that relate to hydrologic principles is features of fluvial origin, where locations of current stream systems are assumed to reflect previous fluvial erosion and deposition (Miliaresis and Argialas, 2000). Due to the high geomorphic variability of mountain terrain, a fluvial modeling approach independent of stream order and elevation was required. Deterministic modeling of geomorphic sedimentation processes have proven to be widely scale-dependent (Slaymaker, 2006), and require extensive knowledge on geomorphic properties to apply available physical theory (Dietrich et al., 2003). The simplified method of predicting sedimentation changes based on erosion presented here is less dependent on scaling, and provides a somewhat stochastic method of representing otherwise complex processes (Chase, 1992). Surface roughness is a useful indicator of debris flows, which often exhibit a higher degree of surface roughness with respect to intact landmass, and is used to delineate of mass failures in the lower reaches of valleys (McKean and Roering, 2004). With the downslope progression of the sediment cascade, the boundaries of features become increasingly enigmatic as the slope evolution processes become highly heterogeneous (Wichmann et al., 2009). Talus slopes, which are significantly less steep than bedrock outcrop source regions (Becht et al., 2005), exhibit a 32 higher slope than colluvium further down the sediment cascade regime. These boundaries require increased complexity to delineate, which will require additional morphometric tools. 3.5.2. Hydrogeologic Characterization Ideally, information on vertical anisotropy and hydrogeologic facies thicknesses would be useful for subsurface characterization. These parameters require a detailed hydrostratigraphic characterization, which is beyond the scope of this study. Instead, we characterize hydrogeology in a way that satisfies the scalability of the information from the input dataset (a DEM). The result is suited to assist in the design and calibration of a groundwater model, and potentially reflect variability in net recharge, and groundwater flow beneath the soil interface. These results may also be directly suitable to a stochastic modeling approach (Freeze, 1975). The distribution of heterogeneity was modeled based on the following key principles associated with each landform. Dividing landforms in the upper reaches of the watershed into respective units is crucial for groundwater flow, as landforms have shown to host potentially multiple groundwater flow regimes (Roy and Hayashi, 2009). The modeled distribution of heterogeneity demonstrates that groundwater is likely to recharge and flow quickly, and have a high probability of discharging above lower conductivity sediments. The bulk hydrostratigraphic units delineated within mass- wasted material are meant to reflect large representative elemental volumes, due to the variability in depositional processes and positions on glacial valley slopes (Roy and Hayashi, 2009). The modeled hydrogeologic characterization reflects the tendency for groundwater to rapidly drain within talus due to its coarse-grained sediment and high topographic gradients (Caballero et al., 2002). Groundwater information at a finer scale is inherently more useful further down the sediment cascade, as heterogeneity likely increases. Colluvium in mountainous watersheds is often paraglacial (related to time after glacial retreat), leaving highly variable sediment distributions that reflect both glacial and mass wasting processes, with subsequent interstitial ice playing a 33 role in their distribution (Millar et al., Sass, 2006). Colluvium and glacial depositional features host multiple groundwater flow regimes (McClymont et al., 2011), which will be difficult to differentiate due to similar hydrostratigraphic properties. As such, the hydrogeologic characterization within colluvium and alpine till is broad and requires a finer classification scale that is primarily represented by surface slope in this study. Within till features, the genetic origins are used as the primary indicator of heterogeneity. Alpine till exhibits hydraulic conductivities several orders of magnitude greater than continental till (Ronayne et al., 2012), which is found outside of high relief glaciated valleys. Despite these differences, both types of till are highly heterogeneous and commonly contain secondary hydraulic conductivities (Cuthbert et al., 2010). The heterogeneity depicted by the model represents bulk units of till that should be at an appropriate scale to encompass inter-unit differences in hydraulic conductivity and anisotropy. Local groundwater flow is likely to converge and have longer residence times in fluvial sediments in the lower reaches of the basin, while the majority of subsurface drainage will take place in channelized alluvium. The degree of flow attenuation is reflected in the geomorphometric classification through alluvial fan and outwash subsurface classification, as well as classification of channelized regions. The size of alluvial fans has an effect on the overall transmissivity, as the volume of sediment shares a 1:1 proportional power law relationship with the area (Giles, 2010). As such, calculations of transmissivity must take into account the position on the alluvial fan. Groundwater flow through alluvial fans is likely to be variable through the proximal, medial and distal regions due to the strong correlation between depositional facies (i.e., coarse well-sorted grain-supported channel deposits versus indurated poorly-sorted finer grained overbank deposits) and the distribution of hydrogeologic facies (Fogg et al., 1998). Using elevation as a proxy for lake bed sediment hydraulic conductivity provides a large range of classifications that can be easily changed to adjust for the rate of snowmelt related recharge 34 and the complex groundwater recharge mechanisms of these lakes. Gurrieri and Furniss (2004) have shown that the inclusion of sub-annual transient water exchange is necessary to accurately describe the hydrodynamic nature of alpine lakes. Several lakes studied in Montana exhibited complete water mass balance replacement throughout the course of a year, and lake water levels during freshet are largely buffered by snow melt inflow and hydraulic connections with shallow groundwater. In drier periods, many of these alpine lakes store water and slowly recharge local groundwater systems. Accurate estimates of lake bed vertical transmissivity are crucial for the calculation of water mass balance downgradient. Conversely, Arp et al. (2006) found that large flow-through lakes in lower reaches of the basin may not regularly attenuate basin peak flows during snow melt, but instead play a significant role in discharge attenuation during summer storm events. The implications of this phenomenon for hydrologic simulations are twofold. First, the depression storage of lakes will play a larger role in basin hydrodynamics than groundwater exchange after lake levels have decreased in the summer. Second, with increased lake levels and more exposure to lake shore sediments, the likelihood of groundwater exchange between the lake and surrounding aquifers is increased (Genereux and Bandopadhyay, 2001). As a result, flow-through lakes are better suited to low hydraulic conductivities to better suit large temporal scales in hydrologic simulations, and will be more suited to direct flow through contribution to stream flow while modeling. The hydrogeologic classification of these lakes is still important for the thermal and geochemical properties of these lentic environments, given that their high residence times during drier periods of the year allow for significant groundwater and surface water mixing (Arp et al., 2006). 3.5.3. Assumptions and Limitations This hydrostratigraphic classification is limited to a heterogeneous extrapolation representing the approximate bulk heterogeneity in mountain watersheds. Limitations related to spatial scaling introduce complexity in the application of the results to the calibration of a hydrologic model. Small and linear-scale glacial and fluvial features contain numerous geomorphic 35 attributes not included in the scale of this study. The delineation of streams and channels would typically require resolution finer than 3 m (Benda et al., 2005), and in the order of 1 m to successfully detect changes in small channels in upper reaches of mountain basins (Dietrich et al., 1993). Moreover, some features were in effect misclassified: roche moutonnées were classified as unconsolidated glacial sediment when in effect are bedrock outcrops with low slope. In this particular study area the misclassification was not detrimental; however, in other regions a significant limitation may arise. Other features, such as alluvial fans and alluvial outwash plains, were delineated using the same geomorphometic rule sets. These scale issues neglect sedimentary deposition regimes that are important for groundwater residence times and surface-groundwater exchange. Additionally, temporal effects, such as surface processes that occur after the collection of the DEM, were not included in the classification, but will have implications for long-term hydrologic modeling. Given all segmented geomorphic landforms in the resultant dataset are two-dimensional, any vertical distribution within the subsurface would have to be extrapolated. Overall thickness can only be inferred from the shape of the valley filled. Observed facies in the unconsolidated sediments may also interfinger, superposing units of varying depths. Units will furthermore be hydrodynamically dependent on antecedent conditions and interstitial ice. The applicability of this dataset is therefore conditional to a reasonable method of establishing these boundary conditions. Factors related to subjective geomorphic interpretations and semantics carry the possibility of creating disagreements on the derivation of geomorphometric and surface geology datasets. These factors are no different than typical variability in geomorphologic classifications using indeterminant boundaries (Dehn et al., 2001); however they should be addressed in the context of this study. The boundary transition indices proposed here reflect the qualitative flexibility in model calibration. Finally, certain geomorphic features were not included in our segmentation, e.g., aeolian deposits or anthropogenic alterations to the subsurface. These features may have 36 an unforeseen effect on the subsurface hydrodynamics, which will limit the effectiveness of subsequent simulations. These features were not found to be of particular abundance in the study area, which did not warrant further analysis in this regard. Conclusion and Recommendations The geomorphometric classification and resultant surface geology dataset generated in this study provide estimations of the distribution of hydrogeologic heterogeneity in an alpine periglacial environment using methods that are dynamic, easily adjustable, and suited as a tool to assist with preliminary groundwater model calibration or stochastic groundwater modeling. Not only can the relative distribution of hydraulic conductivity be assessed, indices used in both the geomorphometric and the hydrogeologic facies models can be modified to alter boundary relationships and progressively aid in groundwater model calibration. When borehole data are absent, or their distribution is sparse, as if often the case in mountain watersheds, this will be a useful alternative in gaining an understanding of the hydrostratigraphy. The verification of the modeled distribution of heterogeneity would benefit from its application to a groundwater model calibration. This can be achieved through future groundwater modeling research, which utilizes the methods or results from this study. Using the water mass balance inputs from a hydrometeorological model, a spatially distributed model that routes not only surface runoff but also groundwater recharge and flow, could be generated and verified using the results from this study as input data. The modeling perspectives used in this study were designed to create a platform for novel and useful geomorphometric techniques to assist in the classification of the subsurface. Geomorphometric modeling methods are diverse and are subject to continual development and acceptance (Bishop et al., 2012). It is expected that the logic used in this study provides a platform for improved methods of geomorphometry that reflect hydrostratigraphic factors in even more detail. 37 Chapter 4 Abstract Groundwater Model Discretization and Computation Using Geomorphic-Derived Data for Estimating Unconfined Topographically Driven Groundwater Flow in Mountain Watersheds The discretization of domains for groundwater modelling strongly dictates the performance of simulations and requires a large amount of effort and data to implement. Although this effort can often not be avoided, large-scale modelling in complex and remote terrain could benefit from a coarse-scale first-order estimation. Modelling of shallow groundwater flow at the watershed scale in mountains can be partially addressed by discretization using geomorphic- scale landforms. A method of discretization is presented that utilizes previously generated geomorphic data and topography to automatically assign geometric attributes required to model groundwater flow. Regions treated as equivalent porous media (EMUs) were combined with topographic flow directions to form irregular-shaped spatial units for model discretization. The study employed methods to simulate saturated conditions recharged by outputs of a hydrologic model at the biogenic soil interface. The discretization of the study watershed in two dimensions permitted the Darcy equation to be solved in one dimension for transient unconfined groundwater flow and increased computational efficiency. Groundwater flow throughout the watershed was solved explicitly, and water in excess of the storage capacity was applied to baseflow. The discretization method was tested by applying hypothetical volumes of recharge to the watershed and evaluating the performance, sensitivity of calibration parameters, and explicit solving techniques. The model performed well for the computation of flow throughout the watershed, which required additional iterations to adhere to the von Neumann criterion while computing discharge between EMUs. The model was found to be well suited for the evaluation of basin-wide groundwater discharge and determining the spatial distribution of discharge to baseflow. This research is well suited to be continued with its application with an 38 appropriate hydrologic model to evaluate its potential benefit in the enhancement of the shape of simulated hydrographs. 39 Chapter 4: Groundwater Model Discretization and Computation Using Geomorphic-Derived Data for Estimating Unconfined Topographically Driven Groundwater Flow in Mountain Watersheds Introduction Upstream watersheds, commonly comprised of mountainous terrain, are susceptible to altered hydrodynamics as a result of anthropogenic alterations (Wohl, 2006) and climate change (Green et al., 2011). Much of the available water downstream of these headwaters originates from groundwater (Viviroli and Weingartner, 1999), which is highly susceptible to seasonal climatic variation (Wei et al., 2009). Understanding the timing and amplitude of peak flows, and the contribution to groundwater resources from these watersheds has become increasingly important, especially when considering responses to climate change. To increase our understanding of mountainous watersheds, research and modeling techniques targeting the simulation of surface-groundwater interaction have advanced rapidly in recent years, and are becoming important and useful tools in catchment hydrology (Wilson and Guan, 2004, Manning and Solomon, 2005, Magruder et al., 2009, Gleeson and Manning, 2008, Kao et al., 2012, Ajami et al., 2011). The effectiveness of these studies is often predicated on the availability of data, which require time and resources that are not always available. This issue is especially pervasive in high-mountain catchments that are difficult to access. As such, robust groundwater flow simulations are preferentially used within regions that include a wealth of hydrogeologic information. Simulation of groundwater routing in upstream catchments could benefit from coarser-scale predictive approaches based on sensitivity analyses, especially when research is faced with resource and time limitations (Johnson, 2007). Simulation of groundwater at coarser scales with minimal supporting empirical data require unique discretization, which may be difficult in topographically complex mountain watersheds. These systems are often treated as sequences of 40 bedrock-bounded valleys that exhibit groundwater recharge and flow behaviour respective to the presence of mountain blocks and alluvial basins (Wilson and Guan, 2004). Inclusion of enhanced geometric complexity in basin discretization may be more useful while considering the spatial variability of near-surface groundwater flow in mountain watersheds. The technique described in this study reflects the geologic setting in greater detail and assigns hydrogeologic principles to a computational discretization routine. This approach allows for a greater degree of discretization, which will be useful for identifying regions within the model domain that may require additional discretization and further modelling with greater detail or alternative methods. Geomorphic maps that contain sufficient detail to provide spatial representation of local near- surface geology are a useful source for the estimation of spatial heterogeneity, at least contributing to the delineation of equivalent porous media (Eaton, 2006). These data can also be used to estimate theoretical distributions of hydrogeologic parameters (de Marsily et al., 1998, de Marsily et al., 2005). At this scale, the prediction of groundwater flow is limited to single variables for entire landforms, which are treated as representative elemental volumes. The calculation method chosen to simulate groundwater flow at this scale is limited to methods that represent topographically driven, unconfined groundwater flow, as the presence of nested and deeper confined flow systems cannot be explicitly determined from geomorphic data. The hydrological implications of groundwater simulation at this scale are nevertheless important given the notion that unconsolidated sediments typically have a profound effect on the timing and amplitude of river flow and govern deeper percolation and recharge of deeper groundwater systems (Wittenberg, 1999). A method of automatic discretization using the delineation of geomorphic features and the optimization of dependent variables (i.e., hydrostratigraphic properties and initial conditions) is presented in this study to estimate basin-scale groundwater flow. A constitutive relationship between groundwater storage and discharge in the basin is established using a relationship 41 between basin shape and storage volume. This approach, as opposed to a reverse storage- discharge analysis (Kirchner, 2009), allows spatial representation of groundwater storage. A study watershed is used to apply methods to automatically discretize the area for groundwater flow using the principles of equivalent porous media. These units are delineated using geomorphic data and are further discretized using topographic flow direction. Geomorphic data can be inferred using digital methods (McCleary et al., 2011, Cairns et al., 2014 Submitted, Vannametee et al., 2014), which is opportune for the implementation of this approach. Results from a geomorphometric analysis are used to discretize the basin and provide estimations of distributions of equivalent porous media adequate for the scale of basin-wide topographically driven groundwater flow. This two-dimensional approach to discretization is employed so a one-dimensional explicit application of the Darcy flow equation can be used to compute saturated unconfined groundwater flow. Recharge is applied using estimates of soil water excess from a physically based spatially distributed hydrologic model. 4.1 Study Area The reader is directed to Section 3.2 for a description of the study area used in this research. Methodology This study employed a watershed discretization method that utilizes output data from an appropriate hydrometeorologic model and a recently developed geomorphometric classification method, which characterizes the heterogeneity of subsurface geology. See Chapter 3; (Cairns et al., 2014 Submitted). Results of the hydrometeorologic model provided spatially distributed excess water at the biogenic soil interface to be used as an input from hydrologic response units (HRUs). Results of the subsurface classification method yielded fine-scale spatial distribution of heterogeneity used in watershed discretization, to which topographic data were also appended. Geomorphic data used as inputs in this study represent solely the saturated zone hosted by parent material deposits and do not represent the dynamics of the upper soil domains. The 42 hydrometeorologic model used to provide inputs must have the capability to spatially compute soil water excess within the upper soil domain. Discretization was achieved computationally using gridded surfaces of each respective input layer, which resulted in a series of unstructured units used for computation. Each of these units was delineated using geometric boundary association algorithms to apply two-dimensional geometric attributes (described in further detail below), which allowed for additional simplicity for groundwater flow calculations. The subsequent sections sequentially outline the study area and methods required to achieve the watershed discretization and computation, which were implemented in a Python programming environment using the Numerical Python (Numpy) library. 4.2.1. Spatial Discretization of Domain Several maps depicting local surficial geology were available for the study area, but their spatial representation of landforms did not possess the level of detail required for this study. Rather, a geomorphometric analysis previously completed in the study area (See Chapter 3, and Cairns et al. (2014 Submitted)) was utilized. This method of geomorphic classification was useful in the context of groundwater model discretization, as it incorporates methods and calibration techniques that align with the delineation of heterogeneity of unconsolidated sediments used in the geomorphometric analysis. These sub-units are referred to here as equivalent porous media units (EMUs). EMUs generated through the geomorphometric classification were further refined using terrain flow direction to apply topographic gradient forcing for groundwater flow. This was completed by gridding both the EMU datasets and using a 10 m resolution digital elevation model. Flow direction was computed using an eight direction flow model (Jenson and Domingue, 1988), in which the flow direction consists of the direction of the highest gradient between each grid cell and its eight adjacent cells. Each EMU was subsequently divided into sub-units reflecting variability flow direction. EMUs were assigned unique identifiers throughout the domain using 43 an algorithm to recursively delineate contiguous grid cells with the same properties, which provided the basis for further discretization. Each unique EMU was instantiated into an object and stored in a dictionary to be programmatically accessed. Additional algorithms were developed to compute spatial relationships between each object (EMU) to be used in the calculation of finite-difference groundwater flow. These algorithms included the following functionality:  Connectivity between each EMU was established by calculating the number of grid cells adjacent to downslope EMUs using the flow direction calculated using the DEM. Connected EMUs were stored within the attributes of each EMU object.  Additional attributes related to the geometric connection between each EMU were also added to each object. Using the number of connected cells, the width of EMU connections were calculated. To determine the contributing cross-sectional width between each downslope EMU, the distance to the centre of the contributing EMU was calculated by iteratively counting grid cells in the opposite direction of topographic flow. The contributing area was also calculated for each connecting EMU using the computed width and length upslope.  Placeholders for attributes likely to change in model simulations were created in each unit, such as the effective porosity, hydraulic conductivity and depth calibration parameters, which are described in further detail below. The final addition to the discretization of the model domain included a coupling routine to variably apply recharge from a spatially-distributed hydrologic model. A routine to calculate coincident HRUs with each EMU and the percentage of the total HRU area was developed and added to each EMU object in the program. The HRU attributes in each object provided the means to add recharge due to soil water excess for the groundwater mass balance calculation. 44 Figure 8 Conceptualization of methods used to assign geometric attributes to each EMU object (separated by colour in figure) for computation. Calculated widths between EMUs are denoted by red lines. 4.2.2. Groundwater Flow Simulation and Calibration Saturated flow throughout the watershed was quantified using volumetric storage and discharge potential between EMUs based on connecting geometry, similar to finite volume-based unstructured grid discretization (Panday et al., 2013). Storage was retained within each EMU using absolute volumes, whereby secondary saturated zone calculations were computed using the parameters associated with both storage and the aquifer properties of the EMU. Recharge was integrated into the model through the application of the spatial connection between each EMU and coincident HRUs using outputs generated at the biogenic soil interface from a hydrologic model. Discharge potential (Qt) between EMUs was computed using the Darcy flux (equation 2) in a downward flow direction between each connecting EMUs. 𝜕ℎ 𝑄𝑡 = −𝐾 ∙ 𝐴 ∙ 𝑡 𝜕𝑥 (2) 45 Where 𝑄𝑡 is the discharge potential over a given time step (𝑡), 𝐾 is the hydraulic conductivity between two EMUs calculated using their geometric mean, 𝐴 is the cross-sectional saturated 𝜕ℎ flow area, and is the hydraulic head gradient over the distance between centroids of each 𝜕𝑥 EMU. Of these variables, the hydraulic conductivity and distance between EMUs were the only to be assigned to each EMU prior to calculation. The remainder of the saturated flow variables were computed at each time step using the methods described in the following sections. Cross-Sectional Flow Area The geometric connection between EMUs was calculated in advance during domain discretization, which provided widths to be used to calculate the cross-sectional area of the saturated zone between EMUs. Saturated thicknesses utilized to estimate the cross-sectional flow area were calculated at each time step given the piezometric head calculated in each stress period. The ability to integrate cross-sectional flow area was beneficial for wetting and drying EMUs in simulations. Depths of EMUs were not available using the chosen geomorphometric method to identify heterogeneity. As such, depth to bedrock was estimated using alternative means. The method employed included additional geomorphometric modelling using the relationship between terrain attributes to determine the appropriate bedrock elevation. A relationship between the shape of the glacial valley and depth to bedrock was established using geometric principles identified in several studies. The use of this method stems from the morphology of valleys created by the advance of alpine glaciers, which erodes bedrock in cross- sectional forms that can be estimated by the function 𝑦 = 𝑎𝑥𝑏 (Svensson, 1959), where y is the elevation, x is the lateral distance from the valley bottom, and a and b are constants. Further studies, such as that by Hirano and Aniya (1988), have shown that the shape may exhibit increased similarity to a catenary, which was found to be approximated by the parabolic function 𝑦 = 𝑎𝑥2. These valley shape functions was useful to estimate the depth of sediment that subsequently infilled the valley by subtraction: 46 𝐵 = (𝑦𝑠 + 𝑧) − (𝑦𝑏 + 𝑧) (3) Where 𝐵 is the thickness above bedrock, 𝑦𝑠 + 𝑧 is the surface elevation above a datum, and 𝑦𝑏 + 𝑧 is the bedrock elevation above a datum. Substituting the parabolic catenary estimation described above for 𝑦𝑠 and 𝑦𝑏 yields: 𝐵 = (𝑎 2𝑠𝑥 + 𝑧) − (𝑎 𝑥 2 𝑏 + 𝑧 − 𝑚) (4) Where 𝑚 is a calibrated value of maximum depth of sediment above bedrock at the bottom of the valley. The value of 𝑚 denotes the origin of the catenary that approximates the bedrock surface below the land surface. Equation 4 reduces to: 𝐵 = 𝑥2(𝑎𝑠 − 𝑎𝑏) + 𝑚 (5) Which was the equation utilized to estimate the depth of sediments above an assigned maximum depth. The value of 𝑚 can be optimized using iterative methods such as Monte Carlo to evaluate the sensitivity of the basin to increased or decreased storage. The values of 𝑎 and 𝑏 were also independently optimized using the maximum distance to the valley bottom and elevations obtained from the DEM. Depth to bedrock was generated throughout the watershed by recursively calculating the distance from the valley bottom to each grid cell. This distance was then used in equation 5 to estimate the depth to bedrock throughout the domain. Hydraulic Head Gradient Specific yield from each EMU was estimated using assigned values of porosity, which were reduced or treated as effective porosity in the model due to the unconfined conditions. The piezometric elevation head (ℎ) was calculated using estimates of effective porosity at each EMU and the volume of water storage using equation 6. (2 − 𝑛) ∙ 𝑉 ℎ = (𝑧 − 𝑑) + 𝐴 (6) 47 Where 𝑧 is the elevation of the EMU, 𝑑 is the depth, 𝑛 is the effective porosity, 𝑉 is the total volume of the EMU and 𝐴 is the area. The gradient in piezometric head was recalculated for each day (or any time step) given the changes in groundwater storage at each EMU. Time Step The model was designed to be solved on a specified time step, but this turned out to be problematic, causing numerical instability, a common computation issue inherent to explicit methods. Numeric instability often arises solving boundary problems explicitly, as the time required to apply discharge may not suit the piezometric gradient and transmissivity conditions. As such, the von Neumann criterion was evaluated with respect to the discharge potential between EMUs to determine the appropriate time step (i.e., for numeric stability of linear partial differential equations). The von Neumann criterion was calculated using equation 7, which was iteratively solved to a tolerance of less than 0.5 by continually dividing the time step by a counter (Crank and Nicolson, 1996). 𝑇𝜕𝑡 2 ≤ 0.5, 𝑛𝜕𝑥 (7) where 𝑇 is the transmissivity, 𝑡 is the time step, 𝑛 is the effective porosity and 𝑥 is the distance between EMU centroids. Iterations to determine the appropriate time step were stopped once the von Neumann criterion was satisfied for the entire watershed and the groundwater flow equation was solved for the given time step until equilibrium. Boundary Conditions Recharge added to the shallow groundwater system was applied to saturated flow conditions given the availability of storage volume in each EMU. Throughout the iterative process of groundwater flow calculation, all EMU piezometric elevations that exceeded the land surface elevation were converted to volumes of surface water. Such ‘upward’ discharged groundwater was routed using two different methods to determine the effect on basin-wide runoff. The first 48 method included reintroducing the water into downslope EMUs as additional recharge. In a second method, the discharged water was not reintroduced into the groundwater system, and was assumed to be directly applied to fluvial baseflow exiting the basin in the main river system. Bedrock was treated as an impermeable boundary and recharge to outcropping regions was automatically applied to the surface water routine. Near-surface bedrock at the study area is comprised of Precambrian Belt series, of which the Ravalli and Missoula groups occur near the surface at the southwest (upper) portion of the watershed. The Ravalli group is the oldest unit present, and consists of the Grinnel argillite, Appekunny argillite, and the Altyn limestone. The argillites within this unit contain abundant sand and quartzite, with minor conglomerate. Feldspar and carbonates also occur within the Grinnel Argillite. Precambrian limestone beneath the upper portion of the watershed is largely magnesian, and is suspected of originating in a shallow marine environment due to the abundance of stromatolites. The older limestone is that of the Altyn, while the Siyeh is younger, belonging to the Piegan group. The remainder of the watershed to the northeast is underlain by undifferentiated Cretaceous shales and sandstones (Ross, 1954, Zientek et al., 2005). Aside from secondary permeability in the sedimentary rock, the assumption that bedrock is impermeable is considered adequate. Groundwater discharge at EMUs flowing outside of the watershed (outlet locations) were solved using fictitious piezometric head gradients equal to the contributing EMUs. Output from the model included the discharge of groundwater in the form of baseflow (method 2) and the exiting groundwater discharge. Calibration Sensitivity of groundwater flow simulations is strongly related to hydraulic conductivity and the distribution of representative units of heterogeneity (Sanchez-Vila et al., 2006), which is a source of uncertainty during calibration. Ranges of hydraulic conductivity that coincide with the landform of each EMU were included in the model to be evaluated further by applying additional Monte Carlo optimization techniques. The performance of the groundwater flow 49 simulation could not be appropriately addressed without iterative application of computationally optimized variables for sensitivity analyses and assistance with calibration. Initial conditions of groundwater in the watershed were also a source of uncertainty that required additional optimization. These conditions were addressed by solving the water balance over the same series of stress periods until relative stabilization was achieved. Initial volume conditions were assigned by saturating the shallow aquifer entirely, and monitoring the attenuation of pressure over daily stress periods. The application of storage was also utilized to evaluate the behaviour of the model as a result of the discretization method presented herein. The sensitivity of the basin to preloaded pressure conditions and the numeric stability were checked by the creation of water table elevation maps and maps of discharge at various stages during neutralization of pressure conditions. The application of recharge was tested against simulated hydrologic data generated using the GENESYS model from the St. Mary watershed (MacDonald et al., 2011). Although those data were primarily utilized to evaluate the snow pack, they provided an early estimate for the application of recharge for this study. Model runs included daily recharge applied using the outputs from the GENESYS model, although they were not sufficient to generate a hydrograph due to moderate incompatibility between studies. GENESYS was utilized to evaluate the effects of climate change on the snowpack and did not include additional soil water outflow verification. Results The application of the automatic discretization of the domain converged, which was confirmed by manually calculating the same geometric attributes from several EMUs zooming into a cell- scale extent (see Figure 8 for a visual representation). This was also confirmed by model convergence during the solving of the groundwater flow equation. A total of 699,966 EMUs were generated in the discretization and utilized in model computation. The estimation of depth to bedrock was completed using various maximum depths to confirm the applicability of the method, however only one result is presented in Figure 9. The shape of the depth the bedrock 50 terrain is strongly influenced by the linearity of the valley bottom, which dictated the cross- sectional distance used in calculations. Figure 9 Modelled depth to bedrock using the approximation of the catenary valley shape and a maximum depth of 20 m at the valley bottom. During calibration, the maximum depth may be optimized and the depth re-calculated. The evaluation of the conditions required for calibration included many iterative evaluations regarding both physical and temporal parameters. Many EMUs with lower hydraulic conductivity were highly sensitive to initial conditions, which responded differently to systematic application of artificial recharge. Model sensitivity arose primarily due to variables requiring calibration, such as the adjustment of hydraulic conductivity. Time required to calibrate iteratively also increased significantly due to the calculation of the time steps required in each day. From saturated conditions with variable values of hydraulic conductivity, the number of additional iterations required per daily time step ranged from 12 to 14, depending on 51 the volume applied to storage and the hydraulic conductivity values assigned. Higher values of hydraulic conductivity resulted in increased difficulty in satisfying the von Neumann criterion. The explicit solving method without compliance with the von Neumann criterion caused numerical instability at the daily time step. The degree of instability is illustrated in Figure 10, which provides a comparison between time step methods. For this particular watershed the number of additional time steps per day varied with respect to the amount of basin volume preloading. These additional iterations are likely attributed to the model sensitivity to initial conditions in lower hydraulic conductivity EMUs. Daily time step Time step calculated 1,000,000 100,000 10,000 0 50 100 150 200 250 300 350 Day Figure 10 Comparison of numeric instability prior to compliance with the von Neumann Criterion using applied fully saturated conditions throughout the watershed A range values of hydraulic conductivity and effective porosity for each EMU type were derived from literature (Batu, 1998) and were applied to evaluate the modelled response to groundwater flow. An optimal set of values was used to create the model testing results presented herein are included in Table 5. The chosen values were a result of representative values in the literature, as well as the evaluation of time-step sensitivity (Batu, 1998). Additional calibration could be achieved through model application and verification using objective functions. 52 Simulated Baseflow (m3/day) Effective Porosity Saturated Hydraulic Surficial Geology Type (fraction) Conductivity (m/s) Lacustrine 0.20 10-08 Channel Alluvium 0.30 10-04 Fan and Outwash 0.20 10-06 Alluvium Bedrock Outcrop (0) (0) Talus/Colluvium 0.30 10-03 Mass Waste 0.30 10-04 Till- Thick, Draped 0.20 10-09 Till- Valley Slope 0.25 10-08 Till- High Alpine 0.25 10-07 Periglacial 0.30 10-03 Table 5 Estimates of effective porosity and hydraulic conductivity used in model runs (Batu, 1998) The attenuation of high piezometric gradients can be noted in Table 5, which presents the surface flow results of the basin draining from fully saturated conditions over 15 years of repeat iterations applying the method of groundwater discharge recharging into downslope EMUs. A record of groundwater discharge was also useful in determining the basin response to drainage. Figure 13 presents the spatial distribution of cumulative discharge over the course of 5 years using the same model results that are included in Figure 11 and Figure 12. 53 1.4E+09 Runoff 1.2E+09 1E+09 800000000 600000000 400000000 200000000 0 0 500 1000 1500 2000 Days Figure 11 Attenuation of fully saturated conditions over 2000 days 100000 90000 80000 70000 60000 50000 40000 30000 3000 4000 5000 6000 7000 8000 9000 10000 Days Figure 12 At later stages in attenuation the presence of freshet is apparent in the hydrograph annually. 54 Discharge (m3/day) Discharge (m3/day) Figure 13 Distribution of cumulative groundwater discharge to the surface (not including basin outlet) after 5 years of iterative basin volume attenuation. A satellite image is included for comparison (ESRI, 2014) Groundwater discharging to the land surface was found to be more representative of realistic conditions when repeatedly applied as recharge, as the method applying this volume directly to baseflow yielded exceedingly rapid attenuation of precipitation events. The difference can be noted in Figure 15 and Figure 15. As such, subsequent model runs were completed attempting to recharge runoff in a downslope direction (i.e,. into EMUs downstream of any surface water). 55 Runoff 2000000 1800000 1600000 1400000 1200000 1000000 800000 600000 400000 200000 0 14235 14285 14335 14385 14435 14485 14535 14585 Days Figure 14 Comparison of groundwater discharge methods during iterative attenuation of fully saturated conditions after 14235 days with discharged groundwater being recharged downslope runoff 42000 41000 40000 39000 38000 37000 36000 14235 14285 14335 14385 14435 14485 14535 14585 Days Figure 15 Comparison of groundwater discharge methods during iterative attenuation of fully saturated conditions after 14235 days with no additional recharge The spatial changes in hydraulic head after flow simulations were computed were relatively subtle with respect to the size of the watershed. Mapped hydraulic head under preloaded 56 Discharge (m3/day) Discharge (m3/day) conditions and after 5 years of iteration are presented in Figure 16. Only the valley bottoms and adjacent hillslopes are presented to increase the contrast, whereby variability can be noted. Figure 16 Water table elevation before (left) and after (right) the attenuation of fully saturated pressure conditions throughout valleys Assumptions and Limitations Keeping in mind that the main objective of this study is to obtain first-order estimates of groundwater flow, given the available information from a remote watershed. This should also be noted when considering the approximation of the glacial valley as a catenary, which is an estimation applied at the watershed scale. Many attributes dictate the shape of glacial valleys, such as glacial erosional capacity, which includes glacier ice flux and thermal regime, bedrock type and tectonic activity (Anderson et al., 2006, Yingkui et al., 2001). The modelling methods presented are highly scale-dependent, as the geomorphic dataset chosen to discretize the model domain must represent equivalent porous media that reflect basin-scale modelling. The 57 objective, as described by Eaton (2006), is to simulate groundwater flow on a volumetrically averaged basis. A geomorphic dataset that is too coarse for the basin size may inhibit the success of the overall estimation of hydraulic head distribution and temporal basin discharge. The distribution of heterogeneity is inferred through the macro-scale differences in depositional processes, which is not a source of empirical hydraulic conductivity input. Each EMU is treated as a homogeneous, isotropic unit with effective porosity equal to the specific yield, which is an approximation. Additional scaling limitations arise from the initial conditions and boundary conditions assigned to the watershed. Watershed initial conditions can be iteratively optimized, however the addition of confounding variables to be optimized decreases the efficiency of the model calibration. The evaluation of the additional sensitivity as a result of increased variables was not fully addressed in this study, as it was tailored to address mainly the unique discretization methods. Another noteworthy limitation in boundary conditions include the lack of simulation of groundwater recharging bedrock aquifers (or equivalent aquiclude). The development of a groundwater model using this technique may calibrate well against observed data, however the conditions may inherently include a portion of flow in fractured or conductive weathered bedrock. Watersheds underlain by rock that is resistant to erosion have been shown to host a significant portion of the groundwater that contributes to baseflow, which was found in a steep granitic watershed by Katsura et al (2008). This may introduce uncertainty when evaluating the spatial distribution of groundwater flow, although groundwater flow in bedrock aquifers has not been found to contributed significantly to event-based hillslope runoff and routing (Gabrielli et al., 2012), which indicates this modelling method is suitable for relatively rapid groundwater flow simulations. Recharge does not occur instantaneously, as is assumed in this modelling method. In reality, groundwater movement through the unsaturated zone will create further delay and promote 58 increased variability in storage/discharge potential throughout the basin, especially when considering the distribution of heterogeneity (Harman and Sivapalan, 2009). Regions affected by permafrost, or where the subsurface is seasonally frozen, will also dictate the distribution of recharge and create no flow boundaries. These regions are likely to control groundwater drainage at high elevations, such as the neutralization of pore space in talus slopes (Scapozza et al., 2011), where high hydraulic conductivities are otherwise typically noted. Approximately one square km of permafrost may exist in the study area (Gruber, 2012), which is not expected to contribute significantly to the overall basin hydrodynamics. Discussion Methods utilized in this study include first, the effects of elevation and flow direction on topographic-driven flow; and second, the hydrogeologic variability delineated using geomorphic features. This approach does not follow the typical finite-difference approach, in which the subsurface is discretized using a rectangular mesh. Irregularities in the shapes of geomorphic units warranted an alternative approach that also discounts grid cell to grid cell calculations that would be computationally intensive at the watershed scale. The success of this method of discretization is largely represented by the spatial representation of discharge. This dataset was beneficial in determining the result of convergent/divergent flow conditions based on the automatic discretization method utilized. Regions of highest discharge could be noted in regions where the satellite image is sparse in vegetation, which is likely attributed to significant surface outwash. Using glacial valley shape as a proxy for sediment depth is considered relevant due to results from hillslope flow analyses. Hopp and McDonnell (2009) found that drainage efficiency is reduced at lower slopes, which aligns well with the distribution of higher hydraulic conductivity and lower depths to bedrock that were indicated by the results of the geomorphometric methods utilized. Lower drainage efficiency in valley bottoms is coincident with longer groundwater 59 residence times, thus slower attenuation of pressure conditions and discharge of groundwater, which was noted on the map of discharge (Figure 12). The discretization and computation technique used appear to be opportunistic for increased development as an alternative method of watershed-scale flow simulation in mountain watersheds. Using a detailed geomorphometric analysis, a high degree of variability is demonstrated in the discretization. This becomes apparent when considering a typical finite- difference approach utilizing a rectangular mesh. The lack of discharge in the majority of the watershed confirms groundwater dominance in the study area, which would need to be accounted for when applying a grid for finite-difference computation. It would take abundant work to apply a representative grid resolution to the entire valley system, which is often only evaluated at areas of interest, and requires abundant input data (Smerdon et al., 2009). The study area exhibits groundwater dominance at the upper reaches, which may not have been adequately addressed using alternate discretization techniques. Future work could include further evaluation of the methods utilized in this study by comparison to other discretization techniques employed in the same region and time period. Conclusions This study employed a unique method of automatic watershed discretization in a data-scarce region. Results showed a favourable use of geomorphometric attributes to estimate the distribution of heterogeneity in the subsurface, which was evaluated by a sensitivity analysis and the attenuation of artificial pressure conditions. The application of “virtual” experiments such as this study are beneficial for the evaluation of discretization methods for complex remote terrain such as diverse mountain watersheds. These methods are useful in addressing the complexity of hillslope-discharge relationships, particularly in flow conditions that tend to rapidly drain most of the volume from large precipitation events (Hopp and McDonnell, 2009). The opportunity for advancement of this technique is robust given the structure of implemented programming. The model is coded such that substitutions for methods of the calculation of 60 several calibrated parameters can be augmented. Other means of depth estimation can be easily applied, and optimized to evaluate alternative boundary condition implementation. Similar to depth estimations, the ability to apply recharge is also inherent in the structure of the model programming. Additional research utilizing these methods is recommended using appropriate pre-existing hydrologic models, such as GENESYS. The quantitative performance of the model can then be compared to observed hydrographs from basins of different sizes. Following additional research, it is expected that augmentation can be easily applied to this modeling methodology using the strong base of code. 61 Chapter 5: Thesis Conclusions and Recommendations This study focused on the integration of remote sensing data, geomorphic principles, theoretical distributions of heterogeneity, basin discretization methods, and saturated flow computation to apply a novel technique to understand groundwater behaviour in mountain watersheds. Methods included a geomorphometric analysis to computationally simulate the distribution of geomorphic landforms, which were used to estimate heterogeneity in the shallow subsurface and provide opportunity to evaluate groundwater flow. Morphometric attributes of various landforms were studied and compared to their genetic origin to identify potential landforms. The resulting landforms were subsequently subdivided into equivalent porous media units based on the theoretical distribution of heterogeneity within landform types. EMUs were evaluated as irregular units used to discretize a saturated groundwater flow model. Groundwater flow was calculated using recharge simulated by any hydrometerological model and was routed using Darcian flow from EMU to EMU. An integral part of this research was the enhancement of techniques to help understand groundwater flow in mountain watersheds. This thesis demonstrates the opportunity in bridging multiple research areas (spatial modelling, geomorphology, hydrology), where unique results may be obtained without expensive field programs. This is not to say field programs are not important, which is demonstrated by the need for field verification for portions of this research. These methods are high-level approaches to approximate complicated environmental phenomenon in order to make more informed decisions related to the planning of empirical data gathering and additional detailed modelling. This study presented research that is well suited for additional advancement, particularly in the application of the unique basin discretization technique. Additional methods to augment the geomorphometric analysis are available, and the model is coded such that implementation would be simple and somewhat modular. The object-oriented programming approach utilized is 62 beneficial for this type of additional implementation, and further modification is encouraged. In addition, the groundwater modelling method stands to benefit from additional discretization. This modelling methodology could be augmented to increase numerical robustness and computational efficiency by employing an implicit finite-volume solver. The nature of the discretization methods presented herein are well-suited to be utilized in a finite-volume solver due to the irregular nature of the distribution of EMUs. One caveat of the finite-difference method is the need for unit boundaries to be perpendicular to the direction of unit flow calculation (Panday et al., 2013). The method of discretization presented in this study automatically computes the width of unit connections perpendicular to the flow direction, which leaves it well-suited for the finite-volume approach. A brief description of the finite-volume method of groundwater flow and how it would be beneficial to be implemented into the methods presented in this study are described below. The finite-difference groundwater flow equation was derived using a similar approach to the MODFLOW (Harbaugh, 2005) control volume finite-difference application of the groundwater flow equation: 𝜕ℎ ∇ ∙ (𝐾∇ℎ) = 𝑆𝑠 + 𝑊 𝜕𝑡 (8) Where K is hydraulic conductivity [L/T], h is hydraulic head [L], Ss is specific storage [1/L], t is time [T], and W is the source/sink term [1/T] used to account for input and output from each volume. For a finite-volume approximation the expansion of the groundwater flow equation yields the finite-difference equation for a control volume that accounts for the inflow and outflow of a call (n) based on m neighbors using: 𝐻𝐶𝑂𝐹𝑛(ℎ𝑛) + ∑ 𝐶𝑛𝑚(ℎ𝑚 − ℎ𝑛) 𝑚∈𝜂𝑛 = 𝑅𝐻𝑆𝑛 (9) 63 −𝑆𝑠𝑛𝑉 −𝑆 𝑉 ℎ 𝑡−1 𝐻𝐶𝑂𝐹 = 𝑛 𝑅𝐻𝑆 𝑠𝑛 𝑛 𝑛Where 𝑛 and, 𝑛 = represent the volumetric storage terms (SSn Δ𝑡 Δ𝑡 [1/L] and Vn [L3])over a given time step and initial hydraulic head (Harbough, 2005; MODFLOW-USG manual, 2013). This method is suitable for complex grid geometry that requires only inflow and outflow by means of connecting cells sharing the boundary. The boundaries can be represented by irregular polygons (with several caveats described in further detail below) which allows for irregular discretization. The finite-volume approximation of groundwater flow that could be implemented using the methods presented in this study would be one-dimensional, comprised of inflow and outflow approximated linearly in the downslope direction. Implementation of groundwater mass balance flux computation at each of the grid polygons includes a conductance and gradient balance using: 𝐻𝐶𝑂𝐹𝑛(ℎ𝑛) + 𝐶𝑚 (𝑖𝑛) 𝑛(ℎ𝑚 (𝑖𝑛) ℎ𝑛) + 𝐶𝑚 (𝑜𝑢𝑡) 𝑛(ℎ𝑚 (𝑜𝑢𝑡)ℎ𝑛) = 𝑅𝐻𝑆𝑛 (10) Which is equation 9 expanded to include in and out subscripts for conductance and hydraulic head gradient at incoming and outgoing connecting cells. Equation10 could then be solved implicitly at each EMU, given the current conditions and position in the model domain. 64 References Ajami, H., Troch, P. A., Maddock, T., Meixner, T. and Eastoe, C. (2011) 'Quantifying mountain block recharge by means of catchment-scale storage-discharge relationships', Water Resources Research, 47(4), pp. W04504. Anderson, R. S., Molnar, P. and Kessler, M. A. (2006) 'Features of glacial valley profiles simply explained', Journal of Geophysical Research: Earth Surface, 111(F1), pp. F01004. Arp, C. D., Gooseff, M. N., Baker, M. A. and Wurtsbaugh, W. (2006) 'Surface-water hydrodynamics and regimes of a small mountain stream–lake ecosystem', Journal of Hydrology, 329(3–4), pp. 500-513. Baraer, M., McKenzie, J. M., Mark, B. G., Bury, J. and Knox, S. (2009) 'Characterizing contributions of glacier melt and groundwater during the dry season in a poorly gauged catchment of the Cordillera Blanca (Peru)', Adv. Geosci., 22, pp. 41-49. Batu, V. (1998) Aquifer hydraulics : a comprehensive guide to hydrogeologic data analysis. New York, John Wiley & Sons. Becht, M., Haas, F., Heckmann, T. and Wichmann, V. (2005) 'Investigating sediment cascades using field measurements and spatial modelling', in Walling, D.E. & Horowitz, A.J. (eds.) Sediment Budgets 1 IAHS Publ. 291. IAHS Press, Wallingford, UK., pp. 206-213. Bell, J. E., Autry, C. W., Mollenkopf, D. A. and Thornton, L. M. (2012) 'A Natural Resource Scarcity Typology: Theoretical Foundations and Strategic Implications for Supply Chain Management', Journal of Business Logistics, 33(2), pp. 158-166. Benda, L., Hassan, M. A., Church, M. and May, C. L. (2005) 'Geomorphology of steepland headwaters: The transition from hillslopes to channels', Journal of the American Water Resources Association, 41(4), pp. 835-851. Beniston, M., Diaz, H. F. and Bradley, R. S. (1997) 'Climatic change at high elevation sites: An overview', Climatic Change, 36(3-4), pp. 233-251. Bishop, M. P., James, L. A., Shroder Jr, J. F. and Walsh, S. J. (2012) 'Geospatial technologies and digital geomorphological mapping: Concepts, issues and research', Geomorphology, 137(1), pp. 5-26. Caballero, Y., Jomelli, V., Chevallier, P. and Ribstein, P. (2002) 'Hydrological characteristics of slope deposits in high tropical mountains (Cordillera Real, Bolivia)', Catena, 47(2), pp. 101-116. Cairns, D., Jiskoot, H., Byrne, J., McKenzie, J. and Johnson, D. (2014 Submitted) 'Novel geomorphometric classification of heterogeneity of unconsolidated sediments for hydrogeologic characterization in mountain watersheds', In Preparation (Journal of Geomorphology). Cannon, M. R. (1996) Geology and groundwater resources of the Blackfeet Indian Reservation, northwestern Montana: U.S. Geological Survey Hydrologic Atlas 737, 2 sheets, scale 1:125,000. USGS Numbered Series. Carrara, P. E. (1987) 'Holocene and Latest Pleistocene Glacial Chronology, Glacier National Park, Montana', Canadian Journal of Earth Sciences, 24(3), pp. 387-395. Carrara, P. E. (1990) Surficial geologic map of Glacier National Park, Montana, 1:100,000. U.S. Geological Survey. Chase, C. G. (1992) 'Fluvial landsculpting and the fractal dimension of topography', Geomorphology, 5(1–2), pp. 39-57. Clow, D. W., Schrott, L., Webb, R., Campbell, D. H., Torizzo, A. and Dornblaser, M. (2003) 'Ground Water Occurrence and Contributions to Streamflow in an Alpine Catchment, Colorado Front Range', Ground Water, 41(7), pp. 937-950. 65 Coulthard, T. J. and Van de Wiel, M. J. (2012) 'Modelling river history and evolution', Philosophical Transactions of the Royal Society a-Mathematical Physical and Engineering Sciences, 370(1966), pp. 2123-2142. Crank, J. and Nicolson, P. (1996) 'A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type', Advances in Computational Mathematics, 6(1), pp. 207-226. Cuthbert, M. O., Mackay, R., Tellam, J. H. and Thatcher, K. E. (2010) 'Combining unsaturated and saturated hydraulic observations to understand and estimate groundwater recharge through glacial till', Journal of Hydrology, 391(3/4), pp. 263-276. de Marsily, G., Delay, F., Gonçalvès, J., Renard, P., Teles, V. and Violette, S. (2005) 'Dealing with spatial heterogeneity', Hydrogeology Journal, 13(1), pp. 161-183. de Marsily, G., Delay, F., Teles, V. and Schafmeister, M. T. (1998) 'Some current methods to represent the heterogeneity of natural media in hydrogeology', Hydrogeology Journal, 6(1), pp. 115-130. Dehn, M., Gärtner, H. and Dikau, R. (2001) 'Principles of semantic modeling of landform structures', Computers & Geosciences, 27(8), pp. 1005-1010. Dietrich, W. E., Bellugi, D. G., Sklar, L. S., Stock, J. D., Heimsath, A. M. and Roering, J. J. (2003) 'Geomorphic Transport Laws for Predicting Landscape form and Dynamics', Prediction in Geomorphology: American Geophysical Union, pp. 103-132. Dietrich, W. E., Wilson, C. J., Montgomery, D. R. and McKean, J. (1993) 'Analysis of Erosion Thresholds, Channel Networks, and Landscape Morphology Using a Digital Terrain Model', The Journal of Geology, 101(2), pp. 259-278. Dunlop, P. and Clark, C. D. (2006) 'The morphological characteristics of ribbed moraine', Quaternary Science Reviews, 25(13/14), pp. 1668-1691. Dymond, J. R., Derose, R. C. and Harmsworth, G. R. (1995) 'Automated mapping of land components from digital elevation data', Earth Surface Processes and Landforms, 20(2), pp. 131-137. Eaton, T. T. (2006) 'On the importance of geological heterogeneity for flow simulation', Sedimentary Geology, 184(3–4), pp. 187-201. Fisher, A. T. and Barnhill, M. (1998) 'The relationship between hydrogeologic properties and sedimentary facies: an example from', Ground Water, 36(6), pp. 901. Fogg, G. E., Noyes, C. D. and Carle, S. F. (1998) 'Geologically based model of heterogeneous hydraulic conductivity in an alluvial setting', Hydrogeology Journal, 6(1), pp. 131-143. Freeze, R. A. (1975) 'A stochastic-conceptual analysis of one-dimensional groundwater flow in nonuniform homogeneous media', Water Resources Research, 11(5), pp. 725-741. Gabrielli, C. P., McDonnell, J. J. and Jarvis, W. T. (2012) 'The role of bedrock groundwater in rainfall–runoff response at hillslope and catchment scales', Journal of Hydrology, 450– 451, pp. 117-133. Genereux, D. and Bandopadhyay, I. (2001) 'Numerical investigation of lake bed seepage patterns: effects of porous medium and lake properties', Journal of Hydrology, 241(3–4), pp. 286- 303. Giles, P. T. (2010) 'Investigating the use of alluvial fan volume to represent fan size in morphometric studies', Geomorphology, 121(3–4), pp. 317-328. Gleeson, T. and Manning, A. H. (2008) 'Regional groundwater flow in mountainous terrain: Three-dimensional simulations of topographic and hydrogeologic controls', Water Resources Research, 44(10). Green, T. R., Taniguchi, M., Kooi, H., Gurdak, J. J., Allen, D. M., Hiscock, K. M., Treidel, H. and Aureli, A. (2011) 'Beneath the surface of global change: Impacts of climate change on groundwater', Journal of Hydrology, 405(3–4), pp. 532-560. Gruber, S. (2012) 'Derivation and analysis of a high-resolution estimate of global permafrost zonation', The Cryosphere, 6(1), pp. 221-233. 66 Gurrieri, J. T. and Furniss, G. (2004) 'Estimation of groundwater exchange in alpine lakes using non-steady mass-balance methods', Journal of Hydrology, 297(1–4), pp. 187-208. Hall, K. and Thorn, C. (2011) 'The historical legacy of spatial scales in freeze–thaw weathering: Misrepresentation and resulting misdirection', Geomorphology, 130(1–2), pp. 83-90. Harbaugh, A. W. (2005) 'MODFLOW-2005, The U.S. Geological Survey Modular Ground- Water Model—the Ground-Water Flow Process', U.S. Geological Survey Techniques and Methods 6-A16. Harbor, J. M. (1995) 'Development of glacial-valley cross sections under conditions of spatially variable resistance to erosion', Geomorphology, 14(2), pp. 99-107. Harman, C. and Sivapalan, M. (2009) 'Effects of hydraulic conductivity variability on hillslope- scale shallow subsurface flow response and storage-discharge relations', Water Resources Research, 45. Hashimoto, A., Oguchi, T., Hayakawa, Y., Lin, Z., Saito, K. and Wasklewicz, T. A. (2008) 'GIS analysis of depositional slope change at alluvial-fan toes in Japan and the American Southwest', Geomorphology, 100(1–2), pp. 120-130. Hayashi, M. and van der Kamp, G. (2009) 'Progress in Scientific Studies of Groundwater in the Hydrologic Cycle in Canada, 2003-2007', Canadian Water Resources Journal, 34(2), pp. 177-186. Hevesi, J. A., Flint, A. L. and Flint, L. E. (2002) 'Preliminary Estimates of Spatially Distributed Net Infiltration and Recharge for the Death Valley Region, Nevada- California', U.S. Geological Survey Water Resources Investigations Report, 02-4010(Sacramento, California, 2001). Hirano, M. and Aniya, M. (1988) 'A rational explanation of cross-profile morphology for glacial valleys and of glacial valley development', Earth Surface Processes and Landforms, 13(8), pp. 707-716. Hopp, L. and McDonnell, J. J. (2009) 'Connectivity at the hillslope scale: Identifying interactions between storm size, bedrock permeability, slope angle and soil depth', Journal of Hydrology, 376(3–4), pp. 378-391. Hubbard, S. S. and Linde, N. (2011) '2.15 - Hydrogeophysics', in Editor-in-Chief: Peter, W. (ed.) Treatise on Water Science. Oxford: Elsevier, pp. 401-434. Jenson, S. K. and Domingue, J. O. (1988) 'Extracting Topographic Structure From Digital Elevation Data for Geographic Information System Analysis', Photogrammetric Engineering and Remote Sensing, 54(11), pp. 1593-1600. Johnson, R. H. (2007) 'Ground Water Flow Modeling with Sensitivity Analyses to Guide Field Data Collection in a Mountain Watershed', Ground Water Monitoring & Remediation, 27(1), pp. 75-83. Kao, Y.-H., Liu, C.-W., Wang, S.-W. and Lee, C.-H. (2012) 'Estimating mountain block recharge to downstream alluvial aquifers from standard methods', Journal of Hydrology, 426–427, pp. 93-102. Katsura, S. y., Kosugi, K. i., Mizutani, T., Okunaka, S. and Mizuyama, T. (2008) 'Effects of bedrock groundwater on spatial and temporal variations in soil mantle groundwater in a steep granitic headwater catchment', Water Resources Research, 44(9), pp. W09430. Kirchner, J. W. (2009) 'Catchments as simple dynamical systems: Catchment characterization, rainfall-runoff modeling, and doing hydrology backward', Water Resources Research, 45(2), pp. W02429. Koltermann, C. E. and Gorelick, S. M. (1992) 'Paleoclimatic Signature in Terrestrial Flood Deposits', Science, 256(5065), pp. 1775-1782. Koltermann, C. E. and Gorelick, S. M. (1996) 'Heterogeneity in sedimentary deposits: A review of structure-imitating, process-imitating, and descriptive approaches', Water Resources Research, 32(9), pp. 2617-2658. 67 Lapp, S., Byrne, J., Townshend, I. and Kienzle, S. (2005) 'Climate warming impacts on snowpack accumulation in an alpine watershed', International Journal of Climatology, 25(4), pp. 521-536. Larsen, J. C. (2004) Rectified images of selected geologic maps in the northern Rockies area, Idaho, Montana, Washington, and Wyoming. Menlo Park, CA: U.S. Geological Survey. MacDonald, R. J., Byrne, J. M. and Kienzle, S. W. (2009) 'A Physically Based Daily Hydrometeorological Model for Complex Mountain Terrain', Journal of Hydrometeorology, 10(6), pp. 1430-1446. MacDonald, R. J., Byrne, J. M., Kienzle, S. W. and Larson, R. P. (2011) 'Assessing the potential impacts of climate change on mountain snowpack in the St. Mary River watershed, Montana', Journal of Hydrometeorology, 12(2), pp. 262-273. Maeda, K. and Yuasa, T. (2009) Performance estimation of countermeasures for falling rock using DEM. Prediction and Simulation Methods for Geohazard Mitigation. edited by Fusao Oka, Akira Murakami, Sayuri Kimoto” CRC Press, Taylor & Francis Group, Baton Rouge, FL, USA. Magruder, I. A., Woessner, W. W. and Running, S. W. (2009) 'Ecohydrologic Process Modeling of Mountain Block Groundwater Recharge', Ground Water, 47(6), pp. 774-785. Manning, A. H. and Solomon, D. K. (2005) 'An integrated environmental tracer approach to characterizing groundwater circulation in a mountain block', Water Resources Research, 41(12). Marquinez, J., Duarte, R. M., Farias, P. and Sanchez, M. J. (2003) 'Predictive GIS-based model of rockfall activity in Mountain Cliffs', Natural Hazards, 30(3), pp. 341-360. Marquis, G., Baldassarri, T., Hofer, T., Romeo, R. and Wolter, P. (2012) 'FAO's Current Engagement in Sustainable Mountain Development', Mountain Research and Development, 32(2), pp. 226-230. McCleary, R. J., Hassan, M. A., Miller, D. and Moore, R. D. (2011) 'Spatial organization of process domains in headwater drainage basins of a glaciated foothills region with complex longitudinal profiles', Water Resources Research, 47, W05505. McCloskey, T. F. and Finnemore, E. J. (1996) 'Estimating hydraulic conductivities in an alluvial basin from sediment facies models', Ground Water, 34(6), pp. 1024-1032. McClymont, A. F., Roy, J. W., Hayashi, M., Bentley, L. R., Maurer, H. and Langston, G. (2011) 'Investigating groundwater flow paths within proglacial moraine using multiple geophysical methods', Journal of Hydrology, 399(1–2), pp. 57-69. McKean, J. and Roering, J. (2004) 'Objective landslide detection and surface morphology mapping using high-resolution airborne laser altimetry', Geomorphology, 57(3–4), pp. 331-351. Miliaresis, G. C. and Argialas, D. P. (2000) 'Extraction and delineation of alluvial fans from digital elevation models and landsat thematic mapper images', Photogrammetric Engineering and Remote Sensing, 66(9), pp. 1093-1101. Millar, C. I., Westfall, R. D. and Delany, D. L. 'Thermal and hydrologic attributes of rock glaciers and periglacial talus landforms: Sierra Nevada, California, USA', Quaternary International. Minár, J. and Evans, I. S. (2008) 'Elementary forms for land surface segmentation: The theoretical basis of terrain analysis and geomorphological mapping', Geomorphology, 95(3–4), pp. 236-259. Mitchell, S. G. and Montgomery, D. R. (2006) 'Influence of a glacial buzzsaw on the height and morphology of the Cascade Range in central Washington State, USA', Quaternary Research, 65(1), pp. 96-107. Motaghian, H. R. and Mohammadi, J. (2011) 'Spatial Estimation of Saturated Hydraulic Conductivity from Terrain Attributes Using Regression, Kriging, and Artificial Neural Networks', Pedosphere, 21(2), pp. 170-177. 68 Napieralski, J., Harbor, J. and Li, Y. (2007) 'Glacial geomorphology and geographic information systems', Earth-Science Reviews, 85(1–2), pp. 1-22. Napieralski, J. and Nalepa, N. (2010) 'The application of control charts to determine the effect of grid cell size on landform morphometry', Computers & Geosciences, 36(2), pp. 222-230. Panday, S., Langevin, C. D., Niswonger, R. G., Ibaraki, M. and Hughes, J. D. (2013) 'MODFLOW–USG version 1: An unstructured grid version of MODFLOW for simulating groundwater flow and tightly coupled processes using a control volume finite-difference formulation', U.S. Geological Survey Techniques and Methods, book 6(chap. A45), pp. 66 p. Pike, R. J. (2002) 'A Bibliography of Terrain Modeling (Geomorphometry), the Quantitative Representation of Topography- Supplement 4.0', USGS Open-File Report, (02-465). Roering, J. J., Kirchner, J. W. and Dietrich, W. E. (2001) 'Hillslope evolution by nonlinear, slope- dependent transport: Steady state morphology and equilibrium adjustment timescales', Journal of Geophysical Research: Solid Earth, 106(B8), pp. 16499-16513. Ronayne, M. J., Houghton, T. B. and Stednick, J. D. (2012) 'Field characterization of hydraulic conductivity in a heterogeneous alpine glacial till', Journal of Hydrology, 458–459, pp. 103-109. Ross, C. P. (1954) 'The Geology of Glacier-National-Park and Vicinity in Montana', Science, 119(3097), pp. 658-659. Roy, J. W. and Hayashi, M. (2009) 'Multiple, distinct groundwater flow systems of a single moraine–talus feature in an alpine watershed', Journal of Hydrology, 373(1–2), pp. 139- 150. Saadat, H., Bonnell, R., Sharifi, F., Mehuys, G., Namdar, M. and Ale-Ebrahim, S. (2008) 'Landform classification from a digital elevation model and satellite imagery', Geomorphology, 100(3/4), pp. 453-464. Saito, K. and Oguchi, T. (2005) 'Slope of alluvial fans in humid regions of Japan, Taiwan and the Philippines', Geomorphology, 70(1–2), pp. 147-162. Sanchez-Vila, X., Guadagnini, A. and Carrera, J. (2006) 'Representative hydraulic conductivities in saturated groundwater flow', Reviews of Geophysics, 44(3), pp. RG3002. Sass, O. (2006) 'Determination of the internal structure of alpine talus deposits using different geophysical methods (Lechtaler Alps, Austria)', Geomorphology, 80(1–2), pp. 45-58. Scanlon, B. R., Healy, R. W. and Cook, P. G. (2002) 'Choosing appropriate techniques for quantifying groundwater recharge (vol 10, pg 18, 2002)', Hydrogeology Journal, 10(2), pp. 347-347. Scapozza, C., Lambiel, C., Baron, L., Marescot, L. and Reynard, E. (2011) 'Internal structure and permafrost distribution in two alpine periglacial talus slopes, Valais, Swiss Alps', Geomorphology, 132(3–4), pp. 208-221. Schmidt, J. and Hewitt, A. (2004) 'Fuzzy land element classification from DTMs based on geometry and terrain position', Geoderma, 121(3/4), pp. 243-256. Slaymaker, O. (2006) 'Towards the identification of scaling relations in drainage basin sediment budgets', Geomorphology, 80(1–2), pp. 8-19. Smerdon, B. D., Allen, D. M., Grasby, S. E. and Berg, M. A. (2009) 'An approach for predicting groundwater recharge in mountainous watersheds', Journal of Hydrology, 365(3-4), pp. 156-172. Stone, D. B., Moomaw, C. L. and Davis, A. (2001) 'Estimating recharge distribution by incorporating runoff from mountainous areas in an alluvial basin in the great basin region of the southwestern United States', Ground Water, 39(6), pp. 807-818. Svensson, H. (1959) 'Is the cross-section of a glacial valley a parabola?', Journal of Glaciology, 3(25), pp. 362-363. Tarboton, D. G., Bras, R. L. and Rodriguez-Iturbe, I. (1991) 'On the Extraction of Channel Networks from Digital Elevation Data', Hydrological Processes, 5(1), pp. 81-100. 69 Trevisani, S., Cavalli, M. and Marchi, L. (2012) 'Surface texture analysis of a high-resolution DTM: Interpreting an alpine basin', Geomorphology, 161–162, pp. 26-39. UNEP (2007) 'Fourth Global Environment Outlook: Environment for Development.', Geneva, Switzerland. USGS (2009) 'National Elevation Dataset (NED) Raster Digital Data'. Available at: http://nationalmap.gov (Accessed January 2012). van Leeuwen, C. J., Frijns, J., van Wezel, A. and van de Ven, F. H. M. (2012) 'City Blueprints: 24 Indicators to Assess the Sustainability of the Urban Water Cycle', Water Resources Management, 26(8), pp. 2177-2197. Vannametee, E., Babel, L. V., Hendriks, M. R., Schuur, J., de Jong, S. M., Bierkens, M. F. P. and Karssenberg, D. (2014) 'Semi-automated mapping of landforms using multiple point geostatistics', Journal of Geomorphology, 221 (2014), pp. 298-319. Vienken, T. and Dietrich, P. (2011) 'Field evaluation of methods for determining hydraulic conductivity from grain size data', Journal of Hydrology, 400(1–2), pp. 58-71. Viviroli, D. and Weingartner, R. (1999) 'The hydrological significance of mountains: from regional to global scale', Hydrol. Earth Syst. Sci., 8(6), pp. 1017-1030. Wada, Y., van Beek, L. P. H., van Kempen, C. M., Reckman, J., Vasak, S. and Bierkens, M. F. P. (2010) 'Global depletion of groundwater resources', Geophysical Research Letters, 37 issue L20402. Webster, T. L., Murphy, J. B., Gosse, J. C. and Spooner, I. (2006) 'The application of lidar- derived digital elevation model analysis to geological mapping: an example from the Fundy Basin, Nova Scotia, Canada', Canadian Journal of Remote Sensing, 32(2), pp. 173-193. Wei, M., Allen, D., Kohut, A., Grasby, S., Ronneseth, K. and Turner, B. (2009) 'Understanding the Types of Aquifers in the Canadian Cordillera Hydrogeologic Region to Better Manage and Protect Groundwater', Watershed Management Bulletin, 13(1), pp. 10-18. Whipple, K., Kirby, E. and Brocklehurst, S. (1999) 'Geomorphic limits to climate-induced increases in topographic relief', Nature, 401, pp. 39-43. Wichmann, V., Heckmann, T., Haas, F. and Becht, M. (2009) 'A new modelling approach to delineate the spatial extent of alpine sediment cascades', Geomorphology, 111(1–2), pp. 70-78. Wilson, J. L. and Guan, H. (2004) 'Mountain Block Hydrology and Mountain-Front Recharge', Washington, DC: American Geophysical Union, Groundwater Recharge in a Desert Environment: The Southwestern United States, pp. 113-137. Winter, T. C. (1999) 'Relation of streams, lakes, and wetlands to groundwater flow systems', Hydrogeology Journal, 7(1), pp. 28-45. Wittenberg, H. (1999) 'Baseflow recession and recharge as nonlinear storage processes', Hydrological Processes, 13(5), pp. 715-726. Wohl, E. (2006) 'Human impacts to mountain streams', Geomorphology, 79(3–4), pp. 217-248. Yingkui, L., Gengnian, L. and Zhijiu, C. (2001) 'Longitudinal variations in cross-section morphology along a glacial valley: a case-study from the Tien Shan, China', Journal of Glaciology, 47(157), pp. 243-250. Zientek, M. L., Derkey, P. D., Miller, R. J., Causey, J. D., Bookstrom, A. A., Carlson, M. H., Green, G. N., Frost, T. P., Boleneus, D. E., Evans, K. V., Gosen, B. S. V., Wilson, A. B., Larsen, J. C., Kayser, H. Z., Kelley, W. N. and Assmus, K. C. (2005) Spatial databases for the geology of the Northern Rocky Mountains - Idaho, Montana, and Washington. Zou, K. H., Warfield, S. K., Bharatha, A., Tempany, C. M. C., Kaus, M. R., Haker, S. J., Wells Iii, W. M., Jolesz, F. A. and Kikinis, R. (2004) 'Statistical validation of image segmentation quality based on a spatial overlap index1: scientific reports', Academic Radiology, 11(2), pp. 178-189. 70