UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Landscape organization based on application of the process domain concept for a glaciated foothills region McCleary, Richard James 2011

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
24-ubc_2011_spring_mccleary_richard.pdf [ 6.76MB ]
Metadata
JSON: 24-1.0071632.json
JSON-LD: 24-1.0071632-ld.json
RDF/XML (Pretty): 24-1.0071632-rdf.xml
RDF/JSON: 24-1.0071632-rdf.json
Turtle: 24-1.0071632-turtle.txt
N-Triples: 24-1.0071632-rdf-ntriples.txt
Original Record: 24-1.0071632-source.json
Full Text
24-1.0071632-fulltext.txt
Citation
24-1.0071632.ris

Full Text

  LANDSCAPE ORGANIZATION BASED ON APPLICATION OF THE PROCESS DOMAIN CONCEPT FOR A GLACIATED FOOTHILLS REGION  by  RICHARD JAMES McCLEARY  B.Sc., The University of Montana, 1993 M.Sc., The University of Montana, 1996   A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF  THE REQUIREMENTS FOR THE DEGREE OF  DOCTOR OF PHILOSOPHY  in  THE FACULTY OF GRADUATE STUDIES  (Geography)   THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  March 2011   © Richard James McCleary, 2011 ii  Abstract The organization of a glaciated foothills and interior plain landscape was examined by analyzing the spatial distribution of process domains.  Glaciations distribute patches of surficial materials with varying hydrologic properties across the landscape and create complex longitudinal profiles.  These two factors complicate the use of a traditional method, drainage area – local slope plots, for explaining landscape organization; hence, an empirical approach using logistic regression models was applied.  An automated procedure was developed that enabled extrapolation of model results across entire study areas.  A 57 km 2  and a 10,000 km 2  study area were divided into hydrological response units (HRUs) with an average area of less than 5 km 2 .  The HRUs were classified into categories based on topography, hydraulic connectivity, and surface type.  A preliminary drainage network was extracted from a LiDAR DEM.  The network was sub-divided into reaches, each characterized by a suite of topographic predictors.  Ground surveys entailed assigning a process domain class as hillslope, swale, seepage channel, or fluvial channel.  Various types of logistic regression models, including mixed-effects, binary, ordinal, and multinomial, were used to predict process domain class. Predictors included longitudinal profile anomalies measured at the reach scale using a normalized stream length-gradient (SL/k) index.  The candidate predictors, including SL/k index, soil moisture regime, and mean basin slope, were selected in consideration of groundwater flow systems that operate respectively at the reach, hillslope, and basin scales. Drainage area, the dominant factor controlling domain transitions, was complemented by various predictors, depending on transition and HRU type, including reach slope, and the three groundwater flow related predictors.  The study revealed features of a stabilized landscape including an extensive network of swales.  Headward channel migration into the swale network is not expected to proceed uniformly across the glaciated landscape; rather, minimal migration is expected where channel heads are anchored below confluences and within over-steepened areas, and higher migration frequencies are forecast in sections of the network that are more evenly graded in terms of both drainage area and slope.  iii  Preface A version of Chapter 2 has been published.  McCleary, R. J., and M. A. Hassan (2008), Predictive modelling and spatial mapping of fish distributions in small streams of the Canadian Rocky Mountain foothills, Can. J. Fish. Aquat. Sci., 65, 319-333.  The student‟s contribution entailed obtaining funding, conducting field work, developing the methods, completing the analysis, and writing the paper under the supervision and guidance of Marwan Hassan. In November 2010, a version of Chapter 3 was accepted for publication.  McCleary, R. J., M. A. Hassan, D. Miller, R.D. Moore, Spatial organization of process domains in headwater drainage basins of a glaciated foothills region with complex longitudinal profiles, Water Resour. Res.  The student‟s contribution involved preparing research grants, conducting field work, developing the methods, completing the analyses, and writing the manuscripts under the supervision and guidance of Marwan Hassan and R.D. Moore.  Dan Miller developed specific GIS programs and provided guidance.  Lorne Rothman provided mentorship specific to statistical analysis and programming.  Chapter 4 was prepared as a manuscript for publication.  The co-authors for Chapter 4 include Marwan Hassan, R.D. Moore, and Lorne Rothman.  The student‟s contribution involved preparing research grants, conducting field work, developing the methods, completing the analyses, and writing the manuscripts under the supervision and guidance of Marwan Hassan and R.D. Moore.  Dan Miller developed specific GIS programs.  Lorne Rothman provided mentorship specific to statistical analysis, programming, and interpretation of results. iv  Table of Contents Abstract ..................................................................................................................................... ii Preface  ..................................................................................................................................... iii Table of Contents ...................................................................................................................... iv List of Tables ........................................................................................................................... vii List of Figures ........................................................................................................................... ix Acknowledgements ..................................................................................................................xiv Chapter 1 Introduction ................................................................................................................1 1.1 Overview ...........................................................................................................................1 1.2 Aims, objectives, and organization ....................................................................................3 Chapter 2 Predictive modelling and spatial mapping of fish distributions in small streams of the Canadian Rocky Mountain Foothills ..........................................................................7 2.1 Introduction .......................................................................................................................7 2.2 Methods ............................................................................................................................9 2.2.1 Study area ...................................................................................................................9 2.2.2 Fish species .............................................................................................................. 11 2.2.3 An automated procedure to measure basin and reach characteristics.......................... 12 2.2.4 Accuracy of stream slope measures generated from a DEM ...................................... 13 2.2.5 Deriving predictive models of fish occurrence .......................................................... 13 2.2.6 Model interpretation and validation........................................................................... 16 2.2.7 Extrapolating and mapping occurrence models ......................................................... 17 2.2.8 Sources of error ........................................................................................................ 17 2.3 Results............................................................................................................................. 18 2.3.1 Sampling summary ................................................................................................... 18 2.3.2 Best ecological models ............................................................................................. 19 2.3.3 Predicted spatial distributions ................................................................................... 23 2.3.4 Best ecological / land-use models ............................................................................. 24 2.4 Discussion ....................................................................................................................... 28 v  2.4.1 Occurrence patterns of nonnative and native salmonids............................................. 28 2.4.2 Anticipated benefits from using a high-resolution DEM ............................................ 31 2.4.3 Limitations and opportunities for modelling fish occurrence using remote methods ................................................................................................................. 32 Chapter 3 Spatial organization of process domains in headwater drainage basins of a glaciated foothills region with complex longitudinal profiles ......................................... 35 3.1 Introduction ..................................................................................................................... 35 3.2 Study area........................................................................................................................ 37 3.3 Materials and methods ..................................................................................................... 42 3.3.1 Creating the overdrawn drainage network ................................................................. 42 3.3.2 Process domain field surveys .................................................................................... 43 3.3.3 Variables used to predict process domain class ......................................................... 44 3.3.4 Modelling approach .................................................................................................. 47 3.4 Results and discussion ..................................................................................................... 53 3.4.1 Linkages between geology, DEM-derived topographic indices, and field observations ........................................................................................................... 53 3.4.2 Screening candidate predictors prior to establishing relations .................................... 55 3.4.3 Relations between predictors and all four process domains ....................................... 58 3.4.4 Relations between predictors and hillslope, swale and colluvial channel domains ..... 63 3.4.5 Regional scale mapping ............................................................................................ 66 3.5 Conclusion ...................................................................................................................... 69 Chapter 4 Landscape scale modelling of process domains using LiDAR digital elevation models with factors representing reach, hillslope, and basin scale groundwater flow systems .......................................................................................................................... 70 4.1 Introduction ..................................................................................................................... 70 4.2 Study area........................................................................................................................ 76 4.2.1 Hinton region ............................................................................................................ 76 4.2.2 LiDAR data .............................................................................................................. 79 4.3 Materials and methods ..................................................................................................... 79 4.3.1 Basin similarity analysis ........................................................................................... 80 vi  4.3.2 Creating the preliminary drainage network................................................................ 83 4.3.3 Process domain field surveys .................................................................................... 85 4.3.4 Variables used to predict process domain class ......................................................... 90 4.3.5 Modelling approach .................................................................................................. 95 4.4 Results........................................................................................................................... 100 4.4.1 Basin similarity analysis ......................................................................................... 100 4.4.2 Drainage networks and flowpath obstruction at road crossings ................................ 103 4.4.3 Screening candidate predictors ................................................................................ 105 4.4.4 Ordinal, multinomial and binary logistic regression results for all four process domains ............................................................................................................... 109 4.4.5 Application of best models for mapping .................................................................. 121 4.4.6 Evidence for anchored and transient domain transitions in response to land use and climate shifts ................................................................................................. 123 4.4.7 Effects of reducing landscape variability and spatial extent on accuracy of binary logistic model for channel occurrence .................................................................. 125 4.5 Discussion ..................................................................................................................... 126 4.5.1 Basin similarity analysis ......................................................................................... 126 4.5.2 Preliminary drainage network and flowpath obstruction at road crossings ............... 127 4.5.3 Relations between process domain occurrence and predictors ................................. 127 4.5.4 Interpretations for drainage basin evolution ............................................................ 128 4.5.5 Opportunities to improve the model ........................................................................ 129 Chapter 5 Conclusion .............................................................................................................. 130 References .............................................................................................................................. 137 Appendix A. Review of categorical data modelling procedures ......................................... 156 Appendix B. Field card ..................................................................................................... 162   vii  List of Tables Table 2-1.Reach-scale variables grouped by category. ..............................................................15 Table 2-2. Logistic regression models of ecological variables predicting the probability of occurrence of brook trout, bull trout, rainbow trout and all species present. ......................20 Table 2-3. Classification analysis for high and medium predicted occurrence probability categories for best ecological models by species. ..............................................................22 Table 2-4. Logistic regression models of ecological (drainage area, reach slope, elevation, and mean basin slope) and land-use variables predicting the probability of occurrence of bull trout, rainbow trout and brook trout. .................................................................................25 Table 2-5. Logistic regression parameter estimates from the best ecological / land-use models for each of three fish species and all species present including standard error of the model coefficient (SE), odds ratio, and associated lower and upper 95% confidence intervals (CI).  .........................................................................................................................................26 Table 2-6.  Classification analysis for high occurrence category for best ecological / land-use models by species. ............................................................................................................27 Table 3-1. Profile lengths and process domains. .......................................................................44 Table 3-2. Candidate predictor variables generated for every reach with considerations for delineating process domains. ............................................................................................45 Table 3-3. Example of the conversion of an ordinal variable with four classes into three non- redundant cumulative variables.........................................................................................50 Table 3-4. Ordinal logistic regression results for the best process domain models including model selection procedure, domain classes, coefficient estimates and their 95% confidence intervals (lower and upper), standard errors (SE), significance (P), odds ratio and their 95% confidence intervals. .........................................................................................................60 Table 4-1. Variables selected to represent the three factors that influence water movement with hydrologic considerations and methods for each. ..............................................................81 Table 4-2. Bedload sediment yield May – August, 2006 from three sediment budget watersheds in the study area................................................................................................................82 Table 4-3. Key to seepage and fluvial erosion classes. ..............................................................90 Table 4-4. Reach-scale predictors grouped by category. ...........................................................91 Table 4-5. Correlations between continuous class candidate predictors for complete dataset (n=617) and subset (italics in brackets) of sites with DTW index data (n=441). .............. 106 viii  Table 4-6. Predictors selected by stepwise selection procedure repeated for 10 cross-validation datasets, including model type, predictor, mean Chi-square score (score), and frequency of selection. ........................................................................................................................ 110 Table 4-7. Ordinal logistic regression results for the best process domain model including dataset, coefficient estimates and their 95% confidence intervals (lower and upper), standard errors (SE), significance (P), odds ratio and their 95% confidence intervals. ..... 115 Table 4-8. Multinomial logistic regression results for the best process domain model including dataset, coefficient estimates and their 95% confidence intervals (lower and upper), standard errors (SE), and significance (P), where fluvial channel (4) is reference domain class for hillslope (1), swale (2), and seepage channel (3) domain classes. ...................... 116 Table 4-9. Logistic regression models for predicting the threshold (‖) to swales, seepage channels and fluvial channels including dependent variable, coefficient estimates and their 95% confidence intervals (lower and upper), standard errors (SE), significance (P), odds ratio and their 95% confidence intervals. ........................................................................ 117 Table 4-10. Ordinal logistic regression predictors for the DTW index subset (n=441), selected by stepwise selection procedure repeated for 10 cross-validation datasets, including predictor, mean Chi-square score (score), and frequency of selection.............................. 120 Table 4-11. Ordinal logistic regression results for the DTW index subset (n=441), including dependent variable, coefficient estimates and their 95% confidence intervals (lower and upper), standard errors (SE), significance (P), odds ratio and their 95% confidence intervals.......................................................................................................................... 121 Table 4-12. Predictors selected by stepwise selection procedure repeated for 10 cross-validation datasets for six data subsets including HRU parent material and number, sample size and predictors with selection frequency ≥ 50 %. .................................................................... 125  Table A- 1. Comparison of dependent variable types, appropriate multivariate methods and solution forms. ............................................................................................................... 157 Table A- 3. Relations between probability, odds, and logit. .................................................... 158 Table A- 4. Example of the conversion of an ordinal variable with four classes into three non- redundant cumulative variables....................................................................................... 159  ix  List of Figures Figure 2-1.  Map of the fish sampling locations by project including monitoring watersheds (model training data set), Hinton region inventory (model testing data set) and Berland region inventory (model testing data set). .........................................................................10 Figure 2-2. Predicted probability of occurrence by habitat variables (a-d) for brook trout (Salvelinus fontinalis) (grey dotted), bull trout (Salvelinus confluentus) (black dashed), rainbow trout (Oncorhynchus mykiss) (black solid) and all fishes (grey thick dash-dot).  (a) drainage area, (b) mean basin slope, (c) elevation (m) and (d) reach slope.  The range of the x-axis represents the approximate observed range of each variable.  For predictions, all values for other variables were standardized to mean values. ............................................21 Figure 2-3. Maps of predicted occurrence probability for bull trout (Salvelinus confluentus), rainbow trout (Oncorhynchus mykiss) and brook trout (Salvelinus fontinalis) within representative areas of the Lower Foothills, Upper Foothills and Subalpine natural subregions. .......................................................................................................................23 Figure 2-4. Correlations between mean basin slope and land-use variables (a-c) using bin averaged values of each land-use variable with bins distinguished by a 5% step in mean basin slope.  (a) Percent of basin harvested, (b) road density and (c) percentage of sites with downstream barrier present. ..............................................................................................24 Figure 3-1 (a) Inset map of study area location. (b) Inset map of ice sheet flow direction down North Saskatchewan River valley with resistant conglomerate of the Cadomin Formation and Upper Paleozoic limestone, plus extent rectangle of study area.  (c) Candidate sub- basins of the study area. (d) Faults and resistant formations from scanned geology maps [Crombie and Erdman, 1947; Douglas, 1956] that were geo-referenced and digitized. .....39 Figure 3-2 (a) Map of candidate sub-basins and surficial materials with five research basins numbered including field surveyed longitudinal profiles with process domain boundaries. (b) Map of candidate sub-basins and SL/k index values by stream reaches for the overdrawn drainage network with five research basins numbered.  SL/k index values in the 3-6 class represent moderately over-steepened reaches, and values in the >6 class represent severely over-steepened reaches. ......................................................................................41 x  Figure 3-3. Hack profiles with step SL/k index curves for five sub-basins.  Dotted black vertical line is boundary between hillslope and swale process domains.  Black vertical arrow is channel head location.  Dotted grey vertical line is boundary between colluvial and fluvial process domain.  SL/k = (reach slope × distance to source) / (rise to source / log (distance to source)). .......................................................................................................................47 Figure 3-4. Plot of log-bin averaged drainage area versus slope for the five research basins. The plot shows the mean slope of individual 2 m grid cells for each 0.1 log bin in drainage area. Black lines represent longitudinal profiles of seven transects with field determined domain class assigned. Drainage area and slope values for profiles were obtained from reaches of the overdrawn network.  Dashed vertical lines divide the plot into process domains based on average drainage area at transitions.  Channel head locations correspond to the left-most colluvial channel symbol. .................................................................................................48 Figure 3-5. Plots of cumulative logit vs. predictors: (a) area for all domains, (b) area by transect for swale and lower domains (swale -), (c) reach slope for all domains, and (d) reach slope by transect for swale and lower domains (swale -). ...........................................................57 Figure 3-6. Plots of (a) accuracy and (b) penalty with standard error by random effects type for four models predicting all four process domains.  Random effects type 1 = no random effects, type 2 = random effects with population average predictions (PAP), type 3 = random effects and spatial autocorrelation accounting (SAA) with PAP, type 4 = random effects with transect specific predictions (TSP), type 5 = random effects and SAA with TSP. .................................................................................................................................58 Figure 3-7. Random effect estimates by transect for the four all domains models including: area; area + lithology; area + profile complexity index; and area + profile complexity index + lithology. ..........................................................................................................................61 Figure 3-8. Predicted changes in logit with spatial autocorrelation accounting by source distance by transect for (a) the best all domains model (area + lithology + profile complexity index) and (b) the best 1, 2, 3 domains model. .............................................................................62 Figure 3-9. Plots of accuracy of process domain predictions with standard error by (a) length of planform curvature measurement, and (b) length of profile curvature measurement.  Length of 0 m = “area” only model with no curvature measurement. ............................................64 xi  Figure 3-10. Plots of (a) accuracy and (b) penalty with standard error by random effects type for four models predicting process domains including hillslopes, swales and colluvial channels.  Random effects type 1 = no random effects, type 2 = random effects with population average predictions (PAP), type 3 = random effects and spatial autocorrelation accounting (SAA) with PAP, type 4 = random effects with transect specific predictions (TSP), type 5 = random effects and SAA with TSP. ..........................................................65 Figure 3-11. Maps of Basin 1 with (a) predicted process domains from the best model from the logit preview exercise, and (b) predicted process domains from model identified with the stepwise selection procedure. ............................................................................................68 Figure 4-1. Groundwater flow processes with potential to affect discharge – drainage area relations at process domain transitions operate at: (a) the reach scale where spatial changes in alluvial area down the longitudinal profile can create hyporheic flow [Tonina and Buffington, 2009]; (b) the hillslope scale where the upwelling volume is determined by hillslope configuration [Toth, 1963], reflected by soil moisture status; and (c) where contributions to local flow systems increase with general relief at the basin scale [Toth, 1963]. ...............................................................................................................................75 Figure 4-2. Study area maps including (a) map of Canada, with extent rectangle for (b) west central Alberta, a transition area between Rocky Mountains and boreal plain, with extent rectangle for (c) shaded relief map of study area, including extents of DTW index and moisture regime data, and field sampling locations. ..........................................................77 Figure 4-3. Example photographs of process domain classes. ...................................................88 Figure 4-4. Channel features used to differentiate seepage and fluvial channels. .......................89 Figure 4-5. Landscape patterns from glacial and fluvial erosion processes shown on 1:50,000 scale shaded relief maps with aerial photograph interpreted channels. (a) Athabasca Tower Hill with stream in hanging valley draining northwestward into abandoned glacial meltwater channel showing discharge (l/s) at six locations from June 19, 2008. (b) Wildhay River tributaries flowing southeastward from ridge to fluvial terrace (dashed white line) then to floodplain. ............................................................................................................92 Figure 4-6. Hack profiles with step SL/k index curves for (a) stream flowing northwestward from Athabasca Tower Hill in hanging valley, and (b) Wildhay River tributary flowing southeastward from ridge to intersection with fluvial terrace then to floodplain.  SL/k = (reach slope × distance to source) / (rise to source / log (distance to source)). ...................92 xii  Figure 4-7. Spatial patterns of soil moisture regime versus DTW index at the catchment scale (a) and (c), and the reach scale (b) and (d)..............................................................................95 Figure 4-8. Within group error by number of clusters for till dominated basins and basins with lacustrine/deltaic deposits. .............................................................................................. 101 Figure 4-9. Box plots of cluster by mean basin slope and cluster by percent of wetlands for till dominated basins and basins with lacustrine/deltaic deposits. ......................................... 101 Figure 4-10. Map of study area with basin types. .................................................................... 102 Figure 4-11. Flowpath diversion at railroad – major stream intersections for (a) the preliminary network and (b) DTW index flowpaths, and flowpath diversion at road – headwater stream intersections for (c) the preliminary network and (d) the DTW index flowpaths.   For reference purposes, the provincial hydrography layer (Prov. streams) was included to show watercourse locations determined by aerial photograph interpretation. ............................ 104 Figure 4-12. Logit plots of continuous class predictors. .......................................................... 108 Figure 4-13. Logit plots of interaction terms including drainage area and each continuous class predictor. ........................................................................................................................ 109 Figure 4-14. Evaluation plots of (a) mean accuracy with standard error, and (b) mean penalty with standard error for ordinal logistic models with increasing number of parameters (n=617). ......................................................................................................................... 111 Figure 4-15. Evaluation plots of (a) mean accuracy with standard error, and (b) mean penalty with standard error for multinomial logistic models with increasing number of parameters for the complete dataset (n=617). .................................................................................... 112 Figure 4-16. Average accuracy from the 10 model training datasets for each of six model types.  ....................................................................................................................................... 113 Figure 4-17. Predicted probability of occurrence by predictor for cumulative logits for swale, seepage channels, and fluvial channels: (a) drainage area, (b) reach slope, (c) SL/k index, (d) moisture regime, and (e) is mean basin slope.  The range of the x axis represents the observed range for each predictor.  For predictions, all values for other variables were standardized to median values. ....................................................................................... 119 Figure 4-18. Stream network maps of (a) the preliminary network, (b) the predicted process domains including swales, seepage channels and fluvial channels, and (c) the predicted occurrence probability for channels. ............................................................................... 122 xiii  Figure 4-19. Stream network map with domains with a swale predicted where the probability of a channel is less than the 0.32 cut-off and channel where the probability of a channel is ≥ 0.32.  Segments (a), (b) and (c) are predicted to be in the swale domain with probability of being a channel as 0.31, 0.26, and 0.16 respectively........................................................ 124 Figure 4-20. Evaluation plots of accuracy for binary logistic models for channel presence by basin type.  Reference line represents regional model accuracy of 78.1 %. ...................... 125  xiv  Acknowledgements I wish to thank the people and agencies that supported my PhD education.  Marwan Hassan and Dan Moore, my supervisors, provided constant encouragement, support, and guidance throughout my learning process.  Dan Miller, contributed computer programming, encouragement, advice and a thorough review of Chapter 3.  Lorne Rothman‟s statistical mentorship was invaluable.  Uldis Silins and Lori Daniels, my committee members, provided help at important junctures.  Lee Benda made a visit to study area and provided insightful discussion about stream morphology in the complex glaciated landscape.  John Buffington, Brian Klinkenberg, and John Richardson served as examiners during the defense of the dissertation and provided valuable comments that were considered in a revision.  The detailed review provided by John Buffington was especially helpful. Mark Schoenberger, Chris Spytz, and Keith Campbell, my forestry colleagues in Hinton, shared their field knowledge on factors influencing headwater stream occurrence.  Eric Leinberger prepared many of the figures in Chapter 3.  Kevin Andras produced the preliminary stream networks.  Julie Duval and Katie Yalte provided GIS support.  Field assistance for this endeavor and the related sediment budget/wood budget project was provided by Catherine Brown, Mike Climie, Scott Fiddes, Steve Haslett, Jessus Karst, Donald Lougheed, Shireen Ouellet, Eric Rehwald, and Lyndon Rempel. Jim LeLacheur and Rick Bonar were instrumental in establishing support through the Foothills Research Institute, whose funding partners included the Alberta Conservation Association, Alberta Newsprint Company Ltd., Alberta Sustainable Resource Development, ConocoPhillips, Encana, the Forest Resources Improvement Program, Jasper National Park, Hinton Wood Products – a Division of West Fraser Mills Ltd., Spray Lake Sawmills, Suncor Energy, Sundre Forest Products – a division of West Fraser Mills Ltd., Talisman Energy, and Trout Unlimited Canada.  Natural Sciences and Engineering Research Council of Canada (NSERC) provided funding through an Industrial Postgraduate Scholarship (IPS). Kristina McCleary, my wife, encouraged and supported this endeavor from start to finish. 1  Chapter 1 Introduction 1.1 Overview Theories on the connections between hydrogeomorphic processes and the evolution of drainage basin topography have continually developed since an original model was proposed by Horton [1945].  According to Horton, streams develop from the point where overland flow first exceeds a critical length determined by surface slope, runoff intensity, infiltration capacity, and shear strength, thus creating a process-based division of the landscape into hillslopes and channels.  Horton suggested that valley and stream development occur together, as evidenced by the development of concave topography.  Subsequent field studies of the hydrogeomorphic processes near the hillslope-stream transition have revealed an upper network of unchannelled valleys or swales, and introduced the importance of episodic events, some of which are driven by water moving through subsurface pathways rather than overland flow [see reviews by Dietrich and Dunne, 1993; Benda et al., 2005; Hassan et al., 2005].  In arid New Mexico, the threshold location between hillslopes and channels is not stationary; rather, in response to ground cover modifications, overland flow can form a discontinuous channel that appears as a sequence of headcuts interspersed with vegetated ground [Leopold and Miller, 1956].   Sections of a discontinuous channel coalesce during extreme rainfall events as individual headcuts advance.  Retreat of the channel head in such environments is a much slower process.  In forested hillslopes in Japan, convergence of sub-surface storm water drives shallow debris flows that originate in elongated spoon-shaped features at valley heads, also called zero-order basins [Tsukamoto et al., 1982].  In the Oregon Coast Range, the recurrence interval of landslides was 6,000, 1,500 and 750 years from zero-order, 1 st  order and 2 nd  order basins respectively [Benda and Dunne, 1987].   Streams can be distinguished into colluvial or fluvial channels based on whether they develop in material conveyed from adjacent hillslopes, or in stream-deposited material [Hassan et al., 2005].  Colluvial channels display characteristics including non-alluvial boundaries, complex flow structures, and periodic evacuation of material by landsliding [Benda, 1990; Rice and Church, 1996].  A colluvial channel transitions to a fluvial channel where depth or slope becomes sufficient to exceed the shear stress threshold required to erode stream banks and create regular bed forms.  Thus, evidence gathered during field research of hydrogeomorphic processes indicates four process domains, or regions characterized by distinct erosional rates and mechanisms, as hillslopes, swales, colluvial channels, and fluvial channels. 2  When hydrogeomorphic processes act continuously on a landscape over a period of time beginning with the exposure of a new surface, each process domain leaves a unique erosion signature that can be detected using slope-area plots [Montgomery and Foufoula-Georgiou, 1993].  There are three limitations that complicate the use of this approach for delineating all four process domains across landscapes. First, slope-area plots predicted from digital elevation models (DEMs) can be used to estimate hillslope – valley and colluvial channel – fluvial channel transitions.  However, the dynamic nature of the channel head requires acquisition of field data to establish channel network extent.  Second, specific thresholds will vary between landscapes and watersheds; thus, thresholds can only be developed for individual basins.  Third, previous land-sculpting processes, including glaciations, may control channel network architecture.  In fact, GIS-based slope-area plots do not discriminate field-based process domains where glaciations have altered the surface; rather, contemporary hillslope and fluvial processes are superimposed upon a landscape characterized by a hierarchy of glacial troughs [Brardinoni and Hassan, 2006].  In such regions, the complex evolution of a landscape reflects the climatic history and the sequence of dominant erosional processes that have acted over time. The more a channel deviates from the idealized graded (concave-up) profile, the greater the number of partitions required during the application of downstream scaling relations.  This inherent complexity precludes the application of theoretically-based remote predictions of channel types in glaciated mountain basins. Procedures to discriminate field-based process domains across complex regions have not been well developed; consequently, models of channel initiation remain limited to the basin scale [e.g., Istanbulluoglu et al., 2002; Hancock and Evans, 2006].  Local discharge of water and local slope determine whether diffusive or concentrative erosional forces act upon a landscape and to what degree [Smith and Bretherton, 1972].  Discharge and drainage area are typically related through a power function [Leopold et al., 1964]; therefore, models of erosional processes usually substitute drainage area for local discharge of water.  In heterogeneous landscapes where surface discharge is difficult to explain using models limited to drainage area and local slope, additional factors including wetland extent can improve predictions [Lindsay et al., 2004].  Models that consider such factors are required to improve the prediction of the spatial distribution of erosion processes in complex landscapes. 3  There are two general approaches for explaining landscape organization.  The first relies exclusively upon remotely sensed information such as aerial photographs and digital elevation models (DEMs) to detect convergent terrain and model earth surface processes [e.g., Tarboton et al., 1991].  The second approach recognizes a more complex landscape whereby a sequence of process domains including swales, colluvial channels, and fluvial channels occupy linear depressions with transitions between these domains that change over time [Dietrich and Dunne, 1993].  The later approach emphasizes GIS-derived topographic indices to explain patterns in field-based process domains [e.g., Massong and Montgomery, 2000; Brardinoni and Hassan, 2006].  Both approaches are suitable for creating map-based models of erosion processes, with the second including validation measures from field-based determinations of process domain. The second approach – a reliance upon field-based process domains – is used in all three chapters of this dissertation. 1.2 Aims, objectives, and organization Headward extension of individual process domains can result in incision of the land surface and mobilization of large volumes of sediment for subsequent deposition in downstream areas [Leopold and Miller, 1956; Benda, 1990].  Related erosion processes can remove productive lands, alter water tables to change the productive capacity of agricultural or forested land, and damage infrastructure.  The sediment evacuation process can alter downstream water quality, channel morphology, aquatic habitats [MacDonald et al., 1990], and bury downstream lands and infrastructure [Schuster, 1996].  Land-use activities and climate change that alter surface runoff patterns can cause the transition locations between process domains to shift, triggering the release of large volumes of stored sediment [Tucker and Slingerland, 1997]. Thus, at the regional scale, it is important to know the contemporary spatial distribution of process domains, and to understand how the climatic history has affected the sequence of dominant erosional processes that have acted over time.  This knowledge can support predictions of how process domains will shift in response to land-use and climate change. Despite the research that has been completed, there is no model that can accurately predict the spatial distribution of erosion processes in complex terrain.  Models can predict a limited set of process domains for regions with uniform bedrock resistance [e.g., Perron et al., 2009]; however, the extent of such regions is limited across most continents.  In fact, landscape variability is the norm in most areas.  Models that strive to explain the variation, rather than 4  define central tendencies, are required to understand how specific sites will respond to land-use and climate change.  Process domain models have been developed at a limited spatial scale for individual basins [Brardinoni and Hassan, 2006]; yet, means to extrapolate such models across landscapes have not been developed.  Multiple factors influence channel development, including bedrock resistance, sediment supply, basin-scale relief, local relief, and climate.  None-the-less, statistical models that incorporate such a suite of predictors are not well developed.   Land- sculpting processes, including glaciations, and shifts in climate that alter ground-cover and runoff, affect channel head locations.  As a result of these temporal shifts in process domains, deterministic models which rely entirely upon remotely sensed data are at a disadvantage compared to more empirical models that utilize field measures.  Thus, the aims of this thesis were: to explain the organization of a complex glaciated foothills landscape using remotely sensed information; to present the relations between process domain occurrence and a set of predictors within a series of spatial models; and to use an automated system to extrapolate the models across landscapes.  The foothills presented a challenging area for this exercise: sedimentary rock formations in the transition between Rocky Mountains and boreal plain have varying levels of resistance; repeated glaciations by cordilleran and continental ice sheets often mask the previous fluvial erosion signature; and forest, petroleum, and coal resource development affect aquatic resources to varying degrees across the region. The availability of high-resolution LiDAR data within the foothills region increased over the duration of this dissertation.  At the start of the research in 2005, LiDAR data were not obtainable.  Thus, the models in Chapter 2 were developed with a 25 DEM.  In 2006, LiDAR data became available for limited areas.  A LiDAR-derived DEM was obtained for a 60 km 2  area with relatively pristine watersheds.  Field work was completed in 2007.  Models from this work are presented in Chapter 3.  In 2008, the Province of Alberta acquired LiDAR data for large areas of the Foothills and provided the data for a study area covering more than 10,000 km 2 .  Fieldwork was completed in 2008 and 2009.  The research completed with these data is presented in Chapter 4. In Chapter 2, the factors that influence erosion process thresholds are examined indirectly through regional models of spatial distribution of fish occurrence.  Hydrogeomorphic and fish habitat models share two common dimensions.  Local discharge of water and local slope determine whether diffusive or concentrative erosional forces act upon a landscape and to 5  what degree [Smith and Bretherton, 1972].  Thus, these two factors also determine whether the specific shear strength thresholds will be exceeded at the individual process domain transitions. Models for fish occurrence can be developed by choosing variables that represent three habitat dimensions including stream size, energy, and climate [Bozek and Hubert, 1992].  In one Washington state ecoregion, the minimum mean annual flow required to generate food to sustain stream dwelling fish is from 0.01 to 0.10 m 3 s -1 ; and minimum estimated bankfull flows required to generate shear stress required to create habitat for fish is from 0.16 – 1.09 m3s-1 [Trotter, 2000].  Thus, surface discharge and local slope are common elements determining the thresholds between fish presence-absence and individual process domains.  However, drainage area–surface discharge scaling relations are not constant across the foothills region.  Runoff increases in response to higher precipitation and lower evaporation with elevational increases [Veldman, 1997].  In low relief basins where vertical groundwater flow can dominate over lateral flow, larger groundwater recharge is expected [Devito et al., 2005a].  Furthermore, in the glaciated foothills, substrate depth increases as relief subsides [Stelfox, 1981].  Deep substrates are more conducive to the formation of groundwater flow systems of intermediate and regional size.  Such systems are less connected to headwater stream networks than local topographically driven recharge-discharge groundwater systems that are common in intermediate and high relief basins [Toth, 1963; Devito et al., 2005a].  Thus, basin relief is expected to influence drainage area–surface discharge scaling relations at the landscape scale, thereby also influencing drainage area thresholds for individual fish species.  If this factor improves predictions of fish occurrence, it is expected that similar considerations will improve predictions of the spatial occurrence of process domain transitions.  The specific objectives of Chapter 2 were: (1) to map distributions for individual fish species and for all fishes regardless of species; (2) to determine the most important variables for explaining occurrence; (3) to compare distributions of native and non- native fish; and (4) to identify relations between land-use and occurrence of individual species. In Chapter 3, the thresholds between hillslope, swale, colluvial channel, and fluvial channel domains were modelled using a LiDAR DEM within small basins.  The research was conducted in the foothills near Nordegg, Alberta, within Dutch Creek and adjacent watersheds. The research basins, which contain several bedrock formations with varying levels of resistance and have been subject to repeated glaciations, were characterized by complex longitudinal profiles.  Area-slope plots, a proven approach for identifying domain thresholds, provided information on the average drainage area at each threshold.  However, explaining the substantial 6  variation between individual basins was the main challenge of the chapter.  Metrics used at the continental scale to detect knickpoints along a single watercourse [e.g., Perez-Pena et al., 2008] were adapted to produce a knickpoints map, which covered all branches of a drainage network. When coupled with a bedrock geology map, this information explained some of the threshold variation.  Using multiple regression, the knickpoint index, geology, and other topography- based predictors, were used to explain profile variation.  The model was extrapolated across adjacent basins using an automated approach to assign process domains. In Chapter 4, the thresholds between hillslope, swale, seepage channels, and fluvial channel domains were modelled across a 10,000 km 2  glaciated region of the Foothills and Interior Plain.  The models were applied at the landscape scale, similar to Chapter 2, adapting techniques developed with high-resolution DEM data from Chapter 3.  Hydrologic classification can promote the transfer of information on controls of hydrologic systems from the regional scale up to the national and global scale [Buttle et al., 2009].   To facilitate the transfer of information to other regions and across broader scales, Chapter 4 begins with assigning a classification to catchment units.  Soil moisture is a major determinant of wetland status and an important factor affecting erosion.  It varies with topographic position, soil type and aspect; hence, it is difficult to quantify based on DEM-derived topographic descriptors.  Thus, a soil moisture index from an ecological classification was used within the catchment classification.  It was also combined with predictors used in Chapter 3, within a multiple regression model.  The relief in the region was generally below the threshold for landsliding; hillslope erosion was largely limited to seepage processes.  Thus the colluvial channel domain was substituted with the seepage channel domain.  Seepage and colluvial channels share similar physical characteristics.  However, in seepage channels the periodic evacuation of accumulated material occurs through the advance of headcuts rather than landslides.  The automated approach from Chapter 3 was used to extrapolate the process domain model across a much larger region. 7  Chapter 2 Predictive modelling and spatial mapping of fish distributions in small streams of the Canadian Rocky Mountain Foothills 2.1 Introduction Quantifying fish distributions in streams and extrapolating these occurrence patterns beyond sampled areas can provide knowledge to support fish conservation planning and guide detailed experimental investigations [Rosenfeld, 2003].  Impediments to extrapolation of occurrence patterns include identifying key parameters that can be accurately measured with remote methods and addressing scaling issues related to inherent correlations between specific geomorphic descriptors (e.g., negative correlation between drainage area and local slope). Approaches to modelling fish distributions differ between regions and land-use types. Differences in terrain limit model transferability.  Models of fish occurrence used in headwater streams of the mountainous regions of western North America include descriptors of stream size, reach slope and disruptions in stream connectivity created by obstructions to fish migration [Fransen et al., 2006].  Geomorphic characteristics that dictate thermal regime also influence habitat selection by salmonids [Baxter and Hauer, 2000].  Fish habitat and distribution in the Great Plains of central North America have been impacted by agriculture, reservoir construction and urbanization and as a result, models of fish distribution include a greater number of variables to describe environmental gradients and land-use [Oakes et al., 2005]. Impacts to fish habitat from riparian logging are emphasized in some regions, however these effects are not relevant in the Rocky Mountain Foothills (foothills) where riparian buffer strips have been used since the onset of industrial forest harvest in 1956 [Bott et al., 2003]. Likely impacts include sedimentation from logging roads [Spillios, 1999], bank erosion, and habitat loss due to obstructions at culverts [Scrimgeour et al., 2003].  However, relations between land-use, habitat impacts and fish occurrence have not been studied extensively in the continental rain and snow climate of the foothills. Habitat studies have historically focused on large streams that provide water, power and habitat for sport fish; less is known about small streams that represent most of the watercourses on the landscape [Moore and Richardson, 2003].  In comparison to large streams where fish presence can generally be inferred, only a portion of the small streams on the landscape provide habitat suitable for fish.   In small streams, land-use has a greater possibility of altering fish 8  habitat features created by streamside vegetation and large woody debris [Hassan et al., 2005]. Small streams can also be suitable for low-cost culverts at road crossings, however where fish are present, these structures may obstruct fish passage [Wolter and Arlinghaus, 2003].  Fish presence at such locations should be ruled out by conducting electrofishing inventories during different years and seasons and by considering characteristics of the channel network.  In the foothills, most small streams are inaccessible and located in a region of growing industrial activity.  Environmental assessments for proposed projects ideally include electrofishing inventories.  However, information on fish-bearing status may be required on short notice or outside of the summer season during which conditions are suitable for electro-fishing.  The sensitivities of small streams to a variety of land-use impacts and the increasing rate of industrial development warrants a remote sensing approach to provide fish occurrence information. The main objective of this research was to develop an automated method to predict fish occurrence in small foothills streams.  This knowledge is important to support planning of road locations and design of stream crossing structures.  Development of this approach is timely in Alberta because of the rapid expansion of the petroleum exploration and extraction.  Fish occurrence information is also necessary when designing processed-based geomorphic and aquatic studies.  The specific objectives of the research are: (1) to map distributions for individual fish species and for all fishes regardless of species; (2) to determine the most important variables for explaining occurrence; (3) to compare distributions of native and non- native fish; and (4) to identify relations between land-use and occurrence of individual species. I sought to address a number of common pitfalls when predicting distributions of organisms.  First, the accuracy of variables measured with a Geographic Information System (GIS) is often not validated with field-measured data.  Second, many researchers report a model‟s predictive accuracy, however this measure is systematically affected by the frequency of occurrence of the target organism [Olden et al., 2002].  Third, few models are validated with an independent data set [Olden et al., 2002].  Fourth, it is difficult to design a sampling program to capture multiple environmental gradients when land-use variables are included.  Land-use activities are not ubiquitous across the landscape [Van Sickle et al., 2004]; rather, their locations are based on criteria that may be similar to those of native fishes adapted to specific niches.  In the foothills, harvesting has always been limited to terrain with less than 40% slope due to the limitations of logging equipment [Bott et al., 2003].  Exclusionary distributions of land-use and 9  fish occurrence create a layer of complexity that confounds the interpretation of model outputs. Finding solutions to these four difficulties will advance the use of GIS-based resource selection function models. 2.2 Methods 2.2.1 Study area The study area covers a 12,000 km 2  portion of the foothills in west-central Alberta (Figure 2-1).  The foothills, underlain with sandstone and siltstone, represent a transition zone between the interior boreal plains to the east and the carbonate-dominated front ranges of the Rocky Mountains to the west.  The combination of highly erodible bedrock and extensive glaciation has created a rolling topography.  Elevations range from 1,000 to 1,600 m. Lodgepole pine and spruce forests cover uplands and either black spruce / tamarack forests or shrubs dominate wet areas.  Approximately 1% of the land base (120 km 2 )  is harvested annually and the first harvest rotation will be completed by 2040. 10   Figure 2-1.  Map of the fish sampling locations by project including monitoring watersheds (model training data set), Hinton region inventory (model testing data set) and Berland region inventory (model testing data set). 11  2.2.2 Fish species This study focused on bull trout (Salvelinus confluentus), rainbow trout, (Oncorhynchus mykiss), and introduced brook trout (Salvelinus fontinalis).  Bull trout are designated as a species of special concern in Alberta.  The Athabasca River watershed (Figure 2-1) is one of three basins east of the continental divide that support native rainbow trout [Taylor et al., 2006]. Both non-native and native rainbow trout inhabit the study area, but small streams provide important habitat for genetically pure native populations [Taylor et al., 2006].  Brook trout were first introduced in Alberta in 1922 and widespread introductions in the 1960s included headwater lakes in the Athabasca River watershed [Donald et al., 1980].  Bull trout have a very high sensitivity to angling, road-related sediment, road-related migration barriers and increases in water temperature [Rieman and McIntyre, 1993], while rainbow trout have a high sensitivity [Ford et al., 1995b] and brook trout have a moderate to high sensitivity to these impacts [Ford et al., 1995a].  Brook trout present a threat to the long-term conservation of bull trout through hybridization [Leary et al., 1993] and juvenile competition [Selong et al., 2001].  Knowledge of spatial distributions of these native and non-native fishes may assist conservation efforts. I used fish occurrence data from three different projects (Figure 2-1).  I collected information within the monitoring watersheds (MW) specifically for this study.  I obtained data from a Hinton region inventory (HR) and a Berland River region inventory (BR) for the purposes of model validation.  In the HR inventory, forest planners selected sites where fish presence/absence information was required to conserve fish passage at existing and proposed road/stream crossings.  The BR inventory, completed in 1995 and 1997, was a reconnaissance survey of a largely pristine area scheduled for forest harvest. The study area includes the Subalpine, Upper Foothills and Lower Foothills natural sub- regions [Beckingham et al., 1996].  I selected 15 watersheds for detailed analysis that represent the range of biophysical and land-use characteristics for the study area (Figure 2-1).  For the MW and HR inventories, technicians used a Smith-Root Model 12A backpack electro-fisher to sample 300 m long reaches.  For the BR inventory, a Coffelt Mark 10 was used and reaches in 1 st  and 2 nd  order streams were reduced to 150 m.  The reduction in stream length was likely made because sampling time tends to increase because channels of this size tend to be overgrown and densely vegetated.   This reduction may have resulted in a higher minimal detectable population for these reaches as described in Section 2.2.8.  Electrofishing was deferred when elevated flow or turbidity levels could have affected capture.  A single pass 12  survey with no block nets was used to determine fish presence/absence.  I randomly selected sample sites from the various stream order and stream slope combinations within each of the 15 monitoring watersheds and these sites were sampled between 1999 and 2002.  Additional fish presence/absence data within the 15 watersheds was obtained for a small number of sites that were monitored between 1996 and 2002. 2.2.3 An automated procedure to measure basin and reach characteristics To map predicted probabilities of target fish species in the study area, I needed a system to extrapolate findings from specific sample points across the study area.  This approach required a set of predictor variables that could be calculated through GIS analyses of existing data.  An automated procedure to measure these variables included two stages.  Before completing these stages, I prepared a single line stream network and hydrologically corrected the digital elevation model (DEM) [Duhaime, 2003]. In the first stage, reach break creation, the stream network was divided into a series of distinct sections based on topologic (stream network) and topographic (slope) criteria.  Reach breaks were placed at all confluences and topological features such as waterfalls and lake boundaries.  I then used a two-step process for the creation of topologic breaks.  In the first step, preliminary reach breaks were placed every 100 m in an upstream direction from each topologic reach break.  In the second step, preliminary topologic reaches were amalgamated wherever the change in slope from a downstream segment to the next segment upstream failed to exceed a given tolerance.  To achieve an average reach length of 300 m for all streams in the study area, the tolerances changed with downstream reach slope as follows: slope ≤ 1%, tolerance = 0.5%; slope 1 ≤ 2%, tolerance = 1.0%; slope 3 ≤ 4%, tolerance = 1.5%; slope 4 ≤ 6%, tolerance = 2.0%; slope 6 ≤ 10%, tolerance = 3.0%; slope 10 ≤ 20%, tolerance = 5.0%; slope 20 ≤ 40%, tolerance = 10%; and slope ≥ 40%, tolerance = 20%.  The process proceeded in an upstream direction from each stream confluence.  Reach breaks were retained at topologic features and stream confluences. The second stage consisted of the creation of watershed boundaries representing a drainage area polygon originating at the downstream end of each reach.  The polygons were created using a nested approach to ensure consistent boundaries between adjacent watersheds 13  and to facilitate editing.  The drainage area boundary for each reach was saved to permit calculation of basin scale variables for each reach. 2.2.4 Accuracy of stream slope measures generated from a DEM To evaluate the accuracy of measures of stream slope derived from a 25 m DEM, I measured true stream slope with a total station (Lecia Model TPS700) and prism pole at 23 stream reaches.  Surveys extended for a distance equivalent to 50 – 100 channel widths and points were captured at all bed feature and meander transitions.   I found a meandering channel pattern in most reaches that extended into headwater positions.  Streambanks were well vegetated and root structures promoted meander development.  For each survey, line length was the sum of the horizontal distances between successive points.  Straight-line distance was calculated using the start and end points.  Using GIS (ArcMap Version 9.2), I isolated the mapped stream segments that corresponded to our field surveyed reaches.  Then I calculated the change in line length by dividing the total length of the mapped stream reach by the total length of the surveyed stream reach [McMaster, 1986].  Stream reaches with a drainage area greater than 20 km 2  had small changes in line length, indicating reasonable discernment of the channel pattern during aerial photograph interpretation.  Below this threshold the change in line length averaged 78 %, indicating greater generalization and poor representation of channel meanders. Therefore, this channel length correction factor was applied to all reaches with a drainage area < 20 km 2  (corrected slope = slope × 0.78). 2.2.5 Deriving predictive models of fish occurrence The strategy I used to explore relations between fish presence/absence and habitat descriptors has been used in related studies [e.g., Ripley et al., 2005] and is consistent with the information-theoretic approach to model development and selection.  This approach differs from traditional hypothesis testing in three ways [Burnham and Anderson, 2002].  First, the researchers identify a set of candidate models from patterns described in the literature and their observations.  Second, the best model is identified using the Akaike‟s information criterion that considers both parsimony and the fit of each candidate model to the phenomenon or process of interest.  Third, this approach does not use hypothesis tests or p-values. I calculated the small sample version of Akaike‟s information criterion (AIC  c) and included a preliminary test for overdispersion [Burnham and Anderson, 2002].  Following the calculation of AIC c, I estimated an evidence ratio (i) and Akaike weights (wi) to help select the 14  most parsimonious model and further provide rankings of model importance and uncertainty. The model with the lowest i is the best-fit model out of the candidate set of models and wi provides an approximate representation for the probability of a particular model fitting the data compared to the candidate set of models.  To accomplish our objectives, two separate model selection exercises were undertaken. In the first exercise, the candidate models were limited to those containing predictor variables that describe fish habitat at each stream reach.  These variables, which I call ecological variables, represented key habitat dimensions including size, energy and climate [Bozek and Hubert, 1992] and were measured from readily available digital information sources.  The best model developed from the ecological variables provided distribution information exclusive of land-use and was well suited to extrapolation.  The second set of models included combinations of the best ecological model and three land-use variables.  The land-use variables were intended to represent management activities that could have affected the productive capacity of fish habitat in small streams. I identified six ecological variables and three land-use variables (Table 2-1).  All candidate ecological models included area because of its importance in most fish occurrence models [Rosenfeld, 2003].  The complete set of candidate ecological models included all combinations of area and non-correlated variables.  The set of candidate ecological/land-use models included the variables from the best ecological model with all combinations of the three land-use variables. 15  Table 2-1.Reach-scale variables grouped by category. Variable code Predictor variable, measurement methods and fish habitat considerations Ecological parameters Area Drainage area above downstream end of reach (km 2 ).  Common variable in fish occurrence modelling [Rosenfeld, 2003]. Basin slope Mean basin slope (%). Calculated as the mean slope percentage from all DEM cells in the catchment of each reach.  Common variable in fish occurrence modelling.  Related to annual sediment yield in region [McPherson, 1975]. Basin elevation The average elevation (m) from all DEM cells in the catchment of each reach. Ecosystem productivity decreases and annual precipitation increases with elevation within study area [Beckingham et al., 1996]. Percent wetlands Percent of drainage area occupied by wetlands on forest inventory map.  Muskegs are common in region and typically form along lower side slopes and depressional areas [Dumanski et al., 1972].  Wetlands include black spruce (Picea mariana) / tamarack (Larix laricina) and shrub / herb communities.  I observed low fish occurrence and organic substrate within streams that drain wetlands. Reach elevation Elevation at downstream end of reach (m). Reach slope Common variable in fish occurrence modelling [Rosenfeld, 2003]. Land-use parameters  Harvest  Percent of basin harvested since 1956 measured from timber inventory maps.  Percent harvest was associated with a lower probability of bull trout occurrence within region [Ripley et al., 2005]. Road density  Length of permanent roads within each basin (km/km 2 ) measured on road maps.  Impact water quality and peak flow [MacDonald et al., 1990].  Related to decreased bull trout redd numbers [Baxter et al., 1999]. Downstream barrier  All reaches with a potential barrier located in downstream areas were coded yes (1) and others coded to no (0).  Culverts that retain an uninterrupted natural bottom conserve fish migration and all others may present a barrier depending on hang height, water velocity, fish species and life stage [Jackson, 2003]. Therefore, barriers included all culverts, except those with substrate retained within the culvert and without any hanging outfall.  Fish migration barriers at road crossings represent challenge for fish population conservation in some foothills streams of Alberta [Scrimgeour et al., 2003]. Note: Ecological variables were measured at all reaches across the study area.  Land-use variables were measured for all reaches within the monitoring watersheds. Bull trout, Salvelinus confluentus. 16  Prior to model development, co-linearity between all continuous ecological and land-use variables was assessed using Pearson correlation tests.  To avoid co-linearity, I excluded one of the variables when correlations were greater than |0.6|.  I used a data-splitting type of cross validation [Olden et al., 2002], with 80% of the MW samples randomly selected for model- training and the remaining used for model testing.  A non-linear quadratic relation was expected for ecological variables due to the adaptation of most organisms to a specific niche.  Prior to model development, I used a two-step process to identify such patterns.  For each predictor variable, I determined the average percent occurrence across 10 evenly sized classes.  Next, I plotted relations for percent occurrence by predictor variable class to identify non-linear relations.  When these patterns were apparent, I included a quadratic term within model structures (i.e. basin slope + basin slope 2 ).  Logistic regression was used to compare fish species present sites (1) with absent sites (0) based on ecological covariates of interest and candidate models tested. 2.2.6 Model interpretation and validation Interpretation of regression coefficients was complicated by different measurement units among variables.  To ensure that the value of regression coefficients corresponded to the magnitude of the effect on probability, all continuous predictors were scaled to standardized scores with mean value = 0 and variance =1 [Klienbaum et al., 1988].  Odds ratios (ratio of probability that event will occur to probability that it will not occur) and associated confidence intervals were also calculated.  When a regression coefficient = 0, its odds ratio = 1.  Therefore, when a confidence interval includes 1, it can‟t be concluded that the variable is associated with a change in odds. I assessed model performance based on receiver operator characteristic (ROC) and cross-validation.  For the best ecological models, I used four datasets (MW training and testing, HR and BR inventories).  Land-use data were not available for the HR and BR inventories. The ROC evaluation included two components.  First, I reported the area under the ROC curve where 0.5 - 0.7 represented  „low‟, 0.7 - 0.9 corresponded to „intermediate‟ and > 0.9 „high‟ model accuracy [Manel et al., 2001].  Second, I established high and medium probability cut- offs.  Logistic regression is well suited to studies on rare organisms because the probability cut- off for establishing presence/absence varies with prevalence.  The high probability cut-off corresponded to the point where the model sensitivity (probability of correctly predicting a positive outcome) and model specificity (probability of correctly predicting a negative outcome) 17  were both maximized.  I set the medium probability cut-off at half the value of the high probability cut-off.  Probabilities less than the medium cut-off were assigned a low probability. Based on the high and medium cutoffs, I reported sensitivity (true positives correctly classified), specificity (true negatives correctly classified) and overall accuracy.  The Canada Fisheries Act provides protection of all fish habitat, including stream reaches with infrequent use.  Therefore correct identification of positive outcomes (sensitivity) is more important than correct identification of negative outcomes (specificity).  For the sensitivity evaluation at both cut-offs, < 70 represented low, 70 - 90 represented medium and values > 90 represented high performance. I included Cohen‟s kappa as a fourth validation statistic because, unlike sensitivity, specificity, and overall accuracy, it provides an assessment of correct predictions outside of chance expectations [Manel et al., 2001].  Kappa values of < 0.4 indicate low, 0.4 – 0.6 medium and 0.6 – 0.8 high model performance. 2.2.7 Extrapolating and mapping occurrence models I used a simple procedure to extrapolate the best AIC models and produce distribution maps by species for the study area.  First, I exported a table from the GIS with unique reach numbers and values for all ecological variables for the 40,000+ reaches within the study area. This table was appended to the table containing the model training data and the combined table was imported into the modelling program (SPSS Version 10.0.5). I ran the best model for each species and used the function to save predicted probability values for each reach.  The predicted probabilities for all reaches were linked back to the reach table in the GIS.  Probability of occurrence was mapped using three categories (high, medium and low) based on the ROC determined cut-off value for each model. 2.2.8 Sources of error Errors in fish occurrence can result where the minimum number of fish present in a reach is less than the minimum detectable population for the capture technique.  Paul and Post (2001) used similar capture techniques (single pass electrofishing, 300 m reach and no block nets) in the foothills and they determined a minimum detectable density of 8 individuals/300 m. These findings warranted a conservative interpretation of logistic regression model predictions that considered additional information below the traditional statistical cut-off for presence/absence.  In recognition of this limitation for determining occurrence, I used three 18  levels of probability of occurrence (low, medium and high as described in section 2.2.6) rather than using a definitive yes/no cut-off.  Furthermore, when interpreting model outputs for a specific reach assigned with a low probability of occurrence, it is advisable to consider Errors in the ecological variables (drainage area, basin slope, reach slope and reach elevation) are related to the resolution of the digital streams layer and accuracy of the DEM.  A digital elevation model is a grid with elevations that represent the earth‟s surface.  Raw DEM data typically contain erroneous features that prevent the downhill movement of water across the surface.  A drainage enforcement function uses GIS smoothing algorithms to correct these erroneous cells by filling depressions and cropping barriers.  The use of a 10 m drainage enforced DEM improves the accuracy of watershed boundary identification over a 30 m DEM lacking drainage enforcement [Clarke and Burnett, 2003].  I used a 25 m DEM with drainage enforced. Land-use variables were measured from a variety of sources.  Percent harvest was measured from digital maps of cut-over areas.  I confirmed this information using digital orthophotographs.  Percent harvest was intended to provide a surrogate indicator for changes in the hydrologic regime.  However, hydrologic changes are not permanent and recovery occurs with forest regrowth.   Road density was calculated using existing digital road information that was confirmed using digital orthophotographs.  Sediment inputs at a road crossing are related to characteristics of the road surface and ditches [Brooks et al., 2006], however such factors were not considered in our models.   Our criteria for identifying potential migration barriers are consistent with policy recommendations for maintaining fish passage under the Canada Fisheries Act and a portion of the potential barriers may be passable to certain species and life stages.  Beaver dams, logjams, headcuts and waterfalls are common natural migration barriers that are not mapped and thus were not included in our models.  Despite these limitations, this information provides a valuable context for interpreting fish occurrence patterns. 2.3 Results 2.3.1 Sampling summary Our study was limited to small streams with drainage areas of < 20 km 2 .  Bankfull width for streams with drainage areas between 15 and 20 km 2  averaged 4.9 ± 1.5 m.  The channel length correction factor (slope×0.78) was applied to all GIS measures of reach slope for these streams.  Based on this drainage area cut-off, sample sizes by project area were: monitoring 19  watersheds (MW) – 284; Hinton region inventory (HR) – 542; Berland region inventory (BR) – 66.  Prevalence of brook trout among these three project areas was: MW – 7%; HR – 12%; BR – 0% (mean for brook trout present areas = 10%).  Prevalence of bull trout was: MW – 25%; HR – 8%; BR – 12% (mean = 15%).  Prevalence of rainbow trout was: MW – 50%; HR – 29%; BR – 32% (mean = 37%).  Prevalence of one or more of any species was: MW – 66%; HR – 43%; BR – 33% (mean = 47%).  In addition to the three target species, 14 additional species were encountered as incidental catches and they were also included in the „all species‟ model along with the three target species. In the Pearson correlation test, mean basin slope was correlated with mean basin elevation and percent wetlands. The variables retained included: area (A); mean basin slope + mean basin slope 2  (B); reach slope + reach slope 2  (S); and elevation (E).  The eight models tested included A, A+B, A+S, A+E, A+B+S, A+B+E, A+S+E and A+B+S+E. Three of the four independent ecological variables shared similar means and standard deviations among the three projects.  Drainage area averaged 7.4 km 2 , (MW 8.0 ± 5.0 km 2 , HR 6.2 ± 5.1 km 2 , BR 7.9 ± 5.8 km 2 ), reach slope averaged 2.7 % (MW 2.8 ± 1.8 %, HR 3.2 ± 2.8 %, BR 2.2 ± 1.6 %) and elevation averaged 1,250 m ASL (MW 1,270 ± 150 m ASL, HR 1,290 ± 190 m ASL, BR 1,190 ± 230 m ASL).  Although basin slope had averaged 13.2 % between the projects, the mean values were from a broader range (MW 18.2 ± 11.4 %, HR 12.6 ± 7.6 % and BR 8.9 ± 8.5 %).  Land-use variables were only measured for the monitoring watersheds project.  Percent harvest averaged 20.1 % (0.0 % - 97.6 %), and road density averaged 0.34 km•km-2 (0.0– 2.54 km•km-2).  I inspected 302 crossings that intersected mapped streams, including those in intermittent watercourses.  Average hang height of the culverts rated as potential barriers was 0.27 m.  I identified downstream barriers at 25% of electro-fished reaches. 2.3.2 Best ecological models Based on the AICc weights (wi), global models for both brook trout and bull trout occurrence had the highest probability as the best model (Table 2-2).  For the rainbow trout and all species models, there was less evidence for one best model.  According to Burnham and Anderson (2002), an evidence ratio (i) of less than 2.0 between models does not show strong support for any one model and they advise to choose based on ecological / management sense. For rainbow trout, I selected the model with the highest wi as the best model.  This model was the most parsimonious and based on p Ĉ, the model also fit the data well (Table 2-2).  For all 20  species, I selected the model with the second highest wi (area + basin slope + elevation) as the best ecological model.  There was some evidence based on p Ĉ that the model with the highest wi did not fit the data, whereas there was no evidence that the second ranked model didn‟t fit the data (Table 2-2). Table 2-2. Logistic regression models of ecological variables predicting the probability of occurrence of brook trout, bull trout, rainbow trout and all species present. Species Model Name K AICc  wi p Ĉ Brook trout Area + Basin Slope + Elev + Reach Slope 7 100.34 0.00 0.45 0.97  Area + Elev + Reach Slope 5 100.78 0.44 0.36 0.63  Area + Basin Slope + Elev 5 102.72 2.38 0.14 0.88 Bull trout Area + Basin Slope + Elev + Reach Slope 7 192.54 0.00 0.62 0.32  Area + Basin Slope + Elev 5 193.57 1.03 0.37 0.10 Rainbow trout Area + Basin slope 4 267.45 0.00 0.55 0.69  Area + Basin Slope + Elev 5 269.43 1.98 0.20 0.76  Area + Basin Slope + Reach Slope 6 269.66 2.21 0.18 0.02 All species Area + Basin slope 4 246.86 0.00 0.35 0.08  Area + Basin Slope + Elev 5 247.84 0.98 0.21 0.21  Area + Basin Slope + Reach Slope 5 248.91 2.06 0.13 0.10  Area + Elev 3 249.27 2.41 0.10 0.02 Note: K, the number of model parameters; AICc, the small sample version of Akaike‟s information criteria (AIC); I, AIC difference between a given model and highest ranked model; wi, Akaike model selection weights; and pĈ, the Hosmer and Lemeshow goodness-of-fit chi- square.  Under average conditions (i.e. basin slope and elevation), rainbow trout were more prevalent in smaller streams than bull trout and brook trout (Figure 2-2).  Streams with a drainage area greater than 3 km 2  had a high probability of rainbow trout (cut-off = 0.51).  The high probability threshold for both brook trout (cut-off = 0.08) and bull trout (cut-off = 0.25) was near 10 km 2 .  Different ecological niches among the three species were apparent from the predicted probabilities for basin slope and reach slope (Figure 2-2).  Maximum probability for brook trout, bull trout and rainbow trout occurred in basins with average slopes of 16-28%, 30- 40% and 12-28%, respectively.  Bull trout probability increased with elevation while brook trout displayed the opposite trend. 21   Figure 2-2. Predicted probability of occurrence by habitat variables (a-d) for brook trout (Salvelinus fontinalis) (grey dotted), bull trout (Salvelinus confluentus) (black dashed), rainbow trout (Oncorhynchus mykiss) (black solid) and all fishes (grey thick dash-dot).  (a) drainage area, (b) mean basin slope, (c) elevation (m) and (d) reach slope.  The range of the x-axis represents the approximate observed range of each variable.  For predictions, all values for other variables were standardized to mean values. Model performance varied among species and datasets (Table 2-3).  Based on average area under the ROC curve, all models had moderate performance at high cut-off values.  Model sensitivity improved with lower cut-off values for all models.  For rainbow trout and all species, sensitivity improved to high with reduced cut-off values.  The use of the lower cut-off value resulted in reduced specificity, overall accuracy and kappa.  Based on kappa values, all models performed moderately at the high cut-off value, except for the brook trout model, which had a low performance rating.  Performance of all models (based on area under ROC curve, sensitivity and kappa) decreased with the HR inventory.  Sixty percent of the incorrectly classified bull trout sites were located within 2km of a stream with a drainage area > 20km 2 ; thus the bull trout stream size preferences suggested by the model could be improved. 22  Table 2-3. Classification analysis for high and medium predicted occurrence probability categories for best ecological models by species. Species Cut-off Dataset n Prevalence AUC PCC (1) PCC (0) PCC Kappa Brook trout High  Training 238 9 0.92 81 82 82 0.36  (0.08) HR 540 13 0.63 29 88 80 0.16   Mean   0.77 55 85 81 0.26  Medium  Training 238 9  95 70 72 0.27  (0.04) HR 540 13  33 82 75 0.12   Mean    64 76 74 0.19 Bull trout High Training 238 25 0.88 77 79 78 0.49  (0.27) Testing 46 22 0.69 60 72 70 0.27   HR 540 8 0.68 40 85 81 0.16   BR 66 12 0.96 63 97 92 0.62   Mean   0.80 60 83 80 0.39  Medium Training 238 25  92 60 68 0.38  (0.13) Testing 46 22  80 58 63 0.26   HR 540 8  53 71 70 0.11   BR 66 12  75 93 91 0.62   Mean    75 71 73 0.34 Rainbow trout High Training 238 50 0.79 71 73 72 0.45  (0.51) Testing 46 52 0.78 67 73 70 0.39   HR 540 29 0.62 47 77 68 0.24   BR 66 32 0.81 57 87 77 0.46   Mean   0.75 61 77 72 0.38  Medium Training 238 50  93 39 66 0.32  (0.25) Testing 46 52  92 32 63 0.24   HR 540 29  86 24 42 0.07   BR 66 32  90 51 64 0.33   Mean    90 36 59 0.24 All species High Training 238 68 0.81 76 74 76 0.48  (0.63) Testing 46 54 0.74 76 57 67 0.34   HR 540 43 0.69 58 75 68 0.34   BR 66 33 0.77 82 68 73 0.45   Mean   0.75 73 69 71 0.40  Medium Training 238 68  99 17 73 0.21  (0.31) Testing 46 54  100 14 61 0.15   HR 540 43  92 15 49 0.07   BR 66 33  91 41 58 0.25   Mean    96 22 60 0.17 Note: High cut-off, the threshold of optimal sensitivity and specificity from ROC curve; Medium cut-off, 0.5 × the high cut-off value; AUC, area under ROC curve. PCC (1), model sensitivity defined as percent of fish occurrence sites correctly classified. PCC (0), model specificity defined as percent of no fish occurrence sites correctly classified.  PCC, overall prediction success.  23  2.3.3 Predicted spatial distributions The different niches for the three target species were also apparent from fish species distribution maps (Figure 2-3).  Brook trout had the greatest predicted extent within the streams originating in the Lower Foothills, while rainbow trout had the greatest predicted extent within streams originating in the Upper Foothills and bull trout within streams originating in the Subalpine natural subregion.  Figure 2-3. Maps of predicted occurrence probability for bull trout (Salvelinus confluentus), rainbow trout (Oncorhynchus mykiss) and brook trout (Salvelinus fontinalis) within representative areas of the Lower Foothills, Upper Foothills and Subalpine natural subregions.  24  2.3.4 Best ecological / land-use models Underlying relations between land-use variables and basin slope present a confounding factor to consider when interpreting correlations between land-use variables and fish occurrence (Figure 2-4).  All three land-use activities including harvest (H), road density (R) and downstream barrier (D) were limited to reaches with basin slopes < 25%.  The eight models tested included the best ecological model (E), E+H, E+R, E+D, E+H+R, E+H+D, E+R+D, E+H+R+D.  Figure 2-4. Correlations between mean basin slope and land-use variables (a-c) using bin averaged values of each land-use variable with bins distinguished by a 5% step in mean basin slope.  (a) Percent of basin harvested, (b) road density and (c) percentage of sites with downstream barrier present.  25  The five top ranked brook trout models shared a i < 2.0, indicating little support for one best model (Table 2-4).  In the highest ranked model, area had the greatest odds ratio, followed by harvest.  The regression parameter values indicate a greater probability of capturing brook trout as harvest increases, but a lower probability as roads increase (Table 2-5).  Percent harvest also had a low standard error.  Road density had a high standard error and also included the value 1 within the odds ratio confidence interval to indicate “no effect” is within the interval. The kappa value indicates medium performance of the best ecological / land-use model, while the ecological model had a lower rating (Table 2-6).  All other performance measures indicate similar trends. Table 2-4. Logistic regression models of ecological (drainage area, reach slope, elevation, and mean basin slope) and land-use variables predicting the probability of occurrence of bull trout, rainbow trout and brook trout. Species Model Name K AICc i wi p Ĉ Brook trout Ecological, harvest, roads 9 93.29 0.00 0.37 0.18  Ecological, harvest 8 93.51 0.22 0.33 0.51  Ecological, harvest, barrier 9 95.37 2.08 0.13 0.79  Global 10 95.45 2.16 0.13 0.19 Bull trout Ecological, barrier 8 188.57 0.00 0.26 0.82  Ecological, harvest, barrier 9 189.20 0.63 0.19 0.91  Ecological, roads, barrier 9 189.79 1.22 0.14 0.99  Ecological, harvest 8 189.89 1.32 0.14 0.21  Global 10 190.36 1.79 0.11 0.93  Ecological, harvest, roads 9 190.74 2.17 0.09 0.38 Rainbow trout Ecological, harvest, roads 6 265.85 0.00 0.25 0.00  Ecological, harvest 5 266.29 0.43 0.20 0.56  Ecological, roads 5 266.78 0.93 0.16 0.34  Ecological 4 267.45 1.60 0.11 0.69  Global 7 267.97 2.12 0.09 0.04  Ecological, roads, barrier 6 268.28 2.43 0.07 0.37  Ecological, harvest, barrier 6 268.34 2.48 0.07 0.36 Note: K, the number of model parameters; AICc, the small sample version of Akaike‟s information criteria (AIC); i, AIC difference between a given model and highest ranked model; wi, Akaike model selection weights; and pĈ, the Hosmer and Lemeshow goodness-of-fit chi- square.  26   Table 2-5. Logistic regression parameter estimates from the best ecological / land-use models for each of three fish species and all species present including standard error of the model coefficient (SE), odds ratio, and associated lower and upper 95% confidence intervals (CI). Model Parameter Estimate SE Odds ratio 95% CI           Lower Upper Brook trout Area 1.763 0.421 5.827 2.554 13.294  Basin slope 0.606 0.807 1.833 0.377 8.920  Basin slope 2  -1.234 0.813 0.291 0.059 1.434  Reach slope -2.522 1.712 0.080 0.003 2.301  Reach slope 2  -2.924 1.533 0.054 0.003 1.083  Elevation -2.088 0.713 0.124 0.031 0.501  Percent harvest 1.419 0.501 4.132 1.548 11.031  Road density -0.949 0.707 0.387 0.097 1.547  Constant -3.557 0.782 0.029 Bull trout Area 0.703 0.260 2.019 1.212 3.363  Basin slope 1.442 0.474 4.230 1.672 10.701  Basin slope 2  -0.326 0.186 0.722 0.502 1.039  Reach slope 0.647 0.380 1.911 0.907 4.025  Reach slope 2  -0.385 0.182 0.680 0.476 0.972  Elevation 0.908 0.296 2.480 1.388 4.432  Barrier downstream -1.295 0.570 0.274 0.090 0.836  Constant -0.906 0.309 0.404 Rainbow trout Area 0.751 0.173 2.120 1.511 2.974  Basin slope 0.434 0.228 1.544 0.988 2.413  Basin slope 2  -0.989 0.192 0.372 0.255 0.542  Percent harvest 0.011 0.006 1.011 0.999 1.023   Constant 0.621 0.241 1.861 Note: To assist interpretation of regression coefficients, all independent variables, except the categorical variable „barrier downstream‟, were scaled to standardized scores. 27  Table 2-6.  Classification analysis for high occurrence category for best ecological / land-use models by species. Species Model Cut-off Dataset n Prevalence AUC PCC (1) PCC (0) PCC Kappa Brook trout Ecological 0.08 Training 238 9 0.922 81 82 82 0.359  Ecological, harvest, roads 0.09 Training 238 9 0.94 86 86 86 0.455 Bull trout Ecological 0.27 Training 238 25 0.883 77 79 78 0.488    Testing 46 22 0.685 60 72 70 0.265    Mean   0.784 68 75 74 0.377  Ecological, barriers 0.31 Training 238 25 0.889 82 84 83 0.595    Testing 46 22 0.732 60 69 67 0.235    Mean   0.811 71 77 75 0.415 Rainbow trout Ecological 0.51 Training 238 50 0.793 71 73 72 0.445    Testing 46 52 0.777 67 73 70 0.392    Mean   0.785 69 73 71 0.419  Ecological, harvest 0.51 Training 238 50 0.795 71 72 72 0.437    Testing 46 52 0.714 67 59 63 0.258    Mean   0.756 69 66 67 0.348 Note: High cut-off, the threshold of optimal sensitivity and specificity from ROC curve; Medium cut-off, 0.5 × the high cut-off value; AUC, area under ROC curve. PCC (1), model sensitivity defined as percent of fish occurrence sites correctly classified. PCC (0), model specificity defined as percent of no fish occurrence sites correctly classified.  PCC, overall prediction success.   The five top ranked bull trout models shared a i < 2.0, indicating little support for one best model (Table 2-4).  The highest ranked model included downstream barriers as the only land-use variable.  The ranking from largest odds ratio in this model was basin slope, elevation then drainage area (Table 2-5).  One parameter in each quadratic term (i.e. basin slope 2  and reach slope) had high standard errors and included 1 within the confidence interval.  The presence of a downstream barrier was associated with a reduction in odds by a factor of 0.274 and the odds ratio confidence interval did not include 1.  The performance evaluation of the highest ranked model indicated improved performance over the ecological model (Table 2-6). 28  There was little evidence (i < 2.0) for one best rainbow trout model (Table 2-4). Based on p Ĉ, the second ranked model, which included harvest, fit the data well.  Whereas there was evidence that the model with the highest wi did not.  Therefore I selected the model that included harvest as the only land-use variable as the best model.  Harvest was associated with a negligible increase in odds and the confidence interval did include 1 to indicate the potential of no effect (Table 2-5).  Overall model performance decreased with the addition of the only plausible land-use variable (harvest) into the best ecological model for rainbow trout (Table 2-6). 2.4 Discussion 2.4.1 Occurrence patterns of nonnative and native salmonids Brook trout had higher probability of occurrence in larger, low elevation channels. These findings were similar to another study in the Rocky Mountains where brook trout tended to persist in the main channels at lower elevations after dispersing from upstream stocking locations [Paul and Post, 2001].  In our study, brook trout occurrence increased in basins with higher harvest levels and, although other researchers have attributed these patterns to habitat degradation resulting from land-use activities, such relations cannot be inferred within this study area. Timber harvest within riparian areas can cause habitat changes, including increased water temperature [MacDonald et al., 1990], and these changes may favor brook trout over bull trout [Selong et al., 2001] and native cutthroat trout (Oncorhynchus clarki clarki) [Shepard, 2004]. However, timber harvesting practices in the study area have always included retention of streamside buffers.  Furthermore, land-use activities were most prevalent in the intermediate relief basins that brook trout commonly inhabited and largely absent for higher relief areas where brook trout were less frequently encountered.  Alternate explanations include brook trout stocking at stream crossings rather than habitat modification. The best ecological model for bull trout performed moderately.  A logistic regression model for bull trout occurrence from a nearby foothills watershed had a similar ROC value (0.80), however it included field measures (stream width and percent fines) combined with GIS measures of stream slope and harvest [Ripley et al., 2005].  They found increased bull trout occurrence with stream width and reduced occurrence as percent fines, stream slope and harvest levels increased.  In contrast, I found a non-linear relation with stream slope and little evidence for including percent harvest in our best overall model.  Our model also indicated a non-linear 29  relation with basin slope tending towards steeper terrain and an increase in occurrence with elevation.  Our best model included the downstream barrier as a land-use variable and this contributed to an improvement in the ROC score from 0.78 to 0.81.  All of the sample reaches with barriers had basin slopes < 25%, which is below the optimal relief for bull trout.  Although there may be some potential for a confounding effect from the exclusive occurrences of bull trout and this land-use variable, I did capture bull trout immediately below an impassible railway crossing in Pinto Creek, but not in any of the upstream tributaries.  Similar patterns were observed in Japan where extinction rates of white-spotted charr (Salvelinus leucomaenis) upstream from dams increased at a greater rate in headwater streams than in downstream areas [Morita and Yamamoto, 2002]. The high odds ratio for basin slope illustrates a strong relation with bull trout occurrence. Basin slope was correlated with mean basin elevation, a climatic descriptor and negatively correlated with percent wetlands, a terrain indicator.  This variable has not been included as a candidate variable in other watershed scale models of bull trout distribution [Paul and Post, 2001; Ripley et al., 2005], but a similar variable was incorporated into other watershed scale models of fish occurrence [Porter et al., 2000]. The use of reach elevation as a surrogate for climate proved successful in a bull trout model within the Rocky Mountains [Paul and Post, 2001].  Improvements to the climatic variable could capture the variations in water temperature due to ground water inputs in intermediate relief areas.  For example, mean summer water temperature (June-September) in a study area stream that supports juvenile bull trout averaged 5.4°C, whereas the mean temperature at a site in a similar adjacent watershed that doesn‟t typically support juvenile bull trout averaged 8.1°C.  These differences were attributed to different ground water inputs between basins [Jablonski, 1980].  A number of reaches shared medium or high probability of bull trout and brook trout, including an important bull trout spawning stream.  Knowledge of such locations where the two species overlap is important for identifying where threats to bull trout persistence occur [Rieman et al., 2006].  Given this overlap, future efforts to model bull trout occurrence should include the locations of ground water upwelling areas, especially within lower elevation streams that may otherwise provide marginal bull trout habitat that are well suited for brook trout. 30  The differences between our findings and those of another correlative investigation into bull trout occurrence in the foothills [Ripley et al., 2005] emphasize the need for process-based approaches that consider the confounding relation between land-use and basin relief.  Bull trout are vulnerable to a range of human impacts including illegal angler harvest, disruption of migration at stream crossings and dams, loss of habitat productivity due to sedimentation, and increases in water temperatures.  To support effective conservation efforts, future investigations should endeavor to rank these impacts. In comparison to bull trout, rainbow trout occurred more frequently in reaches that also supported brook trout and rainbow trout were more common in watersheds that supported high levels of land-use.  Rainbow trout tended towards small streams that were located in intermediate relief terrain ideal for timber production and were less prevalent in lower relief areas that support extensive muskegs.  I observed that muskeg streams had sinuous channels but these low gradient reaches often lacked important rainbow trout habitat attributes including gravel substrate and pools.  Given these observations, I were surprised that reach slope was not indicated within the best rainbow trout model.  This absence of reach slope from the model has two divergent explanations.  One previous study found no indication that trout populations were negatively affected by increases in stream slope [Isaak and Hubert, 2000].  However, two studies found that channel gradient constrained the upstream extent of coastal cutthroat and rainbow trout [Latterell et al., 2003; Fransen et al., 2006].   I suggest the use of higher resolution DEM data in future studies that explore the relations between rainbow trout occurrence and stream slope.  There are also important considerations when interpreting the suggested positive relation between rainbow trout occurrence and land-use.  Our findings were similar to those for winter steelhead (Oncorhynchus mykiss) in Oregon, where the proportion of watershed converted to shrub/young forest was positively related to winter steelhead redd abundance [Steel et al., 2004].  This relation may be confounded by a strong affinity that this fish species may have for streams located within terrain that is well suited for land-use.  Also consider that, contrary to our findings, one long-term study within our study area found that rainbow trout fry densities decreased after harvest when water yield exceeded a critical flow during the incubation period [Sterling, 1992].  Additional processed-based investigations are required to identify land-use activities that may negatively decrease survival and growth rates of native rainbow trout across their range. 31  The all species model performed moderately for predicting fish occurrence among all datasets.  These results support the use of the all species model as a tool for fish habitat conservation efforts in the study area.  Maintaining and restoring connectivity between streams is a cornerstone of most fish conservation strategies and also a legal requirement under the Canada Fisheries Act.  Twenty four percent of the 302 stream crossings were identified as potential migration barriers within medium or high probability streams.  Culvert owners could complete a simple GIS exercise to determine the extent of potential fish habitat upstream of each crossing and use this information to identify infrastructure maintenance priorities or identify sites requiring additional inventories. 2.4.2 Anticipated benefits from using a high-resolution DEM There are several anticipated benefits from using a LiDAR DEM for modelling fish occurrence.  First, a network of swales, and seepage channels that are typically represented on a synthetic network predicted from area-slope thresholds [e.g., Montgomery and Foufoula- Georgiou, 1993] are under-represented in the standard drainage network map (provincial streams layer), which was derived from aerial-photograph interpretation.  This issue may not affect the performance of fish occurrence models because streams at the no fish/fish threshold are of sufficient size to be well represented on the provincial streams layer.  However, many of the land-use activities with potential to alter hydrologic and erosion processes occur in swales and seepage channels within the under-represented portion of the drainage network.  Such impacts can be conveyed downstream to the fish bearing portion of the network.  Thus, to assess the effects of land-use on aquatic habitat, fish occurrence models derived from the 25 m DEM / provincial streams layer could be coupled with more complete models of the swale – seepage channel – fluvial channel network derived from a LiDAR DEM (Chapter 4).  A second set of benefits for using a LiDAR DEM relate to higher slope resolution.  The improved resolution should enable more accurate measurement of reach slope and improved detection of expected slope relations.  Furthermore, LiDAR DEMs are expected to improve  the ability to detect localities with high slopes.  For example, the resolution of the 25 m DEM may be inadequate for detecting knickpoints where streams flow over resistant bedrock layers; thus, localized steep reaches may not be apparent on maps.  Waterfalls, chutes, and other natural features associated with knickpoints can prevent upstream fish migration and limit fish occurrence in headwater streams of the steep foothills terrain immediately adjacent to the Rocky Mountains.  Juvenile bull trout (Salvelinus confluentus) are most common in these higher relief areas.  Thus, the use 32  of a LiDAR DEM may improve performance of fish occurrence models for bull trout and other species that inhabit steep regions where natural barriers limit fish distribution in headwater streams.  2.4.3 Limitations and opportunities for modelling fish occurrence using remote methods I identified a number of considerations for related efforts in the future.  First, correlative investigations that quantify fish occurrence from multiple predictor variables are popular approaches used to improve our knowledge of factors affecting fish occurrence, however, they have limitations introduced by non-linear relations and non-measured correlates [Rosenfeld, 2003].  Furthermore, unlike mechanistic models and true experiments, these descriptive approaches cannot indicate cause-effect relations.  Non-linear relations between a species and features of its environment are expected where a species is adapted to a specific ecological niche.  These relations have been identified using artificial neural networks [Oakes et al., 2005] or by including a quadratic-term within a logistic regression approach [Mattingly and Galat, 2002].  I selected the latter approach, which proved useful for indicating habitat partitioning among the target species in headwater streams.  However, the Pearson correlation test was not well suited for identifying relations between non-linear variables.  I identified such relations by plotting predictor variables of interest. Second, the stream slope accuracy assessment (Section 2.2.4) revealed that calibration or application of correction factors may be required to improve the accuracy of remotely-sensed parameters; otherwise, the use of uncorrected values may limit the interpretations of true habitat preferences, particularly for channel slope, which can provide an indication of expected channel unit distribution including pool extent and frequency.  I found that channel sinuosity in small streams is not accurately mapped and meandering channels extended beyond commonly accepted gradient thresholds of between 2 and 4 percent [Church, 1992].  Sinuosity has an important effect on fish habitat because it decreases stream power by reducing stream slope and increasing form roughness within habitat elements including undercut banks and pools [Leopold et al., 1964].  Therefore, fish habitat characteristics are expected to vary between a straight channel with a 3% slope and a meandering channel with 3% slope.  In watershed science there is a trend to use streams automatically generated from high-resolution DEMs rather than air photo 33  interpreted streams [Benda et al., 2007].  Generated streams follow the path from one grid-cell in the DEM to the lowest adjacent grid-cell.  However, line length reduction will remain a factor for these generated streams where the size of the watercourse is smaller than the grid-cell size. Third, the bull trout model could be improved by incorporating additional variables including: measures of the size of downstream habitats such as the Strahler order downstream of the next confluence [Oakes et al., 2005]; the locations of natural migration barriers in downstream areas [Fransen et al., 2006]; and variables that characterize temperature preferenda [Dunham et al., 2003], including indicators of groundwater upwelling areas [Baxter and Hauer, 2000]. Fourth, scaling issues are a key component of resource selection function models [Boyce, 2006].  Future studies on relations between bull trout distribution and land-use should address scaling.  Ripley et al. (2005) included reaches with areas up to 374 km 2 , whereas I limited our study to reaches with a drainage area ≤ 20 km2.  I used this approach because I anticipated that correlations with land-use would be more pronounced due to close linkages between aquatic habitat and terrestrial environment in small channels.  However, study designs that include large streams may be more suited to addressing impacts from angling.  Our automated procedure is well suited to scaling up from the reach-scale to major watersheds > 1,000 km 2  where resource management strategies are applied.  This ability to align a habitat model with the scale suitable to support resource management decisions can increase the utility of models in conservation and management programs [Boyce, 2006]. Fifth, the models presented in this chapter have no dynamic component to indicate potential changes in species occurrence in response to physical or biological changes.  For example, climate shifts may result in modified flow and temperature regimes; however, the implications of such changes cannot be predicted with the given resource selection function models.  Furthermore, biological changes including introductions of exotic species present a threat to native species conservation, the effects of which are difficult to cannot be forecast using the models that have been presented. In conclusion, an automated approach for modelling spatial fish occurrence patterns provides new insight into fish habitat requirements in small foothills streams.  The approach is conducive to validation with external data sets.  Fundamentally, the maps also provide an 34  opportunity for biologists to visually compare their own inherent spatial models of fish occurrence with those developed mathematically and then to generate ideas on how to further improve our understanding of the factors affecting fish distribution.    35  Chapter 3 Spatial organization of process domains in headwater drainage basins of a glaciated foothills region with complex longitudinal profiles 3.1 Introduction Process domains are regions characterized by distinct erosional rates and mechanisms, ultimately creating characteristic hydrogeomorphic patterns in the landscape [Willgoose et al., 1991; Montgomery and Foufoula-Georgiou, 1993].  The basic set of domains includes hillslopes, where diffusive processes dominate to maintain planar or convex slopes originating from drainage divides, and channels, where fluvial transport and debris flows have incised the land surface.  Zero-order basins, also known as hollows or swales, are transitional domains between hillslopes and channels.  They form where convergent terrain collects and discharges storm water and eroded material to create a spoon-shape slope unit that extends downhill in a hollow of varying length [Tsukamoto et al., 1982].  The distinctive forms of the hillslope and channel domains are expressed by a change in slope of the relation between local slope and upslope contributing area as observed in field studies [Montgomery and Dietrich, 1988] and in simulations [Willgoose et al., 1991].  Montgomery and Foufoula-Georgiou [1993] developed a model for delineating process domains in a landscape on the basis of changes in the slope-area relation.  However, use of this model to predict the spatial extent of process domains across a landscape is complicated by the natural migration of process domain transitions, bedrock heterogeneities, and landscape history (e.g., glaciations). When partitioning the landscape based on process domains, it is important to consider that transitions will move both at a sub-region scale according to the spatial extent of extreme rainstorms that cause mass wasting, and at a regional scale over longer time periods in response to climate change [Montgomery and Dietrich, 1988].  The channel head – the point between the downhill end of a zero-order basin and upstream end of the channel – represents the most dynamic transition.  Climate change can induce temporal variation in frequency distribution of transition site characteristics.  For example, unchanneled reaches in zero-order basins are expected to extend during climatic shifts associated with reduced gullying or landsliding [Dietrich and Dunne, 1993]. Models that delineate the landscape into process domains based on slope-area relations apply well to regions where valley spacing is not determined by structural heterogeneities in 36  underlying bedrock [e.g., Perron et al., 2009].  However, slope-area relations may fail to explain channel head occurrence in locations where features such as bedrock springs associated with specific rock formations exert control [Jaeger et al., 2007].  In addition to influencing subsurface flow properties, bedrock structure can disrupt the characteristic negative relation between drainage area and local slope along a stream channel [Hack, 1957], and the corresponding downstream succession of channel types [Montgomery and Buffington, 1997]. For example, knickpoints in a stream long profile are controlled by tectonics/lithology [Troiani and Della Seta, 2008].  However, they can also result from local base level change [Bishop et al., 2005], valley-side mass wasting [Phillips et al., 2010], and climate change [Hayakawa and Oguchi, 2009]. In previously glaciated catchments, glacial landforms largely dictate post-glacial landscape architecture [Montgomery, 2002; Brardinoni and Hassan, 2006].  Contemporary fluvial processes are superimposed upon a modified longitudinal profile characterized by a hierarchy of glacial troughs controlling the spatial distribution of process domains, sediment sources, and deposition locations [Brardinoni and Hassan, 2007; Brardinoni et al., 2009].  The boundaries between contemporary process domains in glaciated topography are thus determined by the glacially formed topography, leading to longitudinal channel profiles that deviate from the expected slope-area relation that is typical of landscapes with near-equilibrium long profiles [Brardinoni and Hassan, 2007]. While Brardinoni and Hassan [2007] clarified the influence of glaciation on the organization of process domains in a coastal mountain environment dominated by intrusive igneous bedrock, their results likely do not apply to glaciated areas with differing geology and physiography. Various approaches have been used to predict the extent of process domains.  Field- based relations can be extrapolated using measures of slope and drainage area derived from a digital elevation model (DEM) to predict channel head locations [e.g., Montgomery and Dietrich, 1992; Miller, 2003].  Where field data are limited, DEM-derived topographic data can be used on their own to estimate channel head locations.  Such approaches combine slope-area relations with measures of local curvature indicative of concave topography to estimate upstream channel extent [e.g., Band, 1986; Clarke et al., 2008].  To predict the headward extent of channels, Heine et al. [2004] developed a logistic regression model with a suite of topographic descriptors that had higher accuracy than using either a constant threshold for area 37  or a slope-dependent area threshold for the same study area.  Similar models have also been applied to predict flow permanence  [Olson and Brouillette, 2006].  Locations of flow permanence and channel heads are related; however flow permanence is more transient than channel heads and harder to predict [Dietrich and Dunne, 1993].  These cases indicate that a hybrid topography-based probabilistic model may prove useful in relating the transitions between process domains to map-derived predictors. Poor DEM resolution once limited the remote detection of channel networks; however, new technologies such as Light Detection and Ranging (LiDAR) have yielded high-resolution DEMs, creating opportunities for a wide range of hydrologic and geomorphic studies [Tarolli et al., 2009].  Applications related to the spatial distribution of process domains include detecting channel head locations and landslide initiation points [Tarolli and Tarboton, 2006; Tarolli and Fontana, 2009], and channel morphology [Cavalli et al., 2008].  High-resolution DEMs provide opportunities to further explore spatial patterns of erosional processes within complex landscapes. Knowledge of the spatial extent of process domains is critical not only to predicting the spatial organization of hillslope and channel characteristics, particularly the locations of process domain transitions, but also to understanding the response of a landscape to climate shifts and land use changes.  The objectives of this paper, therefore, are (1) to explore the relations between the spatial distribution of process domains and a variety of predictors extracted from a high resolution DEM in a glaciated foothills region, and (2) to develop an automated approach to predict the spatial extent of process domains. 3.2 Study area The study area includes Dutch Creek and adjacent unnamed basins which cover 57.4 km 2  of the foothills in West Central Alberta (Figure 3-1).  The area‟s mixed lithology and multiple ice sheet advances contributed to topographic anisotropy characterized by irregular long profiles.  The foothills, typically underlain with Mesozoic sandstones and siltstone, represent a transition zone between the younger Cenozoic formations of interior boreal plain to the east and the older Paleozoic formations of the Rocky Mountains to the west.  Thrust faults and folding have formed parallel rows of foothills ridges that run from the southeast to the northwest, dipping to the southwest and striking to the northeast.  The study area includes three fault regions: one through the upper reaches of Basin 1; a second that follows the lower and mid 38  reaches of Dutch Creek; and a third along the eastern base of the Brazeau Range (Figure 3-1c and d).  The faulting has contributed to an unusually complex lithology for the foothills, with a rare outcrop of Upper Paleozoic limestone located in the northeast part of the study area (Figure 3-1b) [Crombie and Erdman, 1947; Douglas, 1956].  This formation dips strongly to the southwest, where younger Cretaceous marine shales, sandstones and conglomerates, typical of the foothills, are exposed (Figure 3-1d).  The youngest formation comprised of erodible shale of the Blackstone formation underlies Basin 2.  Basins 1 and 3 are underlain by sandstone, shale, and conglomerate from the Mountain Park Formation.  The Cadomin Formation, a thin exposure of highly resistant chert and quartzite conglomerate, closely borders the limestone (Figure 3-1b), and also underlies a plateau in the central study area including portions of Basins 4 and 5 (Figure 3-1d).  Along the southern boundary of the study area, the North Saskatchewan River flows nearly perpendicular to the fault through a feature known as the Gap (Figure 3-1c).  Thus, a cross-section of the Mesozoic and Upper Paleozoic formations are exposed in the small study area drainages on the north flank of the river valley.  Differential erosion is expected due to varying rock cohesions of these formations. 39   Figure 3-1 (a) Inset map of study area location. (b) Inset map of ice sheet flow direction down North Saskatchewan River valley with resistant conglomerate of the Cadomin Formation and Upper Paleozoic limestone, plus extent rectangle of study area.  (c) Candidate sub-basins of the study area. (d) Faults and resistant formations from scanned geology maps [Crombie and Erdman, 1947; Douglas, 1956] that were geo-referenced and digitized. 40  In the western and central thirds of the study area, the combination of highly erodible sandstone/siltstone bedrock and extensive glaciations contribute to rolling foothills topography, while the eastern third, with the isolated outcrop of erosion-resistant carbonate bedrock, supports rugged mountainous topography.  Elevations range from 1100 – 1900 m.  Annual precipitation averages 595 mm, with 70 percent as rain and 30 percent as snow.  Surface erosion of the forest floor is not an active process due to high infiltration rates associated with thick litter and duff layers typical of forested foothills soils [Archibald et al., 1996]. The study focused on 14 lower relief sub-basins within the 57.4 km 2  area (Figure 3-1c). These candidate basins most closely represent the broader foothills where the Cordilleran ice sheets repeatedly advanced in an easterly direction through main valleys before fanning out over the region.  These advances reorganized the fluvial landscape within the foothills as evidenced by common features including buried valleys, where original valleys are filled with glacial drift up to 290 m in depth [Roed, 1968].  Surficial material depths and type are highly variable [Bayrock and Reimchen, 1975] (Figure 3-2a).  Deeply leached pre-Wisconsin till, with depths greater than 2 m, covers a high elevation plateau in the northeast corner of the study area. Lower elevation toe slopes have slightly leached late Wisconsin till with depths less than 3 m. Shallower colluvial deposits extend across mid-elevations.  Within the 14 candidate sub-basins, surficial deposits are stable and the terrain lacks the steepness required for mass wasting.  Mean basin slope averages 13 degrees (24 percent) within these 14 sub-basins, which are relatively pristine with minimal roads and no commercial forest harvest.  Bulldozer trails constructed for seismic-based petroleum exploration are the main anthropogenic disturbance with the potential to alter erosion processes. These trails are ubiquitous across the foothills region of Alberta.  Five basins were selected from the 14 candidates for more detailed study (Figure 3-1c); they represent the cross-section of valley positions, lithology and surficial material types.  41   Figure 3-2 (a) Map of candidate sub-basins and surficial materials with five research basins numbered including field surveyed longitudinal profiles with process domain boundaries. (b) Map of candidate sub-basins and SL/k index values by stream reaches for the overdrawn drainage network with five research basins numbered.  SL/k index values in the 3-6 class represent moderately over-steepened reaches, and values in the >6 class represent severely over- steepened reaches. 42   3.3 Materials and methods This project was completed in five steps. 1) An overdrawn drainage network, whereby flowpaths were intentionally overextended into hillslopes from the swale origins, was extracted from a high-resolution DEM.  The network was sub-divided into reaches, which were characterized by a suite of topographic predictors.  2) Geo-referenced ground surveys were completed in five candidate sub-basins.  Starting at an origin point from the overdrawn network, each transect followed a steepest-descent flowpath through the series of process domains.  3) For each transect, all corresponding reaches of the overdrawn network were assigned to one of four process domain classes, with 1 = hillslope, 2 = swale, 3 = colluvial channel and 4 = fluvial channel.   4) Ordinal logistic regression was used to fit models for predicting the probability of occurrence for each of the four process domain classes based on the reach-scale topographic predictors.  5) The best model was applied to all reaches of overdrawn network. Reaches that were assigned to the hillslope domain were truncated to produce the final drainage network model. 3.3.1 Creating the overdrawn drainage network The LIDAR data were collected during snow-free conditions in May 2005 from a helicopter flying 400 m above ground level at a forward speed of 90 kmhr-1.  A Helix-1 LIDAR system was used with a Riegl model Q140HR laser pulsating at 17 800 Hz.  The scan angle was 60°, with 30° on either side of nadir; LIDAR swath width was 460 m with 30% overlap; and scan line rate was 18 Hz.  A 1.3 m x 1.3 m grid pattern point spacing was used with a density of 0.52 pointsm-2. A Sony TRV900 progressive scan digital video camera was oriented in nadir position. The point LIDAR data were used to build raster DEMs with three resolutions: 1 m, 2 m and 5 m.  The 1 m DEM contained numerous closed depressions, while loss of detail was a concern with the 5 m grid.  The 2 m DEM was the best compromise, consistent with other process domain studies employing LIDAR data [Tarolli and Fontana, 2009]. The overdrawn drainage network and stream reach attributes were created using two FORTRAN computer programs (Bld_grds and Netrace) that were developed for DEM analysis [Miller, 2003; Clarke et al., 2008].  The programs share a set of approximately 70 user-specified 43  parameters designed to account for regional geographic factors that influence stream channels. In the overdrawn drainage network, channels originated at drainage areas between 0.05 and 0.50 ha, and drainage density averaged 20.0 km/km 2  across the 57.4 km 2  study area (Figure 3-2b). The overdrawn drainage network was divided into reaches averaging 13 m in length, which translated to an average of 14.4 bankfull channel widths across nine field sites (8.0 minimum – 22.6 maximum).  A set of 67 sub-basin polygons (mean area of 86 ha) were created. 3.3.2 Process domain field surveys The specific field criteria for each process domain are as follows.  Hillslopes are characterized by convex or planar slopes.  Swales (zero-order basins) originate in areas of topographic convergence evidenced by concave slopes.  Some distance down each swale, sediment transport becomes sufficiently concentrated to form a colluvial channel within definable banks.  Colluvial channels originate at the channel head.  The colluvial/fluvial channel distinction was consistent with Montgomery and Buffington [1997].  In non-landslide-prone foothills terrain, soil creep is the dominant process for moving sediment from adjacent hillslopes to form the bed and banks of colluvial channels.  Field identification criteria include: a wide range of grain sizes in the channel bed; significant amounts of organic material, including stochastically distributed wood and live roots, which exert a strong influence on sediment storage; and an overall high complexity due to lack of hydraulic conditions required to create regular sequences of bed structures and bank features.  The transition to fluvial channels is expected where an increase in stream power and hence competence becomes sufficient to create typical small stream morphologies [e.g., Church, 1992].  In fluvial channels, channel-bed and bank material is primarily transported from upstream channel reaches rather than adjacent hillslopes. Field surveys were completed within five sub-basins, including three within the sandstone/shale bedrock type representing upper, middle and low positions along the valley side, plus two basins with conglomerate exposures bordering the limestone region.  Using GIS, UTM coordinates were obtained for a number of origin locations from the overdrawn drainage network for each candidate sub-basin.  In the field, a GPS was used to navigate to an origin location in each sub-basin.  This location served as the start location for a longitudinal profile survey.  The origin location was found to be within the hillslope process domain in all but one survey (Profile 2).  For that survey, the hillslope-swale transition was located 50 m upslope of 44  the origin location from the overdrawn network.  Each survey followed the fall-line of a hillslope until a hollow was encountered.  At the base of the hollow, the survey continued downslope through the swale, colluvial channel, and fluvial channel sequence.  A GPS location was recorded at the base of the hollow and each subsequent process domain boundary.  All significant features (e.g. headwall of hollow, base of hollow, headcuts, knickpoints, and fans) were also recorded.  Surveys were terminated at intersection points with roads and seismic trails where obvious changes in water and sediment inputs corresponded to changes in process domain.  In these cases, the GPS was used to navigate to another channel origin location within the sub-basin and another survey was started.  Data from the truncated profiles were considered to be representative of natural conditions and were included in the analysis. A total of seven profiles were surveyed including two truncated profiles within one basin (Figure 3-2a).  Transect length varied between sub-basins and all four process domains were found on four of the transects; three transects lacked at least one domain (Table 3-1). Table 3-1. Profile lengths and process domains. Profile Profile length Process domain 1   (km) H S C F 1 1.13     2 1.16     3a 0.47    3b 0.18   3c 1.13     4 0.72    5 2.33     1  H = hillslope, S = swale, C = colluvial channel, F = fluvial channel. 3.3.3 Variables used to predict process domain class All field surveyed reaches (n = 686), along with the remaining reaches of the overdrawn drainage network (n = 82 819), were characterized by a set of predictors (Table 3-2).  Of these, the SL/k index and profile complexity index had not been used in other process domain classifications and are herein described.  45  Table 3-2. Candidate predictor variables generated for every reach with considerations for delineating process domains. Predictor Scale Measurement method 1. Drainage area Basin Drainage area (ha) above the downstream end of each reach calculated with FORTRAN program [Miller, 2003]. 2. Reach slope Reach Channel gradient (%) calculated for each channel pixel, then averaged over the length of the reach using a FORTRAN program [Miller, 2003]. 3. Profile curvature 6, 10, 20, 50, 100 and 200 m A topographic index describing the shape of the line following the maximum slope (typically the flowpath).  Positive values indicate increasing slope, zero indicates consistent slope, and negative values indicate decreasing slope in a downstream direction. The index was calculated from the midpoint of each reach for six different line lengths using a FORTRAN program [Miller, 2003] that incorporated the formulas from Zevenbergen and Thorne [1987]. 4. Planform curvature As above A topographic index describing the shape of the line perpendicular to the line of maximum slope. Positive values indicate divergent (convex) topography, zero values indicate a flat surface, and negative values indicate convergent (concave) topography.  Calculated for six different lengths as above. 5. Mean basin slope  Basin scale using 2 and 6 m grid cells Hillslope gradient (%) averaged over the entire upstream contributing area for each reach using a FORTRAN program [Miller, 2003].  An important predictor of stream channel presence [Heine et al., 2004]. 6. SL/k index Reach Calculated with a FORTRAN program developed for this project. 7. Profile complexity index Entire transect Standard deviation of the SL/k index for all reaches in a transect. 8. Lithology Reach Coded as (0) for sandstone/shale and (1) for limestone based on Crombie and Erdman [1947] and Douglas (1956).  46  A straight line on a semi-logarithmic plot of elevation (linear) versus distance (logarithmic) represents an equilibrium or graded long profile [Hack, 1957].  This plot is hereinafter called the Hack profile.  However, most natural streams do not follow a single logarithmic profile along their entire course, so it is more informative to create a curved Hackprofile with the elevation and distance measures for successive reaches (e.g. Figure 3-3). By comparing the curved Hack profile to the straight line of the graded profile, under-steepened and over-steepened reaches become apparent.  The SL/k index, which is sensitive to stream slope, provides a measure of the deviation of an individual reach from its graded river profile [Hack, 1973].  Long profile anomalies have been quantified using the SL/k index [Perez-Pena et al., 2008] (where S = reach slope, L = distance from reach to source, and k = graded river gradient), computed as ( ) ln( ) s f t h h k L              (1) where hs is the river head elevation (m), hf  is the elevation of the river mouth (m) and Lt is the total length of the river (m).  Because Perez-Pena et al. [2008] used the SL/k index to detect anomalies related to tectonic and lithologic features along entire river courses, k was calculated using the river mouth as the downstream point of reference. In this study, the point of interest is the downstream end of each reach rather than an arbitrary stream mouth and the formula for k was adjusted accordingly, so that the SL/k index was computed as follows: ln( ) ( ) r r s r Sx x SL k h h         (2) where xr is the distance from the stream source to the middle of the reach (m) and hr is the elevation at the mid-point of the reach (m).  It is informative to display the curved Hack profile with height on the primary Y axis and a step curve of the SL/k index using a secondary Y axis (e.g. Figure 3-3), where the SL/k index value is held constant over the length of each reach to create the step-like appearance [Chen et al., 2003]. 47   Figure 3-3. Hack profiles with step SL/k index curves for five sub-basins.  Dotted black vertical line is boundary between hillslope and swale process domains.  Black vertical arrow is channel head location.  Dotted grey vertical line is boundary between colluvial and fluvial process domain.  SL/k = (reach slope × distance to source) / (rise to source / log (distance to source)).  3.3.4 Modelling approach I developed models to augment the information on process domain transitions derived from slope-area plots (Figure 3-4).  Important findings from the slope-area plots are described herein.  A positive slope-area relation within the hillslope domain, typical of non-glaciated drainage basins [Montgomery and Foufoula-Georgiou, 1993] and the glaciated Olympic Mountains [Montgomery, 2001], was not evident from the bin-averaged slope data from the five research watersheds.  However, the consistent negative relation that steepened sharply within the largest drainage area classes was similar to relations from a number of glaciated mountain drainage basins in British Columbia [Brardinoni and Hassan, 2006].  The average drainage area was 1.9 ha (range from 0.5 ha to 4.6 ha) at the hillslope to swale transition, 18.4 ha (7.0 ha to 48  32.0 ha) at the swale to colluvial channel transition, and 56.2 ha (37.9 ha to 79.3 ha) at the colluvial to fluvial channel transition.  The drainage areas of the hillslope-swale and swale- colluvial channel transitions were more than one order of magnitude greater (i.e. 1.9 ha versus 0.1 ha, and 18.4 versus 1.0 ha) than those determined by Montgomery [2001] for the Olympic Mountains, consistent with differences in continental versus coastal precipitation.  The slope- area plots for individual transects reveal complex longitudinal profiles (Figure 3-4).  Explaining the different drainage area thresholds at domain transitions between these complex profiles was the main challenge in modelling spatial distribution of erosion processes.  Thus, factors in addition to drainage area and slope that affect domain transitions were explored using other modelling procedures.  Figure 3-4. Plot of log-bin averaged drainage area versus slope for the five research basins. The plot shows the mean slope of individual 2 m grid cells for each 0.1 log bin in drainage area. Black lines represent longitudinal profiles of seven transects with field determined domain class assigned. Drainage area and slope values for profiles were obtained from reaches of the overdrawn network.  Dashed vertical lines divide the plot into process domains based on average drainage area at transitions.  Channel head locations correspond to the left-most colluvial channel symbol.  49  Process domains have been previously modelled using a variety of variable types. Brummer and Montgomery [2003] used a quantitative response variable of the interval class – median bed surface grain size (D50) – to model a shift in trend from downstream coarsening to downstream fining as an indicator of the transition from colluvial to fluvial channels.  Heine et al. [2004] used a categorical response variable of the binary type – channel/non channel – to predict drainage network extent.  In comparison to categorical type variables, interval variables (e.g., D50) have high information content and the potential to reveal more about relations with predictors [Venables and Ripley, 2002].  When developing models from categorical data, the specific samples near transitions provide more information than samples that are more dispersed [McCullagh and Nelder, 1989].  These facts favor the use of interval over categorical variables; however, there is no single quantitative variable that spans and is diagnostic of all of the process domain classes of interest in this project.  Thus, for this study, the response variable was a categorical variable with four ordered classes, coded as 1 = hillslope; 2 = swale; 3 = colluvial channel; and 4 = fluvial channel. Cumulative logit regression can be used to analyze ordinal response variables.  Prior to model development, cumulative logit plots can be used to explore relations between individual predictors and the process domain class.  The first step in plot preparation is to divide each predictor of interest into a series of bins (i.e. n = 20 in Figure 3-5).  The response variable, process domain class for example, is converted into J – 1 sets of cumulative dichotomous variables (Table 3-3).  The count of each cumulative variable within each bin is tallied and cumulative logits are calculated. The logit, defined as the logarithm of the odds ratio, is   1 p Logit O log p                 (3) The cumulative logit plot is produced by plotting the midpoint value for each bin of the interval variable against its respective cumulative logit.  Each dichotomous variable is plotted as a separate line.  On the plot, a horizontal reference line can be added at cumulative logit = 0 to represent the threshold where the odds of the event of interest = 50/50 (i.e. Figure 3-5).  50  Table 3-3. Example of the conversion of an ordinal variable with four classes into three non- redundant cumulative variables. Process domain class L1 Hillslope L2 Swale and hillslope (Swale -) L3 Seepage channel, swale, and hillslope (Seepage channel -) 1 1 1 1 2 0 1 1 3 0 0 1 4 0 0 0  Cumulative logit regression predicts the cumulative probability, Fij, that individual i is in the j th  category or lower, rather than being in a higher category.  The cumulative logit model is a set of J – 1 equations with different intercepts but common slopes as 1 1 2 2     ... 1 ij j p p ij F log X X X F                          (4) where α1 to αj-1 are the set of category specific intercepts, β1 to βp are the set of shared slopes, and X1 to Xp are the predictors.  Each equation in the set is used to calculate a sequential dichotomous variable.  The proportional odds assumption requires that the set of common regression coefficients function equally well across the series of dichotomous variables. Cumulative logit plots can be examined for parallel trends to help assess this assumption.  For each sample, the category (1...j) with the highest logit is assigned the predicted class for that sample.  Additional details of three aspects of categorical response variable modelling are provided in Appendix A: (1) dependent variable types and options for modelling; (2) the logit model as the base model for categorical dependent variables; and (3) the ordinal logistic model as a special form of the logit model. The analysis proceeded in two steps: (1) exploration of the relations between process domain class and the topographic and geological controls for surveyed reaches, and (2) development of regression relations to predict process domains for unsurveyed reaches. Step (1) involved viewing cumulative logit plots prior to fitting regression models.  To avoid over-fitting and discovery of effects that are spurious, every effort was made to limit the number of predictors in each candidate model: one of the variables was excluded when independent variable correlations were greater than |0.6|; use of quadratic terms (i.e., x1 + x1 2 ) to quantify 51  potential non-linear relations were limited; and interaction terms (i.e. x1 × x2) were not included. Regression parameters were examined for further statistical evidence to confirm suspected relations.  The best model was selected based on highest accuracy and least unexplained variation.  In step (2), the criteria for candidate variables were relaxed to include interaction terms, and quadratic terms in cases where non-linear relations were expected.  Stepwise selection was used to select the best model.  The analysis was complicated by the structure of the sampling scheme, which produced data that included fixed and random effects, and that were spatially autocorrelated. Fixed effects are model parameters associated with an entire population.  For example, in this study, each reach had an associated measure of drainage area and the entire range of interest for this predictor was sampled.  Random effects represent individual experimental units drawn randomly from a population.  For example, in medical research, individual patients that are randomly selected from a population and periodically re-assessed are a random factor effect. In this study, where one transect was located within each randomly selected watershed, transect was a random effect. Each reach was nested as part of a sequence of reaches within a transect, and thus the reaches were not statistically independent.  For example, the lithology of reaches from the same transect was more likely to be equal than the lithology of reaches from separate transects, and drainage area of two reaches that are closer in space along one transect were more likely to have similar values than measurements that were farther apart.  This lack of independence can translate into obvious spatial patterns in the residuals of a model.  If the residuals from a statistical model are correlated, then the assumptions of independence and identical distribution are violated, potentially biasing parameter estimates and inflating Type 1 error rates (incorrectly rejecting the null hypothesis of no effect).  These problems are well recognized in the statistical and ecological literatures [Wang and Goonewardene, 2004; Dormann et al., 2007]. Mixed-effects models were appropriate for this type of hierarchical data structure [Littell et al., 2006], with transect being a random effect and the predictor variables being fixed effects. The mixed-effects model for ordinal logistic regression for the j th  category is the basic form of the cumulative logit model (Equation 4), plus random variable bk  as: 52  1 1 2 2     ... 1 ij j p p ij k F log X bX X F                    (5) where bk  represents the deviation from the population mean for each randomly chosen transect k.  When the goal was to extrapolate the model to un-sampled watersheds within the original population of watersheds, the full solution including the random effects solutions bk, could not be applied and only the intercepts and fixed-effects coefficients were used.  Mixed-effects models can also accommodate spatial autocorrelation within transects, as well as the unequal numbers of reaches among transects and the fact that some transects did not contain all process domains. The mixed-effects ordinal regression analyses were performed with the GLIMMIX procedure in SAS 9.2 [SAS, 2008], which uses a Gaussian quadrature approximation to log- likelihood (quadrature approximation). In terrestrial ecological studies, spatial autocorrelation is typically quantified in relation to a straight-line distance in two dimensions [Peterson et al., 2007].  In stream networks, it is effective to quantify spatial autocorrelation in relation to distance along a stream network from a point of interest to sample location [Gardner and Sullivan, 2004].  For this study, distance from the source for each sample site was calculated using a FORTRAN program.  Herein, the term source is used to refer to the source of the overdrawn network. To address the issues with random effects and spatial autocorrelation accounting (SAA), three different error structures were evaluated for each candidate model: (1) no random effects; (2) random effects; and (3) random effects + SAA.  The models with random effects were evaluated based on population-average predictions and best linear unbiased predictions.  The population-average predictions included only the fixed effects (i.e., the predictor variables) and thus had applications for regional mapping; however, information on transect-specific variation was lost.  In contrast, the best linear unbiased predictions (BLUP) included both the fixed effects plus the random effects estimates for each transect, thereby allowing predictions within an individual transect to vary from the population average estimates [Littell et al., 2006]. Random effects estimates can also be used to assess subject-specific responses to additional predictors.  Hence, they typically have higher accuracy and can be explored for additional information on the processes of interest.  Therefore, in order of complexity, the five random effects configurations were: (1) no random effects, (2) random effects with population average predictions, (3) random effects and SAA with population average predictions, (4) random 53  effects with transect-specific predictions, and (5) random effects and SAA with transect-specific predictions. The first three models could be applied to reaches in un-sampled transects and were thus applicable for regional mapping. Models 4 and 5 could only be applied to the sampled transects, and would be useful in a context where only a sample of reaches was studied on each transect, but where predictions were desired for the entire stream. Cross-validation was used to choose models.  Model training and validation datasets were created using the hold-out approach where the full dataset was portioned into 10 non- overlapping pieces of similar size.   Each of 10 model training datasets (i = 1 …10) was comprised of a unique combination of nine pieces and the remaining piece was used for validation.  All candidate models were executed 10 times, each time with a unique model training dataset i. Total mean accuracy (Atot), for each candidate model was calculated as 2 10 i i tot Atrain Avalid A            (6) where Ai train was the accuracy of Model Yi applied to the training dataset i and Ai valid was the accuracy of Model Yi applied to the validation dataset i.  Total mean penalty for each candidate model was calculated by substituting mean squared error for accuracy in Equation 4, but not dividing by 2.  Model choice was based on maximum accuracy and minimum penalty. The fixed effects covariates were log transformed as required to achieve a more normal distribution.  Curvature descriptors included negative and positive values and were transformed using the Box-Cox transformation in the TRANSREG procedure in SAS 9.2 [SAS, 2008].  To remove scale effects and assist with model parameter and odds ratio interpretation, predictors were standardized to mean = 0 and standard deviation = 1. 3.4 Results and discussion 3.4.1 Linkages between geology, DEM-derived topographic indices, and field observations Irregular long profiles were evidenced by numerous stream reaches with slopes that are steeper than expected based on a graded river profile (Figure 3-2b).  A coarse comparison between the spatial distribution of SL/k values and bedrock geology from 1:63,360 scale maps [Crombie and Erdman, 1947; Douglas, 1956] supports the findings of Crombie and Erdman 54  [1947], that the topography is largely determined by a combination of differential erosion of the various formations and faulting.  Mesozoic formations that lack conglomerates, including the Blackstone (Basin 2) and Luscar Formations (upper plateau in Basin 4 and 5), had frequent under-steepened reaches (SL/k index < 3).   Basin 1 occupied a lower valley position at the transition between mid-slope colluvium and valley bottom late Wisconsin till (Figure 3-2a).  It had a drainage area of 0.65 km 2  where it met a larger stream draining 7.07 km 2 .  The curved Hack longitudinal profile of Transect 1 resembled a hanging valley with increasing gradient towards its outlet (Figure 3-3).  The step curve of the SL/k index highlighted three over-steepened regions.  The first region – closest to the source – was located within 50 m of the nearest fault and 100 m from the origin of swale (Figure 3-2b and Figure 3-3).  The channel head was located in the second region at the confluence of two swales.  The fluvial channel started in the third region as the stream dropped into a gorge associated with the mainstem channel near the Mountain Park - Luscar Formation boundary. Basin 2 occupied a valley side position within the shale of the Blackstone Formation, and had complete coverage of colluvial surficial materials (Figure 3-2a).  Of all five basins, the curved Hack profile for Transect 2 most closely resembled a graded river (Figure 3-3).  The step curve of the SL/k index for Basin 2 showed this same general pattern with the colluvial / fluvial transition located in a moderately over-steepened section. Basin 3 originated at the divide between the North Saskatchewan River valley and the next valley to the north (Figure 3-2a).  Three profiles were surveyed in this basin.  Transect 3a was terminated at 472 m where the colluvial / fluvial transition aligned with an abrupt grade change created by an abandoned mining road crossing.  Profile 3b ended at the intersection of the swale and an all-terrain-vehicle trail that followed a seismic exploration line.  Sediment generated from the eroding trail surface created an abrupt transition to a colluvial channel at the intersection.  Headward extension of the fluvial channel into the colluvial channel along Transect 3a, and the colluvial channel into the swale along 3b, resulted from land-use impacts; therefore, the impacted sections were excluded from the analysis.  Transect 3c was not affected by land-use and extended for 1,134 m from the source of the overdrawn network.  The swale origin corresponds to a  point at 102 m from the source of the overdrawn network as shown on the curved Hack profile, with the channel head at 436 m within a moderately over-steepened 55  section of the profile (Figure 3-3).  The Mountain Park - Luscar Formation boundary, located 700 m from the source, was characterized by a slight reduction in SL/k index.  In a 20 m section of the colluvial channel where a detailed channel assessment was conducted, five headcuts contributed to a stepped profile. In comparison, headcuts were absent in a section of fluvial channel. The resistant conglomerate of the Cadomin Formation forms a cap on the divide of the North Saskatchewan River valley in the vicinity of Basins 4 and 5 (Figure 3-2b).  A thick pre- Wisconsin till covers the plateau above the Cadomin formation (Figure 3-2a).  In the 1:63,360 scale geologic map [Crombie and Erdman, 1947], the Cadomin Formation follows the 5,600 foot contour line through Basin 4; however, the valley associated with Basin 4 is not depicted by the contour lines of the original map.  Thus, the location for the intersection of the longitudinal profile and the Cadomin Formation was corrected based on the LiDAR topography (Figure 3-3). In Basin 4, the step curve of the SL/k index showed the hillslope / swale transition within a moderately over-steepened section, and the channel head within a severely over-steepened section associated with the corrected location of the Cadomin Formation.  Bedrock was frequently exposed within the colluvial channel.  The profile in Basin 5 had the greatest distance between an origin of the overdrawn network and a hillslope/swale transition, with the transition located on the plateau (Figure 3-3).  The channel head was located approximately 150 m downslope of the Cadomin Formation within an over-steepened section of the profile where the SL/k index increased to a value near 6. These results indicate the SL/k index maps developed from LiDAR DEMs may be valuable for refining bedrock geology maps, especially in glaciated regions where glacial drift covers much of the landscape and bedrock exposures are limited.  Field surveys could target over-steepened reaches where bedrock formations may have a higher likelihood of exposure. 3.4.2 Screening candidate predictors prior to establishing relations To minimize the risk of misleading results from the multivariate analyses, which were intended to highlight relations between predictors and process domains, candidate predictors were screened using two criteria.  First, the data from all transects were pooled and each predictor was examined for compliance with the proportional odds assumptions. To test this assumption, cumulative logit plots were examined for parallel trends.  Secondly, transect- specific data for each predictor were examined on cumulative logit plots for evidence of a 56  consistent relation across transects.  Predictors that did not display consistent positive or negative trend across at least 70% of transects were screened out. Drainage area was the only continuous reach-scale predictor that both met the proportional odds assumption (Figure 3-5a) and displayed consistent trend across transects for all four process domain classes (Figure 3-5b).  Two other predictors (lithology and profile complexity index) were also used.  Drainage area thresholds at transitions between domain classes (cumulative logit = 0) occurred near 2, 20 and 70 ha for the hillslope/swale, swale/colluvial channel, and colluvial/fluvial channel respectively (Figure 3-5a).  These thresholds align closely with those from the slope-area plots of 2, 18, and 56 ha respectively (Figure 3-4), confirming congruency between the two different methods.  For individual transects, the drainage area at the channel head (swale/colluvial channel transition) ranged from 7 – 30 ha (Figure 3-5b).  Although reach slope generally met the proportional odds assumption (Figure 3-5c), it failed to meet the requirement for consistent trend across transects (Figure 3-5d).  The negative relation for the pooled data (Figure 3-4 and Figure 3-5c) was obscured by the variation between profiles.  Thus, to facilitate the exploration of other factors influencing process domain transitions, reach slope was excluded from the first modelling exercise intended to highlight relations within individual transects.  These results are contrary to other studies where slope is the dominant variable for controlling process domain transitions [e.g., Hupp, 1986].  Other excluded variables included mean basin slope from 2 and 6 m grids, planform and profile curvature at all lengths.  Several predictors met requirements for a limited set of process domains (hillslope, swales, and colluvial channels), including SL/k index, profile curvature at 50 and 100 m, and planform curvature at 100 and 200 m. 57   Figure 3-5. Plots of cumulative logit vs. predictors: (a) area for all domains, (b) area by transect for swale and lower domains (swale -), (c) reach slope for all domains, and (d) reach slope by transect for swale and lower domains (swale -).  58  3.4.3 Relations between predictors and all four process domains Candidate models for predicting all four process domains included Area (Model 1), Area + lithology (Model 2), Area + profile complexity index (Model 3), and Area + lithology + profile complexity index (Model 4).  Area was expected to perform as a dominant predictor and was coded in the model for a potential non-linear relation (area + area 2 ), whereas all other predictors were coded for linear relations.  For regional mapping applications, Model 4 had the highest accuracy (87%), the lowest penalty (20%), and the highest performance in each of the three types of predictions (Figure 3-6a and b).  For Models 2 and 3, performance across the three regional mapping configurations was inconsistent.  These patterns can be explained by comparing the fixed and random-effects parameter estimates from the four models.  Figure 3-6. Plots of (a) accuracy and (b) penalty with standard error by random effects type for four models predicting all four process domains.  Random effects type 1 = no random effects, type 2 = random effects with population average predictions (PAP), type 3 = random effects and spatial autocorrelation accounting (SAA) with PAP, type 4 = random effects with transect specific predictions (TSP), type 5 = random effects and SAA with TSP.  59  Negative values for a fixed-effect estimate indicate a decreased probability of being in a lower order class as values increase.  The two drainage area parameter estimates (area and area 2 ) showed this negative value (Table 3-4), which was consistent with the trend in the raw data (Figure 3-5a).  The estimate for profile variation index was positive, indicating that as profile complexity increases, there was a greater probability of being in a lower ordered class (i.e. thresholds between process domain classes have larger drainage areas).  This trend towards thresholds at smaller drainage areas in simple longitudinal profiles was apparent in the comparison of the transition locations along the profiles in Basin 2 (profile complexity index = 0.052) and in Basin 4 (profile complexity index = 2.64) (Figure 3-3).  Brardinoni et al. [2009] observed such a shift to lower ordered process domain in response to increased profile complexity in a glaciated mountain drainage basin where a fluvial channel transitioned to a colluvial channel as a hanging valley stepped down to a larger glacial trough. 60  Table 3-4. Ordinal logistic regression results for the best process domain models including model selection procedure, domain classes, coefficient estimates and their 95% confidence intervals (lower and upper), standard errors (SE), significance (P), odds ratio and their 95% confidence intervals. Selection Domain classes Variable Estimate Lower Upper SE P Odds Ratio Lower Upper Logit plot preview All domains Intercept 1 -11.37 -13.01 -9.73 0.84 <0.01 -- -- --   Intercept 2 -3.31 -4.05 -2.57 0.38 <0.01 -- -- --   Intercept 3 3.85 2.99 4.72 0.44 <0.01 -- -- --   Area -7.61 -8.68 -6.54 0.54 <0.01 <0.001 <0.001 <0.001   Area 2  -1.24 -1.52 -0.95 0.15 <0.01 -- -- --   Lithology = limestone 1.82 1.20 2.45 0.32 <0.01 6.20 3.30 11.63   Profile variation index 1.04 0.57 1.51 0.24 <0.01 2.83 1.76 4.54  Domains 1-3 Intercept 1 -13.02 -16.10 -9.94 1.57 <0.01 -- -- --   Intercept 2 -1.64 -2.83 -0.44 0.61 0.01 -- -- --   Area -8.87 -11.06 -6.68 1.11 <0.01 <0.001 <0.001 <0.001   Area 2  -1.57 -2.20 -0.94 0.32 <0.01 -- -- --   Lithology = limestone 0.75 -0.59 2.09 0.68 0.27 1.79 0.51 6.35   Profile curvature (100 m) 1.63 0.98 2.27 0.33 <0.01 4.76 2.57 8.83   Profile variation index 1.82 0.75 2.89 0.54 <0.01 6.76 2.41 18.99 Forward stepwise All domains Intercept 1 -9.94 -11.57 -8.31 0.83 <0.01 -- -- --   Intercept 2 -1.65 -2.24 -1.05 0.31 <0.01 -- -- --   Intercept 3 5.27 4.43 6.12 0.43 <0.01 -- -- -- 1*  Area -7.81 -8.93 -6.70 0.57 <0.01 <0.001 <0.001 <0.001 2  Area 2  -0.93 -1.33 -0.54 0.20 <0.01 -- -- -- 3  Lithology = limestone -0.87 -1.23 -0.51 0.18 <0.01 9.92 4.82 20.43 4  Profile variation index 0.30 0.04 0.57 0.14 0.03 2.73 1.86 3.99 5  SL/k index 2   0.78 0.36 1.21 0.22 0.00 -- -- -- 6  SL/k index -0.23 -0.43 -0.02 0.11 0.03 1.97 1.36 2.86 7  Basin slope 2  -0.36 -0.95 0.22 0.30 0.22 3.52 1.74 7.15 8  Area × basin slope  -0.88 -1.36 -0.39 0.25 0.00 -- -- -- 9  Reach slope 2   2.29 1.57 3.02 0.37 <0.01 0.46 0.26 0.82 10  Area × reach slope 1.00 0.62 1.38 0.19 <0.01 -- -- -- *order in which variables were selected for final model based on score of the chi-square statistic 61  For the lithology predictor, transects in basins with erosion-resistant conglomerate with deep till (Cadomin) were coded to 1, and transects in the younger shale/sandstone formations were coded to 0.  The positive value for the lithology parameter estimate (1.82) indicates a greater probability of being in a lower order class when lithology equals 1 (Table 3-4). Cadomin lithology was associated with an increase in odds of being in a lower ordered domain class by a factor of 6.20.  The 95% confidence interval did not include 1 (which is the value specified by the null hypothesis) and hence the odds are significantly different from 1 (i.e. 50/50).  Along the Cadomin profiles, this translates to transition locations located further down the profile.  In comparison to the Area model, the Area + lithology model reduced the random effects estimate for Transects 4 and 5 located in the Cadomin region (Figure 3-7); the Area + profile complexity index model reduced the absolute value of the random effects estimate for Transect 2 – the profile with the least complexity that most closely resembles a graded profile (Figure 3-7); and the combined model showed reduced absolute value for all random effects.  Figure 3-7. Random effect estimates by transect for the four all domains models including: area; area + lithology; area + profile complexity index; and area + profile complexity index + lithology.  Accounting for spatial autocorrelation improved the accuracy of the best random effects model from 86 to 88% (Figure 3-6a).  There were no differences between fixed effects 62  parameter estimates for the model with no random effects (Table 3-4) and the model with random effects and SAA; however, the significance level of the profile complexity index estimate changed from <0.001 to 0.058.  These consistencies support the use of the no random effects model for regional mapping and indicate low probability of Type 1 error for its individual coefficient estimates.  The slope and magnitude of the predicted change in logit with “distance to source” differed between transects (Figure 3-8a), with transects 3a, 3b and 4 having a negative slope, 1 and 5 having a positive slope, and 2 and 3c having showing no effect.  These findings highlight the influence of individual long profile character on the spatial distribution of process domains.  Figure 3-8. Predicted changes in logit with spatial autocorrelation accounting by source distance by transect for (a) the best all domains model (area + lithology + profile complexity index) and (b) the best 1, 2, 3 domains model.  The values of the intercept estimates in the best model provide statistical evidence for distinct process domain classes.  Specifically, there is no overlap in the confidence intervals for the three intercepts in the best model (Table 3-4).  For example, Intercept 1 (hillslope domain 63  regression equation) has a value of -11.37 and associated confidence intervals of -13.01 and - 9.73.  There is no overlap with the confidence intervals for intercept 2 (swale and hillslope regression equation). 3.4.4 Relations between predictors and hillslope, swale and colluvial channel domains The three steps used to identify the strongest relations between the seven candidate predictors and the occurrence of the process domain classes 1, 2 and 3 were: (1) sensitivity analysis to determine the best length for measuring profile and planform curvature; (2) comparison of models containing only drainage area with models containing drainage area and a single candidate predictor; and (3) comparison of models containing only drainage area with models containing a combination of two or three candidate predictors.  This iterative strategy was intended to highlight the most parsimonious model, thereby avoiding over-fitting which can result in identification of spurious effects [Burnham and Anderson, 2002]. All curvature measures included in the sensitivity analysis had parameter estimates that were different from 0 (P <0.05); however, model performance was variable.  The optimal length for measurement of both planform and profile curvature was 100 m, but increases in prediction accuracy were largely limited to the swale domain (Figure 3-9a and b).  For regional scale mapping, the “planform 100 m + area” model had a lower accuracy than the “area” model, and the “profile 100 m + area” model had no change in accuracy.  Accuracy of transect-scale predictions improved in both cases. 64   Figure 3-9. Plots of accuracy of process domain predictions with standard error by (a) length of planform curvature measurement, and (b) length of profile curvature measurement.  Length of 0 m = “area” only model with no curvature measurement.  All of the fixed effects parameter estimates for the various “area + one predictor” models were different from 0 (P < 0.05); however, for regional scale mapping, this did not translate to improved accuracy over the “area” model.  Accuracy of transect scale predictions improved in all SAA models (complexity = 5) and the “area + profile 100 m” model had the greatest improvement. The “area + lithology + profile 100 m + profile complexity” model had the highest accuracy for regional scale mapping, with a 3% improvement over the “area” model (Figure 3-10a).  Penalty was reduced for the top three combinations (Figure 3-10b).  In a comparison of the global model and the “area” model, accuracy of transect scale predictions increased from 84 to 90% for the random effects model and from 94 to 95% for the random effects + SAA model. These findings highlight the importance of transect scale factors and SAA for explaining spatial distribution of process domains. 65   Figure 3-10. Plots of (a) accuracy and (b) penalty with standard error by random effects type for four models predicting process domains including hillslopes, swales and colluvial channels. Random effects type 1 = no random effects, type 2 = random effects with population average predictions (PAP), type 3 = random effects and spatial autocorrelation accounting (SAA) with PAP, type 4 = random effects with transect specific predictions (TSP), type 5 = random effects and SAA with TSP.  The fixed effects parameters had similar slopes to those that were shared with the best “all domains” model (Table 3-4).  Lithology was included in all top 3 models (Figure 3-10a); however, there was no evidence that the parameter estimate differed from 0 (P = 0.27), and the confidence interval for the odds ratio estimate included “1”, indicating no evidence that the odds associated with a shift to Cadomin lithology is different than 1 (Table 3-4).  Positive values for profile curvature indicate increasing downslope gradient and acceleration of flow [Zevenbergen and Thorne, 1987].  The positive values for the regression coefficient for profile curvature measured over 100 m indicate that higher values of curvature are associated with an increased probability of being in a lower ordered domain class.  These results are consistent with those from the sensitivity analysis where accuracy improvements were limited to the swale domain. All parameter estimates in the “random effects and SAA” version of the best model fell within 66  the confidence intervals associated with the “no random effects” version model; however, the significance level associated with lithology and profile complexity index changed from 0.27 to 0.49 and <0.001 to 0.65, respectively.  These results support the use of the “no random effects” version as the best model remains for regional mapping, but warrant a closer look at the effects of SAA on model performance. Accounting for spatial autocorrelation improved the accuracy of the best model from 90 to 95% (Figure 3-10).  The slope and magnitude of the predicted change in logit with distance to source differed between transects (Figure 3-8b), with Transects 1, 3a, and 4 having a negative slope, and 2, 3c, and 5 having a positive slope.  These findings highlight the influence of individual long profile character on spatial distribution of process domains. 3.4.5 Regional scale mapping The results from the analyses to identify relations between predictors and process domain class indicated that fixed-effects models were more accurate for regional mapping than the mixed-effects models and regression coefficients were similar in both cases.  Therefore, analyses were limited to fixed-effects models in this section.  The 15 variables evaluated using ordinal logistic regression with the stepwise selection procedure included: (1) area, (2) area 2 , (3) reach slope, (4) reach slope 2 , (5) mean basin slope from 2 m grid, (6) mean basin slope 2 from 2 m grid, (7) mean basin slope from 6 m grid, (8) mean basin slope 2 from 6 m grid, (9) SL/k index, (10) SL/k index 2 , (11) area × reach slope, (12) area × mean basin slope, (13) area × SL/k index, (14) lithology, and (15) profile complexity index.  To determine the best grid cell size (2 m versus 6 m) for measuring mean basin slope, separate models were completed for each. There were no differences in model coefficients or performance for the 2 m versus the 6 m grid cell sizes.  The complete model included 10 predictors with only basin slope, reach slope and area × SL/k index removed (Table 3-4).  In the stepwise procedure, predictors are added to the model in the order of importance based on their score Chi-square statistic.  The first four predictors were the same predictors identified through the logit plot preview exercise and their odds ratios were similar.  This congruence corroborates the results in both cases.  When SL/k index is increased by one standard deviation, the odds of being in a lower ordered category increase by a factor of 1.97 (Table 3-4).  This indicates that higher SL/k values are associated with larger drainage areas at thresholds between process domain classes.  Higher SL/k index values are indicative of greater stream power; hence this relation is contrary to what is expected. 67  When mean basin slope is increased by one standard deviation, the odds of being in a lower ordered category increase by a factor of 3.52.  This is also contrary to what is expected.  When reach slope increases by one standard deviation, the odds of being in a lower ordered category decrease by a factor of 0.46.  From a physical process perspective, stream slope is a determinant of stream power, so this pattern is expected but is contrary to the trend from the cumulative logit plot (Figure 3-5c).  These three outcomes are difficult to explain, but may be a function of the inherent relations in the set of geomorphic variables and dominance of drainage area for explaining process domain transitions. The advantage of the more complex model is most apparent when comparing the maps of predicted process domains from the best model from the logit preview exercise (Figure 3-11a) to the completed model from the stepwise selection procedure (Figure 3-11b). The predicted process domain transition locations are as much as 275 m from the true transition locations in the former, while the latter are within 20 m.  The latter also predicts an isolated colluvial channel that transitions back to a swale.  This type of transition to a lower ordered domain was observed on Transect 3a near the apex of a fan where a continuous colluvial channel shifted to a discontinuous channel.  However, the isolated colluvial channel is predicted where SL/k index values indicate under-steepened reaches (Figure 3-2) to indicate that the low reach slope and SL/k index values may be incorrectly influencing the prediction. 68   Figure 3-11. Maps of Basin 1 with (a) predicted process domains from the best model from the logit preview exercise, and (b) predicted process domains from model identified with the stepwise selection procedure. The test scores for the proportional odds assumption from the best relation model (area + area 2  + lithology + profile variation index) and the regional model were both less than 0.01 to indicate rejection of the models.  This test will reject the null hypothesis more often than warranted, especially with many predictors and large sample sizes [Allison, 1999].  Evidence contrary to the test score to suggest that these models are valid includes: high accuracy scores and low penalty scores from training and validation data sets; maps of predicted process domains that resemble true distributions; similarities between trends in cumulative logit plots and regression coefficients; and similar regression coefficient values in different model configurations. 69  3.5 Conclusion This study involved quantitative analyses of the spatial organization of process domains in a landscape with a complex litho-structural setting that has been affected by multiple ice sheet advances.  The complex geomorphology precluded the use of established approaches, such as area-slope plots or slope-dependent area thresholds, to identify process domain boundaries and predict the extent of the channel network.  As an alternative, ordinal logistic regression models were used to (1) establish relations between individual predictors and process domains, and (2) produce regional scale maps of process domain distribution.  Models that contained only drainage area were 82% accurate and this improved by up to 5% when other reach and transect- scale predictors were added.  Reach slope at process domain thresholds was highly variable between individual profiles; additional predictors were required to explain transect specific variations in domain thresholds.  Models for a limited set of domains that incorporated transect- scale character through the use of random effects and SAA had accuracies that were up to10% higher than the drainage area only model.  These results indicate that in complex landscapes, thresholds between process domains can be explained by exploring variations along individual profiles.  Nonlinear landscape-scale relations can be further quantified by using models that include quadratic and interaction terms.  A model of this type, developed for regional scale mapping, was 95% accurate.  Given the statistical importance of data points around transitions, future field surveys should focus on identifying transitions along individual profiles with multiple profiles randomly distributed across the study area.  Ordinal logistic regression models that account for mixed effects and spatial autocorrelation would be well-suited for such a design.   70  Chapter 4 Landscape scale modelling of process domains using LiDAR digital elevation models with factors representing reach, hillslope, and basin scale groundwater flow systems 4.1 Introduction Drainage basins are dynamic complex systems shaped by tectonic and geomorphic processes.  In addition to fluvial erosion, glaciation is another land-sculpting process that has prevailed across mid and high latitudes including most of North America and Europe north of 40°N [Ehlers and Gibbard, 2007].  Thus, the evolution of a landscape reflects the climatic history and the sequence of dominant erosional processes that have acted over time [Montgomery and Foufoula-Georgiou, 1993; Brardinoni and Hassan, 2006].  Theories on the connections between hydrogeomorphic processes and the evolution of drainage basin topography have continually developed since an original model was proposed.  According to Horton [1945], streams develop from the point where infiltration excess overland flow first exceeds a critical length, thus creating a process-based division of the landscape into hillslopes and channels.  Horton suggested that valley and stream development occur together, as evidenced by the development of concave topography.  Subsequent field studies of the hydrogeomorphic processes near the hillslope-stream transition have revealed an upper network of unchannelled valleys or swales, and introduced the importance of episodic events, some of which are driven by water moving through subsurface pathways rather than overland flow [see reviews by Dietrich and Dunne, 1993; Benda et al., 2005; Hassan et al., 2005].  In arid New Mexico, the threshold location between hillslopes and channels is not stationary; rather, in response to ground cover modifications, overland flow can form a discontinuous channel that appears as a sequence of headcuts interspersed with vegetated ground [Leopold and Miller, 1956].  Sections of a discontinuous channel coalesce during extreme rainfall events as individual headcuts advance.  Retreat of the channel head in such environments is a much slower process. In forested hillslopes in Japan, convergence of sub-surface storm water drives shallow debris flows that originate in elongated spoon-shaped features at valley heads, also called zero-order basins [Tsukamoto et al., 1982].  In the Oregon Coast Range, the recurrence interval of landslides was 6,000, 1,500 and 750 years from zero-order, 1 st  order and 2 nd  order basins respectively [Benda and Dunne, 1987].   Streams can be distinguished into colluvial or fluvial channels based on whether they develop in material conveyed from adjacent hillslopes, or in 71  stream deposited material [Hassan et al., 2005].  Colluvial channels display characteristics including non-alluvial boundaries, complex flow structures, and periodic evacuation of material by landsliding [Benda, 1990; Rice and Church, 1996].  A colluvial channel transitions to a fluvial channel where shear stress exceeds the threshold required to erode stream banks and create regular bed forms.  Thus, evidence gathered during field research of hydrogeomorphic processes indicates four process domains, or regions characterized by distinct erosional rates and mechanisms, as hillslopes, swales, colluvial channels, and fluvial channels. Knowledge of the spatial extent of process domains at the landscape scale is critical not only for predicting the spatial organization of hillslope and channel characteristics, particularly the locations of domain transitions, but also for understanding landscape evolution and the future response of a landscape to climate shifts, disturbances including fire, and land use changes.  Domain transitions can be anchored to features associated with jumps in shear stress, such as those that occur in over-steepened reaches at hanging valley outlets [Brardinoni and Hassan, 2006], and confluences.  In other cases, where longitudinal profiles are graded and drainage area increases in small increments, transitions may be more transient and subject to headward extension.  Thus knowledge of the whether a transition is anchored or more transient can contribute towards predicting landscape response. When hydrogeomorphic processes act continuously on a landscape over a period of time, beginning with the exposure of a new surface, each domain leaves a unique erosion signature that can be detected using slope-area plots [Montgomery and Foufoula-Georgiou, 1993].  There are three limitations that complicate the use of this approach for delineating process domains across landscapes.  First, while slope-area plots can be used to estimate hillslope-valley and colluvial channel-fluvial channel transitions, the dynamic nature of the channel head requires acquisition of field data to establish channel network extent.  Second, specific thresholds will vary across landscapes and between watersheds, thus thresholds can only be developed for individual basins.  Third, previous land-sculpting processes, including glaciations, may control channel network architecture.  In fact, GIS-based slope-area plots do not discriminate field- based process domains where glaciations have altered the surface; rather, contemporary hillslope and fluvial processes are superimposed upon a landscape characterized by a hierarchy of glacial troughs [Brardinoni and Hassan, 2006].  The more a channel deviates from the idealized graded (concave-up) profile, the greater the number of partitions required during the 72  application of downstream scaling relations.  This inherent complexity and the lack of resolution in commonly available 25 m DEMs preclude the application of theoretically-based remote predictions of process domains in glaciated mountain basins. The goal of this project was to develop a landscape-scale model to predict the spatial distribution of process domains in a region with complex terrain owing to repeated glaciations and foothills topography.  Strategies were required to address the limitations identified by Montgomery and Foufoula-Georgiou [1993], and complexity and DEM resolution obstacles encountered by Brardinoni and Hassan [2006].  Two general approaches for spatial modelling of erosion processes include the classic deterministic method [e.g., Montgomery and Foufoula- Georgiou, 1993; Tarboton, 1997], and a more empirical approach that assigns domain class using a probabilistic manner [e.g., Heine et al., 2004].  Glacial-induced landscape complexity precluded the use of deterministic models [Brardinoni and Hassan, 2006]; however, identifying relations through such an empirical approach is a reasonable intermediate step in the development of deterministic models [Dietrich and Dunne, 1993].   Given the complex topography where multiple factors influence spatial distribution of erosion processes, the probabilistic approach was well-suited for this study. Hydrologic classification can promote the transfer of information on controls of hydrologic systems from the basin scale, up to landscape, continental and global scales [Buttle et al., 2009].  To address drainage basin variability across a complex landscape and facilitate the transfer of information to other regions and across broader scales, this project included assigning a classification to catchment units.  The classification was also used to test the hypothesis that the accuracy of empirical process-domain models is better when the landscape is stratified into hydrologic units. Local discharge of water and local slope determine whether diffusive or concentrative erosional forces act upon a landscape and to what degree [Smith and Bretherton, 1972]. Discharge and drainage area are typically related through a power function [Leopold et al., 1964]; therefore, models of erosional processes usually substitute drainage area for local discharge of water.  Transitions between process domains are limited to headwater portions of the drainage network where discharge values remain relatively small (e.g. Trotter [2000] estimated bankfull flows required to create pool morphologies is ranges between 0.16 and 1.09 m 3 s-1).  In these headwater positions, landforms and surficial material shifts that promote either 73  groundwater downwelling or upwelling can exert strong local influence on streamflow volumes, which translate to changes in shear stress, potentially causing a change in process domain.   In heterogeneous landscapes where surface discharge is difficult to explain using models limited to drainage area and local slope, additional factors can improve predictions [Lindsay et al., 2004]. Models that consider such factors are required to improve the prediction of the spatial distribution of erosion processes in complex landscapes. Groundwater flow systems operate across the reach, hillslope, and drainage basin scales. Each of these systems has the potential to sufficiently affect discharge–drainage area relations and alter the location of process domain transitions.  At the reach scale, the hyporheic flux is determined by spatial changes in the energy head gradient, the cross-sectional area of alluvium, and changes in hydraulic conductivity [Tonina and Buffington, 2009].  Channel confinement typically decreases down the length of a channel profile [Montgomery and Buffington, 1997]; thus, the corresponding increase in alluvial area will promote downwelling (Figure 4-1a).  This downstream pattern is disrupted at features including knickpoints and hanging valley outlets where local slope increases.  A reduction in alluvial area in the transition to such over-steepened reaches will promote upwelling.   Hyporheic fluxes in headwater positions require further consideration as a result of a more extensive drainage network and its associated alluvial deposits that developed during the period of maximum landscape instability following de- glaciation between 13,000 and 11, 000 years BP [Hickman and Schweger, 1991].   A well- developed system of swales are now super-imposed upon the formative channel network, with channel heads located between 200 and 300 m downslope from the origin point of zero-order basins [McCleary and Hassan, 2008].  As a result, the material that underlies the contemporary swale and seepage channels may include alluvial deposits with a high capacity for hydraulic conductance.  Thus, in the study area, changes in confinement along swale, seepage channels, and fluvial channels may affect the hyporheic flux, and domain transition locations. At the hillslope scale, the degree of upwelling at low topographic positions, including those associated with channel network, will varying with adjacent hillslope configuration [Toth, 1963].  Hillslope configuration is very difficult to quantify and for any given channel segment, may differ from one valley side to the other.  However, groundwater discharge features located along lower side slopes and in depressional areas, as indicated by soil moisture status, can provide information on upwelling quantity [Dumanski et al., 1972], thus providing a surrogate 74  for hillslope scale upwelling (Figure 4-1b).  Wetland features such as regions of saturated flow are often included in spatial hydrologic models for low relief regions including the Boreal Plain [Devito et al., 2005b; Creed et al., 2008; Sass and Creed, 2008; Clark et al., 2009], and Boreal Shield [Lindsay et al., 2004]. At the drainage basin scale, areas with negligible local relief and limited general slope are not conducive to the development of local flow systems [Toth, 1963].  Local flow systems convey most of the water that drives erosion processes in headwater basins.  Thus, contributions to local flow are expected to increase with general slope of the terrain (Figure 4-1c).  75   Figure 4-1. Groundwater flow processes with potential to affect discharge – drainage area relations at process domain transitions operate at: (a) the reach scale where spatial changes in alluvial area down the longitudinal profile can create hyporheic flow [Tonina and Buffington, 2009]; (b) the hillslope scale where the upwelling volume is determined by hillslope configuration [Toth, 1963], reflected by soil moisture status; and (c) where contributions to local flow systems increase with general relief at the basin scale [Toth, 1963].    76  DEM data at a 25 m resolution are not suitable for delineating process domains that can be detected in the field [Brardinoni and Hassan, 2006].  Prior to 2008, the availability of high resolution DEM data within the study area was limited to specific drainage basins [i.e., Chapter 3].  In 2009, LiDAR data became available for landscape scale investigations.  Thus this project investigated the use of high resolution DEM data for delineating process domains at the landscape scale. A landscape-scale model to predict the spatial distribution of process domains in a region with complex terrain was developed in a sequence.  First, the landscape was classified into units based on topography, hydrologic connectivity, and surface type.  Field determinations of process domain were obtained across the range of slope and drainage area combinations within each type of landscape unit.  A preliminary drainage network was developed from a LiDAR DEM.  A number of reach-scale topographic variables obtained from the network were combined with remotely sensed measures of soil moisture to predict the probability of occurrence of each process domain category. 4.2 Study area 4.2.1 Hinton region The study area, covering approximately 10,000 km 2  of west-central Alberta (Figure 4-2), can be divided into two regions based on bedrock structure: the Interior Plains and the Foothills [Dumanski et al., 1972].  The Interior Plains are underlain by weakly consolidated beds of sandstone and siltstone from the Paskapoo formation of the Tertiary age, which are mainly flat lying.  The Foothills region is composed of a thick sequence of sedimentary rocks mainly belonging to the Brazeau formation of the Cretaceous and Tertiary ages.  Thrust faults have formed parallel rows of foothill ridges that run from southeast to northwest, dipping to the southwest and striking to the northeast.  Erosion has removed the shale and siltstone leaving erosion-resistant sandstones and conglomerates to form the ridge crests.  Valleys are worn into belts of more erodible bedrock and along fault zones to form a trellis drainage pattern. 77   Figure 4-2. Study area maps including (a) map of Canada, with extent rectangle for (b) west central Alberta, a transition area between Rocky Mountains and boreal plain, with extent rectangle for (c) shaded relief map of study area, including extents of DTW index and moisture regime data, and field sampling locations. Cordilleran ice sheets repeatedly advanced in an easterly direction through main valleys of the Rocky Mountains before fanning out over the region.  Continental ice sheets advanced into the eastern portion of the study area.  These advances reorganized the fluvial landscape 78  within the foothills as evidenced by common features including buried valleys, where original valleys are filled with glacial drift up to 290 m in depth [Roed, 1968].  Surficial material depths and type are highly variable.  Glacial lakes, evidenced by lacustrine and deltaic deposits, formed along the western edge of the Continental ice sheet.  Calcareous tills define the eastern extent of the Cordilleran ice sheet, while non-calcareous tills define the western limit of the Continental ice sheet.  Fluvial terraces are common along major rivers.  Organic soils have developed in low-relief areas and toe slope positions.  Elevations range from 1,000 to 1,600 m.  The diversity of bedrock and surficial materials in the study area presents challenges for modelling erosion processes that are typically dependent upon consistent local discharge – drainage area relations. In early nonglacial times, the drainage network in the glaciated Rocky Mountain Foothills region produced more sediment and had a greater extent.  Using an analysis of lake sediments, Hickman and Schweger [1991] found that surface erosion rates were highest in the post-glacial Holocene between 13 000 and 11 000 years BP, when the arid cold climate limited vegetation to a discontinuous cover of Artemisia and Poaceae.  Landscape stabilization initiated near 11,000 years BP with increasing moisture and temperature that allowed the establishment of arboreal taxa including Populus, Betula and Salix.  Woodland and forest vegetation have maintained complete ground cover through to the present.  Surface erosion is no longer an active process due to high infiltration rates associated with thick litter and duff layers typical of forested foothills soils [Archibald et al., 1996].  In addition to a climate shift that promoted stabilization through increased ground cover, sediment sources – in the form of  landforms left unstable following ice retreat – were also exhausted during early nonglacial times between 13 000 and 10 500 BP [Jackson et al., 1982].  At another Foothills location approximately 100 km south of the study area, the estimated drainage density including well-formed unchannelled valleys, or swales, was 9.16 kmkm-2 in comparison to a density of 2.37 kmkm-2 for active channels [Chapter 2].  The system of swales, typical of the glaciated Foothills, may represent the extent of the drainage system associated either with the previously described period of maximum landscape instability or the period of elevated runoff associated with glacial recession. Consistent drainage area-discharge relations form the foundation of slope-area plots and other remote methods to delineate process domains.  For the study area, data to confirm these relations within small watersheds are limited – only two of the 13 stations with more than nine 79  years of data are located in basins with a drainage area less than 10 km 2 [Veldman, 1997].  This scarcity of streamflow data for drainage basins near the seepage-fluvial channel transition is typical of most regions.  The predictions of peak flow (a function of drainage area) for Whiskeyjack Creek (drainage area = 3.13 km 2 ) were significant lower than the larger drainage basins in the region, and these differences were attributed to specific topographic conditions [Veldman, 1997].   Of the total precipitation falling within the adjacent Cache Percotte basin (drainage area = 7.15 km 2 ), 66 percent was lost to evapotranspiration, 16-25 percent formed streamflow, and 9-18 percent was lost to basin leakage and underflow [unpublished Research Council of Alberta data  presented in Dumanski et al., 1972].  The water budget allocation is expected to change with topography; thus, drainage basin morphology is an important consideration when developing models of erosion process domains drainage area-local discharge relations. Lodgepole pine and spruce forests cover uplands and either black spruce-tamarack forests or shrubs dominate wet areas.  Approximately 1% of the land base (120 km 2 )  is harvested annually and the first harvest rotation will be completed by 2040.  An extensive transportation network including roads, railways, pipelines, and petroleum exploration trails has been developed.  Coal mining and petroleum extraction also occur in the region. 4.2.2 LiDAR data The LIDAR data were collected in August 2007 by Airborne Imaging Inc., for Alberta Sustainable Resource Development, Resource Information Management Branch. A helicopter flying 1,200 m above ground level at a forward speed of 160 knots was equipped with an Optech ALTM 3100, 543 LIDAR system.  A 1 m grid pattern point spacing was used with a density of 1.5 pointsm-2.  Horizontal accuracy was 0.45 m and vertical accuracy was 0.30 m. 4.3 Materials and methods This project was completed in five steps, which are detailed in the following sections. As an overview, the five steps were: 1) The 10, 000 km 2  study area was divided into 2,584 hydrological response units (HRUs) with an average area of 4.3 km 2 . Using a fuzzy clustering procedure, the HRUs were classified into six categories based on lacustrine soils (present/absent), average basin relief, and extent of wetlands.  2) A preliminary drainage network was extracted from a high-resolution DEM. The network was sub-divided into reaches, which were characterized by a suite of topographic predictors.  3) The study area was stratified 80  using HRU type, reach drainage area and reach slope.  Geo-referenced ground surveys were completed at 718 reaches across 54 strata.  Each surveyed reach was assigned a process domain class as hillslope, swale, seepage channel, or fluvial channel.  4) Three different types of logistic regression were used to fit models for predicting process domain class based on the reach-scale topographic and soil wetness predictors.  5) The best models were applied to all reaches within the preliminary network. 4.3.1 Basin similarity analysis Basin classification is prerequisite for the study of spatial patterns of erosion processes across landscapes.  For this project, maximum basin size was set to 10 km 2 , or approximately 0.1 % of the 10,000 km 2  study area.  This resolution overlaps well with size of stream that provide important fish habitat [McCleary and Hassan, 2008], and was detailed enough to recognize emerging patterns of the landscape.  More detailed modelling of erosion processes using LiDAR data occurred at the reach scale within these basins.  Polygons for classification included 2 nd  order watersheds, plus 1 st  order watersheds and other slope units that flow directly into larger (i.e. 3 rd  order +) watercourses.   With the inclusion of slope units with no outlet streams, not all polygons were true drainage basins; therefore, all polygonal features were more accurately named hydrologic response units (HRUs), consistent with Devito et al. [2005a]. In comparison to other modelling procedures where relations are established between one dependent variable and a number of independent variables, classification procedures generally involve grouping samples based on similarities among a set of independent variables; there is no dependent variable involved.  For watershed classification, a scheme has been suggested to select the variables important for classification based on the factors that influence water fluxes.  The three main factors influencing water movement are topography, hydrologic connectivity and surface type, hence meaningful basin classification systems should be developed with a set of parameters that are carefully selected to represent each of these factors [Buttle, 2006].   Mean basin slope, wetland extent, and presence of lacustrine deposits were used to represent these three factors respectively (Table 4-1). 81   Table 4-1. Variables selected to represent the three factors that influence water movement with hydrologic considerations and methods for each. Factor and variable Hydrologic considerations and methods 1. Topography: mean basin slope Topography is incorporated into hydrologic and geomorphic models by specific catchment area (drainage area per unit contour length) [Tarboton, 1997], or mean basin slope (slope averaged for flow-accumulation area) [Heine et al., 2004].  Most drainage basin morphometrics that describe topography are highly correlated [Lindsay et al., 2004]. Therefore, mean basin slope was selected because it could be rapidly calculated across the numerous HRUs using standard raster GIS software.  High resolution DEM data were not warranted for this basin-scale summary statistic; hence the 25 m DEM was used.  Specifically, slope (%) was derived for each cell in the DEM for the study area. Then, the average slope of all cells within each HRU was calculated. 2. Hydrologic connectivity: wetlands extent An ecological land classification [Downing, 1998], which was completed between 1994 and 2004, contained the most accurate information on distribution of wetlands for the study area.  The dataset was provided by Hinton Wood Products, a division of West Fraser Mills Ltd.  The relevant details of this classification are as follows.  Aerial photographs at 1:15,000 scale were acquired for the study area.  This resolution exceeded the 1:20,000 scale recommended for vegetation-based inventories [Alberta Environmental Protection, 1996].  The study area was manually divided into 285,096 ecological unit polygons (average area = 0.037 km 2 ) that shared similar plant communities, site characteristics (i.e., moisture regime, nutrient regime, topographic position, slope, and aspect), and soil characteristics (i.e., organic thickness, surface texture, depth to mottles, drainage, parent material).   The ecosite type – ecological units that develop under similar climate, moisture and nutrient regime – was assigned to each polygon based on Beckingham et al. [1996].  Ground truthing was completed at 22 % of the polygons using a visual survey of indicator plant species to confirm ecosite type (n = 33,162), or an ecological sample that also included soil moisture regime and texture determination (n = 29,273).  Each ecosite type is characterized by a moisture-nutrient regime [Beckingham et al., 1996].  The eight candidate moisture regime classes are xeric (2), subxeric (3), submesic (4), mesic (5), subhygric (6), hygric (8), or subhydric (9).  There are five nutrient regime classes ranging from very poor to very rich.  Using the moisture and nutrient regime from the ecosite classification, ecological units can also be classified according to wetland type, as meadows, swamps, bogs, fens and marshes [National Wetlands Working Group, 1997].  For the basin similarity analysis, ecosite polygons that met the specific criteria for these wetland types were used to calculate the wetland extent.  The Pearson correlation coefficient of -0.59 between mean basin slope and wetland extent confirmed a negative relation.  Topographic position is an important factor determining wetland occurrence [Toth, 1963; Dumanski et al., 1972], and may explain additional variation. 82  Factor and variable Hydrologic considerations and methods 3. Surface type: presence of lacustrine deposits Depth and type of surficial deposits are controlling factors in the hydrology of boreal plain catchments [Devito et al., 2005b].  In the foothills, the depth of surficial glacial deposits in many locations is negatively correlated with surface slope [Stelfox, 1981]; therefore, depth is implicitly represented by the topographic descriptor mean basin slope. Surficial materials in the study area, originally mapped by Dumanski et al. [1972], were grouped into three types, including soils developed on till materials, soils with peat accumulations of various thickness, and soils developed on lacustrine or water-laid material.  Soils developed on till material in upland areas do not display any unique properties with respect to water movement.  Soils with peat accumulations were described by the wetland extent.  Soils developed on lacustrine or deltaic materials have high surface erosion potential [Dumanski et al., 1972].  These materials are also landslide prone because of high cohesion and textural shifts at interfaces between layers that retard downward movement of water [Chatwin et al., 1994].  In a sediment budget study of three small catchments in the study area (Figure 4-2), unpublished data indicate that bedload sediment yield was at least an order of magnitude higher in a small basin where lacustrine and deltaic silts were the locus of landslides, than in two till dominated basins with higher relief where no landslides were observed (Table 4-2).  Landslides in lacustrine silts were a primary source of sediment in major drainage basins of the glaciated regions on the east side of the Rocky Mountains, including the Peace and Liard Rivers [Church et al., 1989] located immediately north of the study area.  Thus, basins that contain lacustrine deposits warrant special consideration because they are expected to yield a disproportionately high quantity of sediment in comparison to basins that do not contain such deposits.  For each HRU, lacustrine soils were noted as present or absent.  Table 4-2. Bedload sediment yield May – August, 2006 from three sediment budget watersheds in the study area. Basin Name Drainage area (km 2 ) Relief (% slope) Parent material Bedload yield 1 (t•km-2) Athabasca Tower 3.7 19 Till / weathered bedrock 0.01 Moberly Tower 1.1 24 Till / weathered bedrock 0.32 Wildhay Oxbow 1.7 10 Lacustrine / deltaic deposits 2.29 2  1 Total mass of material collected with correction factor based on ratio of channel width to width of bedload trap opening. 2 This trap was filled to capacity on five occasions, hence actual yield is much higher. 83  Fuzzy clustering techniques are appropriate for assembling patterns into groups when working with systems that lack discrete boundaries [Zadeh, 1965].   Rather than assigning each object in a collection to a specific group, a membership function associated with each group is used to assign a membership grade from 0 to 1 for each object.  Magnitude of the membership grade provides an indication of the strength of the relation between each object and its respective group.  Fuzzy clustering algorithms have been used for spatial classifications in disciplines including climatology [McBratney and Moore, 1985], biogeography [Dragicevic, 2010],  and geomorphology [Cheong, 1998; Gorsevski et al., 2009].  Separate clustering exercises were completed for HRUs with lacustrine soils present versus absent.  Mean slope and wetland extent were standardized prior to analysis.  Fuzzy c-means clustering was completed using the cmeans function in the e1071 package of Dimitriadou et al. [2009].  This function requires that the user input the target number of clusters for the analysis.  To determine the optimal number of clusters, the cmeans function was executed five times, starting with two clusters and increasing by one cluster to a maximum of six clusters.  Outputs from these iterations include within-groups sum of squares, which were plotted as a function of the number of clusters.  The graph was examined to determine where the gain in information associated each additional group was optimized. 4.3.2 Creating the preliminary drainage network The point LIDAR data were used to build raster DEMs with a 2 m resolution, which was consistent with other process domain studies employing LIDAR data [Tarolli and Fontana, 2009].  The preliminary drainage network and stream reach attributes were created using two FORTRAN computer programs (Bld_grds and Netrace) that were developed for DEM analysis [Miller, 2003; Clarke et al., 2008].  The programs share a set of approximately 70 user specified parameters designed to account for regional geographic factors that influence stream channels. Intersection points between the LiDAR-derived drainage network and elevated road beds posed a problem when trying to create a digital representation of the true drainage system.  Specific details on the problem and the steps that were implemented to resolve it are described herein. Road beds and ditches are ubiquitous features in managed landscapes that affect the downslope movement of water.  They are readily detectable with LiDAR DEMs and can influence flowpath routing at intersection locations [Schiess et al., 2004; Mouton, 2005; Murphy et al., 2008].  There are two basic approaches to predicting flowpath locations in landscapes 84  with road grades.  The first is to alter the LiDAR DEM to create „natural‟, pre-disturbance flowpaths [e.g., Mouton, 2005].  The second is to alter the DEM only at locations where drainage structures are present to represent actual flowpaths on the altered landscape [e.g., Murphy et al., 2008].  When developing a stream network from a LiDAR DEM in an area with existing road grades, it is important to state which approach will be used so that the resulting network can be properly verified in the field.  The differences between the pre-disturbance and contemporary flowpath maps can provide a measure of hydrological consequences of road grade construction [Murphy et al., 2008].   For this research, the first approach to predicting flowpath locations was used and anthropogenic alteration of flow routing was ignored.  Some domain boundaries may be influenced by road-relating flowpath diversion; however, considerable work would be required to accurately represent such features.  Water routing through swales and seepage channels will reflect the type of drainage structures at each road crossing, and such information is difficult to obtain.  Furthermore, drainage structure performance can be dependent upon maintenance, especially for small culverts that can be subject to debris and ice blockages.  An analysis of soil moisture conditions upstream and downstream of individual crossings with synthetic aperture radar imagery [e.g., Sass and Creed, 2008] may help to confirm the locations of actual flowpath diversion. In some cases, at the intersection point of the preliminary network and a road or railway grade, the flowpath was diverted into the drainage ditch on the upstream side of the grade rather than crossing the hump in the terrain and continuing down the natural swale.  To create a drainage system that more closely resembled the natural network, subroutines were created in the Bld_grds and Netrace FORTRAN programs to smooth the DEM and force the flowpath across the grade.  Specifically, an uncorrected channel network was created and intersected with the road and railroad network to create a point coverage representing crossing locations.  The 2 m DEM was smoothed for a 20 m radius around each crossing location and a corrected network was produced using the smoothed DEM. The study area boundary corresponded to the extent of the land where both LiDAR and wetland data were available.  To ensure that accurate measures of drainage area were available for all stream reaches within the study area, the preliminary drainage network extended beyond the original study area boundary to include the complete drainage basins of all tributaries, except those of the large rivers originating in Jasper National Park (Figure 4-2).  This expanded 85  area was divided into 36 individual catchments, each with an associated preliminary drainage network.  Average drainage density of the preliminary network across all catchments was 5.4 km/km 2 .  The extent of the network was purposefully over-predicted to provide DEM-traced channels that extended into the hillslope domain.  This overly-dense drainage network was well suited for subsequent filtering based on predicted process domain class.  Specifically, any stream reaches in the preliminary network that were assigned to hillslope domain during the modelling process were removed from the final regional drainage network.  The preliminary drainage network was comprised of 2,350,000 reaches averaging 36 m in length.  Reach length increased with drainage area. 4.3.3 Process domain field surveys Accessibility was a major consideration when establishing the sampling program.  Safety issues and costs associated with accessing remote regions of the study area precluded the use of a random sampling approach.  Therefore, to work within budgetary limitations and to meet safety requirements, a systematic approach that utilized the expansive road network was developed.  Stratified sampling was employed to capture the spatial variation in process domain occurrence.  Each of the six HRU types were divided into three area classes (<0.1 km 2 , 0.1 – 1.0 km 2 , and 1.0 to 10.0 km 2 ), and three reach slope classes, leaving 54 potential strata.  The reach slope class cut-offs in each HRU were adjusted so that three similar sized groups were represented.  The preliminary stream network was intersected with the road network to create point coverage of GPS locations for field surveys.  Surveys were conducted by travelling all roads in a target region that were navigable by a four wheel drive vehicle and assigning process domain classes at every road crossing encountered.  Headwater streams are disturbed by the clearing of vegetation from the road right-of-way, and also by the construction of the road crossing.  Therefore, all field surveys were conducted upstream from the section of the watercourse influenced by road construction.  Field surveys for process domain analysis were limited to lands suitable for industrial forestry with ecological land classification data where the road network is well developed.   In comparison, McCleary and Hassan [2008], also inventoried streams in steeper watersheds outside of lands suitable for forestry with lower levels of road development.  Given that the object of this project was to explore spatial distribution of process domains within the pre-disturbance channel network, each field survey location was examined for road related influence in upstream areas.  Sites with upstream flowpath alteration were flagged for consideration during the statistical analysis.  The limited extent of pristine 86  watersheds within the study area and the high costs associated with accessing them led to the exclusion of sites in road-less areas from the study design.  A tally by strata type for of all completed surveys was maintained and regions were selected and targeted on a systematic basis to provide a broad representation.  Over the course of two field seasons, data were obtained for 718 locations that were suitable for analysis. The sequence of process domain classes (i.e., hillslope, swale, seepage channel, and fluvial channel) was consistent with the classifications by Montgomery and Foufoula-Georgiou [1993], and Church [2002] with one exception – seepage channels were substituted for colluvial channels.  In glaciated areas in the Foothills and Interior Plains, landslides are limited to basins with lacustrine soils and sufficient relief [e.g., Hartman, 2005].  In such moderate and low relief terrain, seepage rather than colluvial processes transport hillslope sediment into channels [e.g., Terajima et al., 1997; Sayer et al., 2006; McCleary et al., 2008], and periodic upward migration of headcuts remove sediment accumulations [Chapter 2].  The specific field criteria for each class were as follows.  Hillslopes are areas of topographic divergence characterized by convex and planar slopes (Figure 4-3a).  Swales originate in areas of topographic convergence evidenced by concave slopes.  A swale may be dry or wet.  In a dry swale, the soil moisture regime in the depressions is similar to adjacent hillslopes (Figure 4-3b).  In other cases, the soil is wet due to seepage for at least a significant portion of the year (Figure 4-3c).  Discontinuous channels can occur in the transition zone between a swale and a seepage channel (Figure 4-3d). In steep terrain, they appear as a swale with a channel that is either migrating upstream through headward extension or in the recovery process from this type of event [Leopold et al., 1964]. Erosion typically initiates at a headcut and sediment is transported a short distance downstream to a deposition fan. Therefore, the bed of the channel may be interspersed with patches of vegetated ground.  Some distance below each swale, the competence of flowing water becomes sufficient to form a continuous seepage channel within definable banks (Figure 4-3e).  Field identification criteria for seepage channels include: a wide range of grain sizes in the channel bed; significant amounts of organic material, including randomly distributed wood and live roots, which exert a strong influence on sediment storage; and an overall high complexity due to lack of shear force required to create regular sequences of channel features.  The transition to fluvial channels is expected where an increase in competence of flowing water becomes sufficient to create typical small stream morphologies as described by Church [1992] (Figure 4-3f).  Example features that were used in the field to differentiate seepage and fluvial channels 87  included presence of head cuts, variable channel width, presence of undercut banks and organic bridges for seepage channels, and riffle-pool or cascade-pool sequences for fluvial channels (Figure 4-4). A key with eight criteria was used to ensure consistency when assigning domain as seepage versus fluvial (Table 4-3) and this was included as part of the field form (Appendix 2). A total of 57 fluvial channels were located on portions of the drainage network without obvious flowpath diversions at road crossings.  Gradients were less than 6% in all of these cases and averaged 2.6 % (18 locations < 2%,  31 locations 2–4%, and 8 locations 4–6% slope).  Steeper alluvial channel-reach morphologies, including the step-pool as described by Montgomery and Buffington [1997], were not encountered.  88   (a) Hillslope (b) Swale - dry  (c) Swale – wet (d) Discontinuous channel  (e) Seepage channel (f) Fluvial channel Figure 4-3. Example photographs of process domain classes.   89   a. Head cut  b. Variable width channel   d. Riffle – pool sequence c. Undercut banks and organic bridge Figure 4-4. Channel features used to differentiate seepage and fluvial channels.  90  Table 4-3. Key to seepage and fluvial erosion classes. Feature Number Seepage Channel features Fluvial channel features 1 Texture of fine material on channel bed is predominantly silt-sized or smaller (< 50% sand), with asilt or loam or silty loam texture a . Texture of fine material on channel bed is >50% sand (sand or silty sand or loamy sand) a . 2 Unconsolidated channel bed indicated by: thalweg boot penetration > 10cm b . Consolidated channel bed indicated by: thalweg boot penetration < 10 cm b . 3 No steps / riffles created by recently mobile gravel or cobbles c . Steps / riffles with regular spacing created by mobile gravel or cobbles c . 4 No pools present c . Pools present with regular spacing c . 5 Organic bridges present c . No organic bridges present c . 6 Head cuts present c . No head cuts present c . 7 Channel maximum width >3x the minimum width Channel maximum width <3x  the minimum width 8 Total undercut width > bank full width Total undercut width < bank    full width a  Use the “ribbon test” [Beckingham et al., 1996] to texture the fine channel bed material. A ribbon will form with < 50 % sand.  No ribbon or short thick ribbon will form with > 50 % sand. b  Sink boot into thalweg of watercourse and measure how far your heel sinks into the channel bed. c  See reference photos (Figure 4-4). 4.3.4 Variables used to predict process domain class All reaches within the preliminary drainage network were characterized by a set of topographic and soil descriptors (Table 4-4).  Of these, the SL/k index (detailed in Chapter 3), has had limited use  in other process domain classifications.  91  Table 4-4. Reach-scale predictors grouped by category. Predictor Scale Measurement method Topographic parameters 1. Drainage area Basin Drainage area (ha) above the downstream end of each reach calculated with FORTRAN program [Miller, 2003]. 2. Reach slope Reach Channel gradient (%) calculated for each channel pixel, then averaged over the length of the reach using a FORTRAN program [Miller, 2003]. 3. Mean basin slope  Basin scale using 2 m grid cells Hillslope gradient (%) averaged over the entire upstream contributing area for each reach using a FORTRAN program [Miller, 2003].  A measure of drainage-basin scale groundwater contributions to streamflow [Toth, 1963].  An important predictor of stream channel presence [Heine et al., 2004]. 4. SL/k index Reach Calculated with a FORTRAN program developed for this project. An indicator of reach-scale hyporheic flow, with increasing SL/k index values providing a surrogate for upwelling due to decreasing spatial extent of alluvial deposits. 5. Valley width index Reach  Calculated with FORTRAN program [Miller, 2003] as width of valley measured at a specified height above the channel equal to 5 times the predicted bankfull width, where bankfull width is a function of drainage area. Soil parameters 6. Soil moisture regime Reach The average soil moisture regime within a 10 m buffer centered on each reach.  Soil moisture regime was determined from a soil and vegetation based ecological land classification [Beckingham et al., 1996] completed for the study area in 2000. An indicator of hillslope scale groundwater contributions to streamflow [Toth, 1963]. 7. DTW index Reach Mean depth-to-water (DTW) index [Murphy et al., 2009] within a 10 m buffer around each channel segment.  The DTW index dataset included an associated flowpath layer.  This flowpath layer was divided into 20 m long channel segments and field sites were snapped to the nearest segment. 8. Lacustrine soil present Reach Presence or absence of lacustrine soil at location of field survey based on Dumanski et al. [1972].  It can be informative to plot the curved Hack profile (height on primary Y axis) and a step curve of the SL/k index (secondary Y axis).  To create the step-like appearance, the SL/k index value is held constant over the length of each reach [Chen et al., 2003].  For example, the 92  changes in gradient along the profiles of the hanging valley near Athabasca Tower (Figure 4-5a), and the intersection of the headwater stream with a fluvial terrace (Figure 4-5b) become more obvious (Figure 4-6).  Figure 4-5. Landscape patterns from glacial and fluvial erosion processes shown on 1:50,000 scale shaded relief maps with aerial photograph interpreted channels. (a) Athabasca Tower Hill with stream in hanging valley draining northwestward into abandoned glacial meltwater channel showing discharge (l/s) at six locations from June 19, 2008. (b) Wildhay River tributaries flowing southeastward from ridge to fluvial terrace (dashed white line) then to floodplain.  Figure 4-6. Hack profiles with step SL/k index curves for (a) stream flowing northwestward from Athabasca Tower Hill in hanging valley, and (b) Wildhay River tributary flowing southeastward from ridge to intersection with fluvial terrace then to floodplain.  SL/k = (reach slope × distance to source) / (rise to source / log (distance to source)). 93  Drainage area and reach slope were included as standard outputs from the FORTRAN programs [Miller, 2003; Clarke et al., 2008] and additional code was added these programs to calculate the mean basin slope and the SL/k index.  Two soil moisture descriptors that were available for the study area were used: the soil moisture regime from an ecological land classification [Beckingham et al., 1996], and the depth-to-water index derived from LiDAR data [Murphy et al., 2009]. The limitations of traditional photo-interpretation techniques for establishing soil moisture include the inferring of wetland status based on community structure of the forest canopy rather than hydrologic conditions, and the small scale of photographs (1:65,000 to 1:85,000) that are typically used [Clark et al., 2009].  These limitations were addressed by the extensive ground truthing that was completed during the ecological land classification [Downing, 1998], and the large scale of the photography (1:15,000) that was used.  To obtain the soil moisture regime data for process domain spatial analysis, each ecosite polygon was coded according to its moisture regime, then each reach in the preliminary network was buffered by 10 m and the average moisture regime was calculated for the buffered area associated with each reach. The depth-to-water (DTW) index provides a measure of how likely soils are to be saturated at their surface [Murphy et al., 2009].  DTW index data based on digitally processed LiDAR surface images of 1 m resolution were provided for a portion of the study area (Figure 4-2) by the Province of Alberta.  The DTW index dataset was produced in two steps: production of a surface water feature grid, and then combination of the surface water grid with a slope grid to produce the DTW index [Murphy et al., 2009].  Relevant elements in the process were as follows. A digital layer of predicted stream channels was produced using a 4 ha minimum upslope area for channel initiation.  Drainage was enforced across road grades where culverts and bridges were known to be located.   Flowpath diversion along ditches was permitted at stream crossings where no culvert was indicated.  This approach was intended to incorporate the hydrologic effects where drainage structures were not installed under graded surfaces, consistent with Murphy et al. [2008].   This predicted stream channel layer was then combined with a layer of known surface water features to produce a water features grid.  The water features grid and slope grid were inputs into the program to compute the DTW index.  A wetness index was generated as the raster output from the program, which was then simplified to show values of 94  consequence, specifically DTW of 1 = most wet and DTW of 4 = least wet, with a value of 0 assigned to all cells where soils are not susceptible to disturbances related to wetness. Limitations of the DTW index include limited ground truthing completed for the study area.  To obtain the DTW index for process domain spatial analysis, the average DTW index was calculated within a 10 m buffer associated with each reach. The soil moisture regime and DTW index data displayed different spatial patterns at the catchment and reach scales.  Measures of soil moisture regime were available for the entire landscape, whereas the DTW index was limited to topographic positions adjacent to flowpaths with all other locations being considered dry (Figure 4-7a and c).   Measures of soil moisture regime were defined for polygons with an average area of 3.7 ha, and hence typically did not vary with the 10 m buffer associated with each reach (Figure 4-7b); however, DTW index typically changed within the buffer as relief increased away from the flowpath line (Figure 4-7d). 95   Figure 4-7. Spatial patterns of soil moisture regime versus DTW index at the catchment scale (a) and (c), and the reach scale (b) and (d). 4.3.5 Modelling approach The response variable for modelling – process domain class – was a categorical variable of the ordinal type.  Three options for modelling categorical data include ordinal logistic regression, nominal logistic regression, and binary logistic regression.  Each modelling option in the sequence reduces the amount of information contained in the dependent variable, but provides increasing flexibility in the model parameters that are used.  The purpose of using three methods was to explore the tradeoffs associated with information content versus model parameter flexibility.  Models were developed using each of these procedures and results were 96  compared. Five aspects of the modelling are detailed in Appendix 2, including: (1) response variable types and options for modelling; (2) the logit model as the base model for categorical response variables; (3) cumulative logit plots for assessing relations between a response variable of ordinal class and continuous type predictors; (4) the ordinal logistic model as a special form of the logit model; and (5) the multinomial logistic regression model as a special form of the logit model. The analysis proceeded in four steps: (1) graphical exploration of the relations between process domain class and the topographic and soil-related controls for surveyed reaches; (2) development of regression relations using the three different logistic regression models; (3) evaluation of model performance with the addition of DTW index for a portion of the study area; and (4) evaluation of the binary logistic model for channel presence/absence developed for subsets of data from each basin type defined in the basin similarity analysis . Step (1) involved viewing cumulative logit plots prior to fitting regression models to screen out any parameters that did not meet the proportional odds assumptions for ordinal logistic regression.  In step (2), ordinal, multinomial, and binary logistic regression models were developed where candidate variables included all predictors, interaction terms, and quadratic terms.  Stepwise selection has been used in related hydrologic models to select the best model from the set of predictors [e.g., Lindsay et al., 2004], and was used in this study.  Ordinal, multinomial, and binary logistic regression analyses with stepwise selection were performed LOGISTIC procedure in SAS [Release 9.2; SAS Institute Inc., Cary NC]. The fixed effects covariates were log transformed as required to achieve a more normal distribution.  To remove scale effects and assist with model parameter and odds ratio interpretation, predictors were standardized to mean = 0 and standard deviation = 1.  This mean centering exercise was also performed because it stabilizes polynomial models. Cumulative logit (ordinal logistic) regression was the first type of model used to predict process domain class.  This type of model predicts the cumulative probability, Fij, that individual i is in the j th  category or lower, rather than being in a higher category.  The cumulative logit model comprises a set of J – 1 equations with different intercepts but common regression coefficients: 97  1 1 2 2     ... 1 ij i p p ij F log X X X F                     (4-1) where αi is the intercept for class i, β1 to βp are the fixed-effects logistic regression coefficients, and X1 to Xp are the predictor variables.  Each equation in the set is used to calculate a sequential dichotomous variable.  The use of different intercepts but identical regression coefficients is consistent with the assumption of proportional odds across each of the cumulative logit classes, which allows the information from the ordered sequence of process domain classes to be incorporated into the model outputs. The second type of model, multinomial logistic regression, does not assume proportional odds, thereby ignoring the information contained in the sequential relation of the process domain classes.  Nominal variables with more than two categories can be estimated by running a set of binary logit models.  Specifically, for the i th  category, where the baseline category is J, the multinomial logistic model is 1 1 2 2 ( )      ... ( ) i i i i ip p J p category log X X X p category                   (4-2) The number of non-redundant sets of solution equations is J – 1.  Each model in the solution set shares the same set of predictors (X1…Xp), but a different intercept (αi) and regression coefficients (β1…βp). In the third type of model, binary logistic regression, cumulative logits were modelled individually.  This approach considers the sequential information inherent in the ordered classes, but does not assume proportional odds, and hence allows for different intercepts, regression coefficients, and predictors in the models for each cumulative logit.  However, in contrast to ordinal and multinomial logistic regression, there is no objective method of assembling the results from a set of binary logistic models of cumulative logits into a single variable representing predicted process domain class.  The basic form of the logit model is: 1 1 2 2     ... 1 i p p i p log X X X p                 (4-3) 98  where α is the intercept.  Logistic regression is well suited to studies on rare events because the probability cut-off for establishing presence/absence varies with prevalence.  The cut-off point was empirically chosen such that the probability cut-off for presence/absence corresponded to the point where the model sensitivity (probability of correctly predicting a positive outcome) and model specificity (probability of correctly predicting a negative outcome) were both maximized [Santner and Duffy, 1989].  Using the classification table option in LOGISTIC procedure in SAS [Release 9.2; SAS Institute Inc., Cary NC], sensitivity and specificity were calculated at all values between (prevalence - 0.20) to (prevalence + 0.20) at every 0.01 interval. For management applications, correct identification of positive outcomes (i.e. channel present) is more important than correct identification of negative outcomes (i.e. channel not present). Therefore, for mapping purposes, the cut-off for presence/absence served as the threshold for high probability of occurrence.  A medium probability cut-off was set where model sensitivity was 95%.  Probabilities less than the medium cut-off were assigned a low probability.  Maps of high, medium, and low probability of occurrence were created using these criteria. Cross-validation was used to choose models.  Model training and validation datasets were created using the hold-out approach where the full dataset was portioned into 10 non- overlapping subsets of similar size.   Each of 10 model training datasets (i = 1, 2 …10) was comprised of a unique combination of nine subsets and the remaining subset was used for validation.  The stepwise selection procedure was executed 10 times, each time with a unique model training dataset i.  Stepwise selection was used because it is a commonly used procedure for identifying important geomorphic variables from a set of predictors [e.g., Murphey et al., 1977; Lindsay et al., 2004].  Accuracies (ACC) for the training and validation datasets were calculated separately, in both cases as 10 iACC ACC               (4-4) where ACCi was the percentage of correctly classified sites when Model Yi was applied to the dataset i. For the ordinal and nominal logistic regression models, mean square error was calculated as follows.  The set of regression equations for the ordinal and nominal solutions were applied to every sample in each dataset to calculate the four separate predicted probabilities: the 99  probability of being in the hillslope category; the probability of being in the swale category; the probability of being in the seepage channel category; and the probability of being in the fluvial category.  For each sample, the predicted category was the one with the highest probability.  The residual (resid) for each sample was calculated as:  1 true classresid prob       (4-5) where probtrue class is the predicted probability of being in the true class for that sample.  The mean square error was calculated for each dataset by averaging the squared values of the residuals.  Total mean penalty for each candidate model was calculated by substituting mean squared error for accuracy in Equation 4.  Criticisms of the stepwise selection procedure include model selection bias resulting from inclusion of redundant and irrelevant predictors, and instability due to influence of outliers [Burnham and Anderson, 2002].  In this study, the potential for model selection bias was minimized by the careful selection and screening of candidate predictors.  The potential for instability caused by outlier influence was reduced by previewing predictor distributions and removing samples that contained erroneous values.  During this process, the size of the dataset was reduced from 718 to 617.  Instability can be expressed through the selection different sets of predictor variables among datasets.  In this study, the results of the stepwise variable selection displayed instability; selected predictor variables varied between cross-validation datasets. Therefore, an additional procedure was applied to determine the best overall model.  For each model run with 10 datasets, a summary table was prepared that listed all of the selected variables ranked by the number of times each was selected.  Changes in model accuracy and penalty were graphically explored for variables that were selected in 5, 6 ... 10 of the cross- validation datasets.  The graph was examined to determine where the gain in information associated each additional variable was optimized, and this corresponded to the model of choice for extrapolation across the region.  100  4.4 Results 4.4.1 Basin similarity analysis For HRUs with lacustrine or deltaic soils and HRUs dominated by till deposits, the maximum improvement in within-groups error occurred when the number of clusters increased from two to three (Figure 4-8); therefore, three clusters were used in the classification for each parent material type.  The box plots confirm the negative relation between mean basin slope and wetland extent for both material types, and also reveal some of the complexities across the cluster types (Figure 4-9).  Cluster 1 (till) and Cluster 4 (lacustrine) have the lowest average slope, but are uniquely characterized by an average of 50 % wetland extent.  The intermediate clusters show the most overlap with adjacent clusters.  Cluster 3 (till) and Cluster 6 (lacustrine) have minimal slope overlap with adjacent clusters but share similar wetland extent.  The extrapolated results of the similarity analysis show an emerging pattern of the landscape with three characteristics: a patchy distribution of lacustrine deposits; high relief till-dominated basins along the west boundary; and medium-relief till dominated basins in the northeast corner (Figure 4-10). 101    Figure 4-8. Within group error by number of clusters for till dominated basins and basins with lacustrine/deltaic deposits.    Figure 4-9. Box plots of cluster by mean basin slope and cluster by percent of wetlands for till dominated basins and basins with lacustrine/deltaic deposits.  102    Figure 4-10. Map of study area with basin types.   103  4.4.2 Drainage networks and flowpath obstruction at road crossings This project used two drainage networks developed from LiDAR DEMs: the first, the preliminary network, was developed from a 2 m DEM with a smoothing algorithm applied to force flowpaths into natural drainage features below road crossings; the second, the DTW index flowpath, was created from a 1 m DEM with an algorithm applied to force drainage across roads at known culvert locations.  Results from each approach have implications for process domain modelling at the landscape scale.  Neither network showed correct drainage enforcement at stream – railroad intersections where valleys were filled to meet grade requirements (Figure 4-11a and b).  This was a specific concern when modelling process domain classes for terrain on the downslope side such crossings because the drainage area information was incorrect.  At locations where headwater flowpaths crossed roads, the smoothing algorithm used for the preliminary network effectively forced the flowpath across the dam created by the roadbed at most crossings (Figure 4-11c); whereas, for the DTW index flowpaths, the algorithm that forced drainage at culvert locations may have incorrectly predicted the flowpath route in cases where culvert data was not available (Figure 4-11d). 104   Figure 4-11. Flowpath diversion at railroad – major stream intersections for (a) the preliminary network and (b) DTW index flowpaths, and flowpath diversion at road – headwater stream intersections for (c) the preliminary network and (d) the DTW index flowpaths.   For reference purposes, the provincial hydrography layer (Prov. streams) was included to show watercourse locations determined by aerial photograph interpretation. 105  4.4.3 Screening candidate predictors Data were collected at 718 sites that were linked to the preliminary network.  The preliminary network was prepared in two iterations.  The network from the first iteration showed numerous flowpath diversions at road crossings; hence, a second iteration was completed in an effort to more accurately represent natural flowpaths at such locations.  After removal of sites associated only with the first iteration, the sample size was 617 (99 hillslopes, 301 swales, 148 seepage channels, and 69 fluvial channels).  The depth-to-water (DTW) index was available for a subset of 441 sites (42 hillslopes, 218 swales, 130 seepage channels, and 51 fluvial channels).  Lindsay et al. [2004] used a correlation coefficient > 0.60  to indicate significant correlations between drainage basin morphometrics.  There was some indication of correlation between a number of predictors (i.e., Pearson‟s correlation coefficient > 0.30), but not enough to suggest a critical lack of statistical independence (i.e., Pearson‟s correlation coefficient > 0.60) (Table 4-5).  There was some indication of a trend between drainage area and the other four topographic descriptors. The trend was positive with SL/k index, and negative with reach slope, mean basin slope, and valley width index.  Reach slope and mean basin slope shared a positive trend.  SL/k index had a negative trend with valley width.  Moisture regime displayed the highest overall independence with the other predictors.  DTW index shared a positive trend with reach slope, mean basin slope and a negative trend with valley width.  106  Table 4-5. Correlations between continuous class candidate predictors for complete dataset (n=617) and subset (italics in brackets) of sites with DTW index data (n=441).  Drainage area Reach slope SL/k index Mean basin slope Valley width Moisture regime DTW index Drainage area 1.00 Reach slope -0.36 (-0.27) 1.00 SL/k index 0.45 (0.37) 0.14 (0.33) 1.00 Mean basin slope -0.31 (-0.25) 0.47 (0.48) -0.08 (-0.03) 1.00 Valley width -0.46 (-0.37) -0.26 (-0.38) -0.36 (-0.33) -0.15 (-0.27) 1.00 Moisture regime 0.26 (0.17) -0.24 (-0.22) 0.02 (-0.04) -0.14 (-0.12) 0.06 (0.13) 1.00 DTW index (-0.16) (0.43) (-0.01) (0.33) (-0.43) (-0.23) 1.00  There was a decreasing probability of being in a lower ordered class as drainage area increased, and the pattern appeared non-linear (Figure 4-12a).  Drainage area thresholds at transitions between domain classes (cumulative logit = 0) were not readily discernable for the hillslope/swale, and occurred near 100 and 300 ha for the swale/colluvial channel, and colluvial/fluvial channel respectively.  This drainage area trend was consistent with the depth- slope product for predicting bed shear stress, whereby increased drainage area translates to greater depth, hence higher shear stress.  Thus, the visual trend indicating a shift to seepage and fluvial channels at specific drainage area thresholds can be explained by changes in erosion potential related to discharge increases.  There were no obvious trends between domain occurrence and stream slopes below 0.03%; however, there was an increasing probability of being in a lower order class when stream slope increased beyond 0.03% (Figure 4-12b).  This trend was contrary to the depth-slope product for predicting bed shear stress where depth is held constant and slope is increased, but generally consistent with the trend of decreasing stream gradient along the longitudinal profile.  Increases in the SL/k and DTW indices were both 107  associated with decreasing probability of being in a lower ordered class as values increased (Figure 4-12c and g).  The negative relation for SL/k index is consistent with the depth-slope product for predicting bed shear stress on two accounts.  First, the SL/k index provides a measure of the degree that a reach is over-steepened.  Second, knickpoints are associated with increased gradient and valley constraint, and groundwater is typically forced to the surface prior to entering such features [Stanford and Ward, 1993]; therefore, increased water depth, may also occur with higher SL/k index values.  Thus, shifts to higher order process domains are expected with increases in SL/k index. For the DTW index, the higher probability of being in a lower ordered site is associated with increased site dryness.  This pattern supports the theory that soil moisture status, an indicator of hillslope-scale upwelling, provides information on locations where drainage area-surface discharge relations are sufficiently modified to induce shifts in process domain transitions.  There were no obvious trends in the cumulative logit plot for mean basin slope (Figure 4-12d); however, there was no obvious violation of the proportional odds assumption.  There was a possible non-linear trend in probability of being in a lower ordered class with increasing moisture regime (Figure 4-12e).  There was an increasing probability of being in a lower ordered class as valley width index increased (Figure 4-12f).   This trend was difficult to explain as it was opposite from the conceptual model of channel development from confined to unconfined valleys, but may be related to the occurrence of planar hillslopes and un- channeled wetlands within the headwater areas.  The cumulative logit plots for the interaction term “drainage area  reach slope” displayed an increasing probability of being in a lower ordered class as values increased (Figure 4-13a).  The plots for “drainage area  SL/k index”, “drainage area  moisture regime”, and “drainage area  DTW index” showed negative trends (Figure 4-13b, d, and f), consistent with the patterns in their basic logit plots (Figure 4-12c, e, and g).  The plot for “drainage area  mean basin slope” showed increasing probability of being in a lower ordered class as values increased (Figure 4-13c).  When the “drainage area  valley width” increased, there was a decreasing probability of being in a lower ordered class (Figure 4-13e), which fit better with the conceptual model for channel development. 108        Figure 4-12. Logit plots of continuous class predictors. 109       Figure 4-13. Logit plots of interaction terms including drainage area and each continuous class predictor.  4.4.4 Ordinal, multinomial and binary logistic regression results for all four process domains The set of 18 candidate predictors included the continuous-class predictors, their quadratic terms, the interaction terms between drainage area and each continuous-class predictor, and the presence/absence of lacustrine soils.  For the ordinal model and multinomial 110  models respectively, 10 and 9 predictors were selected during the stepwise process with at least 50% frequency (Table 4-6).  For both models, five of the seven base predictors (i.e., predictors without quadratic or interaction terms) were selected; only valley width index and lacustrine soils were omitted.   Drainage area was the most significant predictor. Table 4-6. Predictors selected by stepwise selection procedure repeated for 10 cross-validation datasets, including model type, predictor, mean Chi-square score (score), and frequency of selection. Model type Predictor Score Frequency (%) Ordinal 1. Area 187.32 100  2. Reach slope 9.12 100  3. Moisture regime 7.54 100  4. SL/k index 7.07 90  5. Moisture regime 2  5.19 90  6. Mean basin slope 23.20 70  7. Area 2  10.03 70  8. Area  reach slope 5.60 70  9. Area  valley width 24.18 60  10. Area  SL/k index 4.45 50 Multinomial 1. Area 189.89 100  2. Moisture regime 19.57 100  3. Area 2  18.46 100  4. Reach slope 18.04 100  5. SL/k index 2  14.67 100  6. Area  reach slope 20.22 90  7. SL/k index 10.47 90  8. Area  valley width index 11.57 80  9. Mean basin slope 26.73 70 Binary - swale 1. Area 72.47 100  2. Moisture regime 12.91 100  3. Mean basin slope 9.17 100 Binary - seepage 1. Area 153.44 100  2. Reach slope 19.95 100  3. Area 2  10.92 100  4. Area  reach slope 9.37 100  5. SL/k index 9.00 100  6. SL/k index 2  6.53 100  7. Moisture regime 6.15 100 Binary - fluvial 1. Area 84.44 100  2. Mean basin slope 6.54 90  3. Area  valley width index 18.83 80  4. SL/k index 2  5.27 80  5. Area  reach slope 7.19 60 Note: Predictors with selection frequency less than 50% were excluded.  111  For the ordinal and multinomial models, a model selection procedure was used to maximize model performance while minimizing the number of predictors.  Candidate models were identified from the various frequencies of selection (Table 4-6).  For the ordinal logistic model, candidate models were as follows: Model 1 included predictors with 100 % frequency of selection (drainage area, reach slope, and moisture regime); Model 2 with 90% + frequency (Model 1 + SL/k index, and moisture regime 2 ); Model 3 with 80% + frequency (Model 2 + mean basin slope, drainage area 2 , and drainage area × reach slope); Model 4 with 60% + frequency (Model 3 + drainage area × valley width); and Model 5 with 50% + frequency.  For the multinomial logistic model, candidate models were as follows: Model 1 included predictors with 100 % frequency of selection (drainage area, moisture regime, drainage area 2 , reach slope, SL/k index 2 ); Model 2 with 90% + frequency (Model 1 + drainage area × reach slope, SL/k index); Model 3 with 80% + frequency (Model 2 + drainage area × valley width); and Model 4 with 70% + frequency (Model 3 + mean basin slope). From the ordinal logistic model, accuracy with the training dataset ranged from 56.2 to 59.7%, and the penalty ranged from 35.3 to 32.3 (Figure 4-14a and b).  Both statistics generally improved as additional parameters were added; therefore, Model 5 was identified as the best model.  For the multinomial model, accuracy with the training dataset ranged from 59.4 to 60.3%, and the penalty ranged from 32.9 to 31.0 (Figure 4-15a and b).  Both statistics generally improved as predictor variables were added; therefore, Model 4 was identified as the best model.   Figure 4-14. Evaluation plots of (a) mean accuracy with standard error, and (b) mean penalty with standard error for ordinal logistic models with increasing number of parameters (n=617). 112      Figure 4-15. Evaluation plots of (a) mean accuracy with standard error, and (b) mean penalty with standard error for multinomial logistic models with increasing number of parameters for the complete dataset (n=617).  For the three binary logistic models for predicting presence/absence of swales, seepage channels, and fluvial channels, the stepwise procedure selected three, seven, and five variables respectively (Table 4-6).  For the swale, seepage channel, and fluvial channel models, cut-offs were located at predicted probabilities of 0.83, 0.32, and 0.11 respectively.  Accuracies at these cut-offs were 73.6, 78.1, and 82.2 % (Figure 4-16).  The outputs from the three binary models were compiled into a single variable with 45.7 % accuracy, which was lower than both the ordinal and multinomial logistic models. 113   Figure 4-16. Average accuracy from the 10 model training datasets for each of six model types.  For the best ordinal model, the score test for the proportional odds assumption was < 0.0001 for all 10 cross-validation datasets, indicating that logit trends for individual predictors were not parallel, and exploration of models that allow flexibility for individual parameter estimates was warranted.  The best ordinal model shared seven predictors with the best multinomial model and eight predictors across the best binary models (Table 4-7, Table 4-8, and Table 4-9).  The shared selection provides evidence for relations between these predictors and process domain class occurrence.  For the best ordinal model, there was no overlap between the three intercept estimates and their confidence intervals, to indicate a statistical separation between the four classes (Table 4-7).  The odds ratios for drainage area, reach slope, mean basin slope, and moisture regime indicated a decreasing probability of being in a lower ordered class as values of these predictors increase.  The odds ratio for SL/k index indicated the opposite trend.  For all odds ratios, the confidence interval did not include 1 (null hypothesis) and hence the odds were significantly different from 1(i.e. 50/50 odds).  In the set of three equations from the best multinomial model (Table 4-8), there was also no overlap between the three intercept 114  estimates.  The coefficient estimates from the three sets of models showed varying degrees of overlap – SL/k index, mean basin slope, and drainage area × valley width, having the greatest, and reach slope and moisture regime having the least. 115   Table 4-7. Ordinal logistic regression results for the best process domain model including dataset, coefficient estimates and their 95% confidence intervals (lower and upper), standard errors (SE), significance (P), odds ratio and their 95% confidence intervals. Variable Estimate Lower Upper SE P Odds Ratio Lower Upper Intercept 1 -2.37 -2.75 -2.00 0.19 <.0001 -- -- -- Intercept 2 0.98 0.68 1.29 0.15 <.0001 -- -- -- Intercept 3 3.33 2.86 3.79 0.24 <.0001 -- -- -- Area -2.17 -2.48 -1.87 0.16 <.0001 0.06 0.04 0.10 Reach slope -0.61 -0.87 -0.36 0.13 <.0001 0.54 0.42 0.70 SL/k index 0.47 0.23 0.70 0.12 <.0001 1.60 1.26 2.01 Mean basin slope -0.22 -0.42 -0.01 0.11 0.04 0.81 0.66 0.99 Moisture regime -0.47 -0.73 -0.20 0.14 0.00 0.78 0.65 0.94 Area 2  -0.59 -0.88 -0.31 0.15 <.0001 -- -- -- Moisture regime 2  0.22 0.04 0.40 0.09 0.02 -- -- -- Area × reach slope -0.50 -0.77 -0.23 0.14 0.00 -- -- -- Area × valley width 0.22 -0.03 0.46 0.13 0.09 -- -- -- Area × SL/k index 0.29 0.06 0.51 0.11 0.01 -- -- -- 116  Table 4-8. Multinomial logistic regression results for the best process domain model including dataset, coefficient estimates and their 95% confidence intervals (lower and upper), standard errors (SE), and significance (P), where fluvial channel (4) is reference domain class for hillslope (1), swale (2), and seepage channel (3) domain classes. Variable Domain class Estimate Lower Upper SE P Intercept 1 1.44 0.59 2.28 0.43 0.00 Intercept 2 3.41 2.65 4.16 0.38 <.0001 Intercept 3 2.33 1.57 3.09 0.39 <.0001 Area 1 -3.75 -4.78 -2.72 0.52 <.0001 Area 2 -2.91 -3.80 -2.02 0.45 <.0001 Area 3 -1.12 -2.01 -0.23 0.45 0.01 Reach slope 1 -1.61 -2.65 -0.56 0.53 0.00 Reach slope 2 -1.19 -2.18 -0.20 0.51 0.02 Reach slope 3 -0.19 -1.18 0.81 0.51 0.71 SL/k index 1 0.51 -0.39 1.41 0.46 0.27 SL/k index 2 0.23 -0.62 1.09 0.43 0.59 SL/k index 3 -0.21 -1.05 0.63 0.43 0.62 Moisture regime 1 -0.81 -1.36 -0.26 0.28 0.00 Moisture regime 2 -0.01 -0.41 0.39 0.20 0.96 Moisture regime 3 0.24 -0.12 0.60 0.19 0.20 Mean basin slope 1 -0.57 -1.17 0.04 0.31 0.07 Mean basin slope 2 -0.20 -0.75 0.35 0.28 0.48 Mean basin slope 3 -0.34 -0.86 0.19 0.27 0.21 SL/k index 2  1 0.47 0.07 0.87 0.20 0.02 SL/k index 2  2 0.55 0.18 0.91 0.19 0.00 SL/k index 2  3 0.35 0.00 0.70 0.18 0.05 Area 2  1 -0.79 -1.55 -0.03 0.39 0.04 Area 2  2 -1.21 -1.90 -0.52 0.35 0.00 Area 2  3 -0.41 -1.03 0.20 0.31 0.19 Area × reach slope 1 -0.83 -1.72 0.06 0.45 0.07 Area × reach slope 2 -0.64 -1.49 0.21 0.43 0.14 Area × reach slope 3 -0.03 -0.81 0.75 0.40 0.95 Area × valley width 1 0.67 0.02 1.33 0.33 0.04 Area × valley width 2 0.23 -0.31 0.77 0.27 0.40 Area × valley width 3 0.56 0.10 1.02 0.23 0.02 117  Table 4-9. Logistic regression models for predicting the threshold (‖) to swales, seepage channels and fluvial channels including dependent variable, coefficient estimates and their 95% confidence intervals (lower and upper), standard errors (SE), significance (P), odds ratio and their 95% confidence intervals. Threshold  Variable Estimate Lower Upper SE P Odds Ratio Lower Upper Hillslope ‖ swale Intercept -2.22 -2.55 -1.89 0.17 <.0001 -- -- -- Area -1.15 -1.44 -0.86 0.15 <.0001 0.32 0.24 0.43  Moisture regime -0.70 -1.07 -0.33 0.19 0.00 0.50 0.34 0.72  Mean basin slope -0.40 -0.65 -0.15 0.13 0.00 0.67 0.52 0.86 Swale ‖ seepage channel Intercept 1.11 0.82 1.41 0.15 <.0001 -- -- -- Area -2.11 -2.49 -1.73 0.19 <.0001 0.05 0.03 0.09  Reach slope -0.90 -1.23 -0.58 0.17 <.0001 0.41 0.29 0.56  SL/k index 0.50 0.18 0.83 0.17 0.00 2.14 1.47 3.10  Moisture regime -0.32 -0.54 -0.09 0.11 0.01 0.73 0.58 0.91  SL/k index 2  0.26 0.07 0.44 0.09 0.01 -- -- --  Area 2  -0.82 -1.11 -0.52 0.15 <.0001 -- -- --  Area × reach slope -0.83 -1.17 -0.50 0.17 <.0001 -- -- -- Seepage channel ‖ fluvial channel Intercept 3.42 2.81 4.04 0.31 <.0001 -- -- -- Area -2.37 -3.08 -1.66 0.36 <.0001 0.093 0.046 0.19  Mean basin slope -0.47 -0.86 -0.07 0.20 0.02 0.627 0.421 0.933  SL/k index 2  0.37 0.13 0.61 0.12 0.00 1.45 1.138 1.848  Area × valley width index -0.57 -1.03 -0.11 0.23 0.01 -- -- --  Area × reach slope 0.63 0.23 1.04 0.21 0.00 118  Across the three best binary logistic models, the odds ratios for drainage area, reach slope, mean basin slope, and moisture regime indicated a decreasing probability of being in a lower ordered class as values of these predictors increase (Table 4-9 and Figure 4-17).  These trends were consistent with those indicated by the ordinal model.  The odds ratio for SL/k index for both seepage channels and fluvial channels (Figure 4-17) showed non-linear trends, with decreasing probability at lower SL/k index values shifting to increased probability at higher SL/k values.  For the three binary models (Table 4-9), the odds ratio confidence intervals did not include 1.  There was no evidence for the multinomial model or combined binary model as a better overall model than the ordinal model (Figure 4-16).  119       Figure 4-17. Predicted probability of occurrence by predictor for cumulative logits for swale, seepage channels, and fluvial channels: (a) drainage area, (b) reach slope, (c) SL/k index, (d) moisture regime, and (e) is mean basin slope.  The range of the x axis represents the observed range for each predictor.  For predictions, all values for other variables were standardized to median values.   120  When the model selection procedure was repeated with the subset of data containing the DTW index (n= 441), the DTW index was selected during the stepwise process for the ordinal model (Table 4-10), but not for the multinomial or the binary logistic models.  The odds ratios for DTW index and moisture regime were similar, and indicated a decreasing probability of being in a lower ordered class as values of these predictors increase (Table 4-11). Table 4-10. Ordinal logistic regression predictors for the DTW index subset (n=441), selected by stepwise selection procedure repeated for 10 cross-validation datasets, including predictor, mean Chi-square score (score), and frequency of selection. Predictor Score Frequency (%) Area 113.52 100 Reach slope 26.47 100 Area 2  22.68 100 Moisture regime 6.17 100 DTW index 8.67 90 Moisture regime 2  5.87 80 SL/k index 4.96 80 Area × reach slope 5.30 70   121  Table 4-11. Ordinal logistic regression results for the DTW index subset (n=441), including dependent variable, coefficient estimates and their 95% confidence intervals (lower and upper), standard errors (SE), significance (P), odds ratio and their 95% confidence intervals. Effect Estimate Lower Upper SE P Odds ratio Lower Upper Intercept 1 -3.01 -3.52 -2.51 0.26 <.0001 -- -- -- Intercept 2 0.53 0.17 0.88 0.18 0.00 -- -- -- Intercept 3 3.14 2.60 3.68 0.27 <.0001 -- -- -- Area -1.82 -2.13 -1.50 0.16 <.0001 0.09 0.06 0.14 Reach slope -0.66 -0.93 -0.38 0.14 <.0001 0.52 0.40 0.68 Moisture regime -0.64 -0.94 -0.33 0.16 <.0001 0.75 0.60 0.94 DTW index -0.32 -0.55 -0.09 0.12 0.01 0.73 0.58 0.91 SL/k index 0.28 0.02 0.54 0.13 0.03 1.32 1.02 1.71 Area 2  -0.61 -0.83 -0.40 0.11 <.0001 -- -- -- Moisture regime 2  0.35 0.11 0.59 0.12 0.00 -- -- -- Area × reach slope_ -0.42 -0.70 -0.14 0.14 0.00 -- -- --  4.4.5 Application of best models for mapping All reaches from the preliminary network (Figure 4-18a) were assigned the predicted process domain class from the best ordinal model (Figure 4-18b), and predicted channel occurrence from the binary logistic model for predicting the swale – seepage channel threshold (Figure 4-18c).  The improvement in accuracy from 59.7 % for the best process domain class model to 78.1 % for the best channel occurrence model was readily apparent. 122    Figure 4-18. Stream network maps of (a) the preliminary network, (b) the predicted process domains including swales, seepage channels and fluvial channels, and (c) the predicted occurrence probability for channels. 123  4.4.6 Evidence for anchored and transient domain transitions in response to land use and climate shifts A seepage channel is expected to replace a swale when a discharge increase elevates the shear stress beyond the shear strength associated with the ground cover and soil.  The magnitude of the discharge increase required to trigger such a shift will vary across the landscape.  For any given swale segment in the modelled network, the difference between the predicted probability of a channel and the probability threshold for a channel can provide some information on the relative increase in discharge required for a shift.  For example in Figure 4-21, the magnitude of additional discharge required to trigger a shift increases through segments (a), (b), and (c).  The channel head below segment (a) is predicted to be transient, with small changes in discharge resulting in considerable headward advancement.  The channel head below segment (b) is predicted to be less transient into segment (b).  Channel segment (c) has the highest degree of anchoring within the swale domain; hence, it requires the greatest amount of additional discharge to trigger a shift.  These interpretations provide information on specific landscape locations where evolution through headward advancement of the channel network is most likely to occur. 124   Figure 4-19. Stream network map with domains with a swale predicted where the probability of a channel is less than the 0.32 cut-off and channel where the probability of a channel is ≥ 0.32. Segments (a), (b) and (c) are predicted to be in the swale domain with probability of being a channel as 0.31, 0.26, and 0.16 respectively. 125  4.4.7 Effects of reducing landscape variability and spatial extent on accuracy of binary logistic model for channel occurrence The dataset (n = 617) was divided into six subsets using the HRU type boundaries (Figure 4-10).  Based on the number of variables selected, the models for the low relief HRU classes (Type 1 – till dominated, and Type 1 – lacustrine) showed higher complexity than the models for high relief HRUs (Type 3 – till dominated, and Type 3 – lacustrine) that included drainage area as the only predictor (Table 4-12).   The model for low relief till-dominated HRUs had lower model accuracy in comparison to regional model (Figure 4-20).  The model for high relief till-dominated HRU (Basin type 3) showed improved model accuracy.  The models for lacustrine HRU showed the opposite trend, with accuracy highest in low relief HRU (Basin type 1), and decreasing to a similar level as the regional model with increasing relief. Table 4-12. Predictors selected by stepwise selection procedure repeated for 10 cross-validation datasets for six data subsets including HRU parent material and number, sample size and predictors with selection frequency ≥ 50 %. Parent material HRU# n Predictors Till dominated 1 87 Area + basin slope + moisture regime + area × valley width  2 231 Area +  area 2  + SL/k index  3 55 Area Lacustrine 1 64 Area + reach slope + SL/k index  2 104 Area + basin slope + moisture regime  3 76 Area   Figure 4-20. Evaluation plots of accuracy for binary logistic models for channel presence by basin type.  Reference line represents regional model accuracy of 78.1 %. 126   4.5 Discussion 4.5.1 Basin similarity analysis HRUs were split into two groups – till dominated and lacustrine units – because the specific types and rates of erosion within each group translate to different trends in landscape stability and basin-scale sediment yield.  There was no evidence from the three logistic modelling exercises to indicate that presence/absence of lacustrine soils in a drainage basin influenced the location of transitions between erosion processes.  Nonetheless, Quaternary lacustrine deposits that occur on valley sides in close proximity to stream channels are an important sediment source at the drainage basin scale in the Foothills and Interior Plains [Church et al., 1989].  The study area is under intense development pressure from forest and energy industries and this development is associated with changes in the occurrence and structure of native fish assemblages [Scrimgeour et al., 2008].  Headwater streams in the region provide habitat for native bull trout (Salvelinus confluentus) – a species listed as sensitive in the province of Alberta [Alberta Sustainable Resource Development and Alberta Conservation Association, 2009b], and Athabasca rainbow trout – a species considered May Be At Risk [Alberta Sustainable Resource Development and Alberta Conservation Association, 2009a]. While trends in development are related to declines in these two species, the specific hydrogeomorphic effects of the land-use activities on fish habitat remain poorly understood [Ripley et al., 2005; Scrimgeour et al., 2008].   Procedures to assess the effects of land-use typically involve monitoring sediment yield and bed material characteristics [MacDonald et al., 1990].  Thus, because of the high sediment yield associated with landslides and erosion of lacustrine soils, their spatial distribution, as captured with the HRU classification system, should serve as a primary consideration when designing projects to assess the hydrogeomorphic effects of land use on headwater streams. When HRU type was used to partition the landscape into patches for process domain modelling, an interesting pattern emerged, whereby low relief terrain with a high wetland proportion had more complex channel head threshold models (i.e., greater number of predictors) than high relief terrain with a low wetland proportion.   This pattern supports the recommendations of Devito et al. [2005a] to consider the range of factors that influence groundwater flow processes in addition to local slope and drainage area in the Boreal Plain.  In 127  comparison to the regional model, HRU-based models for predicting the transitions to seepage and fluvial channels had fewer parameters, and models for three of the six HRU types showed improved accuracy.  These findings support an alternative approach to using a single complex model at the regional scale to explain process domain transition.  Specifically, to better understand drainage basin evolution in regions with complex terrain, the landscape can be classified into patches with similar topography, hydrologic connectivity, and surface type, consistent with the recommendations of Buttle [2006], and then models for the sequential process domains can be developed for each patch type. 4.5.2 Preliminary drainage network and flowpath obstruction at road crossings Two drainage networks were used in this study, each with different methodology to address flowpath routing at road grade intersections.  The first strategy was to simulate „natural‟, pre-disturbance flowpaths, and the second was to simulate flow associated with the altered land surface.  The first strategy could be improved by permitting the user to alter thresholds, specifically the DEM smoothing radius around each crossing location, used in the smoothing algorithm.  Mouton [2005] found increasing the LiDAR DEM grid cell size to 6 m resolution was required to avoid re-routing flowpaths along ditches.  The range in fill heights at crossings in the study area warrants a flexible approach to addressing this problem.  The second strategy is dependent on information on culvert and bridge locations provided prior to processing and could be improved by ensuring a comprehensive dataset is provided. 4.5.3 Relations between process domain occurrence and predictors The strategy for developing a regional model for process domain occurrence was to select predictors that represented the range of factors influencing local discharge, in addition to drainage area and local slope.  Features of the landscape that complicate the prediction of local discharge include topographic complexity that is related to changes in alluvial area along the stream profile [Tonina and Buffington, 2009], and site-specific variations in soil moisture that are related to local topography, topographic position, and aspect [Beckingham et al., 1996].  The scale-independent SL/k index allowed for rapid map-based exploration and identification of profile anomalies potentially associated with changes in hyporheic flow.  The SL/k index improved predictions of process domain transitions in small drainage basins.  The soil moisture index incorporated the effects of topography, topographic position, and aspect.  Although the 128  resolution was much less than the resolution of LiDAR-derived predictors, the aerial photograph interpretation and extensive ground truthing ensured that actual field conditions were captured. Soil moisture index was an important predictor in almost all process domain models that were developed. The three types of categorical data models used to explain the spatial distribution of erosion processes revealed three main points.  First, ordinal and multinomial logistic regression models performed similarly because the flexibility gained in the multinomial procedure was offset by the information loss in its dependent variable.  Of these two types of models, ordinal logistic regression models were favored because they provided more information about the relations between individual predictors and the erosion process domain distribution.  Second, drainage area is the dominant factor explaining process domain occurrence.  The addition of other variables serves to refine model predictions.  Drainage area was determined from a single cell flow direction algorithm.  Algorithms that allow dispersal of flow in multiple directions have proven more accurate in some cases [Kopecký and Čížková, 2010] and could be evaluated for modelling process domains.  Third, a binary logistic model can provide more information about a single erosion process transition than an ordinal model that simultaneously models three transitions.   Different predictors with specific regression coefficient estimates are permitted in the binary model.  Specific procedures need to be developed to combine the results from a series of binary models into a map showing spatial distribution of all classes. 4.5.4 Interpretations for drainage basin evolution Aside from glaciation, shifts in drainage area-discharge relations over time may be the dominant factor affecting landscape evolution trends.  Such shifts may result from increased storminess or runoff due to climate shifts, or to increased surface runoff in response to specific land use activities that reduce evapotranspiration (timber harvesting) or infiltration (road construction).   The outputs from spatial process domain models can be examined for information related to drainage network confluences and over-steepened reaches that serve to anchor specific domain transition through relatively large increases in discharge.  In contrast, segments lacking such features may display transient domain boundaries and higher erosion rates.  129  4.5.5 Opportunities to improve the model When using presence/absence data for modelling, most information about the process of interest is contained around transition locations.  To maximize the amount of information collected, random effects should be incorporated into the study design [e.g., Chapter 2]. Specifically, rather than determining the process domain for a single reach at each location, each point should be expanded as a transect along the stream profile to include enough reaches to capture the transition between two process domains. Each transect would serve as a random effect during the analysis. Even though soil moisture regime was captured at a lower resolution than other predictors, it contributed to improved model performance.  Thus, additional remotely-sensed predictors that have been validated through ground truthing could be explored.  With extensive ground truthing and improvements to its associated flowpath network, the DTW index may be more useful.  Synthetic aperture radar data have been used for a variety of related applications [Creed et al., 2008; Sass and Creed, 2008; Clark et al., 2009] and may have application. Wetland extent based on probability of wetland occurrence [Lindsay et al., 2004] can also be derived from LiDAR data and could be considered.  In addition to providing elevation information, LiDAR also provides an distinctive intensity measure when the pulse hits water, which has proven useful for differentiating terrestrial and aquatic ecosystems [e.g., Chust et al., 2008].  Therefore, there may be an opportunity to explore the remote detection of wet areas using LiDAR intensity measures. 130  Chapter 5 Conclusion This study demonstrated that organization of complex glaciated landscapes can be explained using a suite of remotely-sensed descriptors to represent the multiple factors that affect normal scaling relations.  Four obstacles to modelling the spatial distribution of process domains at the landscape scale are as follows: (1) the dynamic nature of the swale – colluvial channel transition necessitates the acquisition of field data rather than visual inspection of slope- area plots; (2) thresholds change across landscapes and between basins; (3) previous land- sculpting processes, including glaciations, control channel network organization; and (4) high- resolution DEMs are required to discern the first transition in the hillslope – swale – colluvial channel – fluvial channel sequence [Montgomery and Foufoula-Georgiou, 1993].  These limitations were overcome by specific strategies that varied according to the scale of the investigation.  Important findings from the study of fish presence/absence [Chapter 2] are presented first, followed by findings from Chapters 3 and 4, where spatial process domain models were developed from LiDAR DEMs.  As was discussed in Chapter 1, a congruence of the factors that are related to the fish presence-absence threshold [Chapter 2] and the fluvial channel presence-absence threshold [Chapter 4] was expected. All of the models presented in this dissertation are static in nature.  This is problematic because all thresholds that were modeled in Chapters 2, 3, and 4 are expected to respond to stochastic events such as major storms.  Land-use induced changes in flow regime or sediment source stability are also expected to influence the spatial distribution of domains.  Stochastic models can predict the headward expansion of colluvial channels through landsliding in response to major storms and subsequently route this material through the network [e.g., Benda and Dunne, 1997]; however, accurate representations of the drainage network, such as those presented in Chapters 3 and 4 are required as a starting point for such models.   Furthermore, the static models can provide some information on the probability of a threshold shift occurring for any specific channel segment, and this information could be incorporated into a stochastic model. The research revealed two types of domain transitions with implications for forecasting migration rates in response to climate shifts or land-use change.  Anchored transitions occur at confluences and gradient breaks.  As such, a substantial increase in discharge would be required to trigger headward migration.  Transient domain transitions are located along sections of the 131  network where changes in drainage area or slope occur more gradually.  These transitions would be more responsive to climate or land-use increases in local discharge increases related to extreme events, land-use, or climate shifts. The fish presence-absence models provided a preview of the importance of basin scale relief for predicting process domain boundaries in the transitional landscape between the Rocky Mountain and the boreal plain.  Mean basin slope (basin slope) was the only predictor, other than drainage area, that was included in the three individual species models, and the “all species” model (Figure 2-2).  Bull trout (Salvelinus confluentus) preferred the highest relief basins, while rainbow trout (Oncorhynchus mykiss) and brook trout (Salvelinus fontinalis) preferred intermediate relief basins.  The influence of basin-scale slope has been recognized in studies of species distributions in similar landscapes [e.g., Porter et al., 2000]; nevertheless, this factor was not considered in recovery plans for native rainbow trout within the study area [Alberta Sustainable Resource Development and Alberta Conservation Association, 2009a]. There are several anticipated benefits from using a LiDAR DEM for modelling fish occurrence.  First, the network of swales and seepage channels [Chapter 4] is under-represented in the standard drainage network map (provincial streams layer), which was derived from aerial- photograph interpretation.  This issue may not affect the performance of fish occurrence models because streams at the no fish/fish threshold are of sufficient size to be well represented on the provincial streams layer.  However, many of the land-use activities with potential to alter hydrologic and erosion processes occur in swales and seepage channels within the under- represented portion of the drainage network.  Such impacts can be conveyed downstream to the fish bearing portion of the network.  Thus, to assess the effects of land-use on aquatic habitat, fish occurrence models derived from the 25 m DEM / provincial streams layer could be coupled with more complete models of the swale – seepage channel – fluvial channel network derived from a LiDAR DEM [Chapter 4].  A second benefit of using a LiDAR DEM is the ability to detect localities with high slopes.  For example, the resolution of the 25 m DEM may be inadequate for detecting knickpoints where streams flow over resistant bedrock layers; thus, localized steep reaches may not be apparent on maps.  Waterfalls, chutes, and other natural features associated with knickpoints can prevent upstream fish migration and limit fish occurrence in headwater streams of the steep foothills terrain immediately adjacent to the Rocky Mountains.  Juvenile bull trout (Salvelinus confluentus) are most common in these higher relief 132  areas.  Thus, the use of a LiDAR DEM may improve performance of fish occurrence models for bull trout and other species that inhabit steep regions where natural barriers limit fish distribution in headwater streams.   Slope-area plots are well suited for delineating process domains in landscapes where drainage area-discharge relations hold across the area of interest.  In such cases, domain transitions can be linked to specific drainage area thresholds  [Montgomery, 2002].  Drainage area-discharge relations can be affected by groundwater movement at the reach, hillslope, and basin scale.  To address these processes in complex landscapes, process domains models can be modified by expanding the set of predictors to include factors influencing groundwater flow across multiple scales.  The models developed in Chapter 3 including factors influencing groundwater movement at the reach and basin scale.  For Chapter 4, groundwater movement was considered at the reach, hillslope, and basin scales. Contrasting procedures were used to collect field data and study the spatial variability associated with the sequence of process domains in Chapter 3 versus Chapter 4,.  For the Dutch Creek study, which focused on a 57 km 2  area [Chapter 3], process domain transitions were identified within individual basins by surveying longitudinal profiles (Figure 3-2).  For the Hinton landscape scale project, which covered 10,000 km 2 [Chapter 4], process domains were identified for specific reaches distributed across the study area (Figure 4-2).  Both projects utilized categorical data analyses, which involve making probability-based predictions around transitions.  The data points surrounding these locations contain the most valuable information [Santner and Duffy, 1989], and should be targeted in future studies.  Mixed-effects models were required to analyze the grouped longitudinal data from the Dutch Creek study.  These models contained information about the degree to which individual basins deviated from the central tendency (Figure 3-7), and the degree to which transitions vary as the distance from source increases along each profile (Figure 3-8).  Collection of data along longitudinal profiles is a common approach for geomorphic field studies.  Study designs for such projects involving multiple basins require mixed-effects analyses to properly represent the covariance structure associated with grouped data [Pinheiro and Bates, 2004]; however, such statistical procedures are not commonly applied in hydrogeomorphic research [e.g., Pike et al., 2010].  The field data collection procedure for the Hinton region study [Chapter 4], whereby individual reaches were dispersed across the landscape, was not suited for mixed-effects procedures for analyzing spatial 133  variation of process domains; thus, there were no available means to evaluate basin-scale variability.  In conclusion, longitudinal surveys coupled with mixed-effects models [Chapter 3] provided the best approach for explaining the considerable variability associated hydrogeomorphic processes across multiple drainage basins. Two distinct procedures were employed to account for thresholds that change across landscapes and between basins.  To address landscape diversity, basin scale descriptors were included as candidate predictors within multiple regression models.  In the Dutch Creek study [Chapter 3], basin slope and area × basin slope were selected using a forward stepwise procedure (Table 3-4).  In the Hinton region project [Chapter 4], basin slope was selected in four out of five model types that were used to explore the spatial distribution of process domains (Table 4-6).  These results are consistent with Heine et al. [2004], and indicate that basin slope can help modify drainage-area-based scaling relations within hydrogeomorphic models across multiple basins.  Specifically, increases in both drainage area and basin slope were associated with decreasing probability of being in a lower-ordered process domain category (Figure 4-17a and e).  In the second procedure, a basin similarity analysis was completed by clustering basins with similar topography, hydrologic connectivity, and surface type [Chapter 4].  This exercise revealed a patchy landscape pattern owing to the spatial distribution of lacustrine deposits and an isolated medium relief region in the NE corner of the study area (Figure 4-10).  When models were developed for each of the six basin types, model complexity increased as basin relief decreased (Table 4-12).  Basin classification exercises have proved useful in other studies on landscape organization.  For example, accuracy of models to predict thresholds between ephemeral, intermittent, and perennial streams in New Zealand improved when field sites were grouped according to hydrogeological area [Storey and Wadhwa, 2009]. Glacial-induced deviations from the idealized concave-up graded profile prevent the application of standard scaling relations for predicting process domain transitions [Brardinoni and Hassan, 2006].  A normalized stream length-gradient (SL/k) index was converted to a dimensionless form whereby the index calculation for all segments in the dataset was no longer tied to an arbitrary downstream reference point such as a river mouth.  Rather, each segment‟s downstream point served as the reference.  The use of this form of the SL/k index was found to improve scaling relations in landscapes where profile deviations occur.  In the Dutch Creek study, SL/k index was selected during the forward stepwise procedure (Table 3-4).  In the 134  Hinton region project, SL/k index was selected in four out of five model types that were used to explore the spatial distribution of process domains (Table 4-6). High-resolution DEMs are required to discern the first transition in the hillslope – swale – colluvial channel – fluvial channel sequence [Montgomery and Foufoula-Georgiou, 1993]. High-resolution data were not yet available to predict fish presence/absence, so a 25 m DEM was used [Chapter 2].   To improve reach slope accuracy a line length correction factor [McMaster, 1986] was applied.  Such line length corrections may be required where the rise/run formula is used to determine  stream slope regardless of DEM resolution because the line length is typically simplified and shortened during the interpretation procedures (i.e., aerial photograph interpretation or DEM point interpolation), thereby elevating the slope estimate.  High resolution data permit the accurate calculation of curvature statistics that are used in the detection of the hillslope – swale transition.  In the Dutch Creek study, profile curvature measured over a 100 m distance improved the accuracy of predictions of the swale domain category (Figure 3-9).  For the Hinton region project, the high-resolution DEMs were important for predicting flowpath locations near source areas.  Two drainage networks were used in Chapter 4, each with a different methodology to address flowpath routing at road grade intersections.  The first strategy was to simulate „natural‟, pre-disturbance flowpaths, and the second was to simulate flow associated with the altered land surface.  Mouton [2005] found increasing the LiDAR DEM grid cell size to 6 m resolution was required to avoid re-routing flowpaths along ditches.  These results highlight the importance of matching the DEM resolution to the process of interest and completing accuracy evaluations using field measured data. This research identified a scarceness of hydrometric data for small catchments (i.e., drainage areas less than 10 km 2 ), within which the three process domain transitions occur.  One of the underlying premises of the models that were developed was that predictions of shear stress (proportional to the product of local slope and water depth) at any point within the drainage network, can be improved by considering factors that influence local drainage area – discharge relations.  Additional work is required to confirm this premise.  Specifically, the data from additional long-term hydrometric stations could be used to confirm dominant factors.   The contributions from local, hillslope, and basin scale groundwater systems are also expected to vary temporally.  Investigations into the relative contributions of these three systems through the 135  rising and falling limbs of an event could consider the direct relation between groundwater quality (total dissolved solids) and residence time. Surface water supply is an important resource for the petroleum extraction and processing industries, which operate throughout the region and in downstream areas with extensive oil sands deposits.  Opportunities to augment winter baseflow in major rivers through vegetation management may be explored in the future.  This dissertation has identified a potential sedimentation-related impact that may arise if drainage area – discharge relations are intentionally modified to promote water yield.  Specifically, increases in discharge have the potential to trigger channel head migration into the extensive swale network, thereby releasing stored sediments into the drainage network; however, the degree to which headward migration occurs depends on whether channel heads are anchored at specific features including channel heads, and over-steepened areas.  Outputs from the process domain occurrence models presented in Chapters 3 and 4 could be used to identify specific sites that are most susceptible to headward migration. In complex landscapes, including the Rocky Mountain foothills, additional work is also required to establish linkages between specific land-use activities, impacts to aquatic habitats, and effects of such changes on populations and persistence of aquatic organisms.  Such knowledge is required to make informed conservation decisions [Rosenfeld, 2003].  For example, sedimentation from gravel road run-off is a major impact to streams in forested regions of North America [Waters, 1995]; however, the effects of high traffic levels (e.g., more than 1,000 trips per day for some road sections) on the different types of channels and aquatic organisms that occur in the study area remains poorly understood.  To address such gaps, the fish occurrence [Chapter 2] and process domain models [Chapters 3 and 4], could have application in space-for-time and flume studies.  Space-for-time studies remain an important tool in geomorphic research where processes occur at time scales between 10 and 100 years [e.g., Cover et al., 2010].   A variety of factors determine the spatial distribution of process domains (e.g., drainage area, basin and local relief, bedrock type), and the spatial models were developed based on these criteria; hence these models can assist during study site selection in space-for-time experiments.  The effects of disturbances, including sedimentation, on streambeds is very difficult to study in a natural setting.  As a result, much of the knowledge on how streambeds respond to disturbance has been developed with flume studies.  Such studies 136  have typically focused on low gradient channels with homogeneous bed material.  Flume studies where the specific physical conditions of headwater streams are mimicked, while sedimentation disturbances of varying intensity and duration are applied, provide much promise towards advancing knowledge of these systems [e.g., Zimmermann et al., 2010].  For any region, landscape-scale process domain models could indicate the spectrum of channel types to be mimicked in controlled flume experiments.    137  References Alberta Environmental Protection (1996), Alberta Vegetation Inventory Standards Manual, edited by Resource Data Division - Data Acquisition Branch, p. 42 + App., Alberta Environmental Protection. Alberta Sustainable Resource Development, and Alberta Conservation Association (2009a), Status of the Athabasca Rainbow Trout (Oncorhynchus mykiss) in Alberta, edited by Alberta Sustainable Resource Development, p. 32, Edmonton, Alberta. Alberta Sustainable Resource Development, and Alberta Conservation Association (2009b), Status of the Bull Trout (Salvelinus confluentus) in Alberta: update 2009, edited by Alberta Sustainable Resource Development, p. 48, Edmonton, Alberta. Allison, P. D. (1999), Logistic Regression Using the SAS System: Theory and Application, 304 pp., SAS Institute Inc., Cary, NC. Archibald, J. H., G. D. Klappstein, and I. G. W. Corns (1996), Field guide to the ecosites of southwestern Alberta, 540 pp., Natural Resources Canada, Canadian Forest Service, Northwest Region, Northern Forestry Centre, Edmonton, Alberta. Band, L. E. (1986), Topographic partition of watersheds with digital elevation models, Water Resour. Res., 22(1), 15-24. Baxter, C. V., C. A. Frissell, and F. R. Hauer (1999), Geomorphology, logging roads, and the distribution of bull trout spawning in a forested river basin: implications for management and conservation, Trans. Am. Fish. Soc., 128, 854-867. Baxter, C. V., and F. R. Hauer (2000), Geomorphology, hyporheic exchange, and selection of spawning habitat by bull trout (Salvelinus confluentus), Can. J. Fish. Aquat. Sci., 57(7), 1470-1481. 138  Bayrock, L. A., and T. H. F. Reimchen (1975), Surficial geology Alberta Foothills and Rocky Mountains, Map Sheet No. 3, Alberta Environment and Alberta Research Council, Edmonton, Alberta. Beckingham, J. D., I. G. W. Corns, and J. H. Archibald (1996), Field guide to the ecosites of west-central Alberta, The University of British Columbia Press, Vancouver. Benda, L., and T. Dunne (1987), Sediment routing by debris flows, paper presented at Erosion and sedimentation in the Pacific Rim, IAHS Publ. 165, Corvallis, Oregon, August, 1987. Benda, L. (1990), The influence of debris flows on channels and valley floors in the Oregon Coast Range, USA, Earth Surf. Processes Landforms, 15(5), 457-466. Benda, L., and T. Dunne (1997), Stochastic forcing of sediment routing and storage in channel networks, Water Resour. Res., 33(12), 2865-2880. Benda, L., M. A. Hassan, M. Church, and C. L. May (2005), Geomorphology of steepland headwaters: the transition from hillslopes to channels, J. Am. Water Resourc. Ass., 41, 835-851. Benda, L., D. Miller, K. Andras, P. Bigelow, G. H. Reeves, and D. Michael (2007), NetMap: a new tool in support of watershed science and resource management, For. Sci., 53(2), 206-219. Bishop, P., T. B. Hoey, J. D. Jansen, and I. L. Artza (2005), Knickpoint recession rate and catchment area: the case of uplifted rivers in Eastern Scotland, Earth Surf. Processes Landforms, 30(6), 767-778. Bott, R., P. Murphy, and R. Udell (2003), Learning from the forest: a fifty-year journey towards sustainable forest management, 242 pp., Fifth House Ltd., Calgary, AB. Boyce, M. S. (2006), Scale for resource selection models, Divers. Distrib., 12, 269-276. 139  Bozek, M. A., and W. A. Hubert (1992), Segregation of resident trout in streams as predicted by three habitat dimensions, Can. J. Zool., 70, 886-890. Brardinoni, F., and M. A. Hassan (2006), Glacial erosion, evolution of river long profiles, and the organization of process domains in mountain drainage basins of coastal British Columbia, J. Geophys. Res., 111(F01013). Brardinoni, F., and M. A. Hassan (2007), Glacially induced organization of channel-reach morphology in mountain streams, J. Geophys. Res., 112(F03013). Brardinoni, F., M. A. Hassan, T. Rollerson, and D. Maynard (2009), Colluvial sediment dynamics in mountain drainage basins, Earth Planet. Sc. Lett., 284(3-4), 310-319. Brooks, E. S., J. Boll, W. J. Elliot, and T. Dechert (2006), Global positioning system/GIS-based approach for modeling erosion from large road networks, J. Hydrol. Eng., 11(5), 418- 426. Brummer, C. J., and D. R. Montgomery (2003), Downstream coarsening in headwater channels, Water Resour. Res., 39(10). Burnham, K. P., and D. R. Anderson (2002), Model Selection and Multimodel Inference: a Practical Information-Theoretic Approach, 2nd ed., 488 pp., Springer-Verlag, New York. Buttle, J. (2006), Mapping first-order controls on streamflow from drainage basins: the T-3 template, Hydrol. Proc., 20(15), 3415-3422. Buttle, J. M., I. F. Creed, and R. D. Moore (2009), Advances in Canadian forest hydrology, Can. Water Resour. J., 34(2), 113-126. Cavalli, M., P. Tarolli, L. Marchi, and G. D. Fontana (2008), The effectiveness of airborne LiDAR data in the recognition of channel-bed morphology, CATENA, 73(3), 249-260. 140  Chatwin, S., D. E. Howes, J. W. Schwab, and D. N. Swanston (1994), A guide for management of landslide-prone terrain in the Pacific Northwest, 2nd ed., 220 pp., Research Branch, Ministry of Forests, Victoria, British Columbia. Chen, Y.-C., Q. Sung, and K.-Y. Cheng (2003), Along-strike variations of morphotectonic features in the Western Foothills of Taiwan: tectonic implications based on stream- gradient and hypsometric analysis, Geomorphology, 56(1-2), 109-137. Cheong, A. (1998), Hydrological assessment of basin types, 40 pp, British Columbia Ministry of Fisheries. Church, M., R. Kellerhals, and T. J. Day (1989), Regional clastic sediment yield in British Columbia, Can. J. Earth Sci., 26(1), 31-45. Church, M. (1992), Channel morphology and typology, in The Rivers Handbook : Hydrological and Ecological Principles, edited by P. Calow and G. E. Petts, pp. 126-143, Blackwell Scientific Publications, Oxford. Church, M. (2002), Geomorphic thresholds in riverine landscapes, Freshw. Biol., 47(4), 541- 557. Chust, G., I. Galparsoro, Á. Borja, J. Franco, and A. Uriarte (2008), Coastal and estuarine habitat mapping, using LIDAR height and intensity and multi-spectral imagery, Estuar. Coast. Shelf Sci., 78(4), 633-643. Clark, R. B., I. F. Creed, and G. Z. Sass (2009), Mapping hydrologically sensit ive areas on the Boreal Plain: a multitemporal analysis of ERS synthetic aperture radar data, Int. J. Remote Sens., 30(10), 2619 - 2635. Clarke, S., and K. M. Burnett (2003), Comparison of digital elevation models for aquatic data development, Photogramm. Eng. Rem. S., 69(12), 1367-1375. 141  Clarke, S. E., K. M. Burnett, and D. J. Miller (2008), Modeling streams and hydrogeomorphic attributes in Oregon from digital and field data, J. Am. Water Resourc. Ass., 44(2), 459- 477. Cover, M. R., J. A. de la Fuente, and V. H. Resh (2010), Catastrophic disturbances in headwater streams: the long-term ecological effects of debris flows and debris floods in the Klamath Mountains, northern California, Can. J. Fish. Aquat. Sci., 67(10), 1596-1610. Creed, I. F., G. Z. Sass, M. B. Wolniewicz, and K. J. Devito (2008), Incorporating hydrologic dynamics into buffer strip design on the sub-humid Boreal Plain of Alberta, For. Ecol. Manag., 256(11), 1984-1994. Crombie, G. P., and O. A. Erdman (1947), Geology, Alexo West of Fifth Meridian Alberta, NTS Mapsheets 83B5NW and 83B5SW, Geological Survey of Canada. Devito, K., I. Creed, T. Gan, C. Mendoza, R. Petrone, U. Silins, and B. Smerdon (2005a), A framework for broad-scale classification of hydrologic response units on the Boreal Plain: is topography the last thing to consider?, Hydrol. Proc., 19(8), 1705-1714. Devito, K., I. F. Creed, and C. Fraser (2005b), Controls on runoff from a partially harvested aspen forested headwater catchment, Boreal Plain, Canada, Hydrol. Proc., 19, 3-25. Dietrich, W. E., and T. Dunne (1993), The channel head, in Channel Network Hydrology, edited by K. Beven and M. J. Kirkby, pp. 175-219, John Wiley and Sons, New York. Dimitriadou, E., K. Hornik, F. Leisch, D. Meyer, and A. Weingessel (2009), e1071: Misc functions of the department of statistics (e1071), in R package version 1.5-19, edited, TU Wien, Austria. Donald, D. B., R. S. Anderson, and D. W. Mayhood (1980), Correlations between brook trout growth and environmental variables for mountain lakes in Alberta, Trans. Am. Fish. Soc., 109, 603-610. 142  Dormann, C. F., J. M. McPherson, M. B. Araujo, R. Bivand, J. Bolliger, G. Carl, R. G. Davies, A. Hirzel, W. Jetz, W. D. Kissling, I. Kuhn, R. Ohlemuller, P. R. Peres-Neto, B. Reineking, B. Schroder, F. M. Schurr, and R. Wilson (2007), Methods to account for spatial autocorrelation in the analysis of species distributional data: a review, Ecography, 30(5), 609-628. Douglas, R. J. W. (1956), Geology, Nordegg Alberta, NTS Mapsheet 83C8, Geological Survey of Canada. Downing, D. (1998), Weldwood Canada ecological land classification procedures manual, 57 pp, Timberline Forest Inventory Consultants, Edmonton, Alberta. Dragicevic, S. (2010), Modeling the dynamics of complex spatial systems using GIS, cellular automata and fuzzy sets applied to invasive plant species propagation, Geography Compass, 4(6), 599-615. Duhaime, K. (2003), Watershed delineation and stream classification for the Northern East Slopes Region, 32 pp, Facet Decision Systems Inc., Vancouver, BC. Dumanski, J., T. M. Macyk, C. F. Veauvy, and J. D. Lindsay (1972), Soil survey and land evaluation of the Hinton-Edson area, Alberta, 119 pp, Soils Division, Research Council of Alberta, Edmonton, AB. Dunham, J. B., B. E. Rieman, and G. Chandler (2003), Influences of temperature and environmental variables on the distribution of bull trout within streams at the southern margin of its range, N. Am. J. Fish. Manag., 23, 894-904. Ehlers, J., and P. L. Gibbard (2007), The extent and chronology of Cenozoic Global Glaciation, Quat. Int., 164-165, 6-20. 143  Ford, B., P. Higgins, A. Lewis, K. Cooper, T. Watson, C. Gee, G. Ennis, and R. Sweeting (1995a), Brook Trout (Salvelinus fontinalus), in Literature Reviews of the Life History, Habitat Requirements and Mitigation/Compensation Strategies for Selected Fish Species in the Peace, Liard and Columbia River Drainages of British Columbia  edited by G. Ennis and D. Rowland, pp. 1-22, Triton Environmental Consultants for Department of Fisheries and Oceans and the Province of British Columbia Vancouver, BC. Ford, B., P. Higgins, A. Lewis, K. Cooper, T. Watson, C. Gee, G. Ennis, and R. Sweeting (1995b), Rainbow trout (Oncorhynchus mykiss), in Literature Reviews of the Life History, Habitat Requirements and Mitigation/Compensation Strategies for Selected Fish Species in the Peace, Liard and Columbia River Drainages of British Columbia edited by G. Ennis and D. Rowland, pp. 1-25, Triton Environmental Consultants for Department of Fisheries and Oceans and the Province of British Columbia, Vancouver, BC. Fransen, B. R., S. D. Duke, G. McWethy, J. K. Walter, and R. E. Bilby (2006), A logistic regression model for predicting the upstream extent of fish occurrence based on geographical information systems data, N. Am. J. Fish. Manag., 26, 960-975. Gardner, B., and P. J. Sullivan (2004), Spatial and temporal stream temperature prediction: Modeling nonstationary temporal covariance structures, Water Resour. Res., 40(1), W01102. Gorsevski, P. V., P. E. Gessler, and P. Jankowski (2009), Integrating a fuzzy k-means classification and a Bayesian approach for spatial prediction of landslide hazard, in Handbook of Applied Spatial Analysis, edited by M. M. Fischer and A. Getis, pp. 653- 684. 144  Hack, J. T. (1957), Studies of longitudinal profiles in Virginia and Maryland, U.S. Geol. Survey Prof. Paper, 294-B, 45-97. Hack, J. T. (1973), Stream-profile analysis and stream-gradient index, U.S. Geol. Survey J. Res., 1, 421-429. Hancock, G. R., and K. G. Evans (2006), Channel head location and characteristics using digital elevation models, Earth Surf. Processes Landforms, 31, 809-824. Hartman, G. M. D. (2005), Quaternary stratigraphy and geologic history of the Charlie Lake (NTS 94A) map-area, British Columbia, M.Sc. thesis, 165 pp, Simon Fraser University, Burnaby, British Columbia. Hassan, M. A., M. Church, T. E. Lisle, F. Brardinoni, L. Benda, and G. E. Grant (2005), Sediment transport and channel morphology of small, forested streams, J. Am. Water Resourc. Ass., 41(4), 853-876. Hayakawa, Y. S., and T. Oguchi (2009), GIS analysis of fluvial knickzone distribution in Japanese mountain watersheds, Geomorphology, 111(1-2), 27-37. Heine, R. A., C. L. Lant, and R. R. Sengupta (2004), Development and comparison of approaches for automated mapping of stream channel networks, Ann. Ass. Am. Geogr., 94(3), 477-490. Hickman, M., and C. E. Schweger (1991), A paleoenvironmental study of Fairfax Lake, a small lake situated in the Rocky Mountain Foothils of west-central Alberta, J. Paleolimnol., 6, 1-15. Horton, R. E. (1945), Erosional development of streams and their drainage basins; hydrophysical approach to quantitative morphology, Geol. Soc. Am. Bull., 56(3), 275- 370. 145  Hupp, C. R. (1986), The headward extent of fluvial landforms and associated vegetation on Massanutten Mountain, Virginia, Earth Surf. Processes Landforms, 11(5), 545-555. Isaak, D. J., and W. A. Hubert (2000), Are trout populations affected by reach-scale stream slope?, Can. J. Fish. Aquat. Sci., 57, 468-477. Istanbulluoglu, E., D. G. Tarboton, R. T. Pack, and C. Luce (2002), A probabilistic approach for channel initiation, Water Resour. Res., 38(12), 1325. Jablonski, P. D. (1980), Pretreatment water quality of the Tri-Creeks experimental watershed, 56 pp, Alberta Energy and Natural Resources, Forest Land Use Branch, Watershed Management, Edmonton, Alberta. Jackson, L. E., G. M. MacDonald, and M. C. Wilson (1982), Paraglacial origin for terraced river sediments in Bow Valley, Alberta, Can. J. Earth Sci., 19(12), 2219-2231. Jackson, S. D. (2003), Ecological considerations in the design of river and stream crossings, paper presented at International conference on ecology and transportation, Center for Transportation and the Environment, North Carolina State University, Raleigh, North Carolina. Jaeger, K. L., D. R. Montgomery, and S. M. Bolton (2007), Channel and perrenial flow initiation in headwater streams: management implications of variability in source-area size, Environ. Manag., 40, 775-786. Klienbaum, D. G., L. L. Kupper, and K. E. Muller (1988), Applied regression analysis and other multivariate methods, 2nd ed., 718 pp., PWS-Kent Publishing Company, Boston. Kopecký, M., and Š. Čížková (2010), Using topographic wetness index in vegetation ecology: does the algorithm matter?, Appl. Veg. Sci., no-no. 146  Latterell, J. J., R. J. Naiman, B. R. Fransen, and P. A. Bisson (2003), Physical constraints on trout (Oncorhynchus spp.) distribution in the Cascade Mountains: a comparison of logged and unlogged streams, Can. J. Fish. Aquat. Sci., 60, 1007-1017. Leary, R. F., F. W. Allendorf, and S. H. Forbes (1993), Conservation genetics of bull trout in the Columbia and Klamath River drainages, Conserv. Biol., 7, 856-865. Leopold, L. B., and J. P. Miller (1956), Ephemeral streams - hydraulic factors and their relation to the drainage net, USGS Prof. Pap., 282-A, 37. Leopold, L. B., M. G. Wolman, and J. P. Miller (1964), Fluival processes in geomorphology, 522 pp., Dover Publications Inc., Mineola, N.Y. Lindsay, J. B., I. F. Creed, and F. D. Beall (2004), Drainage basin morphometrics for depressional landscapes, Water Resour. Res., 40(9), W09307. Littell, R. C., G. A. Milliken, W. W. Stroup, R. D. Wolfinger, and O. Schabenberger (2006), SAS for mixed models, 2nd ed., 814 pp., SAS Institute, Cary, NC. MacDonald, L., A. Smart, and R. Wissmar (1990), Monitoring guidelines to evaluate effects of forestry activities on streams in the Pacific Northwest, 166 pp, U.S. Environmental Protection Agency, Region 10, Seattle, Wash. Manel, S., H. Williams, and S. Ormerod (2001), Evaluating presence-absence models in ecology: the need to account for prevalence, J. Appl. Ecol., 38, 921-931. Massong, T. M., and D. R. Montgomery (2000), Influence of sediment supply, lithology, and wood debris on the distribution of bedrock and alluvial channels, Geol. Soc. Am. Bull., 112(4), 591-599. Mattingly, H. T., and D. L. Galat (2002), Distributional patterns of the threatened Niangua darter Etheostoma nianguae, at three spatial scales with implications for species conservation, Copeia, 3, 573-585. 147  McBratney, A. B., and A. W. Moore (1985), Application of fuzzy sets to climatic classification, Agric. For. Meteorol., 35(1-4), 165-185. McCleary, R. J., and M. A. Hassan (2008), Predictive modeling and spatial mapping of fish distributions in small streams of the Canadian Rocky Mountain foothills, Can. J. Fish. Aquat. Sci., 65, 319-333. McCleary, R. J., M. A. Hassan, and R. D. Moore (2008), Partial sediment budgets from two headwater catchments in the Rocky Mountain Foothills, paper presented at 34th Annual Meeting of the Canadian Geophysical Union, Banff, Alberta, Canada. McCullagh, P., and J. A. Nelder (1989), Generalized Linear Models, 2nd ed., 511 pp., Chapman and Hall, London. McMaster, R. B. (1986), A statistical-analysis of mathematical measures for linear simplification, Am. Cartographer, 13(2), 103-116. McPherson, H. J. (1975), Sediment yields from intermediate-sized stream basins in southern Alberta, J. Hydrol., 25, 243-257. Miller, D. (2003), Programs for DEM analysis, 38 pp, Earth Systems Institute, Seattle, Wash. Montgomery, D. R., and W. E. Dietrich (1988), Where do channels begin?, Nature, 336, 232- 234. Montgomery, D. R., and W. E. Dietrich (1992), Channel initiation and the problem of landscape scale, Science, 255(5046), 826-830. Montgomery, D. R., and E. Foufoula-Georgiou (1993), Channel network source representation using digital elevation models, Water Resour. Res., 29, 3925-3934. Montgomery, D. R., and J. M. Buffington (1997), Channel-reach morphology in mountain drainage basins, Geol. Soc. Am. Bull., 109(5), 596-611. 148  Montgomery, D. R. (2001), Slope distributions, threshold hillslopes, and steady-state topography, Am. J. Sci., 301(4-5), 432-454. Montgomery, D. R. (2002), Valley formation by fluvial and glacial erosion, Geology, 30(11), 1047-1050. Moore, R. D., and J. S. Richardson (2003), Progress towards understanding the structure, function, and ecological significance of small stream channels and their riparian zones, Can. J. For. Res., 33, 1349-1351. Morita, K., and S. Yamamoto (2002), Effects of habitat fragmentation by damming on the persistence of stream-dwelling charr populations, Conserv. Biol., 16(5), 1318-1323. Mouton, A. (2005), Generated stream maps using LiDAR derived digital elevation models and 10-m USGS DEM, M.Sc. thesis, 78 pp, University of Washington, Seattle, Washington. Murphey, J. B., D. E. Wallace, and L. J. Lane (1977), Geomorphic parameters predict hydrograph characteristics in the southwest, J. Am. Water Resourc. Ass., 13, 25-37. Murphy, P. N. C., J. Ogilvie, F. R. Meng, and P. Arp (2008), Stream network modelling using lidar and photogrammetric digital elevation models: a comparison and field verification, Hydrol. Proc., 22(12), 1747-1754. Murphy, P. N. C., J. Ogilvie, and P. Arp (2009), Topographic modelling of soil moisture conditions: a comparison and verification of two models, Eur. J. Soil Sci., 60(1), 94-109. National Wetlands Working Group (1997), The Canadian wetland classification system, 2nd ed., 76 pp., Wetlands Research Centre, University of Waterloo, Waterloo, Ontario. Oakes, R. M., K. B. Gido, J. A. Falke, J. D. Olden, and B. L. Brock (2005), Modelling of stream fishes in the Great Plains, USA, Ecol. Freshw. Fish., 14, 361-374. 149  Olden, J. D., D. A. Jackson, and P. R. Peres-Neto (2002), Predicitive models of fish species distributions: a note on proper validation and chance predictions, Trans. Am. Fish. Soc., 131, 329-336. Olson, S. A., and M. C. Brouillette (2006), A logistic regression equation for estimating the probability of a stream in Vermont having intermittent flow, 22 pp, USGS Scientific Investigations Report 2006-5217. Paul, A. J., and J. R. Post (2001), Spatial distribution of native and nonnative salmonids in streams of the Eastern Slope of the Canadian Rocky Mountains, Trans. Am. Fish. Soc., 130, 417-430. Perez-Pena, J. V., J. M. Azanon, A. Azor, J. Delgado, and F. Gonzalez-Lodeiro (2008), Spatial analysis of stream power using GIS: SLk anomaly maps, Earth Surf. Processes Landforms, 34(1), 16-25. Perron, J. T., J. W. Kirchner, and W. E. Dietrich (2009), Formation of evenly spaced ridges and valleys, Nature, 460(7254), 502-505. Peterson, E. E., D. M. Theobald, and J. M. Ver Hoef (2007), Geostatistical modelling on stream networks: developing valid covariance matrices based on hydrologic distance and stream flow, Freshw. Biol., 52(2), 267-279. Phillips, J. D., S. McCormack, J. Duan, J. P. Russo, A. M. Schumacher, G. N. Tripathi, R. B. Brockman, A. B. Mays, and S. Pulugurtha (2010), Origin and interpretation of knickpoints in the Big South Fork River basin, Kentucky-Tennessee, Geomorphology, 114(3), 188-198. Pike, A. S., F. N. Scatena, and E. E. Wohl (2010), Lithological and fluvial controls on the geomorphology of tropical montane stream channels in Puerto Rico, Earth Surf. Processes Landforms, 35(12), 1402-1417. 150  Pinheiro, J. C., and D. M. Bates (2004), Mixed-Effects Models in S and S-Plus, Springer, New York. Porter, M., J. Rosenfeld, and E. Parkinson (2000), Predictive models of fish species distribution in the Blackwater drainage, British Columbia, N. Am. J. Fish. Manag., 20, 349-359. Rice, S., and M. Church (1996), Bed material texture in low order streams on the Queen Charlotte Islands, British Columbia, Earth Surf. Processes Landforms, 21, 1-18. Rieman, B. E., and J. D. McIntyre (1993), Demographic and habitat requirements for conservation of bull trout, 38 pp, U.S. Dep. Agric. For. Serv. Intermountain Research Station, Ogden, UT. Rieman, B. E., J. T. Peterson, and D. L. Myers (2006), Have brook trout (Salvelinus fontinalis) displaced bull trout (Salvelinus confluentus) along longitudinal gradients in central Idaho streams?, Can. J. Fish. Aquat. Sci., 63, 63-78. Ripley, T., G. J. Scrimgeour, and M. S. Boyce (2005), Bull trout (Salvelinus confluentus) occurrence and abundance influenced by cumulative industrial developments in a Canadian boreal forest watershed, Can. J. Fish. Aquat. Sci., 62(11), 2431-2442. Roed, M. A. (1968), Surficial geology of the Edson-Hinton area, Alberta, Ph.D. thesis, 279 pp, University of Alberta, Edmonton. Rosenfeld, J. (2003), Assessing the habitat requirements of stream fishes: an overview and evaluation of different approaches, Trans. Am. Fish. Soc., 132, 953-968. Santner, D. E., and T. J. Duffy (1989), The Statistical Analysis of Discrete Data, 367 pp., Springer-Verlag, New York. SAS (2008), SAS/STAT Software: Release 9.2, edited, SAS Institute Inc., Cary, N.C., USA. Sass, G. Z., and I. F. Creed (2008), Characterizing hydrodynamics on boreal landscapes using archived synthetic aperture radar imagery, Hydrol. Proc., 22, 1687-1699. 151  Sayer, A. M., R. P. D. Walsh, and K. Bidin (2006), Pipeflow suspended sediment dynamics and their contribution to stream sediment budgets in small rainforest catchments, Sabah, Malaysia, For. Ecol. Manag., 224(1-2), 119-130. Schiess, P., F. Krogstad, and F. Damian (2004), Locating ditch relief culverts to reduce sediment delivery to streams - an interactive design tool, paper presented at Joint Conference of IUFRO 3.06 Forest Operations under Mountainous Conditions and the 12th International Mountain Logging Conference, Vancouver, B.C., Canada. Schuster, R. L. (1996), Socioeconomic impact of landslides, in Landslides: Investigation and Mitigation, edited by A. K. Turner and R. L. Schuster, pp. 12-35, Transportation Research Board, Washington, D.C. Scrimgeour, G. J., P. Hvenegaard, J. Tchir, S. Kendall, and A. Wildeman (2003), Cumulative effects of watershed disturbance on stream fish communities in the Kakwa and Simonette River Basins, Alberta., 126 pp, Alberta Conservation Association, Edmonton, Alta. Scrimgeour, G. J., P. J. Hvenegaard, and J. Tchir (2008), Cumulative industrial activity alters lotic fish assemblages in two boreal forest watersheds of Alberta, Canada, Environ. Manag., 42(6), 957-970. Selong, J. H., T. E. McMahon, A. V. Zale, and F. T. Barrows (2001), Effect of temperature on growth and survival of bull trout, with application of an improved method for determining thermal tolerance in fishes, Trans. Am. Fish. Soc., 130(6), 1026-1037. Shepard, B. B. (2004), Factors that may be influencing nonnative brook trout invasion and their displacement of native westslope cutthroat trout in three adjacent southwestern Montana streams, N. Am. J. Fish. Manag., 24, 1088-1100. 152  Smith, T. R., and F. P. Bretherton (1972), Stability and conservation of mass in drainage basin evolution, Water Resour. Res., 8(6), 1506-1529. Spillios, L. C. (1999), Sediment intrusion and deposition near road crossings in small foothill streams in west central Alberta, M.Sc. thesis, 60 pp, University of Alberta, Edmonton. Stanford, J. A., and J. V. Ward (1993), An ecosystem perspective of alluvial rivers - connectivity and the hyporheic corridor, J. North. Am. Benthol. Soc., 12(1), 48-60. Steel, E. A., B. E. Feist, D. W. Jensen, G. R. Pess, M. B. Sheer, J. B. Brauner, and R. E. Bilby (2004), Landscape models to understand steelhead (Oncorhynchus mykiss) distribution and help prioritize barrier removals in the Willamette basin, Oregon, USA, Can. J. Fish. Aquat. Sci., 61, 999-1011. Stelfox, H. A. (1981), Ecological Land Classification Red Deer - James, 115 pages pp, Alberta Energy and Natural Resources, Resource Evaluation and Planning Division, Edmonton, Alberta. Sterling, G. S. (1992), Fry emergence survival of rainbow trout, Oncorhyncus mykiss (Walbaum), following timber harvest in two foothills streams of west-central Alberta, M.Sc. thesis, 116 pp, University of Alberta, Edmonton, Alberta. Storey, R., and S. Wadhwa (2009), An assessment of the lengths of permanent, intermittent and ephemeral streams in the Auckland region, Auckland Regional Council Technical Report No. 028, 25 pp, National Institute of Water and Atmospheric Research, Auckland, New Zealand. Tarboton, D. G., R. L. Bras, and I. Rodrigueziturbe (1991), On the extraction of channel networks from digital elevation data, Hydrol. Proc., 5(1), 81-100. Tarboton, D. G. (1997), A new method for the determination of flow directions and upslope areas in grid digital elevation models, Water Resour. Res., 33(2), 309-319. 153  Tarolli, P., and D. G. Tarboton (2006), A new method for determination of most likely landslide initiation points and the evaluation of digital terrain model scale in terrain stability mapping, Hydrol. Earth Syst. Sci., 10(5), 663-677. Tarolli, P., J. R. Arrowsmith, and E. R. Vivoni (2009), Understanding earth surface processes from remotely sensed digital terrain models, Geomorphology, 113(1-2), 1-3. Tarolli, P., and G. D. Fontana (2009), Hillslope-to-valley transition morphology: new opportunities from high resolution DTMs, Geomorphology, 113(1-2), 47-56. Taylor, E. B., P. Tamkee, G. S. Sterling, and W. Hughson (2006), Microsatellite DNA analysis of rainbow trout (Oncorhynchus mykiss) from western Alberta, Canada: native status and evolutionary distinctiveness of "Athabasca" rainbow trout, Conserv. Genet., 8, 1-15. Terajima, T., T. Sakamoto, Y. Nakai, and K. Kitamura (1997), Suspended sediment discharge in subsurface flow from the head hollow of a small forested watershed, northern Japan, Earth Surf. Processes Landforms, 22(11), 987-1000. Tonina, D., and J. M. Buffington (2009), Hyporheic exchange in mountain rivers I: mechanics and environmental effects, Geography Compass, 3, 1063-1086. Toth, J. (1963), A theoretical analsysis of groundwater flow in small drainage basins, J. Geophys. Res., 68, 4795-4812. Troiani, F., and M. Della Seta (2008), The use of the stream length-gradient index in morphotectonic analysis of small catchments: a case study from central Italy, Geomorphology, 102(1), 159-168. Trotter, P. C. (2000), Headwater fishes and their uppermost habitats: a review as background for stream typing, 69 pp, Cooperative monitoring, evaluation and research committee and Washington Department of Natural Resources, Forest Practices Division, Seattle, Washington. 154  Tsukamoto, Y., T. Ohta, and H. Noguchi (1982), Hydrological and geomorphological studies of debris slides on forested hillslopes in Japan, paper presented at Recent Developments in the Explanation and Prediction of Erosion and Sediment Yield., IAHS Publ. 137. Tucker, G. E., and R. Slingerland (1997), Drainage basin responses to climate change, Water Resour. Res., 33(8), 2031-2047. Van Sickle, J., J. Baker, A. Herlihy, P. Bayley, S. V. Gregory, P. Haggerty, L. Ashkenas, and J. Li (2004), Projecting the biological condition of streams under alternative scenarios of human land use, Ecol. Appl., 14(2), 368-380. Veldman, W. (1997), Hydrologic operational manual, 26 pp, Hydroconsult EN3 Services Ltd., Calgary, Alberta. Venables, W. N., and B. D. Ripley (2002), Modern Applied Statistics with S, 4th ed., 495 pp., Springer, New York. Wang, Z., and L. A. Goonewardene (2004), The use of MIXED models in the analysis of animal experiments with repeated measures data, Can. J. Anim. Sci., 84(1), 1-11. Waters, T. F. (1995), Sediment in streams: sources, biological effects and control, 251 pp., American Fisheries Society, Bethesda, Maryland. Willgoose, G., R. L. Bras, and I. Rodriguez-Iturbe (1991), A physical explanation of an observed link area-slope relationship, Water Resour. Res., 27(7), 1697-1702. Wolter, C., and R. Arlinghaus (2003), Navigation impacts on freshwater fish assemblages: the ecological relevance of swimming performance, Rev. Fish Biol. Fish., 13, 63-89. Zadeh, J. (1965), Fuzzy sets, Inform. Control, 8, 338-353. Zevenbergen, L. W., and C. R. Thorne (1987), Quantitative analysis of land surface topography, Earth Surf. Processes Landforms, 12(1), 47-56. 155  Zimmermann, A., M. Church, and M. A. Hassan (2010), Step-pool stability: Testing the jammed state hypothesis, J. Geophys. Res., 115(F2), F02008.    156  Appendix A. Review of categorical data modelling procedures This supplement aims to make the modelling of categorical data more accessible to watershed scientists.  Three topics are reviewed: (1) response variable types and options for modelling; (2) the logit model as the base model for categorical response variables; and (3) the ordinal logistic model as a special form of the logit model.  Cumulative logit plots, an important tool to assess relations between a response variable of ordinal class and continuous type predictors, are also discussed. The four classes of response variables that can be modelled using multiple regression techniques include interval, ordinal, nominal, and binary (Table A- 1).  Classes with higher information content can be transformed to lower content classes and analyzed accordingly; however, it is not possible to work in the opposite direction (i.e., transform a lower-information variable type to one with higher information). 157   Table A- 1. Comparison of dependent variable types, appropriate multivariate methods and solution forms. Variable type Variable class Example Level of information in variable Appropriate Multivariate method Multivariate type Solution form Assumptions Quantitative Interval -3.2, 1.6, 6.4 high linear regression using least squares  single equation standard Categorical Ordinal low, medium, high low logistic regression using maximum likelihood ordinal logistic regression single equation with N 1  – 1 intercepts standard + proportional odds  Nominal red, green, blue very low  multinomial logistic regression N 1  – 1 equations standard  Binary present / absent very low  binary logistic regression single equation standard 1 N = number of categories 158  The logit (logistic) model serves as a base model for regression analysis for all the three categorical variable classes – ordinal, nominal, and binary – with different logit model forms for each.  To understand the logit model, it is important understand the relationship between probabilities, odds and logits.  Where p = the probability of an event occurring and 1-p = the probability of no event, the odds of an event (O) is  1 p O p        (A1) The logit, defined as the log of the odds, is            1 p Logit O log p           (A2) Using these relations, it is easy to translate between probabilities, odds, and logits (Table A- 2). Table A- 2. Relations between probability, odds, and logit. Probability Odds Logit 0.10 0.11 -2.21 0.20 0.25 -1.39 0.30 0.43 -0.84 0.40 0.67 -0.40 0.50 1.00 0.00 0.60 1.50 0.41 0.70 2.33 0.85 0.80 4.00 1.39 0.90 9.00 2.20  The basic form of the logit model is:  1 1 2 2     ... 1 i p p i p log X X X p                 (A3) where α is the intercept, β‟s are the logistic regression coefficients, and X1 to Xp are the predictors.  The logit model also can be solved for pi where the solution will always be between 0 and 1 as: 159    1 1 2 2   ... 1  1  p p i X X X p e                    (A4) There are different forms of the logit model to suit binary, nominal and ordinal categorical data (Table A 1).  Binary data can be estimated with the basic form of the logit model.  Nominal variables with more than two categories can be estimated by running a set of binary logit models.  Specifically, for the ith category, where the baseline category is J, the model is  1 1 2 2 ( )      ... ( ) i i i i ip p J p category log X X X p category                 (A5) The number of non-redundant sets of solution equations is J – 1. Ordinal dependent variables can be analyzed using the cumulative logit model.  For this model, categories are ordered from j = 1, …J.   This model uses a procedure that transforms the dataset containing J ordered classes into J-1 sets of dichotomous variables for analysis using the binary logit model.  This approach allows the model to force the logistic regression coefficients to be the same for all categories, but allows a different intercept for each of the J-1 sets of variables.    The cumulative logit model will be explained in three steps: (1) the production of cumulative logit plots; (2) the model equation; and (3) explanation of the proportional odds assumption. Cumulative logit plots are used to explore relations between an individual predictor and a response variable of the ordinal class.  The first step is to divide the predictor of interest into a number of groups (n).  In the example plot (Figure 3 5a), n = 20.  The lower and upper values from each of the defined groups are used to form a series of bins.  The next step is to convert the ordered class variables into J – 1 sets of dichotomous variables (Table A- 3). Table A- 3. Example of the conversion of an ordinal variable with four classes into three non- redundant cumulative variables. Process domain class L1 L2 L3 1 1 1 1 2 0 1 1 3 0 0 1 4 0 0 0  160   Cumulative variables (L1 … Lj-1), defined for J – 1 categories (Table A- 3), are created as 1 11,  1,  0L where j elseL     2 2 1,   2,  0L where j elseL       1 11,  1 ,  0 J JL where j J elseL         (A6) and the count of each cumulative variable within each bin is tallied.  Within each bin, cumulative logits are calculated from the counts of each cumulative variable.  Logits are undefined where outcome rate is 0% or 100%; therefore, Santner and Duffy [1989] recommend the addition of a small constant to both the numerator and denominator of the logit formula as  1 1 i i i L L L m log M m            (A7) where m = the number of cases and M = the number of events for each cumulative variable Li.  The cumulative logit plot is produced by plotting the midpoint value for each bin of the interval variable against its respective cumulative logit.  Each dichotomous variable (L1 … Lj-1), is plotted as a separate line.  The plots show associations only with individual predictors. Cumulative logit regression predicts the cumulative probability, Fij, that individual i is in the jth category or lower, rather than being in a higher category.  The cumulative logit model is a set of J – 1 equations with different intercepts but common regression coefficients as             1 1 2 2     ... 1 j ij p p ij F log X X X F                   (A8) Each equation in the set is used to calculate a sequential dichotomous variable.  The proportional odds assumption requires that the set of common regression coefficients function equally well across the series of dichotomous variables.  Cumulative logit plots can be examined for parallel trends to help assess this assumption.  Once the model equations are produced, the 161  predicted class for each sample in a dataset can be assigned.  First, the cumulative probabilities are converted into probabilities for individual classes as 11 L p p  2 12 L L p p p  1 1 jj L p p    Then for each sample, the category (1...j) with the highest probability (p1…pj) becomes the predicted class for that sample.  162  Appendix B. Field card 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0071632/manifest

Comment

Related Items