Forest Investment Account Forest Science Program

FOREST INVESTMENT ACCOUNT FOREST SCIENCE PROGRAM - 2004/2005 ANNUAL TECHNICAL REPORT - Project Y051293 HYDROLOGIC DECISION MAKING TOOLS FOR SUSTAINABLE FOREST MANAGEMENT IN RAIN DOMINATED COASTAL BC WATERSHEDS April 2005 Prepared By: Younes Alila, Ph.D., P.Eng.1 E-mail: [email protected] Website: faculty.forestry.ubc.ca/alila/main.html Markus Schnorbus, M.A.Sc.1 E-mail: [email protected] 1 Respectively, Associate Professor and Research Scientist, Department of Forest Resources Management, University of British Columbia, #2030 - 2424 Main Mall Vancouver, BC, Canada V6T 1Z4 ABSTRACT A physically-based distributed parameter model was used to investigate the impacts of forest management on streamflow from the rainfall-dominated 9.8 km2 Carnation Creek watershed located on the west coast of Vancouver Island, British Columbia. The model used accounts explicitly for the effects of vegetation upon runoff production and is also unique in its explicit portrayal of preferential subsurface runoff processes common to forested watersheds. Flood frequency analysis was used to quantify the impact upon the peak discharge regime as a result of both forest removal and road construction. The relative change in peak discharge is scale invariant and ranges from -4 to 30% for cut rates ranging from 10 to 100%, drainage areas ranging from 0.88 to 9.8 km2 and return periods of 0.17 to 20 years. The relative change in peak discharge is strongly correlated to both rate of cut and return period, increasing with increasing rate of cut and decreasing return period. However, statistically significant (α = 0.05) increases in peak discharge are restricted to cut rates greater than 30% and return periods less than 2 years. The effects of roads alone is to generate relative increases in peak discharge that range from -2 to 11%, an effect that is analogous to changes attributed to forest harvesting at cut rates ranging from 10 to 40%, depending upon return period. Relative changes in peak discharge increase with decreasing return period; however, changes are statistically insignificant (α = 0.05) at all return periods. The combined effect of roads and forest harvesting is to increase peak discharge more so than would have occurred with either forest harvesting or roads in isolation. However, the effect of roads and forest harvesting is only additive at the lowest return period (0.17 years); the combined effect is increasingly less than additive with increasing return period. Changes to peak discharge due to forest removal at return periods of 0.17 years and less are related to increased runoff efficiency arising from higher soil moisture whereas impacts at return periods higher that 0.17 years are affected by decreasing rainfall interception efficiency with increasing storm size. The weakening response of peak flow change to roads with increasing return period is linked to the higher subsurface flow rates associated with preferential runoff and the upper limit on the water table response to precipitation once preferential flow is activated. i Table of Contents ABSTRACT......................................................................................................................... i Table of Contents................................................................................................................ ii List of Figures .................................................................................................................... iii List of Tables ..................................................................................................................... iv 1.0 INTRODUCTION .................................................................................................. 1 2.0 STUDY AREA ....................................................................................................... 3 3.0 METHOD ............................................................................................................... 3 3.1 Distributed Hydrology Soil Vegetation Model (DHSVM)................................. 3 3.2 DHSVM Input Data and Calibration .................................................................. 6 3.3 Simulation Experiment ..................................................................................... 12 3.4 Scenarios ........................................................................................................... 13 3.4.1 Cut Rate .................................................................................................... 13 3.4.2 Scale.......................................................................................................... 14 3.4.3 Roads......................................................................................................... 14 3.5 Analysis............................................................................................................. 15 4.0 RESULTS ............................................................................................................. 21 4.1 Cut Rate Effects ................................................................................................ 22 4.2 Scale Effects...................................................................................................... 24 4.3 Road Effects...................................................................................................... 28 4.4 Scale Invariance ................................................................................................ 31 5.0 DISCUSSION ....................................................................................................... 32 6.0 CONCLUSION AND IMPLICATIONS.............................................................. 37 7.0 REFERENCES ..................................................................................................... 40 ii List of Figures Figure 1. Carnation Creek study area. ................................................................................ 4 Figure 2. Carnation Creek DHSVM application input.. ..................................................... 8 Figure 3. Model stream network with stream classes ....................................................... 10 Figure 4. Model road and stream network with culvert locations. ................................... 12 Figure 5. Hypothetical cut blocks for rate-of-cut scenarios.............................................. 14 Figure 6. Rate-of-cut scenarios with hypothetical cut block layouts and cut rates........... 15 Figure 7. Scale scenarios................................................................................................... 17 Figure 8. Road scenarios................................................................................................... 18 Figure 9. Plotting positions of PDS peak flow sample for basin-wide control scenario .. 22 Figure 10. Plotting positions of PDS peak flow sample for basin-wide control scenario. 23 Figure 11. Sample autocorrelation of peak flow PDS for basin-wide control scenario.. . 24 Figure 12. Flood frequency comparison of management scenario RC100 to Control. .... 26 Figure 13. Relative change in flood quantiles with rate-of-cut for Tp = 0.17, 0.5, 2, 10, and 20 years.. ............................................................................................................ 27 Figure 14. Relative change in flood quantiles with rate-of-cut showing scale effects. .... 30 Figure 15. Diagnostic graphs of linear regression relationship between ∆QTp dependent upon ROC for Tp = 0.17 years................................................................................... 33 Figure 16. Relationship between ∆QTp and ROC for all drainage areas showing individual values and fitted linear regressions ........................................................................... 34 iii List of Tables Table 1. Pre-management Structural Characteristics of DHSVM Forest Classes. ............. 7 Table 2. Constant Vegetation Parameters........................................................................... 7 Table 3. Soil Model Parameters.......................................................................................... 9 Table 4. Channel Drainage Properties. ............................................................................ 10 Table 5. Road Drainage Properties. ................................................................................. 11 Table 6. Description of Management Scenarios ............................................................... 16 Table 7. Control Quantile Magnitude and Relative Quantile Change for Cut Rate Management Scenarios. ............................................................................................ 25 Table 8a. Control Quantile Magnitudes for Scale Effects Management Scenarios with Harvesting in C Tributary ......................................................................................... 28 Table 9a. Control Quantile Magnitudes for Scale Effects Management Scenarios with Harvesting in E Tributary. ........................................................................................ 29 Table 10. Control Quantile Magnitude and Relative Quantile Change for Road Effects Management Scenarios. ............................................................................................ 31 Table 11. Parameters of Linear Regression relationship between ∆QTp and ROC for scenarios 1 through 30 .............................................................................................. 33 iv 1.0 INTRODUCTION Forest harvesting practices in the coastal mountains of British Columbia (BC) have the potential to alter the streamflow regime by increasing the magnitude, duration, and frequency of flows. Such shifts in the regime of what could be geomorphically and ecologically significant events can have substantial consequences arising from changes to the distribution of water, energy, sediment, and nutrients within the channel and its associated riparian area [Ziemer and Lisle 1998]. Over time, the altered trajectory of channel development may seriously jeopardize the availability and quality of streamassociated resources. Therefore, forest planners require knowledge of how forest operations may affect streamflow regime of a watershed. Despite the forest hydrology community’s long history with this issue, considerable uncertainties associated with the effects of forest management on streamflow remain. The considerable conjecture and heated debate in the forest hydrology community regarding the effect of forest management on peak flows [Jones and Grant 1996; Thomas and Megahan 1998; Beschta et al. 2000] is perhaps the best manifestation of the extent of uncertainties associated with the effect of forest management on streamflows. It has been recognized that a better understanding of the effect of forest management on streamflow characteristics can be accomplished by synergistically supplementing experimental results with numerical modelling [Ziemer et al. 1991; Dunne 2001]. It has been demonstrated that hydrology models can alleviate some of the problems associated with traditional paired watershed studies by acting as a control to filter out effects of climate variability [Bowling et al. 2000] and are useful for linking forest management effects to controlling physical processes [Whitaker et al. 2002; Beckers et al. 2004; Schnorbus and Alila, 2004]. Models can also be used to assess how the ability of roads to affect catchment scale storm response is related to the arrangement of road segments relative to hillslopes and stream network [Bowling and Lettenmaier 2001; Tague and Band 2001; Wemple and Jones 2003]. However, physically based spatially distributed hydrologic models such as DHSVM can only be used to quantify the effects of forest management after careful testing of the various components of the hydrologic cycle within the model using carefully designed field experiments [Beven 1989; Grayson et al. 1992]. Carnation Creek is a site of an ongoing multi-disciplinary study regarding the impact of forestry activities on small stream salmonid ecosystems since 1971 [Hartman and Scrivener 1990] and has a long-term hydro-meteorological database that offers a rare opportunity for evaluating the worth of numerical modelling as a means of quantifying and assessing forest management effects on hydrology in the maritime coastal region of BC. In this project, we used the Distributed Hydrology Soil Vegetation Model (or DHSVM) developed originally by the University of Washington in collaboration with the US Pacific Northwest National Laboratory (BATELLE) [Wigmosta et al. 1994]. Since then the model has been further developed and tested by the University of Washington, 1 BATELLE, and researchers at the University of British Columbia. DHSVM is a physical parameter model that explicitly estimates the spatial distribution of moisture, energy fluxes, and runoff generation by subdividing the model domain into small computational grid elements using the spatial resolution of an underlying digital elevation model (DEM). Additionally, using surface and subsurface flow measurements at Carnation Creek combined with general conceptual ideas of macropore flow development, a numerical representation of preferential hillslope runoff has been incorporated in a modified version of DHSVM. Unlike the original matrix flow only model, the modified DHSVM is able to simultaneously explain observed rapid runoff and water table response at Carnation Creek under a range of antecedent conditions [Beckers and Alila 2004]. We believe the rapid preferential runoff response captured by this modified DHSVM is a common characteristic of watersheds throughout coastal British Columbia. We subsequently used this modified DHSVM model to reconstruct the effects of the complex management history of Carnation Creek. In the process, we tested the performance of the modified DHSVM model under both disturbed and undisturbed conditions. This work provides insight into one of the most heated and controversial debates in forest hydrology (namely the effect of timber harvesting and roads on small vs. large peak flows). Our results indicate that the percent peak flow increases due to harvesting and roads decrease as storm event size increases and thus concur with conclusions drawn by Thomas and Megahan [1998; 2001] while contradicting those of Jones and Grant [1996; 2001]. This work has already been submitted and is currently under review by Water Resources Research [Beckers et al. 2004]. In this project we used the DHSVM hydrology model of Carnation Creek to advance the work of Beckers et al. 2004 to quantify the effect of forest management scenarios other than the one implemented in the actual experiment itself. We have addressed the following operational relevant issues: i) ii) iii) Historic logging at Carnation amounts to ~50% of the 9.8 km2 drainage harvested over a 30-year period. How would the hydrology have been affected if the watershed had been subject to more aggressive management (i.e. logging 60%, 70%, etc)? While in practice we may not log this much out of any drainage this investigation allow us to stress the watershed to the limit and quantify worst case scenarios. Historic cut blocks at Carnation Creek have occurred throughout the basin. What if logging covered only headwater basins? How would the magnitude of impact propagate downstream? The headwater basins may be affected but how far does the impact persist downstream? In paired watershed studies, usually only the combined effect of forest removal/roads is quantified. We will use the DHSVM to separate the effects of roads and forest removal. This project extends the worth of the millions of dollars that have historically been invested to study the effects of only one management scenario at Carnation Creek. 2 2.0 STUDY AREA Carnation Creek watershed is located on the west coast of Vancouver Island and drains into Barkley sound, near the town of Bamfield, British Columbia (see Inset, Figure 1). The basin area of 9.8 km2 features rugged topography with elevation ranging from sea level to 900 m and basin slopes of 40 to over 80%; although the lower 3 km of Carnation Creek flows through a relatively wide (50 to 200 m) valley bottom (Figure 1). Slope soils, which are underlain by watertight bedrock of volcanic origin, are derived from colluvial materials and are of gravely sandy-loam and loamy-sand texture and classified primarily as Orthic Ferro-Humic podzols. Regosols occur on steep rocky areas and around the more active alluvial channels [Oswald, 1982]. Soils are generally shallow, with a mean depth of 0.75 m, and are highly permeable [Hetherington, 1982]. The watershed lies entirely within the Coastal Western Hemlock (CWH) biogeoclimatic zone and prior to management was covered by climax forest composed primarily of western hemlock (Tsuga heterophylla), western red cedar (Thuja plicata), and amabilis fir (Abies amabilis), with some Douglas fir (Pseudotsuga menziesii) in drier sites, Sitka spruce (Picea sitchensis) in the valley bottom, and red alder (Alnus rubra) along the stream margins [Oswald, 1982]. About 50% of the basin was clearcut logged in various stages between 1971 and 1990 [Hartman and Scrivener, 1990]. The climate of the study site is strongly influenced by the Pacific Ocean, with mild, wet winters and mild summers. Annual precipitation ranges from 2100 mm at sea level to over 4800 mm at high elevation, with 95% falling as rain. Over 75% of precipitation falls during October through March during frequent frontal storms. Some snow occurs at high elevations in most years, but it disappears quickly. For the most part, Carnation Creek lies below the transient snow zone and streamflow closely follows rainfall patterns and exhibits flashy runoff [Hetherington, 1982]. 3.0 METHOD 3.1 Distributed Hydrology Soil Vegetation Model (DHSVM) The original DHSVM consists of a two-layer canopy module for interception and evapotranspiration, a multi-layer unsaturated soil module, a saturated subsurface flow module, an energy balance module for canopy and ground snow accumulation and ablation, a surface overland flow module and a channel flow routing module. Onedimensional energy and water balances are solved individually over a matrix of 10-m grid cells in the model domain at an hourly time step, considered the best representation of the diurnal meteorological fluctuations and rapid streamflow response that occurs within the study area. Each cell exchanges water with its adjacent neighbours, resulting in a three-dimensional redistribution of surface and subsurface water across the watershed. 3 Figure 1. Carnation Creek study area. 4 A complete description of model structure and governing equations is given by Wigmosta et al. [2002], although the canopy, unsaturated soil, saturated soil, and channel routing modules are summarized in the following paragraphs. Each grid cell may contain an overstory canopy and either an understory or bare soil. The overstory may cover all or some fraction of the grid cell, whereas the understory, if present, covers the entire cell area. Both the overstory and understory may contain wet (from rainfall interception) and dry fractions. The partitioning of the canopy into wet and dry fractions is based on the volume of intercepted water and the leaf-area-index (LAI) [Wigmosta et al., 2002]. All rainfall is automatically stored on each canopy until its maximum interception capacity is reached, at which point any excess precipitation passes as throughfall with no attenuation. If present, intercepted moisture is evaporated at the potential rate while transpiration from the dry fraction is modelled using a PenmanMonteith approach. Every grid is specified to have one or more root soil layers and each vegetation layer removes water from one or more soil layers based on specified root zone fractions. The number of soil layers in the model equals the number of root soil layers plus one and the bottom layer extends from the root zone to bedrock. The properties of the bottom layer are set equal to that of the deepest root soil layer, but evapotranspiration cannot remove water from this layer. Direct evaporation from the top soil layer is only simulated in the absence of an understory, otherwise soil evaporation is ignored. Moisture movement through unsaturated and saturated soil in the original DHSVM was represented strictly as matrix flow [Wigmosta et al., 1994]. This module has been subsequently modified by Beckers and Alila [2004] to account for the rapid preferential runoff that occurs at Carnation Creek. This modified model uses multiple soil moisture stores to better account for non-linear runoff response to rainfall events by representing both matrix and preferential subsurface runoff. The description of the preferential flow modification that follows for the remainder of the paragraph is summarized from Beckers and Alila [2004]. During a rainfall event infiltration into the soil column is initially assumed as matrix flow and macropore channels do not become active until the infiltration capacity of the soil is exceeded. The infiltration capacity is estimated using the Green-Ampt equation, which accounts for soil moisture at the start of the precipitation event, cumulative infiltration since the start of the event, and wetting front suction. Initial soil moisture and cumulative infiltration are reset for rainfall events separated by 8 hours or more. Any moisture in excess of the infiltration capacity is assumed to collect in the macropore system, where it vertically by-passes the soil matrix and is routed instantaneously to the soil bedrock interface. By-pass infiltration is only activated during rain or rain-on-snow events as typical snowmelt rates are well below the infiltration capacity of the soil; hence, all snowmelt infiltrates into the soil matrix. In the initial stages of a rainstorm the preferential flow network is assumed to be spatially disconnected and is routed down slope at a velocity v1. As the water table rapidly rises, faster preferential flow at a velocity v2 >> v1 is initiated at a point when the cumulative rainfall depth exceeds a threshold value R*, halting the rapid water table rise. The preferential flow velocities v1 and v2 are assumed to vary spatially according to slope gradients. Preferential flow routing is performed at a shorter time step than all other 5 model components and is set based on the maximum v2 estimated for the model domain. Lateral subsurface flow in the saturated soil matrix occurs as per the original model, which uses the quasi-three dimensional routing of Wigmosta and Lettenmaier [1999]. In the soil matrix the rate of lateral saturated flow is governed by the soil transmissivity, which decreases by a power law relationship with soil depth, the water table slope, and the flow width, which is controlled by the grid cell resolution and aspect of the water table. In the modified DHSVM there is no direct mixing of water between the matrix and preferential stores, the only feedback being through variations in the water table depth. Flow in stream channels and road drainage ditches is routed using a cascade of linear reservoirs. The stream and, if present, road channel networks are constructed of individual segments, each with its own hydraulic parameters. Each segment is treated as a reservoir of constant width with outflow linearly related to storage [Wigmosta et al., 2002]. The segments of the stream channels and road drainage networks receive water laterally from the cells through which they pass from overland flow and intercepted subsurface flow (i.e. when the channel or ditch bottom intersects the water table). Outflow from a segment may be routed into another segment or exit the watershed. Outflow from road drainage segments may also input back into a watershed cell at a culvert location, where it is added to the surface water of the cell and will be either infiltrated or routed as overland flow. Road segments may also drain into stream segments at river crossings. 3.2 DHSVM Input Data and Calibration The preparation of input data for the Carnation Creek application of DHSVM is described fully by Beckers and Alila [2004] and Beckers et al. [2004] and is only summarized herein. The pre-management vegetation cover, based on vegetation classification and mapping conducted by Oswald [1982], consists of four vegetation classes: channel communities, lower slopes forest, middle slopes forest, and upper slopes forest (Figure 2). The physical characteristics of these forest classes are given in Table 1 and Table 2; the properties in Table 1 vary between forest classes while those of Table 2 are constant between classes. A spatially uniform single root zone layer with a depth of 0.5 m has been assumed for the watershed. However, the depth to bedrock, D (i.e. total soil depth), is highly variable and is represented using a deterministic-stochastic approach. The soil depth is initially divided into three classes based on mapping by Oswald [1982], areas of shallow soil and exposed bedrock (D = 0.5 m, root zone depth), areas of shallow to moderately deep soils (D = 0.75 m, average soil depth in the watershed), and areas of moderately deep soils (D = 1.2 m). Upon the three-way classification was superimposed a Gaussian random soil depth field, with mean of zero and variance s2 = 0.075 m2 (based on soil depth variations between 12 piezometer located in small tributary sub-basin) (see Beckers and Alila [2004] for details). A minimum soil depth of 0.5 m was maintained. The distribution of soil depth for the watershed is given in Figure 2. The soil model parameters are given in Table 3. 6 Table 1. Pre-management Structural Characteristics of DHSVM Forest Classesa. Heightb Crown LAId (m) Closurec (m2/m2) Class Description 1 2 3 4 Channel Communities (Western Hemlock/Cedar/Spruce) Lower Slopes (primarily Western Hemlock) Middle Slopes (Western Hemlock/Douglas Fir) Upper Slopes (Western Hemlock/Douglas Fir) 43 0.75 6.4 43 37 25 0.75 0.60 0.50 6.4 6.4 6.4 a From Beckers et al. [2004] Height of dominant overstory trees c Proportion of land area covered by overstory (vertical projection) d Leaf area index b Table 2. Constant Vegetation Parametersa. Item Parameter Fractional trunk space height (i.e. relative height of 1 canopy bottom above ground) (-) Snow interception efficiency as proportion of 2 precipitation (-) 3 Maximum snow interception capacity (m SWE) Minimum melt needed for mass release of intercepted 4 snow (m SWE) 5 Snow mass release/drip ratio (-) LAI multiplier for maximum rain interception capacity 6 (m) 7 Aerodynamic attenuation coefficient (-) 8 Albedo (-) 9 Solar radiation attenuation coefficient (-) 10 Maximum stomatal resistance (s/m) 11 Minimum stomatal resistance (s/m) Soil moisture level below which transpiration is 12 restricted (-) 13 Vapour pressure deficit forcing stomatal closure (Pa) 14 Critical light level limiting photosynthesis (W/m2) 15 Root fraction in soil layer 1 16 Depth of root zone layer a Overstory Understory 0.2 NAb 0.6 NAb 0.04 NAb 0.002 NAb 0.4 NAb 0.0015 0.0015 3.0 0.1 1.7 5000 400 NAc 0.15 NAd 5000 70 0.1 0.1 4000 30 1.0 0.5 4000 30 1.0 0.5 All items except 9 described in Wigmosta et al. [2004]; item 9 described by Thyer et al. [2004]; Table and all values taken from Beckers et al. [2004] b Understory assumed to be completely covered when snowpack is present c Logarithmic wind profile assumed below overstory canopy d Energy exchange between soil and understory not modelled; soil temperature set equal to air temperature 7 a) b) c) Figure 2. Carnation Creek DHSVM application input showing a) digital elevation model, b) soil depth, and c) vegetation classes. All three panes show derived stream network. (solid blue lines). 8 Table 3. Soil Model Parametersa. Itemb 1 2 3 4 5 6 7 8 9 10 11 12 a b Parameter Porosity (-) Field capacity (-) Lateral hydraulic conductivity of soil matrix (m/s) Vertical hydraulic conductivity of soil matrix (m/s) Power law exponent (-) Effective surface soil hydraulic conductivity for infiltration into soil matrix (m/s) Threshold rainfall depth for initiating fast preferential flow, R* (m) Slow preferential flow velocity, v1 (m/s) Fast preferential flow velocity, v2 (m/s) Pore size distribution index (-) Air bubbling pressure (m) Wilting point (-) Value 0.43 0.29 8.0 x 10-4 1.9 x 10-6 3.2 2.5 x 10-7 0.05 6.0 x 10-3 6.0 x 10-2 0.38 0.15 0.14 From Beckers et al. [2004] Items 1 through 7 derived via calibration The model stream network was derived directly from the 10-m DEM and is shown in Figure 2. The derived stream network has a total length of 60 km and a channel density of 6.1 km/km2; the channel has been divided into the six classes shown in Figure 3, with drainage parameters given in Table 4. The modelled road network was derived from the mapped road layout (provided by MacMillan Bloedel Ltd., now Weyerhaueser Company Ltd.) and road slopes and cutbank heights estimated using ARC/INFO analysis tools. The road network in Carnation Creek has a total length of 34 km and road density is 3.5 km/km2. The road network is composed of 27 classes, with drainage parameters given in Table 5. Stream channel and road drainage geometry (depth and width) were estimated from field surveys while hydraulic parameters were taken from the literature [Shen and Julien, 1993]. The model has a total of 254 culverts; 104 are stream crossings, the remainder are ditch relief culverts. Of the total culverts, 137 were surveyed; the remaining 117 culverts were inferred through GIS analysis and imposed at road/stream crossings and areas of topographic convergence where water runs the risk of overtopping the road (Figure 4). Direct connectivity between the road and stream network occurs only at stream/road crossings; instances of discharge from relief culverts draining directly to stream channels via hillslope gullies (i.e. Jones et al. [2000]) was not incorporated into the analysis. DHSVM requires time series input, at a temporal resolution of the model time step, of precipitation, air temperature, relative humidity, wind speed, solar beam and diffuse radiation, longwave radiation, temperature lapse rate, and precipitation gradient at one or more station locations within the study area. The digital elevation data are used to model 9 Figure 3. Model stream network with stream classes. Parameters for each class are given in Table 4. topographic controls on solar radiation (terrain shadowing and skyview), precipitation, and air temperature. The preparation of the meteorological input data is given by Beckers and Alila [2004]. The DHSVM model was forced with hourly meteorological data collected over the period October 1972 to December 1990 from six climate stations. Precipitation measurements are available from all six stations (A, C, D, EH, F, and L) (Figure 1) and distributed spatially using monthly-average precipitation gradients estimated from monthly precipitation totals at each station. Air temperature and relative humidity were only observed at four stations (A, C, D, and EH). Hourly temperature lapse rates for all four stations were calculated using the A-D, C-D, and EH-D station pairs and temperature records were estimated for stations F and L using the temperature and lapse rates from stations C and EH, respectively. Relative humidity for stations F and L was estimated by directly transferring relative humidity observations from stations C and EH, respectively. Solar radiation and wind speed, observed only at station A, are assumed to be spatially uniform across the watershed; with solar radiation adjusted Table 4. Channel Drainage Properties. Channel Class 1 2 3 4 5 6 Channel Width (m) 1.0 1.0 2.0 4.0 4.0 8.0 Bank Height (m) 0.25 0.25 0.25 0.5 0.5 1.0 Manning’s Total Length Roughness (km) 0.05 47.1 0.14 3.1 0.04 4.4 0.03 1.2 0.06 1.0 0.02 2.9 10 Table 5. Road Drainage Properties. Road Class 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 Ditch Width (m) 1.0 1.0 1.5 1.5 2.0 1.0 1.0 1.5 1.5 2.0 1.0 2.0 1.0 1.0 1.5 1.0 1.5 2.0 1.5 2.0 1.0 1.0 1.5 1.0 1.5 1.0 1.5 Ditch Depth (m) 0.25 0.5 0.5 1.0 1.0 0.25 0.5 0.5 0.75 0.75 1.0 1.0 0.25 0.5 0.5 0.75 0.75 0.75 1.0 1.0 0.25 0.5 0.5 0.75 0.75 1.0 1.0 Manning’s Roughness 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 Total Length (m) 255 2432 281 324 100 3817 6586 781 6391 280 253 1251 611 3414 337 583 1250 126 864 913 100 872 622 354 336 438 244 locally for topographic effects. The temperature threshold used to determine the fraction of precipitation that falls as rain or snow was set such that mean snow contribution was 10% of total annual precipitation for 1972 – 1990. The reader is referred to Beckers and Alila [2004] for details regarding the calibration and validation of the Carnation Creek DHSVM application; what follows is only a summary of that work. The soil subsurface flow model parameters (matrix and preferential flow) were optimized using simultaneous calibration of hourly streamflow measured at B, C, E and H weirs (Figure1). Calibration was conducted using a single storm event on 2-3 February 1974. The parameterization was validated using the entire 1972 to 1990 streamflow record at weirs C and E. Validation was based on three performance measures [Nash and Sutcliffe, 1970]: the volume error δV/V between observed and simulated flow, 11 Figure 4. Model road and stream network with culvert locations. Road network as constructed to access cutblocks harvested from 1969 to 1993. model efficiency E! (which compares both the volume and shape of observed and calculated discharge), and the coefficient of determination D! (which compares the timing of flow)? Values of δV/V close to zero and E! and D! close to unity are indicative of good model performance. For C weir annual performance statistics range from -0.16 to 0.09 for δV/V, 0.69 to 0.86 for E! and 0.84 to 0.93 for D!. For E weir δV/V ranges from 0.18 to 0.13, E! ranges from 0.57 – 0.87 and D! ranges from 0.78 to 0.94. These performance measures suggest that the Carnation Creek DHSVM application is capable of simulating the hydrologic effects of seasonal, inter-annual and spatial variations in weather conditions and runoff processes. 3.3 Simulation Experiment Towards our goal of developing the hydrologic decision making tools for sustainable forest management in rain dominated coastal BC watersheds, we have used long-term computer simulations spanning periods of several decades. Forest management impacts on water quantity were quantified through an analysis of managed (disturbed) and undisturbed (unmanaged) watershed conditions. In general, our study approach adopts the paired-watershed experimental design, where we have used the hydrologic model to simulate both a control (current forest cover) and a treatment (hypothesized harvesting) watershed subject to an ensemble of identical meteorological forcing, with both watersheds based on the same real-world prototype of Carnation Creek. Meteorological forcing utilized the 18-year meteorological record observed at the six climate stations from 1972 to 1990. The experimental steps are as follows: 12 i) ii) iii) iv) Simulate the response of the watershed under fully mature forest across the whole landscape using the 18-year climate data as input to the DHSVM hydrology model; Simulate the response of the watershed to each of the forest management scenarios with the same climate data as in step i above; Compare the long-term characteristics of streamflows of each management scenario to the reference streamflows of the undisturbed conditions of step i using appropriate statistical testing; and Use the results in step iii to establish relations between the hydrologic effects of forest management and design parameters such as cut level, spatial and temporal distributions of cuts, forest road characteristics. The numerical simulation was utilized to generate a sample of peak discharge events by running an 18-member ensemble of annual streamflow simulations for each scenario. From a practical perspective each scenario was conducted as a continuous 18-year simulation and run at an hourly resolution for the period 1 October, 1972 to 31 December, 1990. Initial conditions (soil moisture, interception storage, snow water equivalent, and channel network storage) for 1 October 1972 were derived individually for each scenario by initializing the model from 1 October 1971 to 30 September 1972 by repeating the meteorology of 1 October, 1972 to 30 September 1973. For each successive year within a scenario initial conditions for year n (00:00 hours on 1 October) were based on simulated conditions at the end of year n – 1 (24:00 hours on 30 September). Such an approach allowed for the initial conditions for year n (n = 1, 2,…, 18) to vary between scenarios as a function of management activities. Between scenarios forest cover and road parameters were adjusted to reflect different forest management activities while all remaining basin parameters were held constant. Forest and road parameters were kept constant during each scenario (i.e. no re-growth, vegetation succession or road deactivation). 3.4 Scenarios The 33 management scenarios were designed to address three major issues: 1) the effect of different rates of clearcut harvesting for the entire 9.8 km2 drainage, 2) the effect of scale (drainage area) on forest removal, and 3) the individual effects of roads versus forest removal. In the current study we are concerned with the detection of the immediate impact of forest harvesting upon the magnitude and frequency (i.e. regime) of annual peak discharge and did not incorporate vegetation recovery or road deactivation. The possible effect of soil compaction due to skidding was also ignored. 3.4.1 Cut Rate The effect of cut rate was assessed using 11 hypothetical clearcut harvest scenarios with cut rates ranging from approximately 10 to 100% of basin area in increments of roughly 10%. Roads were not included in these scenarios. The layout of cut blocks was determined randomly: 1) the watershed area was divided into cutblocks by first using the 13 size, location and layout of the 17 historical cutblocks, and then subdividing the remaining mature forest into roughly 25 additional cutblocks based on forest class (Figure 5), and 2) cutblocks were selected randomly without replacement, such that harvest rates accumulated from 0 to 100% (Figure 6). An additional scenario of approximately 50% cut rate was designed by clearcut harvesting all historical cut blocks (harvested prior to 1991) simultaneously. Assessment of hydrologic impacts was based on discharge sampled at the basin outlet (i.e. drainage area of 9.8 km2). Table 6 summarizes this group of scenarios. 3.4.2 Scale Nineteen scenarios were designed to assess how the effects of harvesting in the headwaters of Carnation Creek propagate downstream, in other words, how far downstream the effects of harvesting persist. These 19 scenarios were based on the clearcut harvesting of headwater sub-basins (i.e. along the drainage divide) in both the C and E tributary watersheds at a cut rate of 10% of total basin area (9.8 km2). Hourly discharge was sampled at various locations along the stream network from the outlet of the harvested sub-basins (i.e. cut rate of 100%) to the watershed outlet (cut rate of 10%). Discharge was sampled for drainage areas ranging from 881 to 9830 ha; producing 7 and 13 scenarios for sub-basins C and E, respectively (Figure 7 and Table 6). Roads were not included in these scenarios. 3.4.3 Roads Three scenarios were constructed to specifically examine the effect of roads upon the peak discharge regime. These scenarios utilized the road network that was constructed to Figure 5. Hypothetical cut blocks for rate-of-cut scenarios. 14 Figure 6. Rate-of-cut scenarios with hypothetical cut block layouts and cut rates of a) 9%, b) 22%, c) 30%, d) 40%, e) 50%, f) 60%, g) 70%, h) 81%, i) 90%, j) 100%. Panel k) shows a cut rate of 52% with all historical cut blocks. access cutblocks harvested between 1969 and 1993. The effect of roads alone was first examined by applying the historical road network without cut blocks (i.e. using only premanagement forest cover). The relative effect of roads versus forest removal was examined by simulating the road network in conjunction with two of the rate-of cut scenarios: 1) 52% rate of cut using the historical (prior to 1991) cutblocks, and 2) a 100% rate of cut (Figure 8 and Table 6). 3.5 Analysis The impact of forest management activities was quantitatively assessed by comparing the flood frequency of each scenario to that of the control, or reference, scenario. As a result of scale invariance, the rate of cut and road scenarios all used the same control scenario (i.e. discharge at outlet of 9.8 km2 Carnation Creek basin with pre-management forest cover), whereas each of the 19 scale scenario required its own control scenario (i.e. premanagement forest cover with streamflow estimated at each sample location). The details of the process are as follows: 15 Table 6. Description of Management Scenarios No. Name Roads (km) 1 2 3 4 5 6 7 8 9 10 11 RC009 RC022 RC030 RC040 RC050 RC060 RC070 RC081 RC090 RC100 RC052Op 0 0 0 0 0 0 0 0 0 0 0 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 SC100C SC071C SC062C SC011C SC010C SC009aC SC009bC SC100E SC085E SC040E SC035E SC028E SC022E SC018E SC016E SC012E SC011E SC010aE SC010bE 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 31 32 33 RD000 RD052Op RD100 34 34 34 a b Basin Area Harvest Area (ha) (ha)a Cut rate scenarios 9830 884 9830 2162 9830 2948 9830 3930 9830 4913 9830 5797 9830 6778 9830 7959 9830 8843 9830 9830 9830 5112 Scale scenarios 881 881 1240 881 1410 881 8050 881 9000 881 9580 881 9830 881 957 957 1130 957 2430 957 2740 957 3430 957 4400 957 5290 957 6170 957 8050 957 9000 957 9580 957 9830 957 Road scenarios 9830 0 9830 5112 9830 9830 Drainage area upstream of sample location As percent of forest area upstream of sample location 16 Rate of Cut (%)b Reference Figure 9 22 30 40 50 60 70 81 90 100 52 Figure 6a Figure 6b Figure 6c Figure 6d Figure 6e Figure 6f Figure 6g Figure 6h Figure 6i Figure 6j Figure 6k 100 71 62 11 10 9 9 100 94 40 35 28 22 18 16 12 11 10 10 Figure 7a Figure 7a Figure 7a Figure 7a Figure 7a Figure 7a Figure 7a Figure 7b Figure 7b Figure 7b Figure 7b Figure 7b Figure 7b Figure 7b Figure 7b Figure 7b Figure 7b Figure 7b Figure 7b 0 52 100 Figure 8a Figure 8b Figure 8c a) b) Figure 7. Scale scenarios showing 10% basin-wide rate of cut in headwaters of a) Ctributary, and b) E-tributary. Harvest areas shown as red and alternating white/shaded areas show accumulated drainage area when moving downstream between consecutive sample points (black dots). 17 a) b) c) Figure 8. Road scenarios showing road network in conjunction with cut rates of a) 0% (pre-management vegetation; scenario RD000), b) 52% using all historical (prior to 1991; scenario RD052Op) cut blocks, and c) 100% (scenario RD100). 18 i) ii) iii) iv) v) vi) The lowest of the 18 annual maximum flood peaks, Qa, min for each respective control scenario was established; A partial duration series (PDS) of flood peaks for each control scenario was extracted using Qa, min as the threshold level, Q*, generating a sample of size N; Q* was adjusted for each management scenario such that the extracted PDS sample size was equivalent to the N for the control scenario (necessary for assessing statistical significance; see below); Non-parametric flood frequency analysis was performed on each scenario and flood quantiles, QTp, estimated for partial duration return periods, Tp, of 0.17, 0.2, 0.4, 0.5, 0.75, 1, 2, 5, 10, 15, and 20 years; The relative change in flood quantile, ∆QTp, was estimated for each management scenario using ∆QTp = (QTpm-QTpc)/QTpc, where superscripts m and c refer to management and control scenarios, respectively, and The ∆QTp were assessed for statistical significance using the sampling variability of the non-parametric control quantile estimator. Note that in this report p and Tp are often used interchangeably; the two values are related by Tp = 1 / [λ (1-p)], where λ is the average arrival rate of peak flows > Q* per year, estimated as (N / 18) yr-1. Analysis of the generated flood frequencies resorted to using the PDS (as opposed to the more common annual maximum series) for two reasons: i) due to the multiple occurrence of flood events within a year in coastal watersheds, the use of the annual maximum flood series would tend to ignore a significant amount of information regarding the flood frequency regime (conversely, the PDS extracts more information out of the rather small 18-year sample period), and ii) prior research in coastal watersheds [i.e. Thomas and Megahan 1998; Beschta et al. 2000] as well as preliminary results from Beckers et al. [2004] suggests that only the most frequent (i.e. Tp < 2 years) events tend to respond to forest harvesting. Prior to the extraction of the PDS, each hourly discharge time series was transformed using a 24-hour moving average smooth, the purpose of which was to reduce fluctuations in the hourly discharge time series that may manifest as local, but spurious, peaks. The identification of flood peaks was keyed to the smoothed time series, although the magnitude of the flood peak was based on the original time series. A second constraint, in the form of a 48-hour PDS bandwidth, was also imposed. The purpose of the bandwidth was to extract only the largest untransformed event (that exceeds Q*) within a 48-hour window, the assumption being that peak discharge events that occur within 48-hours of each other are correlated via antecedent moisture conditions. Preliminary analysis of the generated PDS revealed that the flood peaks did not lend themselves to normal parametric flood frequency analysis, due mainly to the observation that the peak appear to derive from at least two distinct populations. As such, the effect of forest management upon the peak discharge regime was quantified using non-parametric food frequency analysis [i.e. Adamowski, 1985]. Let X1, X2, …, XN be independently and 19 identically distributed flood peaks with distribution function F, and let X(1) ≤ X(2) ··· ≤ X(N) denote the corresponding order statistics. Then the quantile function Q is the leftcontinuous inverse function of F such that Q(p) = inf {x : F(x) ≥ p} where 0 < p < 1. An estimate of Q(p) can be obtained using the kernel quantile estimator KQp, given as [Sheather and Marron, 1990]  i − 12  KQ p = ∑ K h  − p  X (i )  N  i =1   N  j − 12    K p − ∑ h   j =1  N  N (1) where Kh = h-1 K(· / h) where h is the bandwidth and K(·) is a kernel density function symmetric about 0. In effect, (1) estimates Q(p) by applying a weighted kernel smoothing to the sample flood peaks, such that greater weight is placed on values of X(i) where (i0.5)/N is closest to p. It is desirable that K(·) have compact support to minimize the effect of using a bounded domain of p on the quantile estimate, particularly for quantiles corresponding to values of p approaching 0 or 1. Two types of kernels are used, an ‘interior’ kernel when p is beyond a bandwidth of the boundary (for interpolating), and a ‘boundary’ kernel when p is within a bandwidth of the boundary (for extrapolating). In this study both ‘interior’ and ‘boundary’ kernels are based on the original and modified forms of Epanechnikov kernel, respectively [Moon and Lall, 1994; Adamowski et al., 1998]. The bandwidth parameter h determines the smoothness of the estimated function; smaller bandwidths result in a rougher estimator while a larger bandwidth averages over a larger number of samples and produces a smoother estimate. Increasing the bandwidth increases the bias while reducing the variance of the estimator. Several techniques have been proposed for selecting an optimum bandwidth based on the sample values, and several are reviewed by Kim and Heo [2002]. Initial attempts to select an optimum h by least-squares cross-validation (LSCV) [Adamowski, 1985; Kim and Heo, 2002] for several of the generated PDS series produced quantile estimates that tended to degenerate to the sample quantile values for all p. As such, h was set, via trial and error, to 0.01 for all PDS series, which produced adequate flood frequency functions in all cases. The necessity of establishing the statistical significance of ∆QTp is based on the premise that KQp is itself a random variable, and its value will change if a new peak flow sample is generated. As such, any change in QTp (estimated using KQp) due to forest management activities must be assessed in light of the sampling variance associated with QTpc. The sampling variance of each QTpc (this must be done for every control scenario) was estimated using bootstrap techniques, which follows the procedure: i) ii) iii) iv) v) Each original control PDS of sample size N was resampled with replacement, generating a new PDS, also of sample size N; Quantiles were estimated from the new PDS using (1); Steps i) and ii) were repeated 1000 times; For a given return period the KQp were ranked as KQp(1) ≤ KQp(2) ··· ≤ KQp(1000), and The upper and lower 95% confidence bounds on QTpc were estimated as KQp(975) and KQp(25), respectively. 20 As sampling variability is a function of sample size, it was necessary to constrain the PDS sample size N for each management scenario to be equal to that of the given control scenario in order to properly assess the statistical significance of the ∆QTp variable. 4.0 RESULTS The efficacy of non-parametric flood frequency analysis is examined in Figure 9, which compares the fit of quantiles estimated using (1) with those estimated using the commonly used combined Poisson – Generalized Pareto (P-GP) parametric model for PDS data [i.e. Stedinger et al., 1993]. The sample points are those peak flows generated for the basin-wide control scenario (i.e. streamflow at the outlet of Carnation Creek for pre-management forest cover), where the empirical probability for the ith largest event is estimated using the distribution-free plotting position formula (i – 0.25)/(N + 0.5) [Adamowski, 1981]. From the figure it would appear that the generated flood peaks can be split into two groups. At Tp < 1.7 years peak discharge magnitude increases with increasing Tp but tends to reach an asymptotic discharge limit of approximately 23 m3/s. At Tp ≈ 1.7 years discharge jumps to 27 m3/s, and then increases with increasing Tp without any apparent asymptotic limit. The parametric P-GP model is clearly unsuitable for describing the distribution of flood peaks, missing the jump in peak discharge eat Tp ≈ 1.7 years and grossly overestimating the magnitude of large (Tp > 5 years) peak flow events. On the other hand, the non-parametric function stays true to the trend in the data, managing to capture the marked shift in flood response at Tp ≈ 1.7 years and the variability of peak discharge at Tp > 1.7 years. Beckers and Alila [2004] hypothesized that the break in slope at Tp ≈ 1.7 was possibly related to a shift in the precipitation regime. The choice of h = 0.01 for the kernel quantile estimator is validated in Figure 10. Although larger bandwidths tend to capture the general flood frequency function, including the differing response for flood peaks with Tp either greater or less than 1.7 years, they tend to over-smooth the jump in peak discharge at Tp = 1.7 years and substantially underestimate the magnitude of the largest peak discharge event. As we have no information that as yet precludes us from accepting the discontinuity at Tp = 1.7 years and the magnitude of the largest event as anything but valid (and perhaps indicative of a third shift in the flood frequency regime), we felt it necessary to adopt a bandwidth that captures these peculiarities of the data with acceptable resolution and have therefore made our choice of h = 0.01. We recognize that the choice of such a small bandwidth implies that we accept the flood frequency distribution is more or less equivalent to the empirical distribution of flood peaks. In order to carry out flood frequency analysis (parametric or non-parametric), it is important to ensure that the flood peaks are independently distributed. The flood peaks shown in Figure 10 derive from a sample of N = 138, which produces a relatively high λ of 7.56 yr-1. Despite this high arrival rate the sample autocorrelation clearly shows that peak flow events are uncorrelated (Figure 11). The Spearman rank order correlation 21 Figure 9. Plotting positions of PDS peak flow sample for basin-wide control scenario (Carnation Creek outlet with pre-management vegetation), N = 138 and λ = 7.56. Also shown are peak discharge estimated using the non-parametric (NP) kernel quantile estimator and the parametric Poisson – General Pareto (P-GP) distribution. coefficient test [Helsel and Hirsch, 2002] on the data shown in Figure 9 fails to reject the null hypothesis (α = 0.05) of independent sample peaks. Further, the Runs Above-Below Median test for general randomness [Pilon and Harvey, 1992] informs us (at α = 0.05) that the flood peaks of Figure 9 are randomly distributed. Although not shown, similar results are derived for the PDS generated by each of the scenarios of Table 6. 4.1 Cut Rate Effects The analysis approach of comparing flood quantiles estimated from the simulated flood peaks is shown graphically in Figure 12, which compares flood quantiles for scenario RC100 compared to the Control. Due to increased sampling variance for less frequent flood events the confidence interval for QTpc (difference between the 97.5th and 2.5th percentile) increases with increasing Tp (dashed lines in Figure 12) from 1.4 to 11.4 m3/s 22 Figure 10. Plotting positions of PDS peak flow sample for basin-wide control scenario (Carnation Creek outlet with pre-management vegetation), N = 138 and λ = 7.56. Also shown are peak discharge estimated using the non-parametric (NP) kernel quantile estimator with bandwidth h = 0.10, 0.05, and 0.01. for Tp ranging from 0.17 to 20 years, respectively. In this extreme harvesting case the discharge following forest harvesting is consistently higher than the Control discharge for all Tp. However, the absolute difference between QTpc and QTpm generally decreases with increasing Tp and the increase in QTpm is only significant (i.e. where QTpm exceeds upper confidence limit) for very small return periods (Tp < 1 year). The flood frequency results for all the rate of cut scenarios (numbers 1 through 11; Table 6) are summarized in Table 7. From these results two general trends are evident, the magnitude of the changes in peak discharge (∆QTp) increases with increasing rate of cut (ROC) and decreasing return period with values of ∆QTp ranging from 30 to -4%. However, as already indicated in Figure 12, these results are only statistically significant at very low return periods, and then only when ROC exceeds 50%. Statistically significant changes in QTp are in all cases restricted to Tp ≤1 year, and are in fact restricted to Tp < 0.5 years for ROC between 50 and 70%. For all cut rates ∆QTp < 10% for Tp ≥ 2 years. 23 Figure 11. Sample autocorrelation of peak flow PDS for basin-wide control scenario. Upper and lower 95% confidence limits are indicated by the dashed lines. Figure 13 compares ∆QTp for ROC from 10 to 100% for return periods of 0.17, 0.5, 2, 10, and 20 years. The relationship between ∆QTp and ROC is generally linear, although the slope of this relationship varies with return period, tending to decrease as Tp increases from 0.17 to 10 years. At Tp = 20 years, this slope increases slightly from that at Tp = 10 years. Of note is how quickly the ∆QTp versus ROC relationship ‘flattens’ out as Tp increases from 0.17 years. Also of interest is the jump in ∆QTp between ROC of 50 (RC050) and 52% (RC052Op), suggesting that at roughly equivalent ROC, the spatial distribution of cut blocks exerts some influence on flood response. However, this jump is most apparent at Tp = 0.17 years and becomes largely indistinguishable from the general ∆QTp versus ROC relationship at increasing Tp. The anomalous ∆QTp value at ROC = 30% for Tp = 2 years requires notice. Although an explanation for this negative ∆QTp value escapes us, it does serve to highlight the often unpredictable manner in which the flood regime can respond following land use changes. 4.2 Scale Effects The results of the scale effects management scenarios (numbers 12 through 30; Table 6) are summarized in Tables 8a through 9b; in all four tables results are presented by scenario in such a way that rate of cut increases but drainage area decreases when reading 24 Table 7. Control Quantile Magnitude and Relative Quantile Change for Cut Rate Management Scenarios. Tp (years) 0.17 0.50 0.75 1 2 5 10 15 20 a QTpc (m3/s) 7.6 16.2 19.1 21.2 27.0 33.5 35.7 39.0 40.8 RC009 2 3 2 1 1 0 0 0 0 RC022 5 6 5 2 1 1 1 1 1 RC030 7 7 7 2 -4 1 1 2 2 ∆QTp by Management Scenario (%) RC040 RC050 RC052Op RC060 RC070 11 14a 18a 18a 20a 10 12 11 11 12 10 11 11 12 13 4 6 7 7 8 2 2 2 4 7 1 2 2 2 2 2 1 2 1 1 2 2 3 2 3 3 3 3 3 3 Change in QTp is significant at α = 0.05. 25 RC081 24a 14 15a 10 7 2 2 3 4 RC090 27a 16a 17a 11a 8 3 2 4 5 RC100 30a 19a 19a 12a 8 3 3 4 5 Figure 12. Flood frequency comparison of management scenario RC100 to Control. the table from left to right. The control quantile discharge for each sample point is given in Table 8a and Table 9a for harvesting in C- and E-tributaries, respectively. Table 8b and Table 9b show the relative change in flood quantiles (∆QTp) for harvesting in C- and Etributaries, respectively. Results for the scale effects scenarios are analogous to those for the rate-of-cut scenarios described above. Both Table 8b and Table 9b show that, as expected, ∆QTp is highest at ROC = 100% (i.e. when discharge is sampled at the outlet of the harvest area), and decreases with decreasing ROC (i.e. as the sample point is moved further downstream); overall values range from 30 to -1%. Significant increases in QTp are restricted to ROC greater that approximately 60% and 35% for harvesting in C- and E-tributaries, respectively. It should be noted that this disparity in the ROC at which significant ∆QTp occur does not suggest substantial runoff process differences between the two sub-basins, but is more an artefact of the stepped and discontinuous manner in which drainage area increases as one moves downstream (resulting from sudden increase in drainage area at the confluence of tributaries). A more precise determination of the point where significant ∆QTp occurs when harvesting in C-tributary could not be established as the ROC drops substantially from 62% to 11% when C-tributary joins the 26 Figure 13. Relative change in flood quantiles with rate-of-cut for Tp = 0.17, 0.5, 2, 10, and 20 years. Statistically significant increases are identified with grey shading. main channel (between scenarios SC062 and SC011) (Figure 7a). The results of both Table 7 and Table 9b suggest that in general this transition likely occurs for ROC between 30 to 50%. Regardless, for harvesting in both C- and E-tributaries, significant ∆QTp is limited to Tp ≤ 0.5 years; although ∆QTp is significant for ROC = 100% and Tp = 2 years for scenario SC100E. Again, ∆QTp never exceeds 10% for Tp ≥ 2 years. As noted earlier for scenarios RC050 and RC052Op, the magnitude of ∆QTp shows some dependence upon the spatial layout of cutblocks, which is evident when one compares ∆QTp at similar ROC when harvesting either in C- or in E-tributary; for instance ∆QTp is 30 and 22% at Tp = 0.17 years for SC100C and SC100E respectively (compare Table 8b and Table 9b). However, any dependence due to spatial variability in cutblock layout seems only to manifest at very low return periods (Tp ≤ 0.75) and high ROC. When Tp ≤ 0.5 years, ∆QTp also appears sensitive to changes in drainage structure as the sampling point is moved downstream. For instance, when harvesting in E-tributary the relationship between ∆QTp and ROC is non-linear, displaying several marked jumps and changes in slope: as the sampling point moves downstream from E-tributary south and into Etributary main (which now includes all of E-tributary north’s drainage area) ROC changes from 40% to 35% and ∆QTp shows a drop; a second drop in ∆QTp occurs when the sampling point moves downstream of the confluence of C-tributary with the main 27 Table 8a. Control Quantile Magnitudes for Scale Effects Management Scenarios with Harvesting in C Tributary Tp (years) 0.17 0.50 0.75 1 2 5 10 15 20 QTp by Management Scenario (m3/s) SC009b SC009a SC010 SC011 7.6 7.5 7.3 6.7 16.2 16.0 15.5 14.3 19.1 18.7 18.2 17.3 21.2 20.7 19.8 18.4 27.0 26.7 25.7 23.8 33.5 33.0 31.7 28.8 35.6 34.6 32.8 30.4 39.0 38.1 36.3 33.4 40.8 39.9 38.1 35.0 SC062 1.0 2.1 2.4 2.7 3.6 4.2 4.8 5.1 5.3 SC071 0.9 1.9 2.2 2.4 3.2 3.8 4.2 4.4 4.5 SC100 0.6 1.3 1.5 1.7 2.1 2.5 2.8 3.0 3.1 Table 8b. Relative Quantile Change for Scale Effects Management Scenarios with Harvesting in C Tributary. Tp (years) 0.17 0.50 0.75 1 2 5 10 15 20 a ∆QTp by Management Scenario (%)a SC009b SC009a SC010 SC011 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 SC062 16a 10a 9 6 3 3 2 2 1 SC071 18a 12a 12 9 4 3 3 2 1 SC100 30a 18a 15 11 6 3 3 3 3 Change in QTp is significant at α = 0.05. Carnation Creek Channel (ROC changes from 16 to 12%) (Figure 14b). It is unclear if the relationship between ∆QTp and ROC for a given Tp potentially differs when the point of interest (sampling point) progresses down different drainage paths. Although, at first glance, Figure 14a and b may suggest that there are differences (particularly for Tp = 0.17 years), again this apparent result is more likely due to the inability to sample streamflow along the C-tributary drainage path at ROC between 11 and 62%. 4.3 Road Effects The impact of roads alone (scenario RD000) and clearcut harvesting with roads (RD052Op and RD100) is summarized in Table 10. When management activities hypothetically only involve constructing a road network, the response of ∆QTp is similar 28 Table 9a. Control Quantile Magnitudes for Scale Effects Management Scenarios with Harvesting in E Tributary. Tp (years) 0.17 0.50 0.75 1 2 5 10 15 20 SC010b 7.7 16.4 19.2 21.3 27.1 33.6 35.8 39.2 40.9 SC010a 7.6 16.3 18.8 20.8 26.7 33.1 34.8 38.3 40.0 SC011 7.3 15.5 18.2 19.3 25.7 31.8 32.8 36.4 38.1 QTp by Management Scenario (m3/s) SC012 SC016 SC018 SC022 SC028 6.78 5.6 5.1 4.6 3.9 14.4 11.7 10.6 9.4 8.0 17.4 13.9 12.6 10.5 8.8 18.4 15.1 13.8 12.0 10.1 23.8 19.6 17.6 15.4 13.1 28.8 22.8 20.2 17.7 15.3 30.4 25.5 23.2 20.4 17.1 33.5 28.0 25.2 22.2 18.3 35.0 29.1 25.9 22.6 18.5 SC035 3.4 6.6 7.6 8.4 11.3 12.9 14.4 15.0 15.1 SC040 3.1 5.9 6.9 7.8 10.4 12.0 12.8 13.4 13.6 SC085 1.4 2.6 3.0 3.5 4.6 5.3 5.8 6.3 6.5 SC100 1.2 2.2 2.6 3.0 4.0 4.6 5.1 5.7 5.8 Table 9b. Relative Quantile Change for Scale Effects Management Scenarios with Harvesting in E Tributary. Tp (years) 0.17 0.50 0.75 1 2 5 10 15 20 a SC010b 0 -1 0 0 0 0 -1 -1 0 SC010a 0 -1 0 0 1 0 -1 -1 0 SC011 0 0 0 0 0 0 0 0 0 SC012 2 1 0 1 1 -1 0 0 0 ∆QTp by Management Scenario (%) SC016 SC018 SC022 SC028 6 5 4 4 1 1 1 2 0 0 2 3 1 1 2 0 0 1 1 1 -1 -1 0 2 0 0 1 1 0 0 0 1 0 0 0 1 Change in QTp is significant at α = 0.05. 29 SC035 6a 3 5 4 2 1 0 1 1 SC040 11a 4 3 3 3 1 0 1 1 SC085 22a 10 11 11 8 1 2 2 2 SC100 22a 14a 11 9 8a 1 2 2 2 Figure 14. Relative change in flood quantiles with rate-of-cut showing scale effects when harvesting in the headwaters of a) C-tributary and b) E-tributary. Statistically significant increases are identified with grey shading. 30 Table 10. Control Quantile Magnitude and Relative Quantile Change for Road Effects Management Scenarios. Tp (years) 0.17 0.50 0.75 1 2 5 10 15 20 a QTpc (m3/s) 7.6 16.2 19.1 21.2 27.0 33.5 35.7 39.0 40.8 ∆QTp by Management Scenario (%) RD000 RD052Op RD100 11 29a 41a 3 12 25a a 6 14 19a 2 9 16a -1 3 6 -2 0 1 0 2 3 2 4 6 3 6 7a Change in QTp is significant at α = 0.05. to that which would occur in the aftermath of forest removal only. The largest ∆QTp of 11% occurs at the lowest return period (0.17 years) and decreases to a minimum (-2 %) at Tp = 5 years, where after it increases to 3% at Tp = 20 years. However, the impact of the road network on QTp drops off much more rapidly between Tp = 0.17 to 0.75 years than generally occurs with forest harvesting alone. Regardless, in no case is the increase in quantile discharge statistically significant. Depending upon return period, the effect of adding 34 km of road to Carnation Creek affects the flood frequency regime to roughly the same degree as a 10 to 40% ROC (Table 7). In the case of both forest removal and road construction, the effect of the road is to generate higher ∆QTp than occurs with forest removal alone (compare RD052Op to RC052Op and compare RD100 to RC100 in Tables 5 and 8, respectively). At Tp = 0.17 years, the effect of roads and forest removal is directly cumulative: ∆QTp is 11% and 18% for scenarios RD000 and RC052Op, respectively, which add directly to give ∆QTp of 29% for RD052Op. Similar results are obtained for a 100% cut rate for scenarios RD000, RC100, and RD100 (i.e. ∆QTp of 11, 31, and 41%, respectively). As the return period decreases, though, the combined effect of roads and forest removal is less than the direct addition of the respective ∆QTp, with the cumulative impact getting progressively smaller as Tp increases. A final note is that the cumulative effect of roads and forest removal is to slightly increase the statistical significance of flood quantile increases for a cut rate of 52%: ∆QTp becomes statistically significant for Tp ≤ 0.75 years (vice Tp < 0.17 years for RC052Op). 4.4 Scale Invariance It appears generally that the relationship between ∆QTp and ROC is scale invariant. For instance when values of ∆QTp are compared to ROC, regardless of drainage area, it is 31 apparent that ∆QTp values are generally identical for similar ROC (Tables 7, 8b, and 9b). When the results of scenarios 1 to 11 are grouped with those of 12 to 30 (Table 6) in this fashion ∆QTp exhibits a strong linear relationship with ROC. In order to examine in greater detail the nature of this relationship, linear regressions of the form ∆QTp = b0 + b1*ROC (2) have been fit for Tp = 0.17, 0.5, 0.75, 1, 2, 5, 10, 15, and 20 years; the regression results are summarized in Table 11. The diagnostic plots for Tp = 0.17 years indicate that (2) is a sound model for describing the relationship between ∆QTp and ROC: residuals do not depart from normality and the 95% prediction limit of ∆QTp is relatively small compared to the observed range of ∆QTp (Figure 15). Interestingly, all regressions in Table 11 are significant at p < 0.001 and all slopes, b1, are significant at p < 0.001. The slope of (2) decreases with increasing Tp, confirming the results of the flood frequency analysis which indicated that sensitivity of peak discharge to forest harvesting decreases with increasing Tp. The trend of the coefficient of determination R2 suggests that the relationship between ∆QTp and ROC weakens with increasing Tp. The results for return periods of 0.17, 0.5 and 20 years are shown in Figure 16, showing clearly that ∆QTp is strongly related to ROC for drainage areas ranging from 0.88 to 9.8 km2. The scatter about the regression relationship defined by (2) is attributed to variations in the spatial distribution of cut blocks, upstream drainage structure and drainage density. 5.0 DISCUSSION The results of this study are in general agreement with the results of Thomas and Megahan [1998] and Beschta et al. [2000], and show clearly that forest harvesting in rainfall-dominated watersheds has the largest impact on peak discharge at very low return period events, and that the effect of forest harvesting diminishes rapidly as return period increases. Additionally, we see that statistically significant increases in peak discharge only manifest at the lowest return periods (Tp ≤ 1 year in all cases), and then only when ROC exceeds roughly 30%. The results from this study and others provide some insight as to why discharge for small events is affected differently than discharge for large events. Increased peak flow following timber harvesting in the pacific Northwest has been attributed to more efficient hillslope runoff due to higher antecedent soil moisture (brought about by reduced evapotranspiration) and decreased rainfall interception [Rothacher, 1973; Ziemer, 1981]. Beckers and Alila [2004] determined that on an annual basis streamflow in Carnation Creek derives from matrix-only flow 67% of the time. However, during rainfall events, storm flow is generated mainly by preferential flow and matrix flow is negligible. Further, events characterized by over 50% of runoff generated by fast (flow velocity = v2) preferential flow occurs when unit area discharge exceeds roughly 2.8 mm/hour. At the basin scale this discharge is 7.6 m3/s, which has a return period Tp = 0.17 years (Table 7). It has been observed in Carnation Creek that when 32 Table 11. Parameters of Linear Regression relationship between ∆QTp and ROC for scenarios 1 through 30a Tp (years) 0.17 0.50 0.75 1 2 5 10 15 20 Intercept, b0 (%) -1.48 -0.74 -0.88 -0.76 -0.82 -0.32 -0.38 -0.37 -0.20 Slope, b1 (-) 0.30 0.18 0.18 0.12 0.08 0.03 0.03 0.04 0.04 SE (%) R2 2.01 2.05 2.20 1.01 1.48 0.74 0.51 0.65 0.96 0.96 0.88 0.87 0.94 0.76 0.65 0.79 0.79 0.65 a All regressions significant at p < 0.001; all slopes (b1) significant at p < 0.001. Figure 15. Diagnostic graphs of linear regression relationship between ∆QTp dependent upon ROC for Tp = 0.17 years. The ‘Predicted Value-Residual Plot’ and ‘X-Residual Plot’ verify that the residuals are randomly distributed, the ‘Q-Q Plot’ verifies that the residuals are normally distributed, and the ‘Regression Line and 95% Confidence and Prediction Limits’ checks the usefulness of the regression relationship. 33 Figure 16. Relationship between ∆QTp and ROC for all drainage areas showing individual values and fitted linear regressions. Statistically significant increases are identified with grey shading. vertical by-pass and lateral preferential flow is activated, subsurface flow rates and water table response are identical between wet and dry antecedent conditions [Hetherington, 1995; Fannin et al., 2000; Beckers and Alila, 2004]. On the other hand, matrix-only flow rates are very sensitive to antecedent soil moisture. As such, matrix-only dominated discharge is expected to be sensitive to forest harvesting and the transition from matrixdominated to preferential flow-dominated runoff generally accounts for the rapid decline in ∆QTp for Tp > 0.17 years. Some other mechanism is likely responsible for the significant increase in QTp observed for certain scenarios for Tp > 0.17 years. The actual volume of intercepted rainfall for any given storm will depend on storm frequency, intensity and duration during a given time of the year. When storms are generally smaller and less frequent, rainfall is intercepted more efficiently (proportionately more rainfall is intercepted during the storm and the canopy dries out between storms; when storms are more frequent and/or more intense, lower interception efficiency results [Rothacher, 1973]. Within Carnation Creek Spittlehouse [1998] observed that rainfall interception amounts to 50% of total rainfall for small storms events ( < 30 mm), but that this decreases to approximately 10% for large storm events (~300 mm). Within DHSVM, canopy parameters have been set such that maximum 34 interception capacities are in the range of 7.2 to 4.8 mm; however, simulated total interception for rainfall events can reach up to 40 mm due to evaporation from wet vegetation surfaces during the event [Beckers et al., 2004]. It is suspected that the reduction in rainfall interception accounts for the positive ∆QTp generally found for the range of 0.17 < Tp < 2 years and that the reduction in interception efficiency with increasing storm size explains the reduction in ∆QTp with increasing Tp. The processes affecting forest harvest impacts (described above) on peak discharge appear to be invariant of the size of the drainage area for areas ranging from 0.88 to 9.8 km2. As such, the relative change in peak discharge magnitude for a given Tp and ROC and be pooled and a general predictive relationship between ∆QTp and ROC (for a given Tp) established. At low return periods (Tp < 2 years) any scatter about the regression is for the most part overshadowed by the strength of the relationship between ∆QTp and ROC. This variance of the sample ∆QTp about the regression given by (2) derives from two sources. Firstly, the random layout of cutblocks for scenarios one through eleven introduces variance that derives from such issues as harvesting in areas of deep versus shallow soils, steep versus gentle slopes, high versus low drainage density, and from changes in runoff synchronization as various tributaries experience forest harvesting. Secondly, additional variance is also introduced in the scale scenarios due to differences in upstream channel morphology (which effects channel routing and streamflow synchronization) between the different drainage paths in C- and E-tributaries. From the statistical significance of the regression results presented (Table 11) one might be tempted to conclude that, contrary to the flood frequency analysis, the model given by (2) suggests that ∆QTp at all return periods is significant. Such an interpretation is incorrect; the regression results simply verify that ∆QTp increases, in a statistically significant fashion, with ROC. The flood frequency analysis tells us more specifically (and more informatively) whether or not the relative increase in QTp at a given Tp and ROC is significant (in most cases it is not). Additionally, the standard errors shown in Table 11 for Tp ≥ 2 years are of the same order of magnitude as the predicted ∆QTp, such that the predictive uncertainty of ∆QTp is almost as large as the observed range of ∆QTp (not shown). In other words, ∆QTp for Tp ≥ 2 years is so small to begin with that there is little utility in attempting to predict it with linear regression. It is apparent from Figure 16 that the ROC threshold for statistically significant increases in QTp is not clearly defined. For example, for Tp = 0.5 years, ∆QTp for ROC = 71 and 85% is statistically significant, however, ∆QTp for ROC = 80% is not statistically significant. This inconsistency is due to the use of Bootstrap technique to derive confidence limits on the control peak flow quantiles. Confidence limits are derived for each control scenario (one control scenario for all 11 cut rate scenarios; one control scenario for each of 19 scale scenarios) based on the sample peak flow quantiles, as such; the confidence limits on QTpC for each return period are also subject to sampling variability, which generates slightly different confidence limits for each control scenario for a given ROC. Generally, statistically significant increases in QTp occur for ROC > 30 and 60% for return periods of 0.17 and 0.5 years, respectively. 35 The general assumption regarding the effect of roads on peak discharge is that the road cut intercepts subsurface water moving slowly down the hillslope, where it is then collected and routed quickly as surface flow through roadside ditches [Megahan, 1987]. Water collected by a road network has the potential to be routed directly to a stream channel, either at stream crossings or by discharging at culverts which are connected to the channel system via gullies [Jones et al., 2000]. The road network thus effectively reduces the time of concentration of streamflow from the roaded portions of the watershed, which can lead to downstream increases in peak flows, depending on how well the increased flows synchronize with flows from other portions of the watershed. Order of magnitude calculations conducted by Beckers et al. [2004] show that flow velocities for matrix-only runoff range from ~10-4 to 10-3 m/s, whereas flow velocities for fast preferential flow are an order of magnitude higher at ~10-2 m/s. Based on some simplifying assumptions(hydraulic radius of 10-2 to 10-1 m and slope = 10%), flow velocities in roadside ditches for typical rain storms were estimated to be in the order of ~10-1 to 100 m/s (velocity increases with storm size). It is to be expected then that the short-circuiting of slow hillslope runoff by the road network will have the greatest impact on peak discharge response for small events dominated by matrix runoff (Tp ≤ 0.17 years). As event size increases, a greater proportion of hillslope runoff will be as fast preferential flow and the effect of the road network will diminish. If matrix only runoff is assumed for all discharge events, then one could surmise that larger events, with higher water tables, and thus a larger proportion of the water table available for interception by a road cut, would be more sensitive to the presence of a road network [Thomas and Megahan, 1998]. The corollary to this concept is that during small discharge events with a road networks built on deep soil, the low water table may pass below the road cut and no subsurface flow would be intercepted (i.e. peak discharge would be unaffected). This may be the case for peak discharge events with Tp ≤ 0.17 years (i.e. those dominated by matrix runoff). However, once preferential flow is initiated the depth of the water table tends to remain fairly constant [Fannin et al., 2000], such that the proportional depth of water table intercepted by the road network will remain the same, regardless of event magnitude. It has been proposed that the combined effect of forest harvesting and road construction will be synergistic, due to an increased volume of subsurface flow intercepted below cut areas [Jones and Grant, 1996]. Although such synergism has yet to be validated, La March and Lettenmaier [2001], using DHSVM simulations in several sub-basins of the Deschutes drainage in Washington, concluded that the effect of roads and forest removal was still additive, and that the relative increase due to roads and forest harvesting increase and decrease, respectively, with increasing return period. Further, La March and Lettenmaier [2001] determined that the effect of roads alone was to increase the mean annual flood (return period of 2.3 years) by 2 to 10%, and the 10-year flood by 3 to 12%, considerably higher than results reported for this study; however, they considered only matrix flow in their simulations. In the current study, the effects of roads and forest harvesting on the relative change in peak discharge is indeed additive for Tp ≤ 0.17 years, when discharge is dominated by matrix flow. However, when return period increases and 36 discharge becomes dominated by preferential flow, the combined effect of roads and forest harvesting is somewhat less than additive; the effect is certainly not synergistic. DHSVM does not account for all potential connectivities between the road and channel networks; specifically, the model does not account for relief culverts that drain into hillslope gullies that carry water directly to the channel. Although the true extent of this phenomenon in Carnation Creek is unknown, it is possible that the increases in drainage efficiency brought about by the road network may be understated in these simulation results. The impact of such connectivity remains to be investigated. The effects of roads on peak discharge can be quite sensitive to the road layout and design, including the spatial arrangement of the road right-of-way with respect to drainage topography and landscape properties (i.e. slope position, hillslope gradient, soil depth), and road design parameters (i.e. road slope, culvert spacing, and cutbank height) [King and Tennyson, 1984; Piehl et al., 1988; Wemple et al. 1996; Jones et al., 2000]. As such, it is not yet clear how results from this study can be extrapolated to other basins and further research in this area is merited. It must be reiterated that the effects analysed represent the maximum potential impacts that are expected to occur immediately following forest harvesting operations. Increases in peak discharge due to forest removal can be subject to reductions flowing regrowth [Thomas and Megahan, 1998], although the nature and extent of hydrologic recovery will likely depend to some degree upon the physiographic and structural properties of the recovering vegetation. As long as roads remain open, they represent a permanent change to the landscape and, as such, effects due to roads alone can be expected to be more lasting. 6.0 CONCLUSION AND IMPLICATIONS This report has presented the results of the use of a physically-based distributed parameter model to investigate the impacts for forest management on the 9.8 km2 Carnation Creek watershed located on the west coast of Vancouver Island, British Columbia. The model used accounts explicitly for the effects of vegetation upon runoff production and is also unique in its explicit portrayal of preferential subsurface runoff processes common to forested watersheds. The effects of forest harvesting upon the peak flow regime were investigated using 33 scenarios which specifically focussed on the effects of cut rate, scale, and roads by subjecting each scenario to long-term simulation based on the use of 18 years of climate data collected in the study area. Flood frequency analysis was used to quantify the relative change in peak discharge for return periods ranging from 0.17 to 20 years. The model results suggest that, for the range of drainage areas examined, the relative change in peak discharge is scale invariant and ranges from -4 to 30% for cut rates ranging from 10 to 100% of drainage area and return periods of 0.17 to 20 years. The relative change in peak discharge is strongly correlated to both rate of cut and return 37 period, increasing with increasing rate of cut and decreasing return period. Statistically significant (α = 0.05) increases in peak discharge are restricted to cut rates greater than 30% and return periods less than 2 years; statistically significant changes do not occur for return periods greater than 0.17 years until the cut rate exceeds approximately 60%. The simulation results also suggest that the relative change in peak discharge is linearly related to rate of cut and that this relationship is statistically significant at all return periods; implying that the potential impact of forest harvesting upon peak discharge can be confidently estimated for any return period for cut rates ranging from 10 to 100%. However, as the standard error of the regression eventually approaches the same magnitude as relative peak discharge change as return period increases, these relationships have limited utility for return periods greater than two years. It has been shown that the effects of roads alone is to generate relative increases in peak discharge that range from -2 to 11%, an effect that is analogous to changes attributed to forest harvesting at cut rates ranging from 10 to 40%, depending upon return period. Relative changes in peak discharge generally increase with decreasing return period; however, changes are statistically insignificant (α = 0.05) at all return periods. The combined effect of roads and forest harvesting has also been examined and results suggest that the combined effect is to increase peak discharge more so than would have occurred with either forest harvesting or roads in isolation. However, the effect of roads and forest harvesting is only additive at the lowest return period (0.17 years); the combined effect is less than additive at higher return periods. The results from the current study in no way suggest that the combined effect are synergistic (i.e. the combined relative increase in peak discharge is greater than additive). In general, the relative change in peak discharge (as result of forest harvesting) with return period can be related to the physical process governing runoff production in lowelevation maritime watersheds. The strongest impact upon peak discharge occurs at return period equal to (and by extrapolation, less than) 0.17 years. This return period marks the transition from subsurface runoff dominated by matrix flow, which is sensitive to changes in runoff efficiency resulting from changes in antecedent moisture, to subsurface runoff dominated by preferential flow, which is already a highly efficient mechanism of converting precipitation to runoff and is generally insensitive to antecedent soil moisture. At return periods greater than 0.17 years, changes in peak discharge are attributed to the reduction in rainfall interception; the relative increase in peak discharge decreases with increasing return period due to the decrease in interception efficiency with increasing storm size. The weakening response of peak flow change to roads with increasing return period is attributed to the higher subsurface flow rates associated with preferential runoff and the upper limit on the water table response to precipitation once preferential flow is activated. Numerical modeling has helped to quantify changes to the peak discharge regime associated with both forest harvesting and roads (and the combined effect of both). These results serve to highlight the strength of the numerical modelling approach over that of conducting traditional statistical analysis of paired-watershed results. In summary, the results of this study suggest that in rainfall-dominated watersheds of coastal British 38 Columbia with high precipitation, steep topography, and thin, permeable soils, the peak discharge regime is fairly robust to land use changes arising from forest management. However, readers are cautioned that these results deal merely with changes to water quantity that manifest within the channel network. We make no claim as to the impact that forest management has upon other aspect of the hydrologic regime, such as the distribution of hillslope soil moisture and water table depth, which may have indirect consequences on such processes as erosion, terrain stability, and the frequency of mass wasting. In addition, results are assessed in terms of their statistical significance. No attempt has yet been made to link the relative changes in peak discharge reported here to their potential physical significance with respect to channel ecology and morphology. Lastly, impacts reported here are those that are expected immediately following forest harvesting, and as such, represent maximum potential changes. No account has yet been made of the persistence of effects due either to forest removal or road construction. 39 7.0 REFERENCES Adamowski, K. (1981). “Plotting formula for flood frequency.” Water Resources Bulletin 17(2): 197-202. Adamowski, K. (1985). “Nonparametric kernel estimation of flood frequencies.” Water Resources Research 11: 1585-1590. Adamowski, K., G. –C. Liang, and G. G. Patry (1998). “Annual maximum and partial duration flood series analysis by parametric and non-parametric methods.” Hydrological Processes 12: 1685-1699. Beckers, J. and Alila, Y. (2004). "A model of rapid preferential hillslope runoff contributions to peak flow generation in a temperate rainforest watershed." Water Resources Research 40: doi: 10.1029/2003WR002582. Beckers, J., Alila, Y. and Hetherington, E. D. (2004). "Peak flow responses to clearcutting and roads in the maritime regions of the Pacific Northwest: a preferential hillslope runoff perspective." Water Resources Research (submitted and reviewed). Beschta, R. L., Pyles, M. R., Skaugset, A. E. and Surfleet, C. G. (2000). "Peakflow responses to forest practices in the western Cascades of Oregon, USA." Journal of Hydrology 233: 102-120. Beven, K. J. (1989). "Changing ideas in hydrology--the case of physically-based models." Journal of Hydrology 105: 157-172. Bowling, L. C. and Lettenmaier, D. P. (2001). Evaluation of the effects of forest roads on flood flows in a mountainous maritime environment. In: Land use and watersheds: human influence on hydrology and geomorphology in urban and forested areas, Water Sci. Appl. Ser. M. S. Wigmosta and S. Burgess (Eds.). American Geophysical Union. 2. Bowling, L. C., Storck, P. and Lettenmaier, D. P. (2000). "Hydrologic effects of logging in Western Washington." Water Resources Research 36: 3223-3240. Dunne, T. (2001). Introduction to section 2- problems in measuring and modeling the influence of forest management on hydromorphic and geomorphic processes. In: Land use and watersheds: human influence on hydrology and geomorphology in urban and forest areas. In: Land use and watersheds: human influence on hydrology and geomorphology in urban and forested areas, Water Sci. Appl. Ser. M. S. Wigmosta and S. Burgess (Eds.). 2: 77-83. Fannin, R. J., J. Jaakola, J. M. T. Wilkinson and E. D. Hetherington (2000). "Hydrologic response of soils to precipitation at Carnation Creek, British Columbia, Canada." Water Resources Research 36(6): 1481-1494. Grayson, R. B., Moore, I. D. and McMahon, T. A. (1992). "Physically based hydrologic modeling 2: is the concept realistic?" Water Resources Research 28: 2659-2666. Hartman, G. F. and Scrivener, J. C. (1990). "Impacts of forestry practices on a coastal stream ecosystem, Carnation Creek, British Columbia." Canadian Bulletin of Fisheries and Aquatic Sciences 223: 148p. Helsel, D. R., and R.M. Hirsch (2002). "Statistical methods in water resources." in Hydrologic Analysis and Interpretation (Chapter A3), United States Geological Survey, available at http://water.usgs.gov/pubs/twri/twri4a3/ 40 Hetherington, E.D. (1982). "A first look at logging effects on the hydrologic regime of Carnation Creek experimental watershed", in Proceedings of the Carnation Creek Workshop: A Ten-Year Review, edited by G.F. Hartman, pp. 45-62, Pac. Biol. Stn., Nanaimo, British Columbia, Canada. Hetherington, E. D. (1995). "Subsurface flow rates over bedrock on steep slopes in the Carnation Creek experimental watershed", in Mountain Hydrology: Peaks and valleys in research and application, edited by B. T. Guy and J. Barnard, Proceedings of the Canadian Water Resources Association Conference, Vancouver, BC, May 16-19. Jones, J.A., F. J. Swanson, B. C. Wemple, and K. U. Snyder (2000). "Effects of roads on hydrology, geomorphology and disturbance patches in stream networks." Conservation Biology 14(1): 76-85. Jones, J. A. and Grant, G. E. (1996). "Peak flow responses to clear-cutting and roads in small and large basins, western Cascades, Oregon." Water Resources Research 32(4): 959-974. _______ (2001). "Comment on "peak flow responses to clear-cutting and roads in small and large basins, western Cascades, Oregon: A second opinion" by R.B. Thomas and W.F. Megahan." Water Resources Research 37(1): 175-178. Kim, K. -D., and J. -H. Heo (2002). Comparative study of flood quantiles estimation by nonparametric models." Journal of Hydrology 260: 176-193. King, J. G., and L. C. Tennyson (1984). "Alteration of streamflow characteristics following road construction in north central Idaho." Water Resources Research 20(8): 1159-1163. La March, J. L., and D. P. Lettenmaier (2001). "Effects of forest roads on flood flows in the Deschutes River, Washington." Earth Surface Processes and Landforms 26: 115-134. Megahan, W. F. (1987). "Effects of forest roads on watershed function in mountainous areas." in Environmental Geotechnics and Problematic Soils, edited by A. S. Balasubramaniam et al., pp. 335-345, A.A. Balkema, Rotterdam, Netherlands. Moon, Y. -I., and U. Lall (1994). "Kernel quantile function estimator for flood frequency analysis." Water Resources Research 30(11): 3095-3103. Nash, J. E., and J. V. Sutcliffe (1970). "River flow forecasting through conceptual models. Part I - A discussion of principles." Journal of Hydrology 10: 82-290. Oswald, E.T. (1982). "Preharvest vegetation and soils of Carnation Creek watershed", in Proceedings of the Carnation Creek Workshop: A Ten-Year Review, edited by G.F. Hartman, pp. 17-35, Pac. Biol. Stn., Nanaimo, British Columbia, Canada. Piehl, B. T., R. L. Beschta, and M. R. Pyles (1988). "Ditch-relief culverts and lowvolume forest roads in the Oregon Coast Range." Northwest Science 62(3): 91-98. Pilon, P. J., and K. D. Harvey (1992). "Consolidated frequency analysis (CFA), version 3.1: Reference manual." Ecosystem Sciences and Evaluation Directorate, Environment Canada, Ottawa, Canada. Rothacher, J. (1973). "Does harvest in west slope Douglas-fir increase peak flow in small streams?", Research Paper PNW-163, 13 pp. U.S. Department of Agriculture, Forest Service, Portland, Oregon. Schnorbus, M. A., and Y. Alila (2004). "Forest harvesting impacts on the peak flow regime in the Columbia Mountains of southeastern British Columbia: An 41 investigation using long-term numerical modeling." Water Resources Research 40, W05205, doi: 10.1029/2003WR002918. Spittlehouse, D. L. (1998). "Rainfall interception in young and mature conifer forests in British Columbia", in Proceedings 23rd Conference on Agricultural and Forest Meteorology, American Meteorological Society, pp. 171-174. Stedinger, J. R., R. M. Vogel, and E. Foufoula-Georgiou (1993). "Frequency analysis of extreme events", in Handbook of Hydrology (Chapter 18), edited by C.R. Maidment, McGraw-Hill, New York, New York. Sheather, S. J., and J. S. Marron (1990). "Kernel quantile estimators." Journal of the American Statistical Association 85(410): 410-416. Shen, H. W. and P. Y. Julien (1993). "Erosion and sediment transport", in Handbook of Hydrology (Chapter 12), edited by D. R. Maidment, McGraw-Hill, New York, New York. Tague, C. and Band, L. (2001). "Simulating the impact of road construction and forest harvesting on hydrologic response." Earth Surface Processes and Landforms 26: 135-151. Thomas, R. B. and Megahan, W. F. (1998). "Peak flow response to clear-cutting and roads in small and large basins, western Cascades, Oregon: a second opinion." Water Resources Research 34(12): 3393-3400. _______ (2001). "Reply to: "Peak flow response to clear-cutting and roads in small and large basins, western Cascades, Oregon: a second opinion." Water Resources Research 37(1): 181-183. Wemple, B. C. and Jones, J. A. (2003). "Runoff production on forest roads in a steep, mountain catchment." Water Resources Research 39(8): 8-1 - 8-17. Wemple, B. C., J. A. Jones, and G. E. Grant (1996). "Channel network extension by logging roads in two basins, western Cascades, Oregon." Water Resources Bulletin 32(6): 1195-1207. Whitaker, A., Alila, Y., Beckers, J. and Toews, D. A. A. (2002). "Evaluating peak flow sensitivity to clearcutting in different elevation bands of a snowmelt-dominated mountainous catchment." Water Resources Research 38: doi: 10.1029/2001514. Wigmosta, M.S., and D. P. Lettenmaier (1999). "A comparison of simplified methods for routing topographically driven subsurface flow." Water Resource Research 35: 225-264. Wigmosta, M. S., Vail, L. W. and Lettenmaier, D. P. (1994). "A distributed hydrologyvegetation model for complex terrain." Water Resources Research 30: 1665-1679. Wigmosta, M. S., B. Nijssen, P. Storck, (2002). "The Distributed Hydrology Soil Vegetation Model", in Mathematical Models of Small Watershed Hydrology Applications, edited by V.P. Singh, D.K. Frevert, pp. 7-42, Water Resources Publications, Highlands Ranch, Colorado. Ziemer, R. R. (1981). "Storm flow response to road building and partial cutting in small streams in northern California." Water Resources Research 17(4): 907-917. Ziemer, R., Lewis, J., Rice, R. M. and Lisle, T. E. (1991). "Modelling the cumulative watershed effects of forest management strategies." Journal of Environmental Quality 20: 36-42. 42 Ziemer, R. and Lisle, T. (1998). "Hydrology", in River Ecology and Management: Lessons from the Pacific Coastal Ecoregion, edited by R. J. Naiman and R. E. Bilby. Springer-Verlag. New York. 43