QUANTIFYING THE WATER BALANCE OF TWO NORTHEASTERN BOREAL WATERSHEDS, BRITISH COLUMBIA by Sina Abadzadesahraei B.Sc., Azad University (Iran), 2006 M.Sc., Algarve University (Portugal), 2011 DISSERTATION SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY IN NATURAL RESOURCES AND ENVIRONMENTAL STUDIES The University of Northern British Columbia August 2018 © Sina Abadzadesahraei, 2018 ABSTRACT Northeastern British Columbia (BC) is undergoing steady development for oil and gas extraction, mainly due to subsurface hydraulic fracturing (fracking), which requires significant quantities of water. Thus, it is of vital importance to obtain accurate long-term water balance information in the complex wetlands of northeastern BC to assist regulators to balance multiple priorities in a way that will not compromise the long-term sustainability of water resources, while minimizing ecological impacts. At the initial phase of this study, all fluxes of the Coles Lake water balance were measured for the 2013_2014 hydrological year. The total storage change was negative (-8.3 mm), and 2013_2014 was considered a relatively dry year. This study also quantifies the water balance fluxes within two boreal watersheds, the Coles Lake and Tsea Lake watersheds, through a combination of observational data analysis and numerical modelling using the MIKE SHE hydrological model for 1979_2014. MIKE SHE model calibration was performed manually based on snowmelt, pressure head, and streamflow, using a trial-and-error parameter adjustment procedure. Similar trends were observed for the Coles Lake and Tsea Lake watersheds although average of actual evapotranspiration (AET = 472.9 mm year-1) was higher while overland flow (OL = 26.3 mm year-1) was lower at the Coles Lake watershed compared to the Tsea Lake watershed (AET= 405.5 mm year-1 and OL = 48.5 mm year-1). Sensitivity simulations with the MIKE SHE model whereby the leaf area index was modified uniformly across the Coles Lake watershed to represent fully open, mixed and closed canopies provided further insights on the role of vegetation on the water balance. Simulated AET = 515, 529, and 558 mm year-1 and OL = 59, 46, and 11 mm year-1 for open, mixed, and closed canopies, respectively. Further, the Coles Lake forcing data were applied for the Tsea I Lake watershed as a sensitivity test while other parameters remained unchanged. The variability of the vegetation canopies and land cover including wetland distribution were the main contributors for different hydrological responses in these two watersheds. Baseline information generated by this study will support the assessment of the sustainability of current strategies for freshwater extraction. II Table of Contents ABSTRACT .............................................................................................................................................I Table of Contents .....................................................................................................................................I List of Tables ......................................................................................................................................... V List of Figures ....................................................................................................................................... IX Glossary ............................................................................................................................................. XVI ACKNOWLEDGMENTS ................................................................................................................. XIX CHAPTER 1: INTRODUCTION ........................................................................................................... 1 1.1. Overview ................................................................................................................................. 1 1.2. Oil and Gas Exploration in Northeastern BC ......................................................................... 2 1.3. Climate Change in Northeastern BC....................................................................................... 3 1.4. Summary of Past Research ..................................................................................................... 6 1.4.1. Past Studies on Boreal Forests ........................................................................................ 7 1.4.2. Past Studies in Northeastern BC Boreal Forests ........................................................... 11 1.4.3. Water Tools in Northeastern BC Boreal Forests .......................................................... 13 1.5. Thesis Objectives and Outline .............................................................................................. 16 1.5.1. Research Objectives ...................................................................................................... 16 1.5.2. Research Questions ....................................................................................................... 16 1.5.3. Thesis Organization ...................................................................................................... 17 CHAPTER 2: STUDY AREA CHARACTERISTICS AND DATA COLLECTION ......................... 19 2.1. Study Area Characteristics .................................................................................................... 19 2.1.1. Watershed Characteristics ............................................................................................. 19 2.1.2. Vegetation and Land Cover .......................................................................................... 23 2.1.3. Geological Characteristics ............................................................................................ 25 2.1.4. Permafrost ..................................................................................................................... 27 2.2. Data Gathering ...................................................................................................................... 28 2.3. Fieldwork Data Collection .................................................................................................... 30 2.3.1. Weather Station Installation .......................................................................................... 32 2.3.2. Precipitation Measurements at Coles Lake ................................................................... 33 2.3.2.1. Rainfall Measurements at Coles Lake .................................................................. 34 2.3.2.2. Snowfall Measurements at Coles Lake ................................................................. 34 2.3.3. 2.3.3.1. Precipitation Measurements under Open, Mixed, and Closed Canopy Types ............. 36 Identification of Vegetation Types ....................................................................... 36 I 2.3.3.2. Rain Gauge Installation......................................................................................... 40 2.3.3.3. Snow Surveys ....................................................................................................... 43 2.3.3.4. Snow Site Selection .............................................................................................. 43 2.3.3.5. Snow Sampling ..................................................................................................... 44 2.3.3.6. Reduction of Sample Size ..................................................................................... 46 2.3.4. Piezometer Installation.................................................................................................. 47 2.3.5. Shallow Groundwater Measurements ........................................................................... 53 2.3.6. Streamflow Measurements............................................................................................ 55 2.3.6.1. Measurement Method ........................................................................................... 58 CHAPTER 3: DATA PREPARATION AND CALIBRATION RESULTS FOR THE HYDROLOGICAL MODEL................................................................................................................ 61 3.1. MIKE SHE Description ........................................................................................................ 62 3.2. Model Setup .......................................................................................................................... 63 3.3. Model Inputs ......................................................................................................................... 64 3.3.1. Model Domain .............................................................................................................. 65 3.3.2. Digital Elevation Model ................................................................................................ 66 3.3.3. Meteorological Data...................................................................................................... 67 3.4. 3.3.3.1. Coles Lake Historical Meteorological Data .......................................................... 67 3.3.3.2. Tsea Lake Historical Meteorological Data ........................................................... 69 3.3.3.3. Potential Evapotranspiration Estimation............................................................... 70 Land Surface Data................................................................................................................. 72 3.4.1. Leaf Area Index ............................................................................................................ 72 3.4.2. Maximum Root Depth .................................................................................................. 75 3.4.3. Overland Flow .............................................................................................................. 76 3.5. Unsaturated and Saturated Zone ........................................................................................... 79 3.6. Stream Network and Hydrometric Flow Data ...................................................................... 82 3.7. Model Calibration ................................................................................................................. 84 3.8. Calibration Results ................................................................................................................ 87 3.8.1. Snowmelt ...................................................................................................................... 87 3.8.2. Pressure Head in the Soil Zone ..................................................................................... 89 3.8.3. Streamflow .................................................................................................................... 94 3.9. Parameters Sensitivity and Experimental Calibration .......................................................... 98 3.9.1. Sensitivity of Snowmelt Parameters ............................................................................. 99 3.9.2. Sensitivity of Hydraulic Head Parameters .................................................................. 101 II 3.9.3. Sensitivity of Streamflow Parameters ......................................................................... 103 3.9.4. Sensitivity of Vegetation Parameters .......................................................................... 103 3.9.5. Sensitivity to the Meteorological Forcing................................................................... 103 CHAPTER 4: QUANTIFICATION OF THE COLES LAKE WATER BALANCE ........................ 105 4.1. Coles Lake Water Balance .................................................................................................. 105 4.1.1. Water Balance Components ........................................................................................ 106 4.1.1.1. Contribution of Precipitation to Coles Lake ....................................................... 106 4.1.1.2. Estimation of Lake Evaporation ......................................................................... 107 Lake Evaporation Results ................................................................................................... 109 Sublimation ......................................................................................................................... 109 4.1.1.3. Streamflow Measurements at Coles Lake........................................................... 110 Correlation Method ............................................................................................................ 111 4.1.1.4. Shallow Groundwater Results............................................................................. 113 4.1.2. Water Withdrawal (W) ............................................................................................... 117 4.1.3. Measured Water Balance Results ............................................................................... 118 4.2. Simulation of Coles Lake Hydrological Processes Using the MIKE SHE Model ............. 119 4.3. Comparison of Measured and Simulated Results ............................................................... 127 4.4. Discussion ........................................................................................................................... 129 CHAPTER 5: INFLUENCE OF FOREST CANOPY ON WATER BALANCE AT THE COLES LAKE WATERSHED ........................................................................................................................ 134 5.1. Rainfall Analysis ................................................................................................................. 135 5.2. Snowfall Analysis ............................................................................................................... 137 5.3. Influence of Forest Cover on the Coles Lake Water Balance Modelling ........................... 141 5.3.1. Discussion ................................................................................................................... 144 CHAPTER 6: THE WATER BALANCE OF TWO NORTHEASTERN BOREAL WATERSHEDS USING MIKE SHE............................................................................................................................. 146 6.1. Characteristics of the Coles Lake and Tsea Lake Watersheds ........................................... 146 6.1.1. Physiographic Characteristics ..................................................................................... 146 6.1.2. Climatic Characteristics .............................................................................................. 147 6.2. Coles Lake and Tsea Lake Water Balance Simulation ....................................................... 148 6.3. Sensitivity Experiments ...................................................................................................... 153 6.3.1. Discussion ................................................................................................................... 155 CHAPTER 7: SUMMARY AND CONCLUSIONS .......................................................................... 161 7.1. Summary ............................................................................................................................. 161 III 7.2. Conclusions ......................................................................................................................... 163 7.3. Limitations .......................................................................................................................... 164 7.3.1. Modelling Framework Limitations ............................................................................. 164 7.3.2. Data Gaps .................................................................................................................... 165 7.4. Recommendations and Future Opportunities...................................................................... 166 BIBLIOGRAPHY ............................................................................................................................... 169 APPENDIX A ..................................................................................................................................... 194 APPENDIX B ..................................................................................................................................... 195 APPENDIX C ..................................................................................................................................... 196 APPENDIX D ..................................................................................................................................... 200 APPENDIX E ..................................................................................................................................... 201 APPENDIX F ..................................................................................................................................... 203 IV List of Tables Table 2-1. UF weather station components. All elements of the weather station were sourced from Campbell Scientific.................................................................................................................. 33 Table 2-2. Percent area of open, mixed and closed vegetation cover in the Coles Lake watershed. The vegetation cover was identified based on percent crown closure of VRI. ............................... 37 Table 2-3. Percent area of open, mixed and closed vegetation cover in the Tsea Lake watershed. The vegetation cover was identified based on percent crown closure of VRI. ............................... 40 Table 2-4. The GPS points for the open, mixed and closed canopy rain gauges.................................. 42 Table 2-5. Start and snow course GPS points for open, mixed and closed vegetation sites................. 43 Table 2-6. Comparison of snow depth and CV at the three different canopies (4 February 2014). ..... 47 Table 2-7. Piezometer installation details. ............................................................................................ 50 Table 2-8. Estimates of the cross-sectional area for Transects 1, 2, and 3. .......................................... 55 Table 2-9. Global Positioning System (GPS) point of location of the inflow and outflow stations along with that on Emile Creek. ........................................................................................................ 56 Table 3-1. Daily cross-correlation results between the NARR dataset, FNA.FLA, and UF. ............... 69 Table 3-2. Daily cross-correlation results between Tsea Lake dataset and NARR and FNA. ............. 70 Table 3-3. Annual PET value for the Coles Lake and Tsea Lake watersheds. ..................................... 71 Table 3-4. Defined vegetation classes with the assigned LAI value for the Coles Lake watershed. ... 73 Table 3-5. Defined vegetation classes with the assigned LAI value for the Tsea Lake watershed. ..... 73 Table 3-6. Assigned Manning’s M estimates for controlling overland flow within the MIKE SHE model. ...................................................................................................................................... 77 Table 3-7. Initial soil parameters for each layer. .................................................................................. 81 Table 3-8. Observed data used to calibrate the model for the Coles Lake and Tsea Lake watersheds.87 Table 3-9. Final parameters used for the snowmelt module. ................................................................ 87 Table 3-10. Calibration results for the daily SWE (October 1979_September 2014). [Zero values were neglected to avoid inflating the performance of the results.] .................................................. 89 V Table 3-11. Statistics of simulation results of daily pressure head in the SZ for all piezometers (February to September 2014). ................................................................................................................. 93 Table 3-12. Final HD value for the Coles Lake and Tsea Lake cross-sections. More cross-sections are defined for the Tsea Lake as the bathymetric map was not available. .................................... 95 Table 3-13. Calibration results of the MIKE 11/SHE daily channel flow for the Inlet, Outlet, and Emile Creek at the Coles Lake watershed. ......................................................................................... 96 Table 3-14. Calibration results of the MIKE 11/SHE daily channel flow for the Tsea Lake watershed. ................................................................................................................................................. 97 Table 3-15. The range of snowmelt parameters used for the snowmelt module. ............................... 101 Table 3-16. The range of soil parameters for each layer. ................................................................... 102 Table 4-1. The t-values and P-values between Coles Lake and Emile Creek stations with the computed Pearson’s product-moment correlation coefficients. The t-value reflects the daily value of the‘t’ test statistic, n is the number of measurements, and the p-value reflects significance from 15 May to 15 July 2014. ............................................................................................................. 111 Table 4-2. Total groundwater fluxes for Transects 1, 2, and 3 (5 February 2014 to 30 September 2014). ............................................................................................................................................... 116 Table 4-3. Estimated outflow average (groundwater to the lake) per month in 2014. ....................... 117 Table 4-4. Final summary of the observed water balance components for Coles Lake ..................... 118 Table 4-5. Overland flow water balance interpretation when the flood code is used (MIKE SHE User Guide, 2017). ......................................................................................................................... 122 Table 4-6. Monthly simulated water balance of Coles Lake (October 2013_September 2014). ........ 124 Table 4-7. Comparison between measured and simulated water balance for Coles Lake (2013_2014). ............................................................................................................................................... 128 Table 4-8. Seasonal variability of simulated hydrological components of Coles Lake using MIKE SHE (2013_2014). .......................................................................................................................... 131 Table 5-1. Monthly rainfall at each site (October 2013_September 2014). ........................................ 136 VI Table 5-2. Rainfall linear model analysis results between open, mixed, and UF vs. closed canopy for June, July, August, September, and October (2013_2014). ................................................... 136 Table 5-3. Rainfall linear model analysis results between open, mixed, and closed canopies for June and July (2013_2014). ............................................................................................................ 137 Table 5-4. Comparison of snow depth and SWE at the three different canopies (4 February 2014). 138 Table 5-5. Post-hoc analysis results of snow depth between open, mixed, and closed canopies (4 February to 28 April 2014). ................................................................................................... 139 Table 5-6. ANOVA analysis results of snow depth between SR50, closed, mixed, and open canopies (4 February to 28 April 2014). ............................................................................................... 139 Table 5-7. Results of the MIKE 11/SHE channel flow for actual, open, mixed and closed vegetation cover relative to observed runoff. The results are computed based on daily values over the simulation period of 2013_2014. ............................................................................................ 142 Table 5-8. Simulated water balance for various vegetation covers at the Coles Lake watershed (2013_2014). Although the model reports small positive and negative values for canopy and snow storage change, their total sum equals zero. ................................................................. 143 Table 6-1.Comparison of precipitation at the Coles Lake and Tsea Lake watersheds (1979_2014). . 148 Table 6-2. Comparison of PET at the Coles Lake and Tea Lake watersheds (October 1979_ September 2014). ..................................................................................................................................... 148 Table 6-3. Correlation between daily precipitation, AET, overland flow, subsurface storage change, and overland storage change for the Coles Lake and Tsea Lake watersheds over the 1979 to 2014 water year period. ......................................................................................................... 150 Table 6-4. Correlation between AET, overland flow, overland storage change, and subsurface storage change for the Coles Lake and Tsea Lake watersheds using the Coles Lake forcing dataset. ............................................................................................................................................... 154 Table 6-5. Water bodies for the Coles Lake and Tsea Lake watersheds (BC Watershed Atlas, 2005). ............................................................................................................................................... 158 VII Table C.7-1. The possible range of values of aerodynamic resistance ra and surface resistance rs under unstressed conditions (actual evaporation = potential evaporation) for a number of lands uses (Hendriks, 2010). ................................................................................................................... 199 VIII List of Figures Figure 1-1. The Province of British Columbia showing regions (in green) with unconventional gas resources. The location of the Horn River Basin is shown in light brown in northeastern BC (retrieved from Natural Resources Canada http://www.nrcan.gc.ca/, 2017). ............................ 3 Figure 1-2. Observed changes in seasonal air temperatures in BC for the period of 1900_2013 (retrieved from the BC Ministry of Environment, 2016). All statistically significant trends are positive and indicate warming. Results were found to be significant at the 95% level. NS indicates that trend is not statistically significant. ........................................................................................... 4 Figure 1-3. Changes in average annual precipitation in percent per century (1900_2013) for northeastern BC, using data from Environment Canada and the BC Provincial Climate Data Set (retrieved from Environmental Reporting BC, 2015). Seasonal trends are based on averages for spring (March – May), summer (June – August), fall (September – November), and winter (December – February). Results were found to be significant at the 95% level. This means that there is a less than 5% percent probability that the results arose randomly. NS indicates that trend is not statistically significant (Environmental Reporting BC, 2015)................................ 5 Figure 1-4. Observed changes in seasonal precipitation in BC for the period of 1900_2013 (retrieved from Ministry of Environment, 2016). All statistically significant trends are positive and reveal increasing precipitation. Results were found to be significant at the 95% level. NS indicates that the trend is not statistically significant. .............................................................................. 6 Figure 1-5. Map of ecozones of Canada’s Boreal Forest (http://www.borealscience.org/boreal/, 2010). ................................................................................................................................................... 8 Figure 2-1. Location of the Coles Lake and Tsea Lake watersheds in northeastern BC within the Petitot River watershed. Inset map shows BC with the red square corresponding to the study area.. 19 Figure 2-2. Monthly mean temperature and total precipitation for 1981 to 2010 Canadian Climate Normals for the Fort Nelson Airport station (retrieved from Environment and Climate Change Canada http://climate.weather.gc.ca/climate_normals/, 2017). ............................................... 20 IX Figure 2-3. Location of the Coles Lake (top left red outline) and Tsea Lake (bottom right red outline) watersheds in northeastern BC. ............................................................................................... 21 Figure 2-4. Digital Elevation Model image of the Coles Lake watershed............................................ 22 Figure 2-5. Digital Elevation Model image of the Tsea Lake watershed. ............................................ 23 Figure 2-6. Boreal White and Black Spruce Biogeoclimatic Zone in BC https://www.sfu.ca/geog/geog351fall07/Group06/Boreal%20Black%20and%20White%20Spr uce.html, 2007). ....................................................................................................................... 24 Figure 2-7. Location of three surveyed boreholes at the Coles Lake watershed. ................................. 26 Figure 2-8. Distribution of the various types of permafrost in Canada. The Coles Lake and Tsea Lake watersheds are located in discontinuous permafrost. The location of the study area is shown in red (retrieved from Pidwirny http://www.physicalgeography.net/fundamentals/10ag.html, 2006). ....................................................................................................................................... 27 Figure 2-9. Observation sites in the Coles Lake watershed and vicinity including the AWS, rain gauges, snow survey sites, nested piezometers and hydrometric stations. The Emile Creek hydrometric station lies just outside the watershed domain used for the MIKE SHE simulations established using provincial shapefiles....................................................................................................... 29 Figure 2-10. Observation sites in the Tsea Lake watershed and vicinity including the AWS and hydrometric stations. The Tsea Lake hydrometric station lies just outside the watershed domain used for the MIKE SHE simulations established using provincial shapefiles. ........................ 29 Figure 2-11. View (looking westward) of Coles Lake with the location of the automated weather and hydrometric stations, northeastern British Columbia (Quicksilver, 2010). ............................. 31 Figure 2-12. Location of the UF automatic weather station near Coles Lake. ..................................... 32 Figure 2-13. Tipping bucket rain gauge installed next to the UF weather station (August 2013)........ 34 Figure 2-14. View of SR50 at the UF weather station (February 2014)............................................... 35 Figure 2-15.Vegetation cover map in the Coles Lake catchment area based on percent crown closure of VRI (Matsuzaki et al., 2012). .............................................................................................. 38 X Figure 2-16. View of the open, mixed and closed vegetation canopies (top to bottom). ..................... 39 Figure 2-17. Vegetation cover map in the Tsea Lake catchment area based on percent crown closure of VRI (Matsuzaki et al., 2012). .................................................................................................. 40 Figure 2-18. View of the rainfall gauge (Davis Bucket Rain Gauge: 165 mm ´ 241 mm with Odyssey ................................................................................................................................................. 41 Figure 2-19. Location of the open, mixed, and closed vegetation sites. ............................................... 42 Figure 2-20. An example of snow course pattern and direction (all points are 10 m apart). Blue points indicate where the snow depths were measured and red points indicate where snow depths and SWE measurements were both taken....................................................................................... 45 Figure 2-21. Example of graduated snow probes used for the Coles Lake snow surveys (www.blackdiamondequipment.com). ..................................................................................... 46 Figure 2-22. The cumulative coefficient of variation of snow depth for the first snow survey (4 February 2014) at closed, mixed and open sites. .................................................................................... 47 Figure 2-23. End of piezometer pipes with 40 holes and welded drive tip (February 2014). .............. 48 Figure 2-24. Locations of three transects around Coles Lake and view of piezometer installation transects layout (not to scale)................................................................................................... 49 Figure 2-25. Schematic diagram illustrating a profile view of the set of nested piezometers at Transect 1 at Coles Lake. ....................................................................................................................... 50 Figure 2-26. Using post-pounder to install the piezometers (February 2014). ..................................... 51 Figure 2-27. Manually measuring depth to the water using the Water Level Sounder (February 2014). ................................................................................................................................................. 51 Figure 2-28. View of the two damaged loggers at Transect 2 (P2-1 and P2-2) in May 2014. ............. 52 Figure 2-29. Location of the inflow and outflow hydrometric stations along with that at Emile Creek. ................................................................................................................................................. 56 Figure 2-30. View of beaver dams. The top is the outflow channel, and the bottom is the inflow channel (July 2015). .............................................................................................................................. 57 XI Figure 2-31. The pattern of separating the stream into rectangles, at each cross-section. The vertilal lines in each rectangle indicates which area was included in calculating water area in each cross-section............................................................................................................................. 59 Figure 3-1. Land surface distinctions within the Coles Lake watershed. ............................................. 74 Figure 3-2. Land surface distinctions within the Tsea Lake watershed. ............................................... 75 Figure 3-3. Manning’s M value for the Coles Lake watershed. ........................................................... 78 Figure 3-4. Manning’s M value for the Tsea Lake watershed. ............................................................. 79 Figure 3-5. An example of the soil profile representation of the unsaturated zone (UZ) and the saturated zone (SZ).................................................................................................................................. 81 Figure 3-6. Bathymetry of Coles Lake at 0.25 m intervals (Matrix Consulting Group Inc., 2010). .... 83 Figure 3-7. Comparison of daily Fort Nelson measured (dots) and Coles lake and Tsea lake simulated (solid line) SWE for the period October 1979_September 2014. The elevations of the Fort Nelson, Coles Lake, and Tsea Lake stations are 370 m, 480 m and 660 m, respectively. ...... 88 Figure 3-8. Comparison of daily measured (dots) and simulated (solid line) pressure heads at Transect 1 for the period October 2013_September 2014. ..................................................................... 90 Figure 3-9. Comparison of daily measured (dots) and simulated (solid line) pressure heads at Transect 2 for the period October 2013_September 2014. ..................................................................... 92 Figure 3-10. Comparison of daily measured (dots) and simulated (solid line) pressure heads at Transect 3 for the period October 2013_September 2014. ..................................................................... 93 Figure 3-11. Comparison of daily measured (dots) and simulated (solid line) discharge at the Inlet hydrometric station on Coles Lake (upstream) for the period October 2013_September 2014. ................................................................................................................................................. 95 Figure 3-12. Comparison of daily measured (dots) and simulated (solid line) discharge at the Outlet hydrometric station on Coles Lake (downstream) for the period October 2013_September 2014. ........................................................................................................................................ 95 XII Figure 3-13. Comparison of daily measured (dots) and simulated (solid line) discharge at the Emile Creek hydrometric station (downstream) for the period October 2011_September 2014. The 2013 discharge appears more reasonable compared to 2011, 2012, and 2014. ....................... 96 Figure 3-14. Comparison of daily measured (dots) and simulated (solid line) discharge at the Tsea Lake hydrometric station for the period (April 2010_October 2011). ............................................. 97 Figure 3-15. Description of hydrological model input variables, model components, and outputs. .... 98 Figure 3-16. Trends of degree-day melting coefficient. The degree-day coefficient default is set at 0 (mm ̊C-1 day-1); this was adjusted during the model calibration process which ranges from 0.5 to 2.0 (mm ̊C-1 day-1). .......................................................................................................... 101 Figure 4-1. Daily rainfall and snowfall (water equivalent) at the UF weather station (October 2013_September 2014). ......................................................................................................... 106 Figure 4-2. Monthly rainfall and SWE pattern at the UF weather station (October 2013_September 2014). ..................................................................................................................................... 107 Figure 4-3. Regional surface conditions including ice coverage at Coles Lake on 7 November 2014 (left) and 6 May 2014 (right). For both images, false-colour is used to restore the look of Landsat images from the Operational Land Imager on Landsat 8. ........................................ 108 Figure 4-4. Monthly evaporation at Coles Lake for October 2013 to September 2014. .................... 109 Figure 4-5. Daily discharge comparison between the inflow, outflow and Emile Creek stations between 15 May and 30 September 2014. Points are manually measured vs. the daily reconstructed discharge at each location. ..................................................................................................... 112 Figure 4-6. Runoff at the inflow, outflow, and Emile Creek during the ice-free period in 2014. ...... 113 Figure 4-7. Daily hydraulic heads in piezometers along each transect during 2014. ......................... 114 Figure 4-8. Daily vertical gradients for all three transects in 2014. Negative gradients imply that groundwater flow is upward. ................................................................................................. 114 Figure 4-9. Daily horizontal gradients for all three transects in 2014. Negative gradients imply that groundwater flows laterally from the land to the lake. .......................................................... 115 XIII Figure 4-10. Daily vertical fluxes at three transects at Coles Lake in 2014. A positive flux indicates that groundwater flow is upward. The vertical gradient record at Transect 2 was short and the overall range of values rather small and inconsistent compared to the other two transects. . 115 Figure 4-11. Daily horizontal fluxes at three transects at Coles Lake in 2014. A positive flux indicates that groundwater flow is towards the lake. ............................................................................ 116 Figure 4-12. A: The Coles Lake shapefile (polygon). B: The Coles Lake watershed including Coles Lake bathymetry and mainstream. ......................................................................................... 120 Figure 4-13. The 2013-2014 water balance terms for Coles Lake, simulated by the MIKE SHE model. ............................................................................................................................................... 125 Figure 5-1. Daily rainfall at the Coles Lake watershed (October 2013_September 2014)................. 136 Figure 5-2. Comparison of snow depth between open, mixed and closed vegetation canopies (4 February to 28 April 2014). Dots outside of the box and whiskers are outliers. Maximum and minimum values are shown at the end of each whisker. The upper and lower quartiles are the ends of the box, and the median is the horizontal line in the centre of the box. .................... 138 Figure 5-3. Snow depth frequency histograms from snow surveys conducted at the Coles Lake watershed (4 February to 28 April 2014). The survey means are represented by a vertical black line, while the fixed-point depth measurements (SR50) are marked with a vertical red line. ........................................................................................................................................ 140 Figure 5-4. Simulated runoff for actual land cover vs. open, mixed, and closed (October 2013_September 2014). ......................................................................................................... 143 Figure 6-1. Comparison of ET, overland flow, subsurface storage change and overland storage change between the Coles Lake and Tsea Lake watersheds for 1979 to 2014 water years. .............. 149 Figure 6-2. Comparison of mean monthly AET between the Coles Lake and Tsea Lake watersheds (October 1979_September 2014). Dots outside of the boxes and whiskers are outliers. Maximum and minimum values are shown at the end of each whisker. The upper and lower XIV quartiles are the ends of the box, and the median is the horizontal line in the centre of the box. ............................................................................................................................................... 150 Figure 6-3. Comparison of mean monthly precipitation between the Coles Lake and Tsea Lake watersheds (October 1979_September 2014). Dots outside of the boxes and whiskers are outliers. Maximum and minimum values are shown at the end of each whisker. The upper and lower quartiles are the ends of the box, and the median is the horizontal line in the centre of the box.................................................................................................................................... 151 Figure 6-4. Comparison of mean monthly overland flow between the Coles Lake and Tsea Lake watersheds (October 1979_September 2014). Dots outside of the boxes and whiskers are outliers. Maximum and minimum values are shown at the end of each whisker. The upper and lower quartiles are the ends of the box, and the median is the horizontal line in the centre of the box.................................................................................................................................... 152 Figure 6-5. Comparison of mean monthly snow storage between the Coles Lake and Tsea Lake watersheds (October 1979_September 2014). Dots outside of the boxes and whiskers are outliers. Maximum and minimum values are shown at the end of each whisker. The upper and lower quartiles are the ends of the box, and the median is the horizontal line in the centre of the box. The sign convention indicates a negative inflow of precipitation and a positive outflow of melting (MIKE SHE User Guide, 2017). .......................................................................... 153 Figure 6-6. Comparison of AET, overland flow, and subsurface storage change, and overland storage change between the Coles Lake and Tsea Lake watersheds forcing the Coles Lake dataset. ............................................................................................................................................... 155 XV Glossary AET Actual evapotranspiration ANOVA Analysis of Variance ASCII American Standard Code for Information Interchange AVG FLA.FNA Average of Fort Liard Airport and Fort Nelson Airport BC British Columbia BCLCS BC Land Cover Classification Scheme BEC Biogeoclimatic ecosystem classification system BWBS Boreal White and Black Spruce zone AGS Above ground surface BGS Below ground surface CV Coefficient of variation DEM Digital elevation model DHI Danish Hydrologic Institute Ea Evaporation ET Evapotranspiration FAO Food and Agriculture Organization of the United Nations FLA Fort Liard Airport FLNRO Forests, Lands and Natural Resource Operations & Rural Development FNA Fort Nelson Airport FNCN Fort Nelson Climate Normals FNFN Fort Nelson First Nation Fracking Hydraulic fracturing FWA Freshwater Atlas GIS Geographic information system GPS Global positioning system GRRG Groundwater Resources Research Group GW Groundwater GW-SW Groundwater and surface water HD Hydrodynamic Data XVI HHG Horizontal hydraulic gradient HydroSHEDS Hydrological data and maps based on SHuttle Elevation Derivatives at multiple Scales HRB Horn River Basin IBCSP International Boreal Conservation Science Panel IPCC Intergovernmental Panel on Climate Change LAI Leaf area index LiDAR Light Detection and Ranging MAP Mean annual precipitation MASL Meters above sea level (datum) MAE Mean Absolute Error ME Mean Error MODIS Moderate-resolution imaging spectroradiometer MOE Ministry of Environment NARR North American Regional Reanalysis NASA National Aeronautics and Space Administration NCEP National Centers for Environmental Prediction NEWA Network for Environment and Weather Applications NEWT NorthEast Water Tool NEXEN NEXEN Company Incorporated NGO Non-Governmental Organization NHG Northern Hydrometeorology Group NHN National Hydro Network NRCan Natural Resources Canada NSERC Natural Sciences and Engineering Research Council of Canada NTS National Topographic System OGC Oil and Gas Commission OL Overland flow P Annual Precipitation PET Potential evapotranspiration PICS Pacific Institute for Climate Solutions XVII RBA Risk-Based Water Monitoring Assessment RFC River Forecast Centre Quicksilver Quicksilver Resources Canada Incorporated RAMP Regional Aquatics Monitoring Program RH Relative Humidity RMSE Root Mean Square Error SFU Simon Fraser University SOG Snow on ground SR50 Sonic Ranging Sensor SW Surface Water SWE Snow Water Equivalent SZ Saturated zone TIFF Tag Image File Format UF UNBC & FLNRO UNBC University of Northern British Columbia US United States USDA United States Department of Agriculture USDS Ultrasonic snow depth ranging sensors USGS United States Geological Survey UZ Unsaturated zone VHG Vertical hydraulic gradient VRI Vegetation Resources Inventory WMO World Meteorological Organization WNA Western North America WP Water Portal WS Wind Speed WY Water Year XVIII ACKNOWLEDGMENTS During my doctoral research at UNBC, I obtained tremendous support from my supervisory committee, colleagues, friends, and many others. First and foremost, I would like to express my profound gratitude and deep regards to my supervisor Dr. Stephen Déry for his guidance and constant encouragement. His support and mentoring have contributed significantly to the completion of my Ph.D. research. I express my sincere gratitude to my committee members, Dr. Ellen Petticrew (UNBC), Dr. John Rex (FLNRO) and Dr. Diana Allen (SFU) for their support and contributions. I also extend my appreciation to my external examiner Dr. Tricia Stadnyk (University of Manitoba) whose comments significantly improved the quality of this dissertation. I sincerely thank Dr. John Rex whose assistance was exceptionally valuable, as it helped me to overcome many obstacles and to advance throughout this effort. Also, I am grateful to my field assistants, Ben McGrath, James Fraser, Derrick Van Tol, and Kyle Siemens for their cooperation during the fieldwork activities. In addition, I am very grateful to Dr. Stefanie LaZerte who helped me with the "R" programming and statistical package. I am fortunate to have had Michael Allchin as my colleague and friend to assist me with Arc GIS. I wish to thank my friend Quinn Bassek for her assistance with reviewing portions of this document. I thank Dr. Dan Ryan for his input to this research, and Mike D’Aloia, Eiji Matsuzaki, Richard Kabzems and Dave Maloney who supported us to initiate this research. Also, I thank the BC Ministry of Forests, Lands, Natural Resource Operations & Rural Development who provided funding for this project, and the Bulkley Valley Research Centre team who managed administrative aspects of the project and organized the yearly seminar talk in Smithers, BC. I also thank Quicksilver staff members, Phil Langille, Mick Somerwil and Bevin Boynton, for their kind support and access to their operation sites. In addition, I am grateful to Geoscience BC for choosing me as a recipient of their 2015 and 2016 scholarship programs. I thank Environment Canada’s Science Horizons program, UNBC, OGC, and NSERC for their additional financial support. Last, but not least, with XIX heartfelt gratitude, I thank my family for their constant encouragement. I especially owe much to my mother, Shole Rastegar who supported me constantly. Finally, I am thankful to my best friend, Sheila Price for her unconditional love and support. XX CHAPTER 1: INTRODUCTION 1.1. Overview In recent years, northeastern British Columbia (BC) has seen significant growth in the oil and gas sector, and this industry has been a major contributor to the region’s and the province’s economy. Oil and gas development requires significant water resources, primarily for hydraulic fracturing, and while there is some understanding of regional water resources, limited research has been undertaken to study local water balances in detail. Therefore, it is essential that industry pursue its activities with minimal detrimental impacts to the natural environment. The Fort Nelson Wetland Project was initiated by the BC Ministry of FLNRO and UNBC to study the historical water balance of two northeastern BC boreal watersheds. The overarching goal of this research was to quantify the water balance of two northeastern boreal watersheds, the Coles Lake and Tsea Lake watersheds, using a combination of observations and hydrological modelling. These watersheds were selected for two main reasons. Firstly, the Coles Lake and Tsea Lake watersheds contain lakes that are being actively used by the Quicksilver Resources Canada Incorporated (Quicksilver) and NEXEN Company Incorporated (NEXEN), respectively. Secondly, historical hydrometeorological and hydrometric data were available. This study was designed and developed in two main phases. The first phase was focused on data collection and gathering baseline information of the Coles Lake and Tsea Lake watersheds, which were used as input to the hydrological model. The second phase comprised the development of hydrological model setups of each watershed using the baseline information from the first phase to simulate the historical water balance. 1 It is noted that the short-term water balance of Tsea Lake has been investigated separately by FLNRO and published as an internal report (Matsuzaki et al., 2013). However, through this research, a water balance of Coles Lake (2013_2014) was provided to support FLNRO regarding water licensing and approvals to Quicksilver (Abadzadesahraei et al., 2017). In addition, long-term (1979_2014) water balances of the Coles Lake and Tsea Lake watersheds were quantified through this research. The results of this research can be used to assist regulators (e.g., FLNRO and Oil and Gas Commission) better understand the water fluxes within northeastern BC watersheds, and more broadly, balance multiple priorities in a way that will minimize ecological impact but not compromise the long-term sustainability of water resources. 1.2. Oil and Gas Exploration in Northeastern BC Over the past decade, there have been revolutionary advancements in the development of extraction techniques for hydrocarbon reserves, where oil and gas are embedded within subsurface formations (Canadian Water Network, 2015). These innovations, such as horizontal drilling and multi-stage high-pressure hydraulic fracking, make it feasible to access these unconventional reserves (Canadian Water Network, 2015). Within northeastern BC (Figure 1-1), the Montney Formation and Horn River Basin (HRB) have come into prominence in the past few years due to their abundant unconventional natural gas resources (Natural Resource Canada, 2017). Unconventional gas development requires large quantities of freshwater, with the most significant volumes of water used for hydraulic fracking (Khyade, 2016). The water volume requirement, per fracture stage, depends on its magnitude and completion method of fracture stage (Johnson and Johnson, 2012). Kennedy (2011) reported 2 that the total water volume per life of a gas well could vary from <1000 m3 to >70,000 m3 in northeastern BC. It is, therefore, necessary to establish a thorough understanding of water availability in this region. Figure 1-1. The Province of British Columbia showing regions (in green) with unconventional gas resources. The location of the Horn River Basin is shown in light brown in northeastern BC (retrieved from Natural Resources Canada http://www.nrcan.gc.ca/, 2017). 1.3. Climate Change in Northeastern BC In general, the province of BC has warmed an average of 1.4°C from 1900 to 2013, which is higher than the global average rate of 0.85°C (BC Ministry of Environment, 2016). Within the same time frame, northern regions of BC have warmed 1.6°C to 2.0°C, with changes mainly occurring in winter months, which equals or exceeds twice the global average (BC Ministry of Environment, 2016). Specifically, winter air temperatures in northern BC have increased by 3.0°C to 3.8°C from 1900 to 2013 (BC Ministry of Environment, 2016), and there is a province-wide warming trend in the spring and summer. From 1900 to 2013, northeastern BC warmed by 1.6°C in spring, while summer air temperatures have warmed 1.4°C to 1.6°C (Figure 1-2). The rate of increase in daily minimum air temperature is also higher than the rate of increase of daily maximum air temperature (Warren et al., 2004). 3 Figure 1-2. Observed changes in seasonal air temperatures in BC for the period of 1900_2013 (retrieved from the BC Ministry of Environment, 2016). All statistically significant trends are positive and indicate warming. Results were found to be significant at the 95% level. NS indicates that trend is not statistically significant. Average air temperatures in BC are projected to increase in all seasons by 2080 but not equally (BC Ministry of FLNRO, 2016). Warming is predicted to be higher in the winter than in the summer, particularly the winter minimum air temperature. In autumn, winter and spring, higher air temperatures are projected to result in an overall increase in precipitation, with more falling as rain instead of snow. Warmer air temperatures in spring will likely trigger an earlier start to the spring snowmelt freshet (IPCC, 2007). If there is a limited success internationally to control future emissions, BC could experience warming of 3°C to 5°C by the 2080s (BC Ministry of FLNRO, 2016). 4 The long-term annual average precipitation in BC, for the period of 1900_2013, has increased by 12% (BC Ministry of Environment, 2016). Annual average precipitation increases have ranged from 10% to 21% across the province, with an 11% average precipitation growth for northeastern BC (BC Ministry of Environment, 2016). Specifically, precipitation rose in winter, summer, and fall, but not in spring within the same time period (Figure 1-3). Overall recorded trends indicate annual precipitation has been increasing across the province (Figure 1-4). However, future projections suggest southern and central BC are to become drier in the summer whereas northern BC is projected to get wetter (BC Ministry of Environment, 2016). Figure 1-3. Changes in average annual precipitation in percent per century (1900_2013) for northeastern BC, using data from Environment Canada and the BC Provincial Climate Data Set (retrieved from Environmental Reporting BC, 2015). Seasonal trends are based on averages for spring (March – May), summer (June – August), fall (September – November), and winter (December – February). Results were found to be significant at the 95% level. This means that there is a less than 5% percent probability that the results arose randomly. NS indicates that trend is not statistically significant (Environmental Reporting BC, 2015). 5 Figure 1-4. Observed changes in seasonal precipitation in BC for the period of 1900_2013 (retrieved from Ministry of Environment, 2016). All statistically significant trends are positive and reveal increasing precipitation. Results were found to be significant at the 95% level. NS indicates that the trend is not statistically significant. 1.4. Summary of Past Research The Coles Lake and Tsea Lake watersheds are located in the northeastern boreal forest of BC, which is the part of the larger boreal forest in Canada, where due to the harsh winters and remoteness, their hydrological processes are less well-understood compared to other forested ecozones (e.g., temperate forests) (Buttle et al., 2000). In recent years, however, significant industrial activities (e.g., oil and gas, logging) have created opportunities for studies aimed at better understanding hydrological processes in boreal forest landscapes 6 (Buttle et al., 2000). The demand for water by industry makes identifying hydrological components all that much more essential. Unfortunately, few studies have been carried out to date (e.g., Ferone and Devito, 2004; Foote and Krogman, 2006; Kabzems et al., 2012). So we have limited understanding of the influential factors that govern the hydrologic functioning of wetland systems in the boreal forest (note: in this thesis, “boreal forest” refers to the area of the northern boreal forest located in BC and Alberta). The purpose of this section is to review previous studies undertaken in the boreal forest (Devito et al., 2005; Buttle et al., 2009; Devito et al., 2012). In particular, Section 1.4.1 discusses previous studies of the hydrological processes in the boreal forest, whereas Section 1.4.2 focuses on some of the studies undertaken in the boreal forest located in northeastern BC (Johnson, 2010; Johnson and Johnson, 2012; Chapman et al., 2012). Additionally, some tools have been developed to support decision-makers to better understand the hydrological processes and water licence approvals in northeastern BC; these tools are discussed in Section 1.4.3. 1.4.1. Past Studies on Boreal Forests The boreal forest is the world’s second largest area of uninterrupted forest. According to the Regional Aquatics Monitoring Program (RAMP, http://www.ramp-alberta.org), about one-third of the world’s forested land is boreal forest, encircling a significant portion of the northern hemisphere (RAMP, 2017). On a smaller scale, approximately 35% (307 million hectares) of Canada’s total land area includes boreal forest, which stretches from west to east (RAMP, 2017). At the provincial scale, BC contains approximately 6% (32 million hectares) of Canada’s boreal forest (Naturallywood, 2010). 7 Canada comprises seven main boreal forest ecozones where there are significant differences in hydrological processes due to their distinct geography and scale (Figure 1-5; http://www.borealscience.org/boreal/, 2010). Each of these boreal forest ecozones contains its own unique features and ranges from mountainous and alpine habitats (Boreal Cordillera and Taiga Cordillera), densely deciduous forests (Boreal Plains), heavily coniferous forest (Boreal Shield), to sparse and often wetland-dense forest regions (Taiga Plains, Taiga Shield and Hudson Plains) (http://www.borealscience.org/boreal/, 2010). Figure 1-5. Map of ecozones of Canada’s Boreal Forest (http://www.borealscience.org/boreal/, 2010). The Coles Lake and Tsea Lake watersheds contain relatively flat plains, the majority of which are identified as the Taiga Plains ecozone, with a small component identified as Boreal Plains ecozone (https://www.for.gov.bc.ca/, 2017). According to the International Boreal 8 Conservation Science Panel (IBCSP) (http://www.borealscience.org/boreal/), Taiga Plain has a sparse and often wetland-dense forest, whereas Boreal Plain is characterized as densely deciduous. Hence, in this section, the emphasis is on the hydrological processes of these two ecozones. Devito et al. (2012) identified that climate plays the primary role in controlling the water balance in the western boreal forest. The boreal forest in northern BC remains a wetlanddominated landscape where the evaporation from open water bodies is a significant component of the hydrologic cycle and water balances (Granger and Hedstrom, 2011). Annual potential evapotranspiration (PET) often exceeds annual precipitation (P) (Buttle et al., 2009; Johnson, 2010). The amount of available moisture is dictated by the difference between evapotranspiration (ET) and precipitation (Devito et al., 2012). As a result of greater PET relative to P, this area is in a water deficit (Devito et al., 2012) and significant interannual variability in soil water storage capacity and soil water deficits has been detected, as indicated by a 30-year average (Devito et al., 2005). Other studies (Buttle et al., 2009; Johnson, 2010; Devito et al., 2012) confirmed similar results, indicating that on average the amount of ET is slightly higher compared to P in the boreal plains ecozone. In the western boreal forest, the typical PET and P values are 520 mm year and 480 mm year , respectively, which puts forests and wetlands at risk due to -1 -1 relatively dry conditions (Devito et al., 2012). As a result, a better understanding of the longterm annual and seasonal averages of P and ET trends could provide essential information to decision-makers in regards to water storage and redistribution (Devito et al., 2012). According to the long-term seasonal trends, P and ET peak simultaneously, during the growing season and ice-free period (Devito et al., 2012). Based on previous studies (e.g., 9 Devito et al., 2012), the majority of the P falls as rain from May to September (growing season) and returns to the atmosphere through ET. Even though the boreal forest is in a net water deficit overall, a net surplus occurs in the non-growing season (freezing period), which may contribute to lessening the water deficit in the boreal forest (Devito et al., 2012). In contrast to the growing season, in the non-growing season (October-April), a much higher proportion of the precipitation (mostly snowfall) accumulates at the surface. Snowfall is temporarily stored as seasonal snowpacks, and the ET is at the lowest level as plants are dormant (Devito et al., 2012). Although snow is one of the main components of the hydrological cycle in the northeastern boreal forest and is a major source of freshwater for oil and gas exploration, there has been no recent effort to measure the Snow Water Equivalent (SWE). In addition, several other studies were conducted in boreal forests to simulate spatial and temporal patterns of ET using a distributed hydrologic model combined with remote sensing (Wigmosta et al., 1994; Chen et al., 2005; Sanford et al., 2007). During the modelling exercise, significant spatial variation in ET was discovered, despite only minor topographic variability (Buttle et al., 2009). Sanford et al.’s (2007) results indicate the importance of the hydrologic processes (ET distribution) as influenced by minor differences in topography and land cover (Buttle et al., 2009). Additionally, Sanford et al. (2007) applied a distributed hydrologic model to simulate the flow at the Batchawana River Watershed located in central Ontario, Canada, noting that the northeastern portion of its basin is occupied by boreal forest. Sanford et al. (2007) simulated the flow regime of about 100 basins of varying size in a similar biogeoclimatic region over 30 years. The purpose of that modelling exercise was to use a 30-year simulated flow record 10 to analyze natural variability in the flow regime at different spatial scales. Sanford et al. (2007) identified, regardless of scale, in times of higher precipitation, flow variability was similar across basins. Under drier conditions, the flow regime was scale-dependent, with smaller basins (<600 ha) showing an extensive range of variability that converged with increasing basin area (Sanford et al., 2007). The variability of flow regimes in smaller basins was correlated positively with the proportion of depressional wetlands in a basin (Sanford et al., 2007). 1.4.2. Past Studies in Northeastern BC Boreal Forests In the last decade, there were several stakeholder groups (including the public, aboriginal groups such as the Fort Nelson First Nation (FNFN), and Non-Governmental Organizations (NGOs)), who were concerned that water allocation in northeastern boreal watersheds for oil and gas development would negatively influence the function of wetland ecosystems. This concern was exacerbated by the fact that regional climate and glacial history of this area make surface water flows vulnerable to change (Kabzems et al., 2012). Hence, since 2010, several studies have been undertaken to address these concerns (Johnson, 2010; Johnson and Johnson, 2012; Matsuzaki et al., 2012; Kabzems et al., 2012; Chapman et al., 2012). A conceptual water balance model for the HRB near Fort Nelson was developed by Johnson (2010). Johnson suggested that generating a reliable numerical water model for northeastern BC was challenging, owing to the presence of discontinuous permafrost, extensive and patchy muskegs, low relief and, more importantly, a lack of in situ observations. Patchy muskegs— defined by predominant sphagnum moss and black spruce forests—are one of the dominant landscapes of these wetlands (Johnson, 2010). Wetland 11 characteristics vary significantly, thus estimating the wetland water balance can be complex (Johnson, 2010), with significant variability in errors (Grundling et al., 2015). The northeastern region also contains clay-rich landforms that are characterized by a very low hydraulic conductivity (Huntley et al., 2011). As a result of low hydraulic conductivity, water contributes to the dominance of near-surface lateral flow paths (Pelster et al., 2008; Johnson, 2010). The water movement near the surface can result in complex spatial patterns of water flow in a landscape with subdued topography (Kabzems et al., 2012). Thus, water discharges from these wetland systems originate from near-surface flows, rather than groundwater (Pelster et al., 2008). For much of the year, subsurface flows are controlled by the depth, thickness, and duration of ground frost and permafrost. Johnson (2010) also discovered that the distinction between fen and bog was vital since they play a fundamental role in stream discharge mechanisms. Fens are water bodies that interact laterally with flowing groundwater, whereas bogs in this region act as water storage and form regions of potential permafrost (Johnson, 2010). Fens are connected to surface water tables that result in less variable water tables and outflow rates compared to bogs (Pelster et al., 2008). Bogs are hydrologically isolated from surface and groundwater inputs (Zoltai and Pollett, 1983). Bogs are the type of wetlands with poor nutrient soil and comprise peatlands, shrubs and maybe some trees (Amsel, 2017). Peatlands which are common in the northeastern boreal forest generally have highly variable discharge. Peak flows usually occur during the spring melt period when the peatland is unable to regulate its flow due to decreased storage capacity during wet and frozen conditions (Roulet and Woo, 1986). Despite the high peak flow conditions, Bowling and Lettenmaier (2010) observed that up to 80% of snowmelt could go into storage (i.e., soil, 12 groundwater, and surface storage). The amount of snowmelt entering storage can be predicted by maximum SWE and the lake storage deficit of the previous summer (Bowling and Lettenmaier, 2010). 1.4.3. Water Tools in Northeastern BC Boreal Forests Water allocation planning efforts were made by Chapman et al. (2012) to develop a webbased hydrological model, called the NorthEast Water Tool (NEWT). One of the primary intentions of designing the NEWT was to create baseline hydrometric information in northeastern BC (Chapman et al., 2012). To fulfill this purpose, they applied available climatological information along with basic land use data to estimate the water balance (runoff = P − ET). The NEWT was constructed based mainly on the available hydrological information, hydrometric data (e.g., monthly and annual averaged runoff), and available water licence and permit records provided by the BC government or OGC. Although the NEWT is a useful tool to support decision-making for water licence approvals, it has some limitations that originate from the underlying data and information on which the tool was originally constructed (Chapman et al., 2012). Uncertainties are one of the main limitations of the NEWT, as the hydrological fluxes are estimated from empirical modelling. Research by Chapman et al. (2012) reported that the median error of the NEWT modelling was 3.7%, and only 78% of the basins used for the model calibration were within ±20% of their observed flow (Chapman et al., 2012). Another limitation observed by Chapman et al. (2012) was that the modelled hydrometric data report the 30-year average (or “normal”) runoff rather than current conditions (Chapman et al., 2012). That is, the NEWT does not capture a particular year’s water availability, despite knowing that interannual variations in hydrology are significant (Holding et al., 2015). 13 In addition, the NEWT does not take into account the effects of natural phenomena (such as climate change and beaver dams), which can result in shifts in the hydrologic cycle. Furthermore, the NEWT reports licenced total freshwater abstractions according to permit records rather than actual abstractions, which may yield unrealistic results (Holding et al., 2015). Further to these limitations, the temporal distribution of abstractions is not considered in the model simulation (Holding et al., 2015); withdrawing large quantities of water over short periods of time can increase stress on water systems (Holding et al., 2015). The BC Ministry of Environment’s Freshwater Atlas (FWA) map coverage was used as the underlying basin and watershed delineation by Chapman et al. (2012). In some areas, the FWA may not have adequately represented the basin delineation (e.g., a few hundred hectares in area) (Chapman et al., 2012). Especially in areas of low relief and muskegdominated headwaters, river flows can be more stable than in other areas, which can cause further uncertainty (Anderson et al., 2009). Lastly, the tool was limited by only representing the surface water system, and interaction of groundwater-surface water (GW-SW) was neglected. The Water Portal (WP) – an online map-based water information tool – was designed in 2014 to provide a broad range of water-related data and information for northern BC (BC Water Tool, 2016). The WP contains water quantity and quality data wherever available including both surface and groundwater. In contrast to the NEWT, the WP tool can produce various graphs and statistics for interpretation. However, both tools are limited by the availability of data — especially for groundwater — resulting in uncertainties and underrepresentation of data-sparse areas in northeastern BC (Holding et al., 2015). In fact, only 14 seven provincial observation wells were used to monitor the water level within northeastern BC, which is one of the main monitoring limitations (Holding et al., 2015). Recognizing the limited information on groundwater in northeastern BC, FLNRO and Geoscience BC, among other partners, have been carrying out several groundwater-related studies in the southern part of the region. To date, no groundwater studies have been conducted in the boreal forest area. However, Holding and Allen (2015) mapped shallow groundwater intrinsic vulnerability throughout the northeastern region. Despite the lack of available data, the map provides a baseline of information of shallow groundwater vulnerability, which refers to the physical characteristics of the aquifer system that make it more or less vulnerable to groundwater contamination. Previous studies discovered that the combination of climate change, shale gas development activities, along with physical characteristics of northeastern BC watersheds (e.g., gentle topography, extensive wetlands) make hydrologic studies in the region particularly challenging (Johnson, 2010). Although several valuable studies have been undertaken to identify some of the hydrological components in the northeastern region (Johnson, 2010; Matsuzaki et al., 2012; Kabzems et al., 2012), insufficient research has examined the integrated historical, temporal, and regional hydrological dynamics of this area. Therefore, this research examines the hydrology of the northeastern region and simulates hydrological processes and water balances to attain a better understanding of boreal wetland dynamics. 15 1.5. Thesis Objectives and Outline 1.5.1. Research Objectives The overarching objectives of this project are: • To employ observational data to estimate the water balance of Coles Lake (for hydrological year 2013–2014) and to assess the temporal relative contribution of hydrological components within a wetland dominated landscape using a distributed hydrological model (Chapter 4). • To determine the impact of forest canopy on liquid precipitation vs. solid precipitation; to apply the hydrological model to quantify the influence of canopy density on the water balance of northeastern boreal watersheds based on the amount of overland flow (Chapter 5). • To characterize the hydrological components within two northeastern boreal watersheds, Coles Lake and Tsea Lake, through a combination of observational data analysis, and numerical modelling (Chapter 6). 1.5.2. Research Questions 1) Can observations of water balance components be used to verify the performance of a distributed hydrological model (i.e., MIKE SHE) for northeastern BC boreal watersheds? 2) What is the impact of forest canopy on rainfall and snowfall accumulation in northeastern BC boreal watersheds? How does forest canopy impact the water balance of boreal watersheds, as simulated by the MIKE SHE model? 16 3) What are the main contributors of physiographic (e.g., vegetation cover, leaf area index) and climatic (e.g., temperature) contributors to different hydrological responses and trends? 1.5.3. Thesis Organization This thesis is organized into seven chapters that follow sequentially: Chapter 1: This chapter provides the objectives of this research along with an overview of past studies. Chapter 2: This chapter characterizes the Coles Lake and Tsea Lake watersheds including their physiography, climatology, vegetation, surface, and subsurface geologic formations. It also discusses the fieldwork procedures and data collection involved in this research. Chapter 3: This chapter outlines the preparation of datasets for the MIKE SHE hydrological model that was used to simulate the regional water balance of the two northeastern boreal watersheds. The model setup and calibration are also presented. Chapter 4: This chapter quantifies the Coles Lake water balance based on the conceptual and numerical model. The measured and simulated water balances results are validated. A water year from 1 October 2013 to 30 September 2014 is used as the temporal framework within which to estimate the balance, as this period begins and ends when both discharge and storage were at their minimum levels. Chapter 5: This chapter examines the influence of forest canopy changes on the precipitation and overland flow in Coles Lake. To this end, the role of 17 vegetation cover at Coles Lake was examined by adjusting the leaf area index (LAI) value within MIKE SHE to reflect the open, mixed, and closed canopy. Chapter 6: This chapter describes the hydrology of two wetland-dominated basins in northeastern BC, with a focus on the simulation of its water balance using a hydrological model. Here, the development and results of hydrological models of the Coles Lake and Tsea Lake watersheds are discussed in detail. Chapter 7: This chapter summarizes the significant findings of each chapter as well as outlines the outcomes of this research, followed by recommendations (e.g., timing of water extraction), and potential future research for northeastern BC (e.g., data gaps, limitations, etc.). 18 CHAPTER 2: STUDY AREA CHARACTERISTICS AND DATA COLLECTION 2.1. Study Area Characteristics 2.1.1. Watershed Characteristics The Coles Lake and Tsea Lake watersheds are located in northeastern BC (National Topographic System (NTS) map area of 094O) between Fort Nelson and the Northwest Territories border (Figure 2-1). The Coles Lake watershed is located approximately 140 km from Fort Nelson and has a drainage area of 227.9 km2. Coles Lake is part of the Peace River Land District and is situated southeast of the confluence of Emile Creek and the Petitot River (Figure 2-1). Figure 2-1. Location of the Coles Lake and Tsea Lake watersheds in northeastern BC within the Petitot River watershed. Inset map shows BC with the red square corresponding to the study area. 19 The Coles Lake watershed is located approximately 30 km northwest of the Tsea Lake watershed with an area of 84.6 km2, which is located about 70 km northeast of Fort Nelson (Figure 2-1). Both watersheds are located in the moist and cool boreal white and black spruce subzone (Delong et al., 2011). They are characterized by wetland complexes of discontinuous permafrost, fens, bogs, swamps and marshes on a glaciolacustrine plain with extensive organic deposits (Kabzems et al., 2012; Johnson, 2010; Golder Associates, 2010; Huntley et al., 2011). Annual air temperature at the Fort Nelson Airport (FNA) weather station (Station ID: 1192940; Latitude: 58°50'11" N; Longitude: 122°35'50" W; Elevation: 381.9 m) averages 1.7°C and ranges from -2.4°C to 3.6°C (Environment and Climate Change Canada, 2017). The region is relatively dry; with mean annual precipitation of 451.7 mm (Environment and Climate Change Canada, 2017). The annual, monthly maximum snow depth, which occurs in January, February, or March, averages 51.0 cm. Permanent snow cover lasts from early November until early or mid-May, depending on the year. Overall, due to long, cold winters, and short growing seasons, forest productivity is low (Figure 2-2). Figure 2-2. Monthly mean temperature and total precipitation for 1981 to 2010 Canadian Climate Normals for the Fort Nelson Airport station (retrieved from Environment and Climate Change Canada http://climate.weather.gc.ca/climate_normals/, 2017). 20 Coles Lake is a small and shallow water body with a maximum depth of 2.2 m. The southern and western portions of the Coles Lake watershed drain to the west and north through Emile Creek and flow into the Petitot River, whereas the northern and eastern sides of the Coles Lake watershed drain to the northeast through Fortune Creek and flow also into the Petitot River (Figures 2-1 and 2-3). Figure 2-3. Location of the Coles Lake (top left red outline) and Tsea Lake (bottom right red outline) watersheds in northeastern BC. 21 Elevation in the Coles Lake watershed ranges from 485 m to 700 m with an average of 524 m above sea level. The Coles Lake watershed contains an elevated central highland, which acts as a drainage divide (Figure 2-4). Coles Lake is the main body of water situated within the Coles Lake watershed (59°46'57''N and 122°36'27''W) with an area of 1.7 km2 (based on the GIS data file from the Government of BC). Figure 2-4. Digital Elevation Model image of the Coles Lake watershed. 22 The Tsea Lake watershed is located about 70 km northeast of Fort Nelson with an area of 84.6 km2 (Figure 2-1). Elevation ranges from 660 m to 771 m with an average of 689.3 m above sea level. In contrast to the Coles Lake watershed, the Tsea Lake watershed contains several larger lakes including Tsea Lake (59°23'9.28"N and 122° 2'54.16"W) and drains northeast also to the Petitot River (Figure 2-1 and Figure 2-5). Both watersheds have very low relief (Holland, 1976). Figure 2-5. Digital Elevation Model image of the Tsea Lake watershed. The streams within both watersheds flow through extensive low-lying muskeg-type terrain, with frequent meanders and beaver dams. The area is predominantly composed of fens that are affected by mineral soil waters (ground and surface) where groundwater inflows maintain relatively high mineral content (Matrix Solutions Incorporated, 2012). 2.1.2. Vegetation and Land Cover According to the BC biogeoclimatic ecosystem classification (BEC) system, the Coles Lake and Tsea Lake watersheds are located within the Boreal White and Black Spruce zone 23 (BWBS). Most of the BWBS zone forms part of the vast Boreal Forest that stretches in the northeastern corner of BC (Figure 2-6). This zone contains 10% of BC’s total land area, which makes it the largest biogeoclimatic zone in BC (DeLong et al., 1991). The Boreal White and Black Spruce zone covers most of the Plateau in BC’s northeast (DeLong, 1991) as well as lower elevations in the central north. The primary climatic features of this zone are long and cold winters, while the summer growing seasons are warm but short. As a result of the cold climate, trees and most other plants grow slowly, and dead plants decompose slowly. Figure 2-6. Boreal White and Black Spruce Biogeoclimatic Zone in BC https://www.sfu.ca/geog/geog351fall07/Group06/Boreal%20Black%20and%20White%20Spruce.html, 2007). 24 In addition, this zone combines two main ecosystems, that of the upland forests and the muskegs. Mixed stands of trembling aspen and white spruce can be found in upland forest zones, while mixed stands of lodgepole pine and black spruce dominate in this region (DeLong, 1991). In areas with a wetter climate, denser communities of black spruce and moss abound. In contrast, muskeg remains more extensive in the northeastern lowland areas of poor drainage. Muskeg is a peatland combination of bogs and nutrient-poor marshes that cover extensive parts of northeastern BC. Along with muskeg, stunted black spruce and tamarack trees, brown mosses, and boreal grasslands dominate this area (McDowell, 1996). 2.1.3. Geological Characteristics Subdued topography with extensive bog and fen deposits characterize the landscapes of northeastern BC. A thick blanket of clay-rich till deposited by the Laurentide Ice Sheet overlies the bedrock throughout much of the area (Smith et al., 2005). Silt and clay-rich morainal deposits are typically present at the surface in well-drained forested areas but are invariably overlain by organic materials and/or glaciolacustrine sediments in poorly drained areas (Smith et al., 2005). Morainal landforms include low-relief plains, crevasse-squeezed ridges, flutes, and rolling, recessional, and interlobate moraines (Simandl et al., 2005). The clay-rich landforms that dominate landscapes near Fort Nelson (Huntley et al., 2011) have low hydraulic conductivity and contribute to the dominance of near-surface lateral flow paths for water (Pelster et al., 2008; Johnson, 2010). The Fort Nelson area is underlain by gently dipping marine shales and siltstones of the Lower Cretaceous Buckinghorse Formation, Fort St. John Group (Thompson, 1977). This formation has a minor sandstone component (Johnson et al., 2004). 25 To further investigate the Coles Lake watershed, Harder Associates Engineering Consulting Incorporated was retained by Quicksilver to define the subsurface sediment strata. Quicksilver was one of the leading oil and gas companies in northeastern BC, and their operation was located within the Coles Lake watershed. Three boreholes, with depths ranging from 30 m to 120 m were drilled at the Coles Lake watershed with sediment samples collected during drilling (Figure 2-7). Figure 2-7. Location of three surveyed boreholes at the Coles Lake watershed. The subsurface conditions are relatively uniform at all three borehole locations. The borehole stratigraphy consists of, from the surface downward, peat, clay, and clay till. The clay till extends to the maximum exploration depth (30.5 m), which varies slightly amongst the three boreholes. The sediment is relatively uniform with slight variations. Peat is present at the surface of all three boreholes and extends to an average depth of 0.32 m. Clay exists below the peat in all three boreholes at an average depth of 0.32 m below the surface. Clay till appears at an average depth of 1.83 m and extends to the maximum depth in all three boreholes. 26 2.1.4. Permafrost The Coles Lake and Tsea Lake watersheds are located within the sporadic discontinuous permafrost zone (Figure 2-8). Permafrost occurs when the ground remains below 0°C for a minimum of two years consecutively (Miceli and Lewkowicz, 2011). Northeastern BC is part of the sub-arctic region located in the discontinuous permafrost zone and could be strongly affected by the future climatic change (Bonnaventure et al., 2012). Warmer air temperatures can lead to the degradation of discontinuous permafrost, which is relatively warm and potentially thin (Romanovsky et al., 2010; Smith et al., 2010; Lewkowicz et al., 2011). Figure 2-8. Distribution of the various types of permafrost in Canada. The Coles Lake and Tsea Lake watersheds are located in discontinuous permafrost. The location of the study area is shown in red (retrieved from Pidwirny http://www.physicalgeography.net/fundamentals/10ag.html, 2006). 27 Overall, the distribution and characteristics of permafrost are controlled by climate, primarily by air temperature and snow. However, at the local scale, surface and subsurface factors impact the permafrost layer (Smith and Riseborough, 1996). Although microclimate, terrain and soil properties exercise much control over permafrost conditions in the discontinuous zone, changes in permafrost distribution and thickness differ on a regional scale (Kwong and Gan, 1994). The existence of permafrost may have a substantial impact on the local water cyc1e but be of less importance on a regional scale. In watersheds where permafrost is prevalent, it acts as an aquiclude, preventing recharge of the groundwater by precipitation and consequently increasing ET (Hartman and Carlson, 1973). The net result is a concentration of water near the surface and decrease of groundwater storage beneath the permafrost (Johnston and Brown, 1964). 2.2. Data Gathering A field campaign was undertaken to collect high-frequency hydroclimatological data to address the research questions of this study. Details of the frequency of sampling for each parameter are reported in this Chapter. Fieldwork was conducted from June 2013 to September 2014 at Coles Lake, to measure the hydrological components of this watershed in detail (Figure 2-9). Challenges related to the field efforts included the remoteness of the basin and difficulty accessing the area, as well as frequent adverse weather conditions. Although access to the Tsea Lake site was not permitted by Nexen, meteorological datasets (2011_2012) and hydrometric data (2010_2011) were provided by the BC Ministry of FLNRO. According to the report by Matrix Solutions Incorporated (2012), both the automatic weather station (AWS) and hydrometric instrumentation met the BC Ministry of Environment (MOE) exposure standards (Figure 2-10). 28 Figure 2-9. Observation sites in the Coles Lake watershed and vicinity including the AWS, rain gauges, snow survey sites, nested piezometers and hydrometric stations. The Emile Creek hydrometric station lies just outside the watershed domain used for the MIKE SHE simulations established using provincial shapefiles. Figure 2-10. Observation sites in the Tsea Lake watershed and vicinity including the AWS and hydrometric stations. The Tsea Lake hydrometric station lies just outside the watershed domain used for the MIKE SHE simulations established using provincial shapefiles. 29 2.3. Fieldwork Data Collection This section focuses on in situ data collection and the associated fieldwork procedures. Results of the fieldwork analysis and numerical models are discussed in subsequent chapters. A general description of each field procedure and its purpose is outlined below: • A Campbell Scientific, Inc. weather station (controlled by a CR1000 data logger) was installed in a large clearing area ~400 m from Coles Lake (Figures 2-11 and 2-12). Its purpose was to record the hydrometeorological conditions that contribute to the water balance of Coles Lake and its watershed (see Appendix A for the meteorological station equipment configuration and accuracy). Rainfall, snow depth, wind speed (WS) and direction, air and soil temperatures, atmospheric pressure, relative humidity (RH), and solar radiation were recorded at the automatic weather station (see Table 2-1 for details). • In addition to the rain gauge at the weather station, three Davis Instruments tippingbucket rain gauges were deployed, each under a different vegetation canopy (i.e., open, mixed, and closed). This instrumentation was selected to determine the amount of rainfall contributed to the watershed for the various vegetation canopies. • Three snow survey sites were established, one at each of the Davis Instrument rain gauge locations. These sites were selected to capture the contribution of snowfall to the watershed for the different vegetation canopies. • Nine nested piezometers equipped with the Odyssey™ capacitance water level recorders (Dataflow Systems Limited) were installed along three transects at the shoreline of Coles Lake. This was to examine the variability of seasonal shallow groundwater flows and to quantify the contribution of shallow groundwater to Coles Lake (see Appendix B for details). 30 • Two Onset Computer Corporation hydrometric stations (controlled by HOBO® data loggers) were established, one at the inflow and another at the outflow of Coles Lake (Figure 2-11). The hydrometric stations allow measurement of the streamflow discharge at the inflow and the outflow stations. • Quicksilver established a hydrometric station on Emile Creek, 4 km downstream of the lake outflow (Figure 2-11), providing additional streamflow data on this system. Figure 2-11. View (looking westward) of Coles Lake with the location of the automated weather and hydrometric stations, northeastern British Columbia (Quicksilver, 2010). • Amounts of lake water extracted for Quicksilver’s exploration operations were provided by the industrial partner, providing the total amount of water withdrawn from the lake. Each of the above procedures is explained in more detail in the following sections. 31 2.3.1. Weather Station Installation To collect meteorological data at Coles Lake, the UNBC & FLNRO (UF) automatic weather station was installed in June 2013. The automatic weather station was located in a large clearing ~500 m2 in area (Figure 2-12). There was no vegetation near the station where the ground was bulldozed flat approximately three years prior to deployment. The surrounding area of the clearing was a mixed forest composed of white and black spruce with alders measuring nearly 10 m to 15 m in height. Figure 2-12. Location of the UF automatic weather station near Coles Lake. The weather station contained three temperature sensors, one installed at a depth of 50 cm, one at 18 cm underground and another 15 cm above ground. Other instrumentation was also deployed at the weather station to measure WS and direction, air temperature and RH, atmospheric pressure, and snow depth (Table 2-1). The precipitation gauge was set up on a plywood stand at the height of 101 cm above ground. In addition to the UF station, the Tsea Lake station was installed by Nexen to record the local climate data in its watershed. Thus, mean air temperature, rainfall, RH, and WS were 32 recorded at 60 min intervals within the Tsea Lake watershed. Unfortunately, access to the Tsea Lake station was not permitted to UF personnel by Nexen. However, according to the report by Matrix Solutions Incorporated (2012), the weather stations were installed far from any existing road or construct, and instrumentation met the BC Ministry of Environment (MOE) exposure standards (MOE, 2011). Table 2-1. UF weather station components. All elements of the weather station were sourced from Campbell Scientific. Device (Model) Measurement and function Instrument accuracy Instrument height Data Logger (CR1000) Control System ~120 cm Barometer (61205v) T, RH Probe (HMP45C) Pressure Temperature and Relative Humidity (RH) Soil/Air/Snow Temperature (109B) Temperature via Thermistor ±0.12% of reading + offset, -25 ̊C to 50 ̊C ±0.5 hPa -40 ̊C to +60 ̊C, ±2% (0% to 90% RH) ±3% (90% to 100% RH) ±0.1 ̊C (over 0 ̊C to 70 ̊C range increasing to ±0.5 ̊C at -50 ̊C) Sonic Ranger (SR50-45) Snow Depth ±1 cm or 0.4% of distance to target (whichever is greatest) RM Young Wind Monitor (05103-10) Wind Speed (WS) and Direction Wind Speed: ±0.3 m s-1 Wind Direction: ±3° 293 cm above ground Pyranometer Solar Radiation 2.87% 300 cm above ground Tipping Bucket (TE525USW) Rainfall 0.245 mm 101 cm above ground ~120 cm ~ 270 cm above ground 18 cm and 50 cm underground 15 cm above ground 259 cm above ground 2.3.2. Precipitation Measurements at Coles Lake The following section describes the data collection process used to measure the amount of precipitation over Coles Lake itself. The precipitation measurement procedure for various canopies is also described in Section 2.3.2. 33 2.3.2.1. Rainfall Measurements at Coles Lake The tipping bucket rain gauge was installed less than 400 m from Coles Lake in an open area (Figure 2-13). The tipping bucket gauge was only used during the ice-free period, as it would not operate under frozen conditions unless insulated and heated. The rainfall data were recorded at 15 min intervals from 1 October 2013 to 30 September 2014; however, the 15 min values were summed to provide daily totals for analysis purposes. Figure 2-13. Tipping bucket rain gauge installed next to the UF weather station (August 2013). 2.3.2.2. Snowfall Measurements at Coles Lake To measure the amount of snowfall over Coles Lake, an ultrasonic snow depth ranging sensor (USDS) was used. The USDS was employed for taking automated snowfall measurements rather than labour intensive and costly manual observations. Therefore, a Campbell Scientific Sonic Ranging Sensor (SR50) was used to measure the amount of snowfall at Coles Lake (Figure 2-14). Snowfall refers to the amount of new snow that had 34 fallen during the 24-hour measurement period as recorded by the SR50, while the snow depth refers to the total amount of snow on the ground and includes both old and new snow. The sensor is based on a 50 kHz (ultrasonic) electrostatic transducer and can record the amount of snow depth by sending a high-frequency pulse of sound toward the ground (Campbell Scientific, 2013). An SR50 was installed at the UF weather station. One sample was taken per minute, then averaged over 15 minute intervals. The time it takes for the pulse to return to the receiver (after reflecting off a targeted surface) divided by two gives the distance to the target. The deeper the snow on the ground beneath the sensor, the less time it requires for the sound to return to the receiver. Subtracting the reordered number from a fixed reference point creates a Snow on Ground (SOG) measurement. In theory, the change in SOG levels over time provides a snowfall measurement (Fischer and Durocher, 2007). Figure 2-14. View of SR50 at the UF weather station (February 2014). Although the SR50 is a reliable instrument to record snow depth, the ability to derive a snowfall statistic may be affected by some meteorological phenomena (e.g., melting, settling, and wind redistribution of snow). Additionally, some issues have also been identified with the ultrasonic pulse being attenuated due to intense snowfall or snow surface 35 structure (e.g., low snow density), resulting in less reliable return signals (Brazenec and Doesken, 2005). To investigate the reliability of data collected by the USDS, Brazenec et al. (2006) tested the reliability of the USDS by installing 15 of these devices at various locations during 2004_2005 to collect snowfall and snow depth information in USA. Data compared well against 6-hourly and daily manual observations In some cases, however, the USDSs produced occasional spurious data associated with adverse weather conditions including heavy snow, blowing snow, and extreme cold (temperatures <-20°C) (Brazenec et al., 2006). 2.3.3. Precipitation Measurements under Open, Mixed, and Closed Canopy Types The water cycle and vegetation are closely interconnected (Gerten et al., 2004). Water balance is one of the key determinants of the vegetation distribution and its productivity (Gerten et al., 2004). Concomitantly, the distribution and productivity of plant communities is of essential importance for ET and can also regulate runoff generation (e.g., Dunn and Mackay, 1995). In other words, forest vegetation through the interception of rainfall directly affects the availability of water for streamflow, and can consequently alter components of the water balance at the watershed scale. Hence, it is important to capture the variability of precipitation within different vegetation canopies. 2.3.3.1. Identification of Vegetation Types The purpose of identifying vegetation types (open, mixed and closed vegetation) was to assess their variability and influences on interception and ET processes. Therefore, based on crown closure of the most recent vegetation resource inventory (VRI) shapefiles (Ministry of Sustainable Resource Management, 2002), percentage areas of open, mixed and closed 36 vegetation cover within the Coles Lake and Tsea Lake watersheds were identified. The vegetation type information was employed as one of the main inputs for the MIKE SHE hydrological model as discussed in the following chapter. Additionally, this information was utilized in the design of the snow survey sites and rainfall monitoring networks within the Coles Lake watershed. After identifying vegetation types based on VRI crown closure, the percent cover of vegetation types was also measured at two different spatial scales to assess its variability: a LiDAR strip digital image and the entire watershed area. The percent cover of the three vegetation types showed agreement between two spatial scales with the percent differences being less than 10% (Ministry of Sustainable Resource Management, 2002). Finally, the spatial VRI coverage of vegetation types based on crown closure was draped over a vegetation height digital elevation model (DEM; i.e., 3D view of LiDAR) to visually assess the spatial agreement for open, mixed and closed vegetation. The separation of the vegetation types into three cover classes was confirmed visually and the watershed cover distribution being mixed (51.4%) followed by open (43.3%) and closed (5.3%) (Table 2-2 and Figure 2-15). Photos of open, mixed and closed vegetation canopies are found in Figure 2-16. Table 2-2. Percent area of open, mixed and closed vegetation cover in the Coles Lake watershed. The vegetation cover was identified based on percent crown closure of VRI. Canopy Crown closure (%) Area (km2) Percent area (%) Open ≤ 25 104 43.3 Mixed > 25 - < 61 123 51.4 Closed ≥61 13 5.3 Total - 240 100 For the Tsea Lake watershed, 17% qualifies as open, 56% as mixed, and 27% as closed vegetation cover (Table 2-3 and Figure 2-17) (Matsuzaki et al., 2012). Although the 37 vegetation type was identified for both watersheds, the following section focuses on the fieldwork procedures for the Coles Lake Basin, but the complete results of the analyses are presented and discussed in Chapter 5. Figure 2-15.Vegetation cover map in the Coles Lake catchment area based on percent crown closure of VRI (Matsuzaki et al., 2012). 38 Figure 2-16. View of the open, mixed and closed vegetation canopies (top to bottom). 39 Table 2-3. Percent area of open, mixed and closed vegetation cover in the Tsea Lake watershed. The vegetation cover was identified based on percent crown closure of VRI. Canopy Open Mixed Closed Crown closure (%) ≤ 25 > 25 - < 61 ≥61 Area (km2) 14 48 23 Percent area (%) 17 56 27 Total - 85 100 Figure 2-17. Vegetation cover map in the Tsea Lake catchment area based on percent crown closure of VRI (Matsuzaki et al., 2012). 2.3.3.2. Rain Gauge Installation Three additional rain gauges were deployed on tripod stands (101 cm) in various canopies – in addition to the rainfall gauge at the UF weather station – to record rainfall during the summer season (note: the precipitation gauges collect only rainfall and not snowfall) (Figure 2-18). See Appendix A for the accuracy of the rainfall gauge. The locations of the open, mixed and closed canopy stations are displayed in Figure 2-19, and their GPS points are shown in Table 2-4. 40 Figure 2-18. View of the rainfall gauge (Davis Bucket Rain Gauge: 165 mm ´ 241 mm with Odyssey Logger) at the closed canopy site in the Coles Lake watershed (August 2013). Odyssey data loggers were selected and calibrated to record the rainfall data in each of the three rain gauges (one in the open canopy, one in the mixed canopy and one in the closed canopy). All sites were selected based on the following criteria: 1) VRI map and each station represented homogeneous vegetative coverage for each type of canopy and was located within the Coles Lake watershed boundary; and 2) Sites were no closer than 100 m to any existing road or construct. 41 Figure 2-19. Location of the open, mixed, and closed vegetation sites. Table 2-4. The GPS points for the open, mixed and closed canopy rain gauges. Forest Canopy Open Precipitation Gauge GPS Points (Latitude, Longitude) 59°47'6.15" N 122°29'21.83" W Mixed 59° 47' 34.07" N 122°36' 43.43" W 488 m Closed 59° 47' 53.75" N 122° 34' 12.47" W 494 m 42 Elevation 488 m 2.3.3.3. Snow Surveys In addition to automated snowfall measurements collected at the UF meteorological station, a snow sampling procedure was initiated at the open, mixed, and closed canopy sites to measure the SWE and snow accumulation in the Coles Lake watershed. 2.3.3.4. Snow Site Selection Sampling sites were established with the intention of being visited multiple times during winter and spring. The snow courses were chosen, and the GPS points were marked for each station (Table 2-5). Table 2-5. Start and snow course GPS points for open, mixed and closed vegetation sites. Forest Canopy Start points (Latitude, Longitude) Elevation Open 59°47'6.15" N 122°29'21.83" W 488 m Mixed 59° 47' 34.07" N 122°36' 43.43" W 488 m Closed 59° 47' 53.75" N 122° 34' 12.47" W 494 m Snow Course GPS Points (Latitude, Longitude) 59°47'4.14"N 59°47'4.07"N 122°29'34.44"W 122°29'6.54"W 59°47'16.87"N 59°47'17.03"N 59°47'28.87"N 59°47'28.40"N 59°48'1.44"N 59°48'0.71"N 59°47'41.40"N 59°48'13.00"N 59°48'13.85"N 59°47'41.91"N 122°29'6.42"W 122°29'34.49"W 122°36'43.45"W 122°35'38.73"W 122°35'35.57"W 122°36'41.15"W 122°34'55.61"W 122°34'54.60"W 122°33'49.88"W 122°33'49.59"W The sites were selected below the timberline but high enough to be free from premature melting, protected from drifting, and shaded to some extent from the direct rays of the sun. All three sites were within a distance of 2 km of each other, at ~500 m in elevation, and were nearly level. The criteria for establishing the three snow survey courses are listed below: 43 1) Each station represented homogeneous vegetative coverage for each different type of canopy and was located within the Coles Lake watershed boundary; 2) Sites were no closer than 100 m to any existing road or infrastructure; 3) Sites were accessible by foot, skis, and snowshoes; and 4) Site location was selected far away from the influence of any Quicksilver structure since the permanence of the site and the continuity of records were key factors in site selection (Gray and Male, 1981). 2.3.3.5. Snow Sampling Snow surveys were designed with regular intervals at designated stations. Snow courses were established for the open, mixed and closed canopies that included the line of fixed points throughout the winter. Snow courses were permanently marked where the snow surveys were conducted. Periodic sampling was undertaken throughout the winter of 2014; measurements began in early February (4th) and were then conducted bi-weekly throughout the snowmelt period until the end of April. Since it was important to compare the change of snow accumulation in the different canopies, snow measurements needed to be obtained on the same day for all three of the different sites. In the first round of snow surveying, long courses were considered and multiple samples of snow depth were taken – 64 points (8 × 8) at 10 m intervals (Figure 2-20) – although the number of snow depth samples was reduced for the second and subsequent rounds once the data were verified and the spatial variability was fully captured. Gray and Male (1981) used a similar approach in their research study. The reduction of sampling procedure is explained in more detail in Section 2.3.3.6. 44 Figure 2-20. An example of snow course pattern and direction (all points are 10 m apart). Blue points indicate where the snow depths were measured and red points indicate where snow depths and SWE measurements were both taken. A grid around the snow survey site constructed for measurement and recording of snow depth (in centimeters) was performed using graduated snow probes (±1 cm accuracy) (Figure 2-21). The first and last points of a course were permanently marked, and the other points were flagged. The standard snow sampling procedure was used to measure the SWE and density. To measure the SWE, a snow sampling tube was used (±0.5 cm accuracy), and snow density was computed by dividing the SWE by the snow depth. To obtain more accurate results, for the first round of sampling, eight SWE samples were taken from each site and one sample at each line (Figure 2-20). However, after the number of samples was reduced to 36, just six samples were taken from each site and one sample at each line. 45 Figure 2-21. Example of graduated snow probes used for the Coles Lake snow surveys (www.blackdiamondequipment.com). 2.3.3.6. Reduction of Sample Size After snow courses were conducted and the spatial variability was fully captured for three snow sites, the coefficient of variation (CV) of snow depth for each site was assessed (Table 2-6). CV is a useful statistical approach for determining how many samples are required for each site to fully capture the spatial variability in snow depth. The CV is calculated iteratively along the snow course and stabilizes as the number of measurements increases (Neumann et al., 2006). The results as shown in Figure 2-22 indicate that, after 30 46 samples, the CV of snow depth becomes less variable, such that ~30 samples should be sufficient to capture the spatial variability at each site. Table 2-6. Comparison of snow depth and CV at the three different canopies (4 February 2014). Canopy Open Closed Mixed Nº of Samples 64 64 64 Max Snow Depth (cm) 52 67 85 Min Snow Depth (cm) 10 26 25 Average Snow Depth (cm) 30 46 49 Coefficient of Variation (CV) (%) 32 16 19 Figure 2-22. The cumulative coefficient of variation of snow depth for the first snow survey (4 February 2014) at closed, mixed and open sites. Therefore, for the remainder of the measurements, sampling decreased from 64 to 36 (6 × 6) snow depth measurements. Indeed, knowing the number of required samples was especially important since manual snow surveying was labour-intensive and timeconsuming in such a remote and frigid area. 2.3.4. Piezometer Installation As part of this study, piezometers were installed to determine hydraulic head at specific points around Coles Lake. Piezometers are perforated pipes installed vertically on the ground 47 to intercept the groundwater passively and monitor the fluctuation of the water table. Nine piezometers (2 m in length) were constructed out of 1½ inch galvanized steel. The piezometer intake was made by drilling forty ¼ inch diameter holes into the lower 23 cm of each piezometer. The holes were spaced 3 cm apart (Figure 2-23). Figure 2-23. End of piezometer pipes with 40 holes and welded drive tip (February 2014). The hydraulic head can vary both vertically and horizontally. To measure vertical gradients, nested piezometers were used. Nested piezometers are groups of piezometers (here three) of different lengths that are installed side by side to measure hydraulic heads at various depths. Vertical gradients were estimated from the hydraulic head measurements made at the different depths, while horizontal hydraulic gradients were estimated from piezometers of similar depth that were placed laterally apart. In this study, nested piezometers were installed around the lake along three transects, with three piezometers per transect (Figure 2-24). Specifically, Transect 1 was placed in the middle of the lake shoreline, close to the UF weather station, in an area with a 5% grade. Transect 2 was installed on the east side of the lake, and Transect 3 on the west side of the 48 lake. Transects 2 and 3 were approximately 2.1 km apart, both 80 m away from streams in areas of a 6% grade. Figure 2-24. Locations of three transects around Coles Lake and view of piezometer installation transects layout (not to scale). Each transect contained two nested piezometers while the third piezometer was installed 5 m away. Points (A), (B), and (C), respectively represent deep (1.75 m), moderate (1.45 m), and shallow (1.15 m) piezometers. For each transect, two piezometers were installed 30 cm from each other (at two different depths). One piezometer was pounded into the ground to a depth of 115 cm and the other to a depth of 145 cm (Figure 2-25 and Table 2-7). The piezometer nest on each transect was deployed on the shore and in the shallow part of the lake (at 2 m from the lake). The third piezometer was installed at a depth of 175 cm, 5 m upward gradient from the nested piezometers. All depth measurements were made from the ground surface to the middle of the screen. The term “screen” refers to the perforated end of the pipe, and no extra casing was used at the bottom of the piezometers. An example of piezometer placement along transect one is shown in Figure 2-25. 49 Figure 2-25. Schematic diagram illustrating a profile view of the set of nested piezometers at Transect 1 at Coles Lake. Table 2-7. Piezometer installation details. Piezometer Transect Depth of Piezometer (cm bgs) 1 P1-1 (115 cm) P1-2 (145 cm) P1-3 (175 cm) P2-1 (115 cm) Transect 2 P2-2 (145 cm) P2-3 (175 cm) P3-1 (115 cm) Transect 3 P3-2 (145 cm) P3-3 (175 cm) 1 cm bgs – centimeters below ground surface Transect 1 50 The piezometers were installed during the winter (4 to 7 February 2014) when the ground was frozen and covered with ice. Therefore, an auger was used to break through the ice, and then a post-pounder was used to drive the piezometer into the ground (Figure 2-26). Figure 2-26. Using post-pounder to install the piezometers (February 2014). All piezometers were levelled during their installation to ensure that they were vertical. Since the water table had already passed its lowest level (lowest point of the recession), it was essential that the piezometers be installed deep enough to intersect the water table. Therefore, after all of the piezometers were installed, two days were allowed to pass for water to collect and stabilize within the piezometers. After two days, a water interface probe e-type was used to measure the depth to the water. Results indicated that the piezometers had, in fact, reached the water table (Figure 2-27). Figure 2-27. Manually measuring depth to the water using the Water Level Sounder (February 2014). 51 Following installation, each piezometer was equipped with an Odyssey capacitance logger to record water levels using 2 m sensor lengths. The logger has a resolution of ~0.8 mm and a memory capacity of 64 kilobytes. Water levels were recorded at one-hour intervals; this was considered sufficient to detect changes in water level while minimizing the amount of stored data. Prior to installing the loggers inside the piezometers, each logger was calibrated. In addition, data logger water levels were verified against manual measurements of water level, and similar values were obtained. Ultimately, the date, time, and water level data were downloaded manually from each unit to a computer and sensors were cleaned after each download. In three cases (P1-1, P2-2, and P2-3), after the piezometers were pounded into the ground, mud entered the lowest 15 cm of the pipe. Thus, to ensure that the capacitor remained vertical and that it did not bend in the tube, the piezometers were extended. In addition, the elevations of the top of the piezometer casings were surveyed, and other relevant measurements (e.g., piezometer stickup) were made (see Appendix C). Although considerable effort was made to collect continuous data, three piezometers were bent as a result of ice (Figure 2-28). Specifically, two loggers were damaged at Transect 2 (P2-1 and P2-2) and one at Transect 3 (P3-3). The damage was discovered in late May 2014. The damaged piezometers were removed, straightened and driven back into the ground. Figure 2-28. View of the two damaged loggers at Transect 2 (P2-1 and P2-2) in May 2014. 52 2.3.5. Shallow Groundwater Measurements This section describes the methods used to calculate the hydraulic gradients (vertical and horizontal) as well as the interaction of groundwater with Coles Lake. For analysis, hourly data for each piezometer were converted to daily averages for ease of interpretation. The daily hydraulic head (h) for each piezometer was then calculated using Equation (2.1). ℎ" = elev. of the top of casing ˗ depth to water (SA ,SB ,SC ) (2.1) Once the hydraulic head was computed for each piezometer, vertical and horizontal hydraulic gradients were calculated using Equations (2.2) and (2.3), respectively: '( ') '( '. = = (+ – (- (2.2) )+ – )(- – (/ (2.3) .- – ./ where: dh = hydraulic head difference dz = vertical separation between mid-screens of (B) and (C) dx = horizontal separation between mid-screens of (A) and (B) The vertical gradient was computed between the moderate (B) and the shallow gradient (C) for all transects using vertical separation between the intakes (0.3 m). The horizontal gradient was computed between (A) to (B) for all transects using the horizontal separation between them (5 m spacing). Although (A) and (B) were not installed at the same depth, given that there is a slope near the lake, the bottom elevation of (A) is close to the same bottom elevation of (B). Note that the elevations of (A) and (B) in Transect 3 are 0.5 m different, but this should not affect the horizontal gradient calculation. Finally, the components of the groundwater flux (q) in the vertical and horizontal directions were computed using Darcy’s Law (Equations (2.4) and (2.5)). Darcy’s Law is frequently used for quantifying flow between groundwater and surface water, in particular 53 on the scale of an entire lake (Rosenberry et al., 2008). Thus, vertical and horizontal fluxes were computed as the product of the hydraulic gradient and the hydraulic conductivity: qx = - K '( (2.4) '. qz = - K '( (2.5) ') where K= hydraulic conductivity (m s-1) and q = flux (m s-1). The hydraulic conductivity is a proportionality coefficient that describes the ease with which a fluid flows through a porous medium, per unit of flow rate (m s-1) (Post et al., 2007). In this study, the hydraulic conductivity (K= 2.4 × 10-10 m s-1) was determined by Quicksilver from slug tests conducted in the three shallow boreholes in the Coles Lake watershed. To provide a better indication of the range of hydraulic conductivity at a lake or watershed scale, slug tests should be conducted in multiple wells. To determine the average annual exchange of groundwater with the lake, the average annual fluxes in the x and z directions were calculated (this approach ignores the seasonal variability). The average annual components of the flux were then used to calculate the magnitude of the net flux according to Equation (2.6): (2.6) q = 0. 1 + 0) 1 To compute a volumetric inflow/outflow rate, the lake shoreline was divided into segments, based on the location of each transect. Thus, computed fluxes were applied to portions of the lake circumference (7660.53 m). The shoreline was divided into three segments by identifying the mid-points along the lake-shore between each transect and then adding the segment lengths on the left and right of each transect. Table 2-8 shows the shares of the total perimeter for each transect. 54 Table 2-8. Estimates of the cross-sectional area for Transects 1, 2, and 3. Transects Length Between each Transect Total Length of Shoreline Segment Apportioned to each Transect (L) Effective Thickness of Aquifer (b) Area (m2) (L ´ b) From Transect 1 to Transect 2 From Transect 1 to Transect 3 From Transect 2 to Transect 3 2214.73 (m) 2530.85 (m) 2917.96 (m) 2372.79 (m) 2564.84 (m) 2722.90 (m) 1.83 (m) 1.83 (m) 1.83 (m) 4342.21 4693.66 4982.91 The cross-sectional area A of the shoreline segment represents a vertical plane at the shoreline through which water passes to either enter or leave the surface-water body. The cross-sectional area was computed using Equation (2.7): A=L × b (2.7) where A=cross-sectional area (m2), L=shoreline length (m), and b=effective thickness of the aquifer (m). The maximum length of the installed piezometers was 1.75 m, and the thick clay till was encountered at an average depth of 1.83 m. Therefore, the effective thickness of the aquifer was taken as 1.83 m. Using Darcy’s Equation (2.8), the net average annual flow of water through the cross-sectional area associated with each segment is: Q = q ×A (2.8) where Q is water flow (m3 s-1), and q remains the net flux [from the magnitude of the net flux obtained from Equation (2.6)]. 2.3.6. Streamflow Measurements There were limited efforts in measuring streamflow in northeastern BC watersheds including at Coles Lake. Hence, the following section outlines the procedure for flow measurements, which are referred to as the streamflow or discharge. In 2012, Quicksilver established a hydrometric station to measure discharge at Emile Creek, located 4 km 55 downstream from Coles Lake. Given the distance of the hydrometric station from Coles Lake, the recorded data were not representative of the discharge of Coles Lake, although the collected data were used for water balance analysis and hydrological model calibration. Prior to the streamflow measurements, it was essential to establish the network of inflows and outflows to/from the lake. Hence, during fieldwork, two channels were discovered, one on the east shore and one at the west shore of the lake. Therefore, two hydrometric (stream gauging) stations were established, one at the inflow, and another one at the outflow from the lake (Table 2-9 and Figure 2-29). Table 2-9. Global Positioning System (GPS) point of location of the inflow and outflow stations along with that on Emile Creek. Stations Latitude Longitude Elevation (m) Inflow 59°46'59.41"N 122°35'1.98"W 474 Outflow 59°47'0.21"N 122°38'19.71"W 473 Emile Creek 59°48'51.45"N 122°40'39.78"W 463 Figure 2-29. Location of the inflow and outflow hydrometric stations along with that at Emile Creek. 56 The data collection period began in June 2014 and ended in September 2014 because of very low flows (nearly zero). One measurement was taken for each month; however, stream measurements were impacted by unforeseen circumstances after July. Indeed, several beaver dams were built in the inflow/outflow channel (Figure 2-30). Figure 2-30. View of beaver dams. The top is the outflow channel, and the bottom is the inflow channel (July 2015). As a consequence, it was observed that both the inflow and the outflow were almost entirely blocked after July, and no discharge measurement could be made as a result of the low flows. To compute the amount of discharge, the streamflow data from Emile Creek were compared with two Coles Lake hydrometric stations. Although it was not possible to employ the stage-discharge rating curve to estimate the streamflow due to beaver dams, the depth of the cross-section was measured continually, which was later used in building a surface hydrological model (MIKE 11). 57 2.3.6.1. Measurement Method Flow discharge was measured at Coles Lake by establishing two cross-sections, one at the inflow station, and one at the outflow station. The width of a cross-section was 8 m and 20 m at the inflow and outflow sites, respectively. At each site, the stream was relatively straight and uniform for a sufficiently long distance to provide uniform flow through the data collection area. The stream bed was stable and free of large rocks, weeds, and protruding obstructions. After the cross-sections had been established, the velocity profiles were measured by a Swoffer meter Model 2100 at each section (see Appendix A for the equipment configuration and accuracy). The water velocity varies from surface to stream bed. Harrelson et al. (1994) recommended the division of each cross-section into at least 20 segments to accurately obtain flow velocity. However, more measurements were required for broad and complex sites (Harrelson et al., 1994). Therefore, to achieve an accurate measurement of the total discharge through the cross-section, it was first divided into 25 segments. Measurements of width, depth, and velocity were obtained for each section. The required steps were designed according to the BC hydrometric standards and explained below (Resources Information Standards Committee, 2009): • The width of the inflow/outflow was measured at each cross section. • Depth was measured manually at 25 locations on the cross sections, which provides a division of the stream into rectangles (Figure 2-31). A staff gauge was used to get a manual water level reading at a hydrometric station. • Each subsection was at least 0.1 m wide, and for accurate measurement, each subsection should provide 10% or less of total discharge. 58 Figure 2-31. The pattern of separating the stream into rectangles, at each cross-section. The vertilal lines in each rectangle indicates which area was included in calculating water area in each crosssection. • Water velocity (Vi) in each rectangle (i) was estimated using a Swoffer meter (Model 2100), with an accuracy of 0.1 m s-1. The Model 2100 was designed to be easily calibrated in a swimming pool or a quiet water body. The Swoffer meter was calibrated in a stagnant section of Coles Lake. Complete details of the calibration steps are described in the Model 2100 manual. • One or two velocity measurements at each subsection were taken. For water depths that could not be safely waded, velocity measurements were conducted from a canoe. Velocity was measured at 60% of the water depth (d) for water depths up to and including 1 m, and at 20% and 80% of the water depth for depths >1 m. The average values recorded at the 20%, and 80% points were taken as representing the mean velocity for the water column in that cross section. 59 • The area (Ai) in each rectangle was calculated. Once the area was estimated at each section, discharge in each rectangle was estimated = (Area of water in crosssection) × (Water velocity) or (qi = Ai × Vi). • Total stream discharge was computed as the sum of discharges in all sub-section rectangles (Q = Σqi). • In addition to total discharge measurement at the inflow and outflow channels, transducers were installed at each hydrometric station to capture hourly water depths over time. • Further, one transducer was mounted on the shoreline in a nearby tree to measure the barometric pressure every hour. This transducer was used for post-processing to adjust the water level measurements for changes in atmospheric pressure over time. 60 CHAPTER 3: DATA PREPARATION AND CALIBRATION RESULTS FOR THE HYDROLOGICAL MODEL Overall, there are several types of hydrological models available; however, distributed and lumped models are most commonly used in the hydrological community. A spatially lumped model produces low-resolution results using basin-averaged input data, whereas a spatially distributed model provides finer details at higher resolution (Johnson, 2010). Johnson (2010) recommended employing a distributed hydrological model for northeastern BC in areas of complex wetlands. A distributed model can account for heterogeneity of vegetation, soils and land-use characteristics in a watershed (Johnson, 2010). Given the increasing demands for surface water in the complex wetlands of northeastern BC, obtaining a long-term water balance is a critical priority. Hence, the MIKE SHE hydrological model was employed to identify the interactions among the atmospheric, surface, and subsurface components. MIKE SHE was selected for this study because it resolves the principal processes in the hydrologic cycle including ET, overland flow, unsaturated flow, groundwater flow, and channel flow as well as their interactions (Sandu and Virsta, 2015). The water balance tool within the MIKE SHE model was used to estimate the historical water balance of two northeastern BC watersheds (Coles Lake and Tsea Lake) over 35 water years (October 1979_September 2014). Furthermore, the MIKE SHE model was applied to examine the effect of changes in vegetation canopies on surface runoff and the water balance of the Coles Lake watershed. This chapter discusses the data required for the application of the hydrological model. In addition, the detailed calibration process is described at the end of this chapter (Sections 3.7 and 3.8). 61 3.1. MIKE SHE Description The MIKE SHE modelling system is a physically-based, distributed, and integrated hydrological modelling system (Yan and Zhang, 2001). MIKE SHE was produced jointly by the Danish Hydraulic Institute (DHI), the British Institute of Hydrology and the French consulting firm SOGREAH (Abbott et al., 1986). The framework of the MIKE SHE modules consists of interception and evaporation, overland flow, unsaturated zone flow, and saturated flow processes (Foster, 2014). Additionally, river networks, water bodies, and other channel parameters are operated in MIKE 11, which is directly coupled to the MIKE SHE model. A brief description of the various modules is discussed below (Jaber and Shukla, 2012; Foster, 2014): • The actual ET (AET), including the interception and evaporation module, was computed with a user-defined PET, using the Kristensen and Jensen (1975) model. This model requires vegetation dependent parameters such as Leaf Area Index (LAI), vertical root characteristics, and an interception parameter (Foster, 2014). • The hydrological model only calculates the unsaturated flow vertically in 1-D, using a soil moisture retention curve, and the saturated hydraulic conductivity for each soil classification. Within MIKE SHE, the Richards’ equation, gravity flow, and 2-layer unsaturated zone (UZ) options are available to solve for water movement from the UZ to the saturated zone (SZ) or vice versa. Richards’ equation was used in this study since it provides more realistic results compared to other options. In addition, the MIKE SHE model requires PET values as input and can only be selected if the Richards’ equation is used in the UZ. 62 • The overland flow component simulates runoff when rainfall exceeds the infiltration capacity of the soil, when groundwater discharges to the surface, or when stream banks flood (Foster, 2014). Within MIKE SHE, flow direction and runoff rate are controlled by topography and the Manning’s M coefficient (reciprocal of Manning’s n), respectively (Foster, 2014). • The SZ flow component in MIKE SHE is three-dimensional using Darcy’s equation. To simulate realistic results and to produce real-world conditions within the MIKE SHE model, boundary conditions (e.g., fixed head, zero and specified fluxes) can be modified to control the groundwater flow (Foster, 2014). Subsurface conditions are simulated as geological layers and lenses, with assigned hydraulic properties (e.g., horizontal and vertical hydraulic conductivity). • MIKE 11 simulates the surface runoff including the routing of rivers and lakes. Within MIKE 11, the rivers module comprises four main components: the river network, river cross-sections, boundary conditions, and hydrodynamic parameters. MIKE 11 uses the one-dimensional St. Venant equation to solve for channel flow (Thompson et al., 2004). • Once MIKE SHE and MIKE 11 were completed, the coupling between them was made via river links (H-points). 3.2. Model Setup Application of the hydrological model to northeastern BC was established independently for the Coles Lake and Tsea Lake watersheds. Each simulation contains two main phases. Phase one consists of constructing the MIKE SHE module, incorporating ET, vegetation cover, overland flow, and unsaturated and saturated flow components. Phase two consists of 63 the specifics of the lakes and river systems, including the morphology of the river (spatial location, and vertical extent using cross-sections), boundary conditions, and channel bed resistance numbers using the MIKE 11 module. Once MIKE SHE and MIKE 11 were set up for the two watersheds, they were coupled using river links (H-points). This task was accomplished by generating discharge (Q)/stage levels (h) in relation to the cross-section, which is specified in MIKE 11 (MIKE SHE user manual, 2007). That is, MIKE 11 uses the cross-section specifications (i.e., as a minimum upstream and downstream) to compute flow and water levels (MIKE SHE user manual, 2007). The simulation period spanned from 1 October 1979 to 30 September 2014 inclusively, with output time set to 24 hours. Additionally, the vegetation cover within the Coles Lake watershed was modified three times to examine the influence of open, mixed, and closed canopies on the water balance and runoff. During these simulations, only vegetation cover was altered while other parameters remained the same. 3.3. Model Inputs During the hydrological model setup process, an extensive number of parameters must be specified to describe the physical characteristics of the Coles Lake and Tsea Lake catchment areas. MIKE SHE requires many parameters since it operates based on idealized mathematical equations that represent the real phenomenon (Devia et al., 2015). These parameters comprise meteorological/hydrological data, PET, soil moisture content, initial water depth, Leaf Area Index (LAI), maximum root depth, topography, topology, and dimensions of the river network (Abbott et al., 1986). 64 Below is a list the required inputs for both MIKE SHE and MIKE 11 in detail; their specifics are described in the following sections. MIKE SHE A. Model domain B. Digital Elevation Model (DEM) C. Meteorological data (1979_2014) a) Precipitation b) Air temperature c) Potential Evapotranspiration (i.e. Reference Evapotranspiration) D. Land surface data a) LAI b) Maximum root depth c) Overland flow (e.g., land covers, vegetation distribution, and Manning value) E. SZ and UZ flow characteristics MIKE 11 A. Stream Network and Hydrodynamic Data (HD) a) River Network b) Cross-sections c) Boundary data d) Hydrodynamic parameters (e.g., water level and discharge) 3.3.1. Model Domain The model domain was defined for the Coles Lake and Tsea Lake watersheds using the watershed boundary map. The Coles Lake model grid size was defined as 100 m by 100 m, 65 creating a rectangular grid of 185 (x) by 255 (y) points (n=47175). The Tsea Lake model grid size was also defined at the same spatial resolution with a grid of 185 (x) by 105 (y) points (n=19425). Initially lower spatial resolutions of 250 m were considered for these watersheds; however, model results were unsatisfactory. Therefore, higher resolutions of 100 m were applied to improve the accuracy of results. However, using higher resolution caused an increase in computational time. A fine mesh grid was able to better capture the flow interaction within the unsaturated and saturated flow regimes. Additionally, using the higher resolution grid produced a higher correlation between measured and simulated data. Of note, within the MIKE SHE model, the river network was interpolated to the edges of each rectangular grid with MIKE 11. Therefore, the exchange occurs at the edges between grids cells, and more refined grid results represent more accurate outcomes (Graham and Butts, 2006). 3.3.2. Digital Elevation Model The topographic map was assigned to the model domain using the Hydrological data and maps based on SHuttle Elevation Derivatives at multiple Scales (HydroSHEDS) with 30 m horizontal resolution (World Wildlife Fund, 2009). After DEMs were created for both the Coles Lake and Tsea Lake watersheds, sinks and peaks were filled using ArcGIS to optimize them. Sinks and peaks are common errors because of the resolution of the data or rounding of elevations to the most proximate integer value (ArcGIS, 2008). Notably, in low-relief areas where total discharge is low, it is essential to fill sinks and round off peaks to avoid undefined flow directions, to fill sinks to ensure proper delineation of basins, and to define streams. If the sinks are not filled, a derived drainage network may be discontinuous. The 66 boundary conditions assigned to the model consist of a zero-flux lateral boundary, representative of the topographical divides of the watershed with adjoining areas. 3.3.3. Meteorological Data Meteorological datasets (e.g., precipitation, air temperature, and PET) form an essential input to the MIKE SHE model. Therefore, it was important to carefully select a suitable meteorological dataset to best represent the historical climatology of these northeastern BC watersheds. Hence, various datasets were examined, analyzed and compared to the recorded local climate data at the UF and Tsea Lake stations. Data were recorded at the UF and Tsea Lake stations from 1 October 2013 to 30 September 2014, and 1 October 2011 to 8 August 2012, respectively. Long-term historical meteorological data were acquired from the North American Regional Reanalysis (NARR) dataset, produced by the National Centers for Environmental Prediction (NCEP) (Mesinger et al., 2006). The NARR dataset was provided in a grid format of a 32 km by 32 km area. Two separate gridded daily meteorological datasets were extracted for each watershed. Each dataset was re-gridded to 10 km by 10 km by the Northern Hydrometeorology Group (NHG) using the linear interpolation technique. Upon forcing the NARR dataset into the model, the resolution was unified by MIKE SHE. Complete details are discussed in the following section. For analysis purposes, the water year (1 October -30 September) was used as it carries a better seasonal representation for northern environments compared to the calendar year (Greenland, 1994). 3.3.3.1. Coles Lake Historical Meteorological Data Initially, local climatological data from UF were collected for one hydrological year. The collected data were then compared to the Fort Nelson Airport (FNA), and the Fort Liard 67 Airport (FLA) stations. Furthermore, the UF data were compared to the average of the FNA.FLA station, as the UF station is located midway between the FNA and FLA stations. Once the data analyses were performed, the highest correlation, lowest mean absolute error (MAE) and root mean square error (RMSE) were identified between UF and the average of FNA.FLA. The initial UF historical climate data (October 1979_September 2014), including maximum, minimum, and mean air temperature, as well as total rain, snow, and precipitation, were then constructed using an average of the FNA.FLA stations. In addition to these parameters, daily relative humidity (RH) and wind speed (WS) data were required to estimate the total amount of PET. Although both parameters were available at the FNA and FLA stations, the disadvantage was a gap within the historical dataset, especially FLA. The possibility of using high-resolution gridded data was then examined. The NARR dataset was used in conjunction with the northeastern BC meteorological dataset. For that reason, the available parameters from the UF stations along with the constructed historical data were compared with the NARR dataset. The NARR dataset correlates with the daily UF and the constructed dataset from the average of FNA.FLA. Table 3-1 presents the results of the cross-correlation between the NARR dataset and these data. Hence, NARR was considered a suitable choice since it covers 1979 to the near present and provides daily data on a Northern Hemisphere Lambert Conformal Conic grid for a suite of meteorological variables. Although low correlation values between WS and RH were reported, NARR was still considered the most appropriate dataset to force the MIKE SHE simulations. The low correlation of these variables may have derived from the short period of record as well as small-scale (microclimate) effects at the station versus grid scale NARR data. NARR not 68 only provided a higher correlation compared to the average of FNA.FLA but was the only available source that reported continuous daily RH and WS. These two quantities were employed to estimate the daily amount of PET (see Section 3.3.3.3 for details). Table 3-1. Daily cross-correlation results between the NARR dataset, FNA.FLA, and UF. Stations Variables Period MAE RMSE R (Correlation) P-value FNA.FLA Maximum air temperature Minimum air temperature Mean air temperature Snow on ground Total precipitation Mean air temperature 1979-2014 3.0 ̊C 3.8 ̊C 0.97 p < 2.2 × 10-16 1979-2014 2.6 ̊C 3.5 ̊C 0.98 p < 2.2 × 10-16 1979-2014 4.2 ̊C 4.8 ̊C 0.98 p < 2.2 × 10-16 1979-2014 7.1 cm 11.6 cm 0.86 p < 2.2 × 10-16 1979-2014 1.1 mm 2.4 mm 0.60 p < 2.2 × 10-16 2013-2014 2.9 ̊C 3.6 ̊C 0.98 p < 2.2 × 10-16 Snow on ground RH WS 2013-2014 10.1 cm 16.7 cm 0.78 p < 2.2 × 10-16 2013-2014 2013-2014 10.2 % 1.5 m s-1 13.3 % 1.8 m s-1 0.59 0.43 p < 2.2 × 10-16 p < 2.2 × 10-16 FNA.FLA FNA.FLA FNA.FLA FNA.FLA UF UF UF UF 3.3.3.2. Tsea Lake Historical Meteorological Data Limited and insufficient meteorological parameters were recorded by the Tsea Lake station, which remains one of the main limitations of this part of the study. To prepare a historical meteorological dataset, parameters recorded at the Tsea Lake station (i.e., mean air temperature, rainfall, RH, and WS) were compared to the available FNA data (daily mean air temperature and rainfall) since it is the closest station to the Tsea Lake station (~70 km). Despite a significant correlation between mean air temperatures, a low correlation for rainfall was observed (Table 3-2). The mean air temperature, RH, and WS at the Tsea Lake station were compared to the NARR dataset, where a high correlation for the mean air temperature was discovered. 69 Additionally, a significant correlation between RH and an insignificant correlation for WS were identified. Similar to the Coles Lake watershed, the meteorological data from NARR were employed for Tsea Lake since it provided a high correlation of mean air temperature and provided a continuous record of daily RH and WS. Table 3-2. Daily cross-correlation results between Tsea Lake dataset and NARR and FNA. Station Variable Period MAE RMSE FNA Mean air temperature Total rainfall Mean air temperature RH WS 2011_2012 4.2 ̊C 2011_2012 2011_2012 2011_2012 2011_2012 FNA NARR NARR NARR 7.2 ̊C R (Correlation) 0.76 P-value p < 2.2 × 10-16 0.9 mm 4.5 ̊C 4.3 mm 7.4 ̊C 0.10 0.76 0.3584 p < 2.2 × 10-16 17.0 % 2.4 m s-1 21.4 % 3.0 m s-1 0.30 0.15 0.00119 0.03251 Therefore, similar to Prucha et al. (2011) and Wobus et al. (2013), baseline climate data for this study were derived from NARR, which include estimates of daily air temperature, precipitation, RH, and WS over 35 years. Climate parameters were considered as fixed input and were not subject to calibration, while other parameters that control snowmelt were adjusted through the calibration process. For simplification, a spatial distribution of the climate model setup was defined uniformly across the watersheds. 3.3.3.3. Potential Evapotranspiration Estimation PET is the amount of evaporation from soil and other surfaces as well as transpiration from plants that would occur with an available sufficient water source. PET plays a vital role in any hydrological study, specifically as one of the main inputs into the MIKE SHE model. Due to the difficulty of obtaining accurate field measurements, PET is commonly estimated from meteorological data. In past decades, many empirical or semi-empirical equations have been developed for estimating PET from meteorological data. Numerous researchers have analyzed the performance of the various methods for different regions. As a result of an 70 expert consultation held in May 1990, it is recommended to employ the Food and Agriculture Organization (FAO) of the United Nations’ Penman-Monteith method as the standard method for the computation of PET (Raes, 2009). To compute the PET rate, the principal weather parameters including the net radiation, air temperature, RH, and WS were acquired from the NARR dataset. Then, the ETo calculator (http://www.fao.org/land-water/databases-and-software/eto-calculator/en/) was used to compute daily PET (FAO Penman-Monteith) for both the Coles Lake and Tsea Lake watersheds over 35 hydrological years (Table 3-3). The details of PET computation are also described step by step in Appendix C. Table 3-3. Annual PET value for the Coles Lake and Tsea Lake watersheds. Hydrological Year 1979__1980 1980 1981 1981__1982 1982 1983 1983_1984 1984__1985 1985 1986 1986__1987 1987 1988 1988_1989 1989__1990 1990 1991 1991_1992 1992_1993 1993_1994 1994__1995 1995 1996 1996__1997 1997 1998 1998_1999 1999__2000 2000 2001 2001_2002 2002_2003 2003__2004 2004 2005 2005__2006 2006 2007 2007_2008 2008__2009 2009 2010 2010__2011 2011 2012 2012__2013 2013 2014 Annual Average Std. Deviation Coles Lake (mm) Tsea Lake (mm) 603.4 625.7 566.0 556.7 543.5 564.9 575.1 582.1 554.3 588.8 590.0 575.2 556.6 590.2 619.8 633.2 547.4 543.5 640.1 571.3 601.8 615.1 557.0 577.6 617.5 593.7 585.5 540.4 555.4 573.0 572.5 587.3 625.1 580.7 603.3 629.4 647.3 576.9 570.0 566.6 581.6 593.8 603.6 573.8 601.6 614.4 593.6 573.4 607.8 643.1 649.0 560.2 555.8 660.5 586.3 621.7 636.3 576.4 593.9 632.0 610.6 605.5 553.0 572.7 588.1 587.3 605.7 647.9 595.8 620.6 583.2 27.3 601.0 29.2 71 3.4. Land Surface Data The vegetation classes, representing the latest stage of vegetation cover in 2014 (Gonsamo and Chen, 2014), were mapped in ArcGIS and imported into MIKE SHE. Each vegetation class was assigned a representative LAI and rooting depth. Although MIKE SHE allows for changes in vegetation characteristics, to simplify the simulation, vegetation changes (i.e., LAI and rooting depth) were not included and were not subject to calibration during the simulation period. This should be a valid assumption, knowing the cold climate and short growing season prevent tree growth in northeastern BC. A similar approach was used by Foster (2014) and Voeckler et al. (2014). 3.4.1. Leaf Area Index The national LAI map (http://individual.utoronto.ca/gonsamo/index_files/Data.htm), created by the University of Toronto, was used as input for the MIKE SHE model (Gonsamo and Chen, 2014). This map uses an improved LAI algorithm based on the moderateresolution imaging spectroradiometer (MODIS) reflectance data at 250-m resolution (Gonsamo and Chen, 2014). The improvement was attained through pixel-by-pixel consideration of local topography, clumping index, and background reflectance variations (Gonsamo and Chen, 2014). Additionally, the LAI resolution was converted to a finer resolution to match the DEM, as this is an option within MIKE SHE. Overlaying the VRI spatial dataset on the two watersheds revealed that 13 vegetation classes of the BC Land Classification System (LCS) are represented in the Coles Lake basin, while the Tsea Lake basin contains 11 (Table 3-4 and Table 3-5). A Python script was then used to obtain statistics (range, mean and standard deviation) describing the distribution of LAI values associated within each class on the 1st, 72 11th and 21st days of each month from a 250 m pan-Canadian LAI dataset (Gonsamo and Chen, 2014) generated for 2008. The land surface distinctions for each class and their extension are illustrated in Figures 3-1 and 3-2. These values were considered variable seasonally but remained constant throughout the simulation period from 1979 to 2014. Table 3-4. Defined vegetation classes with the assigned LAI value for the Coles Lake watershed. Vegetation Class Jan Feb Mar Apr May June July Aug Sep Oct Nov Dec Water - Lake 0.0 0.0 0.0 0.0 0.1 0.6 0.9 0.6 0.1 0.0 0.0 0.0 Water - Stream 0.0 0.0 0.0 0.0 0.2 1.7 3.3 1.5 0.2 0.0 0.0 0.0 Upland - Bryoid - Open 0.0 0.0 0.0 0.0 0.0 1.3 2.8 1.5 0.2 0.0 0.0 0.0 Upland - Herb - Open 0.0 0.0 0.0 0.0 0.1 1.6 3.1 2.2 0.7 0.0 0.0 0.0 Upland - Herb - Dense Wetland - Herb - Open / Sparse 0.0 0.0 0.0 0.0 0.1 1.2 2.7 1.6 0.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.8 1.6 0.8 0.1 0.0 0.0 0.0 Wetland - Herb - Dense Upland - Shrub - Open / Sparse 0.0 0.0 0.0 0.0 0.2 2.3 3.9 1.8 0.3 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.9 3.9 5.4 3.7 1.2 0.0 0.0 0.0 Upland - Shrub - Dense Wetland - Shrub - Open / Sparse 0.0 0.0 0.0 0.0 0.2 2.0 3.4 1.7 0.3 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.4 2.7 4.2 2.4 0.3 0.0 0.0 0.0 Wetland - Shrub - Dense Upland - Treed - Open / Sparse 0.0 0.0 0.0 0.0 0.4 2.2 3.5 2.0 0.4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.9 3.6 4.5 2.7 0.8 0.0 0.0 0.0 Upland - Treed - Dense 0.0 0.0 0.0 0.0 0.3 2.0 3.3 1.8 0.4 0.0 0.0 0.0 Table 3-5. Defined vegetation classes with the assigned LAI value for the Tsea Lake watershed. Vegetation Class Jan Feb Mar Apr May June July Aug Sep Oct Nov Dec Upland - Exposed Soil 0.0 0.0 0.0 0.0 1.0 4.1 5.5 4.3 1.3 0.0 0.0 0.0 Water - Lake 0.0 0.0 0.0 0.0 0.0 0.3 0.7 0.6 0.2 0.0 0.0 0.0 Upland - Herb - Open 0.0 0.0 0.0 0.0 1.9 6.8 9.1 7.0 1.8 0.0 0.0 0.0 Upland - Herb - Dense Wetland - Herb - Open / Sparse 0.0 0.0 0.0 0.0 0.0 1.3 3.1 2.5 0.7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 1.3 2.8 2.0 0.6 0.0 0.0 0.0 Wetland - Herb - Dense Upland - Shrub - Open / Sparse 0.0 0.0 0.0 0.0 0.1 1.4 2.8 1.5 0.3 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 3.1 4.8 3.4 0.9 0.0 0.0 0.0 Upland - Shrub - Dense 0.0 0.0 0.0 0.0 0.0 1.3 3.2 2.5 0.5 0.0 0.0 0.0 Wetland - Shrub - Dense Upland - Treed - Open / Sparse 0.0 0.0 0.0 0.0 0.4 2.2 3.6 2.7 0.7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.6 2.7 4.0 3.1 1.0 0.0 0.0 0.0 Upland - Treed - Dense 0.0 0.0 0.0 0.0 0.2 1.8 3.2 2.4 0.5 0.0 0.0 0.0 73 Figure 3-1. Land surface distinctions within the Coles Lake watershed. 74 Figure 3-2. Land surface distinctions within the Tsea Lake watershed. 3.4.2. Maximum Root Depth Vegetation in northeastern BC varies considerably due to the transition from the Rockies eastward. According to the BC tree atlas, this area is mixed wood forest, primarily containing tamarack, white and black spruce, trembling aspen, and lodgepole pine in addition to grasslands (Parish et al., 1994). Tamarack typically has a shallow, compact, root system that develops in the upper 30 cm of the soil (Ministry of Forests and Range, 2007). White spruce 75 is usually characterized as shallow rooted (Ministry of Forests and Range, 2007). Although the bulk of the root biomass of black spruce is in the upper 20 cm of the organic horizon, some of its roots may penetrate to a depth of 60 cm in ideal conditions (Ministry of Forests and Range, 2007). Additionally, many other species such as paper birch and lodgepole pine with lateral and shallow root systems exist in northeastern BC. Furthermore, rooting depths depend on soil fertility and structure (Foster, 2014) with most root mass concentrated within the upper 40 cm of soil, regardless of soil quality (Curt et al., 2001). According to the VRI, the root depth varies from a minimum of 25 cm to a maximum of 75 cm in northeastern BC. Based on this available information, a rooting depth value of 40 cm was assigned to all of the tree classes; however, a rooting depth of 0 cm is allocated to the water bodies inside the watersheds. It is assumed that the rooting depth value remains constant during the simulation period. 3.4.3. Overland Flow Within MIKE SHE, overland flow is simulated using a diffusive wave approximation of the St. Venant equation. Overland flow settings require three main parameters: Manning’s M coefficient (reciprocal of Manning’s n) (m1/3 s-1), detention storage (mm), and initial water depth (m). These parameters were subject to calibration. The Manning’s M coefficient is equal to the Strickler roughness coefficient, which controls the magnitude of resistance and rate at which water can move horizontally (MIKE SHE user manual, 2007; Foster, 2014). That is, the resistance is controlled by the “roughness” of the land surface. Higher Manning’s M values lead to the faster routing of overland flow to the channel. In addition to the Manning’s M value, detention storage, and initial water depth values are required by MIKE SHE. Detention values are used to limit the amount of water flow 76 over the ground surface. The depth of ponded water must exceed the detention storage before water flows as sheet flow to the adjacent cell (Voeckler et al., 2014). Low detention storage implies that more water can flow over the surface resulting in higher streamflow (Sandu and Virsta, 2015). The initial water depth determines the available water in the initial period of simulation. Depending on the time constant specified, the initial depth may influence the results during a few days or up to several months of the simulation. During the simulation, both the detention and initial water depth remain as the model default equalling zero. Manning’s M values were inferred from recent VRI, and BC Land Cover Classification Scheme (BCLCS) maps and then assigned to each land classification. However, the initial Manning’s M values were taken from Liu and Smedt (2004) (Table 3-6). Similar to the LAI classification, 13 land surface types were defined for the Coles Lake watershed and 11 classes for the Tsea Lake watershed. A Manning’s M value was then assigned to each class (e.g., upland, wetland) within the MIKE SHE overland flow module. Table 3-6. Assigned Manning’s M estimates for controlling overland flow within the MIKE SHE model. Land Surface Manning’s M estimate (m1/3 s-1) Freshwater Stream beds Upland - Exposed Soil Upland - Bryoid - Open Upland - Herb - Open Upland - Herb - Dense Upland - Shrub - Open / Sparse Upland - Shrub - Dense Upland - Treed - Open / Sparse Upland - Treed - Dense Wetland - Herb - Open / Sparse Wetland - Herb - Dense Wetland - Shrub - Open / Sparse Wetland - Shrub - Dense Wetland - Treed – Open / Sparse 100.0 50.0 50.0 3.4 3.3 2.9 2.5 2.5 2.5 1.3 20.0 18.0 25.0 20.0 33.3 77 The spatially distributed grid files for the Manning’s M coefficients were used for both watersheds (Figures 3-3 and 3-4). Figure 3-3. Manning’s M value for the Coles Lake watershed. 78 Figure 3-4. Manning’s M value for the Tsea Lake watershed. 3.5. Unsaturated and Saturated Zone In MIKE SHE, the unsaturated zone (UZ) is the zone through which the water table rises and falls vertically at each grid point; unsaturated flow in the UZ is modelled using Richards’ equation module. Richards’ equation was employed since it better represents unsaturated flow under dynamic conditions (DHI, 2007). Richards’ equation is based on Darcy’s law and the continuity equation. In Richards’ equation, the driving force for transport of water is the gradient of the hydraulic head, which includes gravitational and pressure components (Graham and Butts, 2006). After precipitation occurs, it infiltrates the UZ, where it will either evaporate or move downward according to the soil’s moisture retention curve and the hydraulic conductivity function (KƟ) (Foster, 2014). The soil moisture retention curve and 79 the hydraulic conductivity function are defined using estimated van Genuchten parameters (3.1). Ɵψ =Ɵr + ƟS -Ɵr [1+(αψ)n ] m (3.1) where Ɵψ is the soil water content (cm3 cm-3), Ɵr is the soil residual water content (cm3 cm-3), Ɵs is the soil saturated water content (cm3 cm-3), ψ is soil water potential (kPa), constant alpha (α) is a scale parameter inversely proportional to mean pore diameter (cm-1), n and m are the shape parameters of the soil water characteristic curve, with m = 1−1/n, 0 < m < 1 (Yang and You, 2013). For saturated flow simulations, the 3-D finite difference method is used. The inputs for the saturated zone (SZ) are geologic properties of the soil profile such as horizontal and vertical hydraulic conductivities, specific yield and specific storage. The soil classes and their depths were determined through fieldwork by Quicksilver at Coles Lake. After that, the soil profile was defined within the model based on the collected soil information for both watersheds. In the model setup, the soil layers and bedrock were assigned to the UZ, while the SZ contains only bedrock. During a simulation, Voeckler et al. (2014) discovered that the water table fell below the depth of the UZ in higher and steeper areas of the Okanagan Basin and the model did not converge. Thus, to overcome this problem, bedrock was also assigned to the UZ (Figure 3-5). A similar approach was used by Foster (2014) at the Cowichan watershed on Vancouver Island. Accordingly, three layers of organic layer muskeg (peatland), glacial till (clay-silt-sand), and bedrock comprised of marine shales and siltstones were defined for the UZ. Layers 1 to 2 comprise the upper “soil” layer in MIKE SHE, and layer 3 is the bedrock. The thicknesses of the first and second layers were set to 0.3 m and 3.0 m, respectively. The bedrock was 80 assigned to a depth of 300 m in both the UZ and SZ domains. Table 3-7 lists the complete soil parameters for each layer. These layers are considered uniform over the entire model domain. Figure 3-5. An example of the soil profile representation of the unsaturated zone (UZ) and the saturated zone (SZ). Table 3-7. Initial soil parameters for each layer. Soil Layer Texture a Depth (m) 1b O 2c G 0.0 to 0.3 0.3 to 3.3 3d B a 3.3 to 300 Saturated hydraulic conductivity Ksat (m s-1) 2.8 × 10-6 Saturated moisture content (Ɵs) (cm3 cm-3) Constant alpha (α) (cm-1) Constant n Specific Yield Bulk density (kg m-3) 0.80 Residual moisture content (Ɵr) (cm3 cm-3) 0.130 0.02 1.60 0.27 200 2.4 × 10-7 0.41 0.095 0.019 1.31 0.10 2100 2.9 × 10-7 0.10 0.050 0.0036 2.75 0.01 1200 O: Organic layer muskeg (peatland), G: Glacial till (clay-silt-sand), B: Bedrock marine shales and siltstones. b Values from Letts et al. (2000). c Values from Leij (1996). d Values from Leij (1996) and Foster (2014). The bedrock layer was discretized as a single computational layer. The hydraulic conductivity (K) values were assigned uniformly over the full 300 m thickness, and the base of the model was assigned as a no flow (zero flux) boundary. The initial K values of 81 horizontal hydraulic conductivity (Kx = 1× 10-7 m s-1), vertical hydraulic conductivity (Ky = 1× 10-7 m s-1), specific yield (Sy = 0.01) and specific storage (Ss = 1× 10-4 m-1) were assigned for both watersheds. Values for horizontal and vertical conductivities were subject to calibration. Specific yield and specific storage values were not subject to calibration given the absence of data. Therefore, these storage parameters remain highly uncertain. In addition, the groundwater module runs several times to reach equilibrium, during the simulation period. Accordingly, the initial groundwater head from the previous simulation was updated as the groundwater system gradually attained equilibrium (DHI, pers. comm., 2016). Upon completion, the UZ and SZ were explicitly coupled (run in parallel). 3.6. Stream Network and Hydrometric Flow Data MIKE 11 uses the HD module to define the channel flow within the domain area in MIKE SHE. The HD module comprises several components; however, only four main components (the river network, river cross-sections, boundary data and HD parameters) were applicable for this study. Within MIKE 11, lakes and rivers are represented as one-dimensional features using the St. Venant equation. For the Coles Lake and Tsea Lake watersheds, the lake and river network maps were created from the National Hydro Network (NHN) (http://open.canada.ca/en, 2017) in ArcGIS. Both Coles Lake and Tsea Lake are small and shallow water bodies. According to bathymetric surveys completed by Matrix Consulting Group Inc. in 2010, the deepest part of Coles Lake is 2 m in depth, with an average depth of 0.5 m near the edge of the lake (Figure 3-6). Unfortunately, the Tsea Lake bathymetric map was not available. Therefore, one of the significant challenges was a lack of lake morphology and runoff data. As a result, in the initial stages of modelling, the runoff was underestimated. To correct this, additional 82 lake cross-sections similar to those along the river network were created to allow more communication between the overland component in MIKE SHE and the stream network in MIKE 11. In fact, the lake is treated similarly to a river whereby a continuation of the river network that widens represents the lake using the cross-section and lake bathymetry. However, it was still possible to define the area of the lake using the flood code inside MIKE SHE and outline the subsurface topography of the lake applying the cross-sections in MIKE 11. Therefore, bed topography and the extent of Coles Lake and Tsea Lake are specified in three dimensions using the flood code in MIKE 11. Figure 3-6. Bathymetry of Coles Lake at 0.25 m intervals (Matrix Consulting Group Inc., 2010). Upon importing the lake boundary and physical specifications into MIKE 11, the river networks were created for both watersheds. The river network shapefile was converted to river nodes (H-points), then discharge (Q) and stage levels (h) were defined for both watersheds. The discharge measurements were extrapolated from points in between input cross sections (Foster, 2014). Stage measurements were calculated at all H-points and were determined by the dynamics of flow within the cross-sections. River cross-sections were assigned to each stream using data obtained through fieldwork. 83 In addition, boundary conditions were assigned for the river network at each stream. A closed boundary was set to the upstream end of the network, whereas a water level (an open boundary) was assigned to the downstream end at Emile Creek and Tsea River. All other HD computation parameters (e.g., wind factor) were set to MIKE 11 default values. Upon completion of the MIKE 11 setup for the Coles Lake and Tsea Lake watersheds, the MIKE 11 models were dynamically coupled with MIKE SHE through river links. 3.7. Model Calibration All physical and climatological datasets were collected, formatted appropriately, and imported into the MIKE SHE model using two different setups for each watershed. The model was calibrated for a simulation period of 35 water years (1979/1980 to 2013/2014). Model calibration provides the perturbation of model parameters within reasonable and realistic ranges. Calibration quantifies the agreement between simulated model predictions and measured data, and precise calibration significantly improves the model outcome (MacLean, 2009). Generally, there are two main methods _ manual or automatic _ used to calibrate the hydrological model. Manual calibration methods include trial-and-error, whereas the autocalibration procedures use an optimization algorithm (MacLean, 2009). For this research, model calibration was completed manually using the trial-and-error approach, with each successive simulation having modified input parameters that were adjusted according to the results of the previous simulation. The manual calibration was chosen rather than autocalibration in that auto-calibration only incorporates streamflow, whereas this research emphasizes utilizing (maximizing) various sources of information (snowmelt and pressure head data). Dekker (2013) employed the same approach for his research for the Lovers Creek 84 sub-watershed near Barrie, Ontario, which uses specific inputs targeted to produce the desired results. For this current research, the snowmelt, pressure head, and streamflow data were used during the trial-and-error calibration (Table 3-8), and the calibration results are described in Section 3.8. During the calibration process, ranges of values were applied to the parameters and the results reviewed using the hydrographical and statistical analysis. After the calibration phase was completed, a simple sensitivity analysis was conducted to examine the influence of each parameter. Analysing parametric sensitivity is vital to identify which parameters have the most influence on the model’s performance (McCuen, 1973). Details of sensitivity analyses for the calibration parameters are discussed in Section 3.9. The trial-and-error approach leads to a stronger understanding of model parameters and can be highly effective when conducted by experienced modellers or modellers with good process-based knowledge of the landscape being modelled, compared to the auto-calibration method (Madsen, 2000). The trial-and-error approach, along with a flexible parameterization, allows researchers to successfully overcome the problems encountered in the calibration phase (Refsgaard and Storm, 1995). The trial-and-error approach was used commonly by various groups (e.g., DHI, GRRG) and researchers (Dekker, 2013; Foster, 2014; Voeckler et al., 2014; De Schepper et al., 2017). Although the trial-and-error technique has been used extensively, this technique remains limited as it does not utilize the full range of parameters considered by the model. Therefore, some researchers use auto-calibration, where parameters are adjusted according to a specified search algorithm, and numerical measures asses the goodness-of-fit for the new parameter set (Madsen, 2000). Though auto-calibration is much faster than trial-and-error 85 calibration (MacLean, 2009), in a study on the Lovers Creek watershed, Dekker (2013) found that an auto-calibration that begins from a simulation that is too far from an acceptable fit with observed results. The simulations may take several weeks to complete and may or may not converge on a solution (Dekker, 2013). Where possible, the manual calibration is first used to bring the model within a practical range, then sensitivity analysis can determine the most sensitive parameters and their expected practical ranges (for this landscape). Autocalibration can then be used to converge on the optimal solution. Furthermore, when Dekker (2013) attempted to set up an auto-calibration, some of the important and sensitive parameters (e.g., Manning’s roughness) could not be used in the calibration. This occurred as a result of the way these parameters had been specified within MIKE SHE (Dekker, 2013). For instance, the hydraulic conductivity is a sensitive parameter within the SZ and thus a good candidate for calibration (Dekker, 2013). Hydraulic conductivity was specified using a DHI two-dimensional grid file (*.dfs2) in MIKE SHE. However, these types of files cannot be modified by the auto-calibration tool. Similar issues apply to the Manning’s roughness and depression storage (Dekker, 2013). Both MIKE SHE and MIKE 11 were calibrated simultaneously since modifications to a calibration parameter within one model could influence results in the other (Thompson et al., 2004). Although the models were set up to simulate 35 years, the model calibration was carried out based on available data (Table 3-8). That is, there were insufficient data available to perform a meaningful split sample model validation. Pressure head data were available for Coles Lake (October 2013_September 2014) but not for Tsea Lake. In fact, long-term concurrent runoff and pressure head data for the two lakes were not available to ensure the 86 reliability of the model to reproduce out-of-sample river flow regimes, a typical situation in ungauged basins (Winsemius et al., 2009). Table 3-8. Observed data used to calibrate the model for the Coles Lake and Tsea Lake watersheds. Observed Data Coles Lake Tsea Lake Snowmelt Pressure Head Streamflow October 1979_September 2014 October 2013_September 2014 Inlet and Outlet October 2013_September 2014 Emile Creek October 2011_September 2014 October 1979_September 2014 Not available Tsea River April 2010_October 2011 3.8. Calibration Results 3.8.1. Snowmelt Snowmelt modelling is an essential first step in the model calibration process as the air temperature affects the snow accumulation during winter as well as the release of frozen water during the freshet seasons (Foster, 2014). Snowmelt was calibrated using SWE data (October 1979_ September 2014), which were measured manually at the Fort Nelson station operated by the BC River Forecast Centre (RFC). Similar to Foster (2014), snow accumulation and melt (timing and release) were calibrated following a step-wise methodology, whereby many parameters were adjusted independently to produce realistic results through the calibration process. The parameters with the largest effect on the snowmelt model included: melting temperature, precipitation, degree-day melting coefficient, and maximum wet snow fraction in snow storage (Table 3-9). Table 3-9. Final parameters used for the snowmelt module. Snowmelt Parameter Melting temperature Degree-day coefficient Minimum snow storage Maximum wet snow fraction Value 0.0 Variable between 0.5 to 2.0 0.0 0.1 Units ̊C mm ̊C -1 day-1 mm Fraction 0.0 0.0 mm Fraction Initial total snow storage Initial wet snow fraction 87 In general, snow accumulation and snowmelt were reasonably well simulated for both the Coles Lake and Tsea Lake watersheds (Figure 3-7), with better agreement for the latter. This is a reasonable result since the Tsea Lake watershed is closer to Fort Nelson than Coles Lake. The calibration fit statistics are given in Table 3-10. Small Nash-Sutcliffe Efficiency (NSE) scores are reported from the MIKE SHE simulations for both watersheds. The closer to 1, the more accurate the model is. When NSE = 1, the modelled data are perfectly matched to the observed data. NSE ranges from -∞ to 1, where a negative score poorer than climatology, and with anything between 0 and 1 having better predictive capacity than the climatology. Although there are several statistical goodness-of-fit criteria (e.g., Kling–Gupta efficiency) that can be used for the calibration process, the NSE was selected to assess the predictive power of the MIKE SHE model (Gupta et al., 2009). According to Gupta et al. (2009), the mean squared error (MSE) and the NSE are the two criteria most widely used for calibration and evaluation of hydrological models using observational data. Figure 3-7. Comparison of daily Fort Nelson measured (dots) and Coles lake and Tsea lake simulated (solid line) SWE for the period October 1979_September 2014. The elevations of the Fort Nelson, Coles Lake, and Tsea Lake stations are 370 m, 480 m and 660 m, respectively. 88 Table 3-10. Calibration results for the daily SWE (October 1979_September 2014). [Zero values were neglected to avoid inflating the performance of the results.] Calibration Station ME (mm) MAE (mm) RMSE (mm) STDEV Residuals (mm) NSE R (Correlation) p-value Fort Nelson/Coles Lake -19.64 23.36 40.69 35.63 0.20 0.78 7.2 × 10-5 Fort Nelson/Tsea Lake -18.85 22.43 39.64 34.87 0.20 0.79 7.8 × 10-5 3.8.2. Pressure Head in the Soil Zone Shallow groundwater levels were measured with nine piezometers installed at three piezometer nests around Coles Lake from 2 February to 30 September 2014. As discussed in the previous chapter, piezometer transects were installed strategically around Coles Lake to capture the overall fluctuation of shallow groundwater near the lake. The complete procedure of the piezometer installation was described in Section 2.1.4. For Transect 1 (P1-1, P1-2, and P1-3), the observed curve is smoother compared to the simulated response (Figure 3-8). In addition, the observed results indicate that the peak flow occurs earlier compared to the simulated while the total observed volume is smaller compared to the simulated volume. Overall, the model simulates a steady pressure head from the beginning of November to May (freezing period) at all three transects. MIKE SHE begins to capture the water fluctuation after the start of the melting period (ice-free conditions). The overall results at Transect 1 indicate a low correlation with the observed data, NSE values are negative, while the model generates pressure heads that are considerably higher compared to the observed values (Figure 3-8). 89 Unlike Transect 1, the simulated pressure heads in Transect 2 (P2-1, P2-2, and P2-3) are well simulated and exhibit higher correlations while the NSE values remain negative (Figure 3-9). However, like Transect 1, the model generates pressure heads that are higher compared to the observed values. P2-1 and P2-2 show a higher correlation compared to P2-3, which may be a result of the shorter comparison period. Overall, MIKE SHE simulated similar patterns, compared to the observed shallow water table at Transect 2. Figure 3-8. Comparison of daily measured (dots) and simulated (solid line) pressure heads at Transect 1 for the period October 2013_September 2014. Simulated trends at Transect 3 (P3-1, P3-2, and P3-3) demonstrate an increase in the water table from February until the end of April (freezing period) relative to observed data (Figure 3-10). After May (melting period), the dynamics of the shallow water table are 90 reasonably well simulated. Similar to Transects 1 and 2, the model produced pressure heads that are higher compared to the observed data. Although the dynamics of the shallow water table are reasonably simulated as above, the overall pressure heads for all nine piezometers are higher compared to their observed values. Similar results were reported by Voeckler et al. (2014) and De Schepper et al. (2017), where MIKE SHE is neither able to simulate the dynamics of observed hydraulic heads in the shallow piezometers, nor able to simulate the differences observed between piezometers. In this study, for the calibration of pressure head, one point is being compared against the grid size of 100 m by 100 m, this could have been improved by reducing the grid size, which would significantly increase the simulation time. Additionally, the existence of a discontinuous permafrost zone within northeastern BC watersheds may influence the calibration results. Permafrost is a factor that might influence the spatial distribution of recharge and the groundwater flow paths. When the ground is frozen (seasonally), the hydraulic conductivity diminishes. But MIKE SHE assigns a single value throughout the year as the model does not support varying hydraulic properties. Thus, if the hydraulic conductivity value is significantly higher, too much water gets into the groundwater system, leading the model to produce higher pressure heads. In this study, MIKE SHE indeed produced higher pressure heads during the calibration due to its inability to resolve frozen ground processes. Frozen ground results in significant modification of groundwater flow systems within this region (Williams, 1970). Frozen ground acts as a confining and impermeable layer that restricts recharge, discharge, and movement of groundwater and limits the volume of unconsolidated deposits and bedrock in which liquid water may be stored (Williams, 1970). Understanding of local 91 variability of permafrost distribution is important as it can impact the hydrology of specific sites (Williams, 1970). Groundwater recharge may occur when the base of the permafrost is wet during a warm period, but may be prevented when the ground is frozen. Figure 3-9. Comparison of daily measured (dots) and simulated (solid line) pressure heads at Transect 2 for the period October 2013_September 2014. Overall, in all transects the simulated pressure head peaks at 486 m as this is the maximum height possible in all locations. Average statistical results for the pressure head include a low RMSE average of 2.86 m, with a correlation average of 0.57, ranging from 0.06 to 0.93. In all piezometers, NSE values are negative (Table 3-11). 92 Figure 3-10. Comparison of daily measured (dots) and simulated (solid line) pressure heads at Transect 3 for the period October 2013_September 2014. Table 3-11. Statistics of simulation results of daily pressure head in the SZ for all piezometers (February to September 2014). Calibration Site ME (m) MAE (m) RMSE (m) STDEV Residuals (m) NSE R (Correlation) p-value Transect 1 (P1-1) -2.08 2.08 2.09 0.27 -59.52 0.34 6.2 × 10-9 Transect 1 (P1-2) -1.61 1.61 1.65 0.28 -29.47 0.40 8.3 × 10-3 Transect 1 (P1-3) -1.34 1.34 1.36 0.23 -25.53 0.51 2.2 ×10-2 Transect 2 (P2-1) -1.10 1.10 1.12 0.17 -18.22 0.88 4.4 × 10-34 Transect 2 (P2-2) -1.27 1.27 1.30 0.26 -11.53 0.93 4.4 × 10-52 Transect 2 (P2-3) -1.04 1.04 1.11 0.36 -2.62 0.91 3.6 × 10-32 Transect 3 (P3-1) -0.77 0.79 0.86 0.39 -4.53 0.06 5.9 × 10-1 Transect 3 (P3-2) -0.84 0.85 0.97 0.49 -3.02 0.13 1.3 × 10-5 Transect 3 (P3-3) -1.22 1.22 1.23 0.18 -32.20 0.93 1.4 × 10-39 93 3.8.3. Streamflow The calibration process for streamflow at the Coles Lake and Tsea Lake watersheds involved comparing the observed daily hydrometric data to the simulated stage and discharge results from MIKE 11. For this purpose, hydrometric data from the Inlet (upstream), Outlet (downstream) and Emile Creek (downstream) stations were employed for the Coles Lake watershed (Figure 2-29). The hydrometric data for the Tsea Lake station (downstream) were used for that watershed. The periods of recorded data vary by the station with lack of historical data and no continuous data recorders. Hydrometric data were recorded daily from 14 May to 30 September 2014 for the Inlet and Outlet at Coles Lake, 8 August 2011 to 30 September 2014 for Emile Creek, and 1 May 2010 to 13 August 2011 for the Tsea Lake Outlet. In both watersheds, the MIKE SHE model provided reasonable results at an early stage. Foster (2014) obtained similar results where simulated water levels matched the observed data points in the initial years of the running period. To achieve a successful calibration, careful attention was given to define the lake bathymetry and stream morphology system using the flood codes classification in MIKE SHE. Additionally, the HD parameters (i.e., Q and h) within MIKE 11 were defined for both Coles Lake and Tsea Lake cross-sections as they contributed significantly to the calibration of the model (Table 3-12). The dynamics of discharge are reasonably well simulated as peaks occur around May and then decline until the flow approaches zero (Figures 3-11 and 3-12). This result was already confirmed during the fieldwork. The simulated results matched the observed data points at the Coles Lake station (Table 3-13). The simulation results at the Emile Creek station (2011_2014) produce a poor correlation compared to the observed discharge (Figure 3-13). 94 Some simulated peak levels are underestimated despite a good correlation for SWE at the Coles Lake watershed (see Section 3.8.1). Table 3-12. Final HD value for the Coles Lake and Tsea Lake cross-sections. More cross-sections are defined for the Tsea Lake as the bathymetric map was not available. Sites Coles Lake Coles Lake Coles Lake Tsea Lake Tsea Lake Tsea Lake Tsea Lake Tsea Lake Tsea Lake Tsea Lake Tsea Lake Tsea Lake Tsea Lake Tsea Lake H (m) 482.0 481.7 481.5 696.8 696.8 662.4 663.3 659.7 659.7 659.7 659.7 659.7 659.7 659.7 Q (m3 s-1) 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 Figure 3-11. Comparison of daily measured (dots) and simulated (solid line) discharge at the Inlet hydrometric station on Coles Lake (upstream) for the period October 2013_September 2014. Figure 3-12. Comparison of daily measured (dots) and simulated (solid line) discharge at the Outlet hydrometric station on Coles Lake (downstream) for the period October 2013_September 2014. 95 Figure 3-13. Comparison of daily measured (dots) and simulated (solid line) discharge at the Emile Creek hydrometric station (downstream) for the period October 2011_September 2014. The 2013 discharge appears more reasonable compared to 2011, 2012, and 2014. Table 3-13. Calibration results of the MIKE 11/SHE daily channel flow for the Inlet, Outlet, and Emile Creek at the Coles Lake watershed. Calibration Station Period ME (m3 s-1) MAE (m3 s-1) RMSE (m3 s-1) STDEV Residuals (m3 s-1) NSE R (Correlation) p-value Inlet Oct 2013_ Sep 2014 1.29 1.29 1.73 1.15 -0.43 0.85 9.49 × 10-22 Outlet Oct 2013_ Sep 2014 1.97 2.31 3.34 2.70 -0.26 0.85 7.98 × 10-40 Emile Creek Oct 2011_ Sep 2014 1.58 2.63 4.48 4.19 0.04 0.40 9.70 × 10-15 Figure 3-14 shows the dynamics of streamflow discharge are well simulated for the Tsea Lake watershed. Although simulated and observed values correlate well, some simulated peak levels are underestimated (Table 3-14). Average statistical results of streamflow discharge include a low RMSE average of 3.18 m3 s-1 and 1.21 m3 s-1 with a correlation average of 0.70 and 0.62, respectively, for the Coles Lake and Tsea Lake watersheds, respectively. 96 Figure 3-14. Comparison of daily measured (dots) and simulated (solid line) discharge at the Tsea Lake hydrometric station for the period (April 2010_October 2011). Table 3-14. Calibration results of the MIKE 11/SHE daily channel flow for the Tsea Lake watershed. Calibration Station Period ME (m3 s-1) MAE (m3 s-1) RMSE (m3 s-1) STDEV Residuals (m3 s-1) NSE R (Correlation) p-value Tsea Lake Apr 2010 _ Oct 2011 0.21 0.84 1.21 1.18 0.37 0.62 9.40 × 10-37 Overall, peak discharge corresponds to enhanced winter precipitation events and the freshet. The timing and magnitude of the peak flows were calibrated by adjusting the overland flow characteristics. Because none of the small streams from the tributaries to Emile Creek (Coles Lake watershed) and Tsea River (Tsea Lake watershed) were included, the peak flows may be underestimated. Flow regime underestimation may also be related to the lack of a frozen ground module in the MIKE SHE model. Further, some storage changes in the sub-surface during winter were observed in the simulation. As mentioned previously, the soil properties cannot be modified in MIKE SHE to account for frozen ground (i.e., permafrost, seasonal frost). This can influence the streamflow calibration results where, at the early freshet stage, the ground is still frozen, but instead of resulting in streamflow peaks, there is less variability in streamflow because water recharge the groundwater system instead of being routed over land. 97 In summary dynamics of the snowmelt, pressure heads, and streamflow discharge are realistically simulated by MIKE SHE for both watersheds. Upon the compilation of model calibration, the MIKE SHE model was employed to simulate the Coles Lake water balance, to quantify the influence of vegetation cover at the Coles Lake watershed scale, and to simulate the historical water balance of the Coles Lake and Tsea Lake watersheds (Figure 3-15). Figure 3-15. Description of hydrological model input variables, model components, and outputs. 3.9. Parameter Sensitivity and Experimental Calibration During the model set up, careful attention was given to the model parameters to simulate the hydrological system with minimum uncertainty. These parameters were classified into three categories: parameters that can be measured, parameters that can be obtained with 98 certainty from the literature, and parameters that need to be calibrated (Foster, 2014). However, the following section emphasizes the calibrated parameters, and discusses how these parameters may affect the model output. According to Jaber and Shukla (2012) and Foster (2014), complete parameterization based on field observations at the watershed scale is nearly impossible. Further, complete parameterization can be especially unrealistic when investigating the simulation of spatial and temporal variability (Jaber and Shukla, 2012), as well as when the study area is located in a remote area. Due to the large number of parameters used in the MIKE SHE model as applied to the Coles Lake and Tsea lake watersheds, much time was assigned to the proper setup and calibration of the hydrological model, rather than a rigorous sensitivity analysis, which would have been time-consuming (Foster, 2014). Both watersheds were calibrated using snowmelt, pressure head, and streamflow values. In general, through the development of the model and calibration process, the following aspects were found to significantly influence the model results: 3.9.1. Sensitivity of Snowmelt Parameters The early stages in the model set up focused on the calibration of the snowmelt model. Due to the remoteness of the study area, historical SWE data were not available specific to the Coles Lake and Tsea Lake watersheds. Therefore, snowmelt was calibrated using SWE data from the Fort Nelson station, which are measured manually and may contain human and instrumental errors. Although the measuring station is less than 140 km from both watersheds, the SWE data from the Fort Nelson station may not entirely be representative of the study area, as the SWE is strongly variable locally (Schmidlin, 1989). 99 Because the Coles Lake and Tsea Lake watersheds are located in snowmelt-dominated hydrological regimes, the degree-day coefficient was a key parameter for calibration. For instance, during simulations, the degree-day coefficient influenced the snowmelt and simulated hydrograph. The spring melt was occurring too slowly which results in simulating less initial spring melt, and a snow melt period that lasted further into the summer. This issue was resolved by increasing the degree-day coefficient. Although the degree-day coefficient value was considered variable (Figure 3-16), the annual patterns remained constant throughout the simulation period from 1979 to 2014. The degree-day coefficient was at the lowest (̴0.5 mm ̊C-1 day-1) at the beginning of winter (November) and highest (2.0 mm ̊C-1 day-1) at the beginning of freshet (May). When the air temperature is above the threshold melting temperature (0 C ̊ ), the snow begins to melt, and consequently, snow storage is reduced (Equation 3.2): qsnow=Degree day factor* Air temp ˗ Threshold temp *∆t - - - - (3.2) where qsnow is snowmelt [mm day-1], Degree- day- factor [mm ̊C-1 day-1] is the amount of snow that melts per day for every 1 ̊C of air temperature above the threshold melting temperature. 3ℎ456ℎ789: ;5<= is the temperature at which the snow starts to melt, taken here as 0 ̊C and ∆t the UZ/ET time step. Snow storage will be reduced to zero if 0>?@A exceeds SWE at the previous time step (MIKE SHE User Guide, 2017). If the air temperature is below the threshold melting temperature, then the ET module will remove water from the snow storage as sublimation before any other ET is removed (MIKE SHE User Guide, 2017). A complete range of snowmelt parameters used during the calibration as well as the final value for the snowmelt module is reported in Table 3-15. 100 Figure 3-16. Trends of degree-day melting coefficient. The degree-day coefficient default is set at 0 (mm ̊C-1 day-1); this was adjusted during the model calibration process which ranges from 0.5 to 2.0 (mm ̊C-1 day-1). Table 3-15. The range of snowmelt parameters used for the snowmelt module. Snowmelt Parameter Melting temperature Degree-day coefficient Minimum snow storage Maximum wet snow fraction Initial total snow storage Initial wet snow fraction Range of Snowmelt Parameter Final Value Units -1 to 1 0 to 3 0 0 to 0.5 0 0 0 0.5 to 2.0 0 0.1 0 0 ̊ C mm ̊C -1 day-1 mm Fraction mm Fraction 3.9.2. Sensitivity of Hydraulic Head Parameters Unfortunately, limited hydraulic head data were available that could be used for the calibration of the SZ. Thus, the model was only calibrated to groundwater levels in a single uniform aquifer unit. Overall during the calibration process, the following parameters were found to be the most sensitive, based on impact on streamflow. 101 § The horizontal and vertical hydraulic conductivity parameters were sensitive within the SZ (see Section 3.5 for complete details of SZ parameters). § When a lower hydraulic conductivity was applied in the SZ, the peak surface water flows increased and the base flows decreased. § Within the UZ the hydraulic conductivity (Ks) and the specific yield (Sy) proved to be more sensitive compared to other parameters. § When a lower hydraulic conductivity was applied in the UZ, the surface water flows increased, and, to a small degree, the groundwater head elevations decreased. § When the lower value was assigned to specific yield (Sy) in the UZ, the surface water flows increased. § The relationship of bulk density and hydraulic conductivity was reciprocal within the UZ parameters. Low bulk density was assigned to the upper layer of the soil, while the bulk density increased within the middle and deeper layers. High bulk density was assigned to explain low soil porosity and soil compaction, which indicates a poor movement of water through the soil layer. § Table 3-16 reports a complete range of soil parameters within the UZ used during the calibration. Table 3-16. The range of soil parameters for each layer. Soil Layer Texture a Depth (m) 1 O 0.0 to 0.3 2 G 0.3 to 3.3 3 B 3.3 to 300 Saturated hydraulic conductivity Ksat (m s-1) 2.8 × 10-6 to 2.8 × 10-6 2.4 × 10-7 to 2.4 × 10-7 2.9 × 10-7 to 2.4 × 10-7 Saturated moisture content (Ɵs) (cm3 cm-3) 0.50 to 0.80 Residual moisture content (Ɵr) (cm3 cm-3) 0.130 Constan t alpha (α) (cm-1) 0.02 Constant n Specific Yield Bulk density (kg m-3) 1.60 0.27 50 to 400 0.31 to 0.41 0.095 0.019 1.31 0.10 1500 to 2400 0.08 to 0.12 0.050 0.0036 2.75 0.01 1000 to 1800 a O: Organic layer muskeg (peatland), G: Glacial till (clay-silt-sand), B: Bedrock marine shales and siltstones. 102 3.9.3. Sensitivity of Streamflow Parameters Overland flow was simulated using a diffusive wave approximation of the St. Venant equation. Overland flow settings require three main parameters: Manning’s M coefficient (reciprocal of Manning’s n) (m1/3 s-1), detention storage (mm), and initial water depth (m). Among these parameters, the Manning’s M coefficient was found to be the most important parameter, as it controls the magnitude of resistance and rate at which water can move horizontally (MIKE SHE user manual, 2007; Foster, 2014). Higher Manning’s M values led to the faster routing of overland flow to the channel. Complete details of Manning’s M values are reported in Section 3.4.3. 3.9.4. Sensitivity of Vegetation Parameters In addition to snowmelt, pressure head, and streamflow parameters, the sensitivity in the vegetation parameters was evaluated by changes in LAI values. Section 5.3 explains how LAI values play a significant role in the water balance of boreal watersheds. AET was elevated in areas where the LAI value is higher (i.e., there is more vegetation), e.g., dense cover compared to mixed or sparse canopies. Complete details of the vegetation sensitivity tests are described in Section 5.3. 3.9.5. Sensitivity to the Meteorological Forcing In addition to the vegetation parameters, meteorological forcing data (e.g., precipitation, air temperature, and PET) form an essential input to the MIKE SHE model. Thus, a meteorological sensitivity experiment was applied to extract indicators from the vegetation cover that influence the Tsea Lake water balance. This was accomplished by removing the influence of atmospheric forcing data on the simulation results. To this end, meteorological forcing data from the Coles Lake watershed were applied to the Tsea Lake watershed 103 (October 1979_September 2014) as a sensitivity test, while other parameters remained unchanged. Complete details of the meteorological sensitivity test are described in Section 6.3. 104 CHAPTER 4: QUANTIFICATION OF THE COLES LAKE WATER BALANCE The primary objective of this chapter is to estimate the short-term water balance (hydrological year 2013_2014) of Coles Lake. Baseline information generated from this study will support the assessment of the sustainability of current and future strategies for freshwater extraction in the region. Upon the estimation of the hydrological components of Coles Lake, the results are evaluated using a hydrological model. Overall, this chapter is divided into three sections, sequentially: 1) Estimation of the Coles Lake water balance components using field observations. 2) Detailed simulation of Coles Lake hydrological processes using MIKE SHE. 3) Evaluation of MIKE SHE simulated water balance with observations. 4.1. Coles Lake Water Balance To minimize errors, careful attention was given during the field campaign as well as the analysis to quantifying the water fluxes and stores within the study area for one water year (October 2013_September 2014). Details of the field campaign were previously discussed in Chapter 2. Wetlands are an essential part of the hydrological cycle within this region, with water balance components including precipitation, evaporation, surface water and groundwater flows (Ingram, 1983). In addition, due to industrial usage, a withdrawal component needs to be considered. Hence, as a first step, an appropriate water balance Equation (4.1) was defined to represent these components for Coles Lake: ∆S = P + Qi ± G - E - Qo - W (4.1) where DS is the change in stored water, P is precipitation (P = R + S, where R denotes rainfall and S represents snowfall), Qi is the mean annual stream inflow to Coles Lake, G is 105 the groundwater exchange with Coles Lake, E is evaporation from the lake, Qo is discharge (outflow) from the lake, and W is the licensed withdrawal of water. All terms are expressed in units of millimeters over the area of the Lake itself for one hydrological year. All water balance components are subject to measurement errors, and these errors are therefore reported in the water balance results (Sokolov and Chapman, 1974). 4.1.1. Water Balance Components 4.1.1.1. Contribution of Precipitation to Coles Lake Local rainfall data were collected at 15-minute intervals; however, the 15-minute values were summed daily to calculate the total contribution of rainfall to the Coles Lake water balance. For all calculations, rainfall from 1 November 2013 to 15 May 2014 was taken as zero because air temperatures were continuously sub-freezing during this period and reported data were thus from snowfall rather than rainfall. In total, 262.1 mm of rainfall was reported for Coles Lake, with maximum rainfall occurring in June (73.9 mm) and July (61.9 mm) (Figures 4-1 and 4-2). Figure 4-1. Daily rainfall and snowfall (water equivalent) at the UF weather station (October 2013_September 2014). Among meteorological parameters, snowfall and snow depth are perhaps the most difficult to measure because of high local variability (e.g., terrain), drifting, and melting between observations (Schmidlin, 1989). To record the amount of snowfall, an SR50 was 106 installed at the UF station next to Coles Lake. Unfortunately, due to the remoteness of the study area and limited budget, it was not possible to monitor the snowfall over the watershed scale. Overall, 1 cm of snow is assumed to correspond to 1 mm of water (Dubé and Rimouski, 2003). Therefore, the total snowfall at the UF weather station was determined to be 167.4 mm (water equivalent) with maximum SWE (56.5 mm) occurring in December based on the ultrasonic ranging sensor. The first snowfall of the 2013_2014 winter occurred near the beginning of October, and the last snowfall occurred at the end of April (Figure 4-2). Figure 4-2. Monthly rainfall and SWE pattern at the UF weather station (October 2013_September 2014). 4.1.1.2. Estimation of Lake Evaporation In general, the estimation of evaporation from lakes is not a simple matter, as there are a number of factors that can affect the evaporation rates, notably the climate and physiography of the water body and its surroundings (Finch and Calver, 2008). Overall, a wide variety of methods for estimating open water evaporation has been reported in different studies, each with different techniques depending on areas of research and types of water bodies. 107 The first step was to identify the ice-free period for the 2013_2014 hydrological year. To this end, Landsat images were downloaded (https://landsat.usgs.gov/, 2014) to determine the start and end time of freezing and melting at Coles Lake. Landsat is a series of Earthobserving satellite missions jointly managed by the United States (US) National Aeronautics and Space Administration (NASA) and the U.S. Geological Survey. Based on the Landsat images, Coles Lake began freezing at the beginning of November 2013 (Figure 4-3). Figure 4-3. Regional surface conditions including ice coverage at Coles Lake on 7 November 2014 (left) and 6 May 2014 (right). For both images, false-colour is used to restore the look of Landsat images from the Operational Land Imager on Landsat 8. The image (left) is relatively cloudy; however, it is still possible to identify the frozen nature of the lake, using the yellow colour of the lake (using false colour). The melting period began around 6 May 2014, and most of the ice had melted off Coles Lake by 15 May 2014, although the wider images indicate that most of the lakes nearby were still ice-covered (Figure 4-3). Based on the Landsat images, the ice-free period was defined from 1 October 2013 to 1 November 2013 and then from 15 May 2014 to 30 September 2014. The Penman-Monteith method (Zotarelli et al., 2010) was selected to estimate evaporation at Coles Lake during the ice-free season. The Penman-Monteith method was selected as it is a standard and accurate method to estimate daily potential evapotranspiration 108 approved by the FAO (Allen et al., 1998). Complete steps of the evaporation calculations are discussed in Appendix C. Lake Evaporation Results The average hourly lake evaporation is plotted during one year. The rate of evaporation began to increase at 07:00 local time while evaporation peaked at 14:00 hours, declining after that. Maximum monthly evaporation occurred in July (97.7 mm), and minimum monthly evaporation occurred in December (2.6 mm). This demonstrates a strong dependence on air temperature, with a total loss of 368.5 mm during the ice-free months (Figure 4-4). Figure 4-4. Monthly evaporation at Coles Lake for October 2013 to September 2014. Sublimation To assess the potential contribution of blowing snow sublimation to the water balance of Coles Lake, the Piektuk blowing snow model (Déry et al., 1998; Déry and Yau, 2001) was run over the frozen surface of the lake using observed meteorological conditions at the UF station. The model inputs include air temperature, RH, WS, atmospheric pressure and snow depth, and it integrates over time four prognostic equations to assess blowing snow 109 frequency and fluxes. The model establishes the occurrence of blowing snow when air temperatures are below freezing, snow is present, and wind speeds exceed a given threshold estimated from (Déry and Yau, 1999): Ut = Ut0 +0.0033 (T + 27.27)2 (4.2) Here T is in degrees Celsius and the minimum value of the threshold wind speed at 10 m above the surface, Ut0 , is ≈7 m s-1 at T = -27.27°C. According to Equation (4.2), resistance to transport increases both near 0°C and at temperatures below -27.27°C (Déry and Yau, 1999). At the Coles Lake meteorological station, however, wind speeds remained low and never exceeded the threshold for transport, suggesting that blowing snow sublimation does not contribute significantly to the water balance of Coles Lake and can be safely neglected in this study. In addition, there was no evidence of snowdrifts along the edges of the lake during the field campaign. 4.1.1.3. Streamflow Measurements at Coles Lake Discharge measurements were conducted at both the inflow and outflow hydrometric stations. Although careful efforts were made, it was not possible to employ the stagedischarge rating curve to estimate the streamflow due to beaver dams. Since these structures blocked water, the water level recorded by transducers in July, August and September levelled off. Consequently, water levels were raised in the upstream stage and therefore rendered the rating curve invalid. To further validate this result, the comparison was made with the Emile Creek station—located about 4 km downstream from Coles Lake— with the inflow and outflow estimates. According to the field measurements, the inflow and outflow values were higher than Emile Creek after mid-summer (July). This result is unrealistic, as Emile Creek formed as a 110 combination of several other streams, including flow from Coles Lake. As a consequence of the construction of beaver dams and remoteness of the site, only four measurements could be quantified during the field work. However, a minimum of 12 to 15 measurements is required to provide a robust stage-discharge relationship. In light of this challenge, an alternative method was employed to obtain representative streamflow data. Correlation Method To compute the amount of discharge, the streamflow data from Emile Creek were compared with two Coles Lake hydrometric stations. Four days of onsite discharge measurements of inflow and outflow, using the mid-point method, correlated highly with data for the same day from the Emile Creek station (Table 4-1). A regression allowed extrapolating the daily discharges (Qi and Qo) at each site for each day of the study period. Table 4-1. The t-values and P-values between Coles Lake and Emile Creek stations with the computed Pearson’s product-moment correlation coefficients. The t-value reflects the daily value of the‘t’ test statistic, n is the number of measurements, and the p-value reflects significance from 15 May to 15 July 2014. Stations Between inflow station and Emile Creek station Between outflow station and Emile Creek station n 4 t-value 9.39 p-value 0.011 Correlation 0.99 4 9.78 0.010 0.99 In general, the estimated values decrease from May to September (Figure 4-5). Prior to the construction of the beaver dams (15 May to 15 July), the stream rose rapidly in response to rainfall events and then declined at a moderate rate as water was shed from the wetland storage. However, after the beaver dams were built (after 15 July), rainfall had minimal effect on the streamflow. The pond behind the dam may have been at a low level, and overflow may have commenced once the water level topped the dam crest. In contrast, if the pond behind the dam was already at a high level, overflow would have occurred 111 immediately, having little effect on modulating channels flows (Woo and Waddington, 1990). Figure 4-5. Daily discharge comparison between the inflow, outflow and Emile Creek stations between 15 May and 30 September 2014. Points are manually measured vs. the daily reconstructed discharge at each location. Although the correlation was not the preferred method, it was the most appropriate approach for this situation. The following equation was used to compute the amount of runoff for each month: R= Q× s m (4.3) A where R denotes runoff, Q represents discharge (m3 s-1), s equals the number of seconds in a month, m is the number of days in a month, and A is the lake area (m2). The ArcGIS flow accumulation tool was used to compute the area at the inflow where the sub-basin fed the eastern end of the lake. Based on the ArcGIS analysis, the total area of Coles Lake is 1.8 km². Based on these calculations, the average annual runoff reaches 106.6 mm, and 198.5 mm for the inflow (Qi), and outflow (Qo), respectively (Figure 4-6). 112 Figure 4-6. Runoff at the inflow, outflow, and Emile Creek during the ice-free period in 2014. 4.1.1.4. Shallow Groundwater Results The hydraulic heads in the piezometers are measured along each transect, and complete details of hydraulic measurements are described in Section 2.3.4. P2-1 data are missing from 15 April to 13 June 2014, and P2-2 and P3-3 data are missing from 5 February to 13 June 2014 (Figure 4-7). Recharge, or downward flow conditions, are characterized by decreasing hydraulic head with depth. Discharge conditions are characterized by increasing hydraulic head with depth, here corresponding to a negative gradient. Therefore, flow was upward for the entire time series in Transect 1. In Transect 2, the flow changed from upward to downward at the beginning of August. For Transect 3, the flow changed from upward to downward in midFebruary and then flipped again in mid-May (Figure 4-8). Thus, groundwater flow directions were seasonally variable at these three transect locations. 113 Figure 4-7. Daily hydraulic heads in piezometers along each transect during 2014. Figure 4-8. Daily vertical gradients for all three transects in 2014. Negative gradients imply that groundwater flow is upward. Additionally, negative horizontal gradients, corresponding to groundwater flow from the land toward the lake, were reported for all three transects (Figure 4-9). 114 Figure 4-9. Daily horizontal gradients for all three transects in 2014. Negative gradients imply that groundwater flows laterally from the land to the lake. Based on the available K estimate, daily components of the flux in the vertical and horizontal directions were graphed (Figures 4-10 and 4-11). Figure 4-10. Daily vertical fluxes at three transects at Coles Lake in 2014. A positive flux indicates that groundwater flow is upward. The vertical gradient record at Transect 2 was short and the overall range of values rather small and inconsistent compared to the other two transects. 115 Figure 4-11. Daily horizontal fluxes at three transects at Coles Lake in 2014. A positive flux indicates that groundwater flow is towards the lake. The magnitude of the net flux at Transect 2 was positive (the square of the flux components are used in Equation (2.6)). However, the direction of the flux was toward the lake and downward at this location. This result is in contrast with the other two transects that have upward and lateral flow into the lake. The vertical gradient recorded at Transect 2 was brief, and the overall range of values was rather small and inconsistent compared to the other two transects. Therefore, the vertical flux at this location was disregarded. The net flux in Table 4-2 only reflects the horizontal flux. Table 4-2. Total groundwater fluxes for Transects 1, 2, and 3 (5 February 2014 to 30 September 2014). Location qx (m s-1) qz (m s-1) q (m s-1) Transect 1 Transect 2 Transect 3 1.29 ×10-11 2.87 ×10-11 6.55 ×10-12 3.47 ×10-10 0.00 5.45 ×10-11 3.47 ×10-10 2.87 ×10-11 9.95 ×10-11 The net average annual flows of groundwater to the lake for Transects 1 through 3 were 1.51 × 10-6 m3 s-1, 1.35 × 10-7 m3 s-1 and 2.74 × 10-7 m3 s-1, respectively (Table 4-3). Therefore, in total, groundwater inflow contributed 4.1 mm to the water balance of the lake 116 each month from the three transects (Flow = Q × (Number of seconds in a given month) / (an area that each transect covers)). Flows to and from the surface-water body were summed to calculate the net groundwater contribution to the entire lake. Table 4-3. Estimated outflow average (groundwater to the lake) per month in 2014. Transects q1 q2 q3 Total Average Annual Flux (q) (m s-1) 3.47 × 10-10 2.87 ×10-11 9.95 ×10-11 --- Area (A) (m2) 4342.21 4693.66 4982.91 14018.78 Water Flow (Q) (m3 s-1) 1.51× 10-6 1.35× 10-7 2.74× 10-7 1.92 ×10-6 Average Outflow per Month (mm) 2.23 0.67 1.17 4.07 Piezometers were installed on 5 February 2014; hence no data were available prior to this date. Groundwater fluxes were only estimated from 5 February 2014 to 30 September 2014. Furthermore, in the initial phase of this research, the ice-free period was identified by Landsat 8, and the groundwater fluxes were only considered during the ice-free period. Although this logical inference is based on sound assumptions, it may have underestimated the flux. Based on these assumptions, 22.6 mm (4.1 × 5.5 months) of groundwater flowed into the lake during the ice-free period. 4.1.2. Water Withdrawal (W) The BC Oil and Gas Commission (OGC) provided Quicksilver a maximum water withdrawal allocation of 38.5 mm from Coles Lake, although this value changes with each permit application based on water availability. The OGC permitted the maximum amount of water withdrawals during the high-flow period, at the beginning of the snowmelt season (May) when lakes typically achieve their highest water levels. The allocated permit decreased gradually towards spring, and the lowest permitted allocation was usually issued for March. However, for this study, the total amount of water withdrawal equals zero, as no 117 water was extracted during this study period due to reduced Quicksilver operations because of low natural gas prices (Quicksilver Resources Canada Inc., pers. comm., 2013). 4.1.3. Measured Water Balance Results Once each water balance component for the selected water year was computed for Coles Lake, it was possible to quantify the total water balance for Coles Lake (Table 4-4). The total computed amount of storage change, computed as the residual of the Coles Lake water balance equation (4.1), equalled -8.3 mm. A negative storage change identifies that 2013_2014 was a relatively dry water year compared to the previous year (see Appendix E for details). Additionally, the amount of precipitation was lower while evaporation was higher compared to the previous water year (2012_2013). Table 4-4. Final summary of the observed water balance components for Coles Lake (1 October 2013 to 30 September 2014). Water Balance Components Rainfall (R) Snowfall (S) Inflow (Qi) Outflow (Qo) *Evaporation (Ea) Groundwater (G) Minimum Water withdrawal (W) ∆S with no water withdrawal Observed value +262.1 mm +167.4 mm +106.6 mm -198.5 mm -368.5 mm +22.6 mm 0.0 mm -8.3 mm Error ± 1.5 mm ± 2.4 mm ± 0.9 mm ± 0.5 mm **± 1.1 mm ± 0.8 mm 0.0 mm ± 7.2 mm *This is only during the ice-free period. ** The error for evaporation is estimated based on the sum of error from temperature, WS, and RH. Although during the study period the total amount of water withdrawal from Coles Lake was zero (Table 4-4), it remained important to identify how water withdrawal may impact the water availability at Coles Lake. In general, water withdrawals may be significant if their rates exceed the replenishment of that water, especially within a shallow water body where water shortages with adverse environmental impacts remain a concern (Ridgway et al., 2016). 118 Identifying the exact impacts of water withdrawal is based on how much water may be extracted at different times of the year, the type of water supply being used (in this case, surface water), and the ability to replace withdrawn water (Ridgway et al., 2016). Water withdrawal from Coles Lake impacts the downstream flow. As a result, it can lead to alterations in wetlands, and affect the aquatic ecosystem health (including the habitat of plants and animals) of a water system (Ridgway et al., 2016). As the water availability is limited during wintertime, a minimum amount of water withdrawal is permitted for oil and gas activities (1 December _ 31 March). The water allocation can increase in spring during the freshet, however, if the ∆S < 0 then the allocation amount must be adjusted accordingly. Therefore, understanding of the hydrologic processes and the temporal and spatial exchange of GW-SW is the foundation of effective and sustainable water resources management. 4.2. Simulation of Coles Lake Hydrological Processes Using the MIKE SHE Model Knowing the complexity of wetlands in northeastern BC, and the extensive water demands by the oil and gas industry, it was essential to quantify the Coles Lake water balance in an integrated fashion. However, due to the remoteness of the area, conducting comprehensive fieldwork remained challenging. Thus, it was valuable to employ a distributed hydrological model to better understand hydrological processes in this region. To this end, the water balance tool within the MIKE SHE model was employed to simulate the Coles Lake hydrological components (i.e., for the lake itself) during the 2013_2014 period. By default, MIKE SHE simulates the water balance of the entire Coles Lake watershed (MIKE SHE User Guide, 2017). Therefore, the Coles Lake shapefile (polygon) was used to extract the water balance for the lake using the sub-catchment option in MIKE SHE 119 (Figure 4-12). First, the model was run for the Coles Lake watershed, and then the water balance utility—a post-processing tool—was used to generate water balance summaries for the lake itself. Figure 4-12. A: The Coles Lake shapefile (polygon). B: The Coles Lake watershed including Coles Lake bathymetry and mainstream. 120 In this study, monthly water balance data were extracted for Coles Lake using the following calculation: P ± ∆Canopy Storage ˗ AET ± ∆Snow Storage ˗ OL Flow to River(Lake) ± ∆OL Storage ˗ OL Boundary Outflow ± ∆Subsurface Storage ˗ Baseflow to River(Lake) + River(Lake) to Baseflow ˗ Subsurface Boundary Outflow ± Model Error (4.3) By default MIKE SHE uses Equation (4.3) to estimate the total water balance for Coles Lake. The term River (Lake) is used to represent Coles Lake because the lake forms as a continuation of the river network, extending through the cross-sections and the lake bathymetry (Figure 4-12). Further, flood codes are used to inundate floodplain areas based on the water level in the lake (MIKE SHE User Guide, 2017). In this situation, therefore, special considerations for water balance calculations are required in active flood code cells as described below (MIKE SHE User Guide, 2017). Overall, MIKE SHE interprets the water on the ground surface as a ponded storage, except in active flood code cells where the cell is flooded, and the water level is controlled by the water level in MIKE 11 (MIKE SHE User Guide, 2017). In general, when the groundwater table is at or above the land surface, water can exchange directly between ponded water and the SZ, and an UZ does not exist (MIKE SHE User Guide, 2017). Thus, if the land surface is an active flood code cell, the water enters into or is removed from the storage available for exchange with MIKE 11 (MIKE SHE User Guide, 2017). Specifically, if the ponded water level is above the ground level in a flood code cell, the water level is set equal to the corresponding water level in MIKE 11 at the beginning of the overland time step (MIKE SHE User Guide, 2017). Then, at the end of the overland flow time step, MIKE SHE computes the change in ponded water level and adds or subtracts this as lateral inflow to MIKE 11 over its next time step(s), covering the period of the MIKE SHE overland time step (MIKE SHE User Guide, 2017). 121 Thus, flow from the SZ is not directly added as lateral inflow, but it is one of the source/sink terms that contributes to the change in storage in flooded cells, which is then added to MIKE 11 as an inflow to the lake (MIKE SHE User Guide, 2017). Accordingly, the flood to lake and lake to flood together represent the net lateral inflow to MIKE 11 from active flood code cells. Summed together, they are the actual exchange between flood code areas and MIKE 11 (MIKE SHE User Guide, 2017). According to the MIKE SHE User Guide (2017), when the flood code is used the sign convention needs to be reconsidered as in Table 4-5. Note, however, that these water balance fluxes are not represented in the overall model water balance discussed below. These exchanges are represented in the detailed results for the individual components of the water balance (e.g., OL flow). Table 4-5. Overland flow water balance interpretation when the flood code is used (MIKE SHE User Guide, 2017). Water Balance Flux SZ to Flood Cell SZ to Flood Cell Flood to Lake (Inflow) Flood to Lake (Outflow) Description Sign Convention in the Water Balance Direct flow upwards from SZ to an active flood code cell. This is a positive upwards flow in the MIKE SHE results files. Inflow - negative (Note sign change in water balance definition) Direct flow downwards from an active flood code cell to SZ. Outflow - positive (Note sign change in water balance definition) Net lateral inflow exchange between active flood code cells and MIKE 11 nodes that are inside the current water balance area Inflow (negative) Outflow (positive) Net lateral inflow exchange between active flood code cells and MIKE 11 nodes that are outside the current water balance area This is always zero unless the water balance is being calculated for a subarea. 122 Inflow (negative) Outflow (positive) The monthly water balance of the Coles Lake watershed for the 2013_2014 water year was simulated in detail by MIKE SHE (Table 4-6 and Figure 4-13). Overall, the total water balance input for the watershed stemmed only from precipitation (i.e., NARR dataset). The other components of the water balance included fluxes out of the system (e.g., evaporation), which are always negative, or storage changes (e.g., snow storage, canopy storage), which can be positive or negative. The water balance results for Coles Lake itself were then extracted. Monthly total precipitation ranged from 19.2 mm in November to 122.7 mm in April totaling 588.0 mm over the duration of the 2013_2014 water year (Table 4-6). Evaporation (Ea in Table 4-6) induces the second component of water losses at Coles Lake. The actual evaporation was 494.9 mm year-1, ranging from 0.2 mm in December to 112.2 mm in July. A total of 76.8% of the evaporation occurred from May to August, while only 29.8% of the precipitation fell in that period. Therefore, only considering precipitation to and evaporation from the lake surface, there would be a net gain of 93.0 mm. However, there are other gains and losses of water to the lake as described below. In Table 4-6, OL boundary flow (in-out) is the net contribution of the OL flow from the area outside the lake (223.0 mm ˗ 7.4 mm = 215.6 mm) (DHI, pers. comm., 2017). This process occurs when the water depth on the boundary increases and the flow discharges across the boundary within the UZ. In this case, the boundary acts as an infinite source of water (MIKE SHE User Guide, 2017). In the SZ, the model simulates the contribution of 527.0 mm of subsurface inflow to Coles Lake from the edges, and 195.0 mm of water lost through the subsurface outflow, with a net gain of water to the lake from the subsurface 123 (527.0 mm – 195.0 mm = 332.0 mm). In fact, the groundwater inflow from outside the lake area to the lake itself is 332.0 mm. Table 4-6. Monthly simulated water balance of Coles Lake (October 2013_September 2014). Period Oct Nov Dec Jan Feb Mar April May June July Aug Sep Total P (mm) 43.2 19.2 22.5 63.5 28.1 64.6 122.7 43.8 54.4 36.1 40.8 49.1 588.0 Canopy storage change (mm) 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Ea (mm) -21.0 -2.3 -0.2 -2.0 -1.9 -10.4 -37.5 -86.0 -103.5 -112.2 -78.3 -39.5 -494.9 Snow storage change (mm) 0.0 31.2 56.6 34.3 13.5 15.3 -95.7 -55.2 0.0 0.0 0.0 0.0 0.0 Lake to OL flow (mm) -45.4 -27.2 -24.7 -21.7 -17.8 -13.4 -160.5 -204.0 -67.2 -4.3 11.6 -22.9 -597.5 0.0 0.0 0.0 0.0 0.0 0.0 70.5 140.5 11.5 0.5 0.0 0.0 223.0 -0.6 0.0 0.0 0.0 0.0 0.0 -1.0 -4.7 -0.5 -0.4 0.0 -0.2 -7.4 42.2 39.7 40.5 40.1 35.9 39.4 39.0 60.5 52.9 49.9 46.4 40.5 527.0 -14.3 -16.5 -19.1 -20.5 -19.3 -22.0 -22.2 -13.1 -7.9 -11.1 -13.6 -15.6 -195.0 -1.1 -1.2 0.0 0.0 0.0 0.0 8.0 -7.0 2.1 -3.3 0.0 0.5 2.0 0.1 -1.0 -1.0 -0.9 -0.9 -1.8 21.2 -4.9 -5.0 -5.9 -3.4 -0.5 -4.0 0.0 0.0 0.0 0.0 0.0 1.7 0.0 0.0 0.3 0.0 0.0 0.0 2.0 Lake to Base flow (mm) -2.0 -2.2 -2.4 -2.4 -2.2 -2.6 -3.2 -2.0 -0.6 -1.0 -1.5 -1.9 -24.0 Total model error (mm) 0.0 0.0 0.0 -0.1 -0.1 -0.1 -0.3 0.6 0.0 0.2 0.1 -0.1 0.2 OL Boundary flow In (mm) OL Boundary flow Out (mm) Subsurface Boundary Inflow (mm) Subsurface Boundary Outflow (mm) OL-storage change (mm) Subsurface storage change (mm) Base flow to lake (mm) P=Precipitation; Ea=Evaporation; OL=Overland Flow. *Ea is the evaporation from the lake’s surface. 124 Figure 4-13. The 2013-2014 water balance terms for Coles Lake, simulated by the MIKE SHE model. Therefore, the net amount of water entering the lake from both OL flow and subsurface flow from the area outside the lake is 215.6 mm + 332.0 mm = 547.6 mm. Once these boundary flows (at the surface or in the subsurface) enter the lake sub-catchment, numerous exchanges can take place. Firstly, there is a net loss of water (24.0 mm ˗ 2.0 mm = 22.0 mm) to the groundwater system at the H points (DHI, pers. comm., 2017). Additionally, there is a series of exchanges between the SZ and the surface zone (Figure 4-13). Of the subsurface 125 inflow (527.0 mm) some is transferred to the surface (497.0 mm), some is transferred back from the surface to the SZ (141.0 mm), and some leaves the lake in the subsurface (195.0 mm) – the net loss is 24.0 mm, which is the lake to baseflow loss. These transfers are not shown in the overall water balance, but are available in the detailed water balance results and are illustrated in Figure 4-13. The total lake to OL flow (597.5 mm) is approximately the same as the net boundary flow into the lake (547.6 mm), accounting for evapotranspiration loss (494.9 mm), suggesting that the lake storage was maintained over the simulation period. The simulated streamflow at the outlet of the lake could not be determined from the water balance. Describing the recharge process between the lake and SZ is complicated due to several factors. For instance, the water table rises as water enters the SZ and this, in turn, affects flow conditions in the storage (MIKE SHE User Guide, 2017). Furthermore, the major challenge in describing the linkage between the surface storage and SZ zones arises from the fact that the two components are explicitly linked and run in parallel in MIKE SHE (MIKE SHE User Guide, 2017). The maximum amount of simulated overland flow occurred in May 2014 during the melting period, and it decreased into September. A similar pattern was found for the Inlet, Outlet, and Emile Creek stations during fieldwork (see Section 4.1.1.3). During this time, some of the overland flow may evaporate, some may infiltrate, and/or some may leave via outlets (DHI, pers. comm., 2017). In addition, some of the overland flow may be retained as overland storage. In this study, overland storage changes approached zero during the freezing period, and its difference between April and May marked the start of the melting phase. 126 Additionally, the MIKE SHE model reports snow storage values, indicating a negative inflow of precipitation and a positive outflow of melting (MIKE SHE User Guide, 2017). In this study, the maximum amount of simulated snow storage change occurred in April during the melting period with the overall snow storage equaling zero. Finally, the water balance tool reports the amount of error that occurs during a simulation due to convergence issues (DHI, 2007). For Coles Lake, the reported error during the simulation period was quite small (0.2 mm). Despite this, even a minor amount could create uncertainties in the water balance results. Factors that play a role in model uncertainties are discussed in detail in Section 7.3.2. Further, calibration error, the absence of long-term measured data for calibration purposes, and uncertainty in the parameters with no calibration data (e.g., ET, the UZ parameters, etc.) may further influence the simulated water balance results. 4.3. Comparison of Measured and Simulated Results The Coles Lake water balance results simulated by MIKE SHE were compared to the measured results for the 2013_2014 water year (Table 4-7). However, it is important to note that two different precipitation datasets were used in the observation (i.e., UF station) and simulation (i.e., NARR) for 2013_2014. Therefore, it is essential to keep this in mind when comparing observed and simulated results. The ratio of each hydrological component to annual precipitation allows easier comparisons between observed and simulated terms (Table 4-7). To compare the results of this study, it was required to define a simplified equation for observed and simulated values (Equation 4.4): ∆S = P + G ˗ Q ˗ Ea 127 (4.4) Table 4-7. Comparison between measured and simulated water balance for Coles Lake (2013_2014). Components Measured +429.5 mm Ratio to Measured P +1.00 Precipitation (P) Simulated +588.0 mm Ratio to Simulated P +1.00 *Evaporation (Ea) -368.5 mm -0.86 -416.5 mm -0.71 ***Net Runoff (Q) -91.9 mm -0.21 -254.1 mm -0.43 ****Groundwater (G) +22.6 mm +0.05 -22.0 mm -0.04 ∆S -8.3 mm -0.02 -104.6 mm -0.18 *This is calculated for one hydrological year based on the observation value during the ice-free period. The total simulated Ea by MIKE SHE is -494.9 mm. **This value is measured the Coles Lake station (Net Runoff = Inflow ˗ Outflow). ***This is taking the total OL flow from the lake and subtracting off the OL boundary flow, the subsurface boundary flow and the exchange with SZ (597.5 mm – 7.4 mm– 195.0 mm– 141.0 mm = 254.1 mm). **** The net exchange between the lake and groundwater at the H points. Measured precipitation (429.5 mm) was less than the NARR value used as inputs to the model (588.0 mm). Although a significant correlation of 0.98 was computed between the NARR and UF dataset (Coles Lake), the reported precipitation by NARR was higher compared to UF for the 2013_2014 water year. A total precipitation of 515.3 mm fell at the FNA station for the 2013_2014 water year, with mean annual precipitation of 451.7 mm for 1981_2010 (Environment and Climate Change Canada, 2017). Although measured precipitation at UF was relatively comparable to the reported mean annual precipitation at the FNA station, unfortunately, a long-term precipitation dataset was not available specific to UF to confirm if the precipitation was underestimated for this site. Therefore, to make the comparison between observed and simulated values possible, the ratio of each hydrological component was reported relative to annual precipitation. In addition, 368.5 mm of evaporation was computed over the lake, which is smaller compared to the simulated value (416.5 mm) during the ice-free period. The evaporation had similar patterns for both measured and simulated values, as the maximum evaporation 128 occurred in July (observed = 97.7 mm, simulated = 122.2 mm, respectively) and the minimum occurred in December (observed = 2.6 mm, simulated = 0.2 mm, respectively). Simulated net runoff is -254.1 mm, although higher than observed runoff (-91.9 mm), it had a seasonal pattern. The maximum amount of runoff occurred in April during the melting period, and it decreased towards October. Similar outcomes were observed during fieldwork at Coles Lake. Differences between simulated and measured net runoff can be explained by the presence of beaver dams that affected the inflow and outflow of Coles Lake during the 2014 summer field campaign. Beaver dams were built around 15 July 2014 leading to upstream water storage and overflow may only have commenced if the water level crested the dam. The representation of beaver dams is not an option in the MIKE SHE modeling framework leading to excess efflux of water from the lake. Additionally, simulated net runoff was higher compared to observations explained by differences between the simulated and measured storage changes (∆S). Additionally, measurements revealed that a small amount of groundwater (22.6 mm) contributed to the Coles Lake water balance due to the region’s low hydraulic conductivity and permeability. However, MIKE SHE simulated a small net loss from the lake to the groundwater (22.0 mm). 4.4. Discussion Due to the lack of historical, in situ data, this section emphasizes the patterns of the hydrological processes rather than the values. Therefore, seasonal variability of hydrological patterns of Coles Lake is discussed below. This approach was also chosen because there were numerous water exchanges that took place in the simulation that were not measured during the fieldwork. Detailed information of the fluxes and patterns provided by the MIKE 129 SHE model can be employed to design a sustainable water management plan for northeastern BC. During the snow accumulation period (November through April), the simulated evaporation rate is low and snow storage change values remain positive. However, during the snow melting period (22 April to 6 May 2014), evaporation increased and snow storage change shifted from positive to negative. Evaporation peaked in summer and the snow storage change equaled zero as the snow had already melted (Table 4-8). Similar trends were observed for snow accumulation and evaporation during the field campaign (2013_2014). The first snowfall occurred near the beginning of October, and the last snowfall occurred at the end of April. Maximum monthly evaporation occurred in July, and minimum monthly evaporation occurred in December. MIKE SHE reports that the flow rate is at a maximum during the melting period. Comparable trends were also identified during the field work, as the maximum flow rate occurs during the melting period and it then falls during the summer. Hence, knowing the maximum flow rate occurs at the beginning of snowmelt (May), I recommend most water withdrawal activities take place during this time frame and decrease in late summer as the water level diminishes. The contribution of the lake to OL flow is small (near zero) during the freezing period, and a similar pattern was observed during the field campaign. This is a realistic result as Coles Lake is a shallow water body and it freezes during winter (freezing period). In contrast, during the melting period, the contribution of the lake to OL flow reaches its highest point, becoming almost equal to OL boundary inflow (Table 4-8). Most of the contribution to the lake is from surface water, rather than groundwater during the melting period (Table 4-8). 130 However, during the summer, as the evaporation rate increases, the contribution of the lake to OL decreases gradually and the amount of lake to OL flow almost equals to the subsurface boundary inflow. Groundwater, rather than surface water, contributes most of the flow of water to the lake during the summer (Table 4-8). This is further confirmed by the estimated groundwater results during the field campaign, as groundwater flowed into the lake during the ice-free period. Table 4-8. Seasonal variability of simulated hydrological components of Coles Lake using MIKE SHE (2013_2014). Period Freezing Period (Winter) Melting Period (Spring) Ice-free Period (Summer) Ea (mm) Low (̴ 0) Medium High Snow storage change (mm) Before 22 April 2014 > 0 22 April to 6 May 2014 < 0 After 6 May 2014 = 0 Low High Medium Lake to OL flow ≃ OL Boundary inflow (Contribution is mostly from Surface) OL Boundary Inflow > OL Boundary Outflow Lake to OL flow ≃ Subsurface Boundary inflow (Contribution is mostly from Subsurface) OL Boundary Inflow < OL Boundary Outflow Lake to OL flow (mm) OL Boundary flow In (mm) Low (̴ 0) OL Boundary flow Out (mm) Low (̴ 0) OL Boundary Outflow < OL Boundary Inflow OL Boundary Inflow < OL Boundary Outflow Subsurface Boundary Inflow (mm) Subsurface Boundary Outflow (mm) OL-storage change (mm) Low (not variable) Low (̴ 0) Subsurface Boundary Inflow = Subsurface Boundary Outflow Subsurface Boundary Inflow = Subsurface Boundary Outflow Medium (start to increase) Subsurface Boundary Inflow > Subsurface Boundary Outflow Subsurface Boundary Inflow < Subsurface Boundary Outflow Low (less compared to spring) Subsurface storage change (mm) Very low Medium Low Base flow to the Lake (mm) Low (̴ 0) Low Base flow to Lake > Lake to Base flow Lake to Base flow (mm) Low (̴ 0) Low Base flow to Lake > Lake to Base flow Low (not variable) 131 Beaver dams were built at both the inflow and outflow of Coles Lake after 15 July. A newly built beaver dam can cause shifts in the capture zone of a wetland pond by increasing hydraulic head behind beaver dams (Feiner and Lowry, 2015). In addition, beaver dams can change the groundwater flux, and the extent of both the capture and discharge zones of the wetland area (Feiner and Lowry, 2015). Although the results from MIKE SHE indicate a local contribution of groundwater to the lake, these results may not be applicable regionally. The clay unit underneath the peat may disconnect this site from regional groundwater flow, and as a result, where there are beaver dams, there are minimal changes to regional groundwater flow paths. However, to examine the extent to which groundwater flow through a peat wetland is affected by a beaver dam, extensive field work is required to identify the local hydrostratigraphy and the connectivity of the wetland to regional groundwater flow paths (Feiner and Lowry, 2015). In addition, a low rate and steady flow from the subsurface boundary inflow/outflow occurs during the freezing period. The amount of subsurface boundary inflow nears the outflow during the melting period, while the subsurface boundary inflow exceeds the subsurface outflow during summer. Thus, groundwater contributes to the lake when the surface water is at the lowest level. These results were further validated during the field campaign using nested piezometers installed around Coles Lake. Further, OL storage change and subsurface storage change approach zero during the freezing period while it increases during the melting period. The simulated OL and subsurface storage changes decrease during summer. Additionally, the amount of baseflow to lake and lake to baseflow nears zero during the freezing period. Although the baseflow to lake and lake to baseflow increase during the melting period, they remain low, as surface 132 water, rather than groundwater, contributes most of the flow from the lake. These patterns shift during summer, as the amount of base flow to the lake is greater compared to the lake to the base flow. 133 CHAPTER 5: INFLUENCE OF FOREST CANOPY ON WATER BALANCE AT THE COLES LAKE WATERSHED Forest vegetation plays a vital role in water availability for watersheds through the interception of precipitation (i.e., rain and snow), transpiration, and evaporation of intercepted water (Winkler et al., 2010a). In fact, when precipitation is intercepted by the forest canopy, and sublimation and evaporation occur, some amount of water returns to the atmosphere, while some becomes surface runoff or just water infiltration to soils (Pike et al., 2010; Winkler et al., 2010a). The amount of precipitation intercepted depends on canopy characteristics (Winkler et al., 2010a). Hence, an increase in forest canopy is anticipated to result in an increase in precipitation interception and therefore a decrease in surface runoff. The results of the Coles Lake water balance (i.e., the lake itself) are discussed in Chapter 4, while the role of vegetation cover was not elaborated upon. In this chapter, the role of vegetation canopies at the watershed scale are explored in detail. For this purpose, the total amount of precipitation (rainfall and snowfall) was quantified within the Coles Lake boreal watershed for three vegetation canopies of open, mixed, and closed. These three sites are located on relatively flat terrain, so a significant difference in air temperature, elevation or slope was not expected. In addition, the close geographic proximity of the sites implies that they are influenced by the same synoptic climate regimes (Figure 2-19). However, the selected sites are different in their canopy characteristics and ground cover. Therefore, forest canopy was considered one of the main factors controlling precipitation within the watershed. Fieldwork procedures were discussed earlier in Chapter 2, Sections 2.3.2 and 2.3.3. Therefore, the following section focuses on measured precipitation under the different vegetation canopies. Additionally, MIKE SHE was employed to further examine the effect of vegetation canopies on surface runoff within the Coles Lake watershed. 134 5.1. Rainfall Analysis To record the total amount of rain, four tipping bucket gauges were deployed at open (two stations), mixed, and closed canopies within the Coles Lake watershed. The UF and open stations were both located at the open sites. At the FNA station (Figure 2-1), the rainfall period typically occurred between the beginning of June to the end of October. Tipping bucket gauges do not function well with frozen precipitation (NEWA, 2016), unless they are insulated and heated. Thus, the analysis was carried out only during the campaign from 1 to 30 October 2013, and 1 June to 30 September 2014 to ensure representative results. Daily rainfall distributions for each canopy are shown in Figure 5-1. Overall, the total amount of rainfall decreased from June to October (Table 5-1). Although the overall patterns are similar, the total amount of rainfall differs between these stations. The highest amounts of total rainfall were recorded at the UF, mixed, open, and closed canopy sites, respectively. A monthly maximum in rainfall occurred in June for the UF, open and mixed canopy sites, while the monthly maximum in rainfall occurred in July at the closed canopy site. As mentioned in Chapter 2, two rainfall gauges — UF and open — were installed in the open canopy. The UF rain gauge was located relatively close to the mixed station (0.8 km) while the other rainfall gauge was installed farther relative to the mixed station (8.6 km). The UF rain gauge received 14.3 mm more rainfall compared to the gauge in the mixed canopy, while the other open gauge recorded 21.6 mm less rainfall than the mixed canopy. This variability amongst the rainfall gauges within the same canopy reveals how precipitation may vary considerably even over relatively short distances. 135 Figure 5-1. Daily rainfall at the Coles Lake watershed (October 2013_September 2014). Table 5-1. Monthly rainfall at each site (October 2013_September 2014). Total Rain (mm) UF (Open) Oct 2013 32.8 June 2014 73.9 July 2014 62.0 Aug 2014 38.6 Sep 2014 36.6 Total (mm) 243.9 Std. Error (mm) ± 1.5 Open 27.8 59.8 59.6 33.0 27.8 208.0 ± 1.2 Mixed Closed 30.4 30.0 64.8 29.8 60.2 41.4 36.6 37.6 37.6 32.2 229.6 171.0 ± 1.4 ± 1.1 Additionally, the linear model analysis (ANOVA) was conducted to explore differences in monthly rainfall among all sites – open, mixed, closed, and UF (open). The ANOVA test was conducted using R, with all sites compared to the closed canopy site, as the analysis demonstrated that closed sites received less rainfall compared to the other sites. Overall, there were no significant differences between sites in monthly rainfall (Table 5-2). However, more rainfall occurred in June and July compared to the other months across all sites relative to the closed site. Thus, to examine differences in rainfall during the rainy season, the analysis is restricted to the period of June and July (Table 5-3). Table 5-2. Rainfall linear model analysis results between open, mixed, and UF vs. closed canopy for June, July, August, September, and October (2013_2014). Closed vs. Mixed Closed vs. Open Closed vs. UF Estimate (mm) Std. Error (mm) t-value P-value 11.72 8.40 14.57 9.25 9.25 9.25 1.27 0.91 1.56 0.22 0.38 0.14 136 Table 5-3. Rainfall linear model analysis results between open, mixed, and closed canopies for June and July (2013_2014). Closed vs. Mixed Closed vs. Open Closed vs. UF Estimate (mm) 26.90 24.10 32.35 Std. Error (mm) 6.11 6.11 6.11 t-value 4.41 3.95 5.29 P-value 0.011 0.016 0.006 5.2. Snowfall Analysis The Coles Lake watershed is located within the Boreal White and Black Spruce zone, and only a portion of the snowfall reaches the ground through gaps in its forest canopy. The remaining snow is intercepted by vegetation or other surfaces and consequently drops or slides from these surfaces onto the ground or is lost through evaporation or sublimation (Winkler et al., 2010a). The amount of snow captured varies between canopy types (Buttle et al., 2000). Consequently, forest canopies influence the variability of SWE on the ground (Buttle et al., 2000). Snow reaching the ground depends on several factors, although the main factors identified are WS, the time between precipitation events, the intensity and duration of precipitation events, and the type of vegetation present (Pike et al., 2010). Once snow reaches the ground, it is referred to as “snow cover.” The snow cover not only varies between stands of different species, stem distribution and canopy density (Winkler and Moore, 2006), but also varies with distance from individual trees (Pomeroy and Goodison 1997; Faria et al., 2000). Forest cover affects both snow accumulation and the rate of snowpack interception. Therefore, it is expected that a decrease in forest cover—from closed to mixed and mixed to open canopy— will increase snow accumulation and accelerate ablation. That is, the amount of forest canopy is inversely related to SWE on the ground; however, the nature of this relationship changes with climatic conditions (Metcalfe and Buttle, 1998). The first round of snow surveys was conducted on 4 February 2014 (Table 5-4). The highest average snow 137 depth was observed within the mixed canopy, followed by closed and open vegetation. Similar results were observed for the remainder of the snow surveys. Table 5-4. Comparison of snow depth and SWE at the three different canopies (4 February 2014). Canopy Number of Samples Open Mixed Closed 8 8 8 Min Snow Depth (cm) 10 25 26 Max Snow Depth (cm) 52 85 67 Average Snow Depth (cm) 30 49 46 Coefficient of Variation of Snow Depth (%) 32 19 16 Average Snow Water Equivalent (SWE) (cm) 13 55 37 Standard Error (cm) 1.2 1.2 1.0 Overall, the snow surveys exhibited differences in snow depth among open, mixed, and closed vegetation canopies (Figure 5-2). An ANOVA confirmed that canopy type affected snow depth (ANOVA P-value = 0.02). Figure 5-2. Comparison of snow depth between open, mixed and closed vegetation canopies (4 February to 28 April 2014). Dots outside of the box and whiskers are outliers. Maximum and minimum values are shown at the end of each whisker. The upper and lower quartiles are the ends of the box, and the median is the horizontal line in the centre of the box. 138 Additionally, a post-hoc analysis was employed to further explore the pair-wise differences in snow depth among mixed, open and closed canopies (Table 5-5). There was no significant difference in snow depth between mixed and closed canopies. However, there was a significant difference between open and mixed and open and closed canopies such that snow depths were greater in closed and mixed canopies compared to open ones. Table 5-5. Post-hoc analysis results of snow depth between open, mixed, and closed canopies (4 February to 28 April 2014). Stations Estimate (cm) Std. Error (cm) t-value P-value Open vs. mixed -30.43 4.89 -6.22 < 0.001 Open vs. closed -24.92 4.89 -5.09 < 0.001 Mixed vs. closed 5.51 4.89 1.13 0.51 The frequency histograms of snow depth at each site were compared to the values recorded by the SR50 (i.e., installed at the UF station) and displays the range of snow depth values at each canopy (Figure 5-3). Perhaps surprisingly, the histograms indicate a better match between the SR50 data (red vertical lines in Figure 5-3) with the survey mean values (black vertical lines) for the closed canopy, followed by the mixed and open vegetation (Table 5-6). Of note, there was more snow in the mixed and closed canopies compared to the open site. In addition, the rainfall analysis demonstrated that closed sites received less rain compared to the other sites. Table 5-6. ANOVA analysis results of snow depth between SR50, closed, mixed, and open canopies (4 February to 28 April 2014). Stations t-value P-value Significant F SR 50 vs. Closed 1.47 0.17 0.48 SR 50 vs. Mixed 0.49 0.64 0.59 SR 50 vs. Open 4.85 < 0.001 0.85 139 Figure 5-3. Snow depth frequency histograms from snow surveys conducted at the Coles Lake watershed (4 February to 28 April 2014). The survey means are represented by a vertical black line, while the fixedpoint depth measurements (SR50) are marked with a vertical red line. 140 5.3. Influence of Forest Cover on the Coles Lake Water Balance Modelling Vegetation parameters are one of the main inputs in hydrological distributed modelling. Vegetation exerts considerable effects on runoff through features such as albedo and interception (Eckhardt et al., 2003), and consequently, plays an important role in the water balance cycle. Growth in forest cover usually results in an increase of ET, and forested catchments have higher ET compared to grass-covered catchments (Zhang et al., 2001). Increasing ET reduces the amount of overland flow (Bosch and Hewlett, 1982). The role of vegetation cover at the Coles Lake watershed was examined by replacing the actual canopy with different land covers and adjusting their LAI values accordingly within MIKE SHE to reflect the open, mixed, and closed canopy. To this end, 13 land cover classes were defined for the Coles Lake watershed (see Section 3.4.1 for details). Further, these 13 land classes were merged into three main categories – open, mixed, and closed canopies – using ArcGIS software, and the BC LCS coding in the VRI dataset. These three categories were represented on the VRI map (Figure 2-15). Then, the average of the existing LAI of 1.4, 1.7, and 2.1 was computed for the open, mixed, and closed canopies, respectively. In these simulations, only the LAI values were altered to represent the actual, open, mixed, and closed canopies and other parameters remained the same during the model simulations. The actual canopy refers to the natural vegetation cover of the Coles Lake watershed containing open, mixed, and closed canopies together (Table 2-2) as was modelled in Chapter 4. The details of the model setup were already discussed in Chapter 3; hence, in this section, only the outcomes of various canopies on the simulation results are discussed. To achieve a dynamically stable state, the model was run from 1979 to 2014; however, only the results of 2013_2014 for the Coles Lake watershed are discussed and evaluated. 141 The model results for the different vegetation covers were compared by assessing the simulated overland flow calibration results. The statistics for the open, mixed and closed canopy were similar compared to the actual canopy (Table 5-7). Low ME, MAE, and RMSE and high correlation were reported for all land covers, with a negative NSE. Overall, for the 2013_2014 hydrological year, the peak flow for the simulation with the actual vegetation cover matched well with that for the open and mixed canopy cases; the peak flow for the open canopy was higher compared to that for the actual canopy, and the peak flow for the actual canopy was slightly higher compared to that for the mixed canopy (Figure 5-4). In contrast, the peak flow for the closed canopy simulation decreased substantially compared to the actual vegetation canopy. Specifically, when the land cover distribution in MIKE SHE was altered to fully open or mixed canopy, the water balance fluxes changed slightly compared to the actual canopy distribution. However, when it was assumed the watershed was a fully closed canopy, the water balance fluxes changed significantly (Table 5-8). These results were reasonable considering the actual land cover contains 43.3% of open and 51.4% of mixed canopy, with a very small portion (5.3%) of the watershed covered by the closed type of canopy. Table 5-7. Results of the MIKE 11/SHE channel flow for actual, open, mixed and closed vegetation cover relative to observed runoff. The results are computed based on daily values over the simulation period of 2013_2014. Vegetation Cover ME (mm) MAE (mm) RMSE (mm) Actual Cover Open Cover Mixed Cover Closed Cover 1.29 1.21 1.23 1.35 1.29 1.21 1.23 1.35 1.73 1.45 1.49 1.66 STDEV Residuals (mm) 1.15 0.80 0.83 0.97 NSE R (Correlation) P-value -0.43 -0.79 -0.89 -1.34 0.85 0.84 0.83 0.83 5.22 × 10-44 7.18 × 10-45 3.94 × 10-42 9.57 × 10-46 The water balance results for the entire watershed were also compared. Overall, as the density of the vegetation cover increased from open to mixed to closed, the simulated ET 142 increased (Table 5-8). This caused a decrease in overland flow, while the overland storage, overland boundary outflow, and subsurface storage change remained almost the same. According to the ANOVA analysis, the differences of ET and overland flow between open, mixed, closed, and actual are significant (p-value between groups = 0.00). Figure 5-4. Simulated runoff for actual land cover vs. open, mixed, and closed (October 2013_September 2014). Table 5-8. Simulated water balance for various vegetation covers at the Coles Lake watershed (2013_2014). Although the model reports small positive and negative values for canopy and snow storage change, their total sum equals zero. Vegetation Cover P (mm) Canopy storage change (mm) ET (mm) Snow storage change (mm) OL flow to river (mm) OL storage change (mm) OL B outflow (mm) Base flow to river (mm) River to Base flow (mm) Subsurface storage change (mm) Total model error (mm) Actual Open Mixed Closed 588 588 588 588 0 0 0 0 -521 -515 -529 -558 0 0 0 0 -53 -59 -46 -11 -8 -9 -8 -4 -3 -4 -3 -1 0 0 0 0 0 0 0 0 -52 -52 -51 -50 3 3 3 3 P=Precipitation; ET= Evapotranspiration; OL=Overland Flow; B=Boundary. 143 5.3.1. Discussion Numerous studies have demonstrated the impact of vegetation on the water balance (Spittlehouse and Winkler 2002; Winkler et al., 2010a; Winkler et al., 2010b). Winkler et al. (2010a) found that altering forest vegetation affected the water balance at the watershed scale. Specifically, Winkler et al. (2010a) identified that forest vegetation affects the amount of water available for streamflow through the interception of precipitation, the evaporation of intercepted water, as well as transpiration. Results from this research also indicated that altering forest vegetation (i.e., altering the LAI value) directly affected the amount of ET and overland flow at each site (i.e., open, mixed, and closed). Spittlehouse and Winkler (2002) reported that forest vegetation plays a significant role in the water balance by reducing the rate of snowmelt, which supports the results obtained from fieldwork at the Coles Lake watershed. The rate of snowmelt was lower at the closed and mixed canopy compared to the open canopy. Figure 5-3 illustrates that the snow at the open site melted by 15 April 2014, while a significant amount of snow remained at the mixed and closed canopies on 28 April 2014. The Coles Lake watershed is located in a boreal forest in northeastern BC where winters are cold, and the existence of dense forest stands is uncommon because of the short growing season (i.e., only 5.3% of the Coles Lake watershed is covered with a dense canopy). However, where forest stands exist, the snow interception and ET are greater compared to the open canopies. Therefore, runoff from forested sites (closed and mixed) is less than compared to the open sites (Table 5-8). In addition, during the fieldwork, the average of snow depth and SWE was lower at the open canopy compared to the mixed and closed canopies. However, the average of rainfall was higher at the open canopy relative to the 144 mixed and closed canopy, which contribute to the total amount of runoff. In particular, ET in the dense forest is higher in comparison with non-forested areas at the Coles Lake watershed. A statistical analysis confirmed that canopy type directly affected the amount of ET (P-value < 0.05) This research did not measure the influence of wind activity (apart from potential snow drift sublimation at the UF weather station), which includes wind speed and direction. However, the results demonstrated that the patterns of precipitation in the form of snow were not as hypothesized, and the snowfall patterns were also dissimilar to the patterns of the observed rainfall. It is possible that wind activity may be present as an additional factor. The first component of the atmospheric hydrological cycle is intercepted precipitation, which is stored in the canopy, but can be released directly back to the atmosphere through canopy evaporation and sublimation (Edwards et al., 1983). High wind speeds impact significantly aerodynamically 'rough' canopies, and the loss of intercepted moisture may be rapid. In areas where the canopy is often wet its moisture is readily available thus promoting evaporation of intercepted water and recycling a significant proportion of the total rainfall (Edwards et al., 1983). Therefore, I recommend further investigation into the role of wind in this area, as wind may play a meaningful role in the factors affecting water balance estimates. In addition, it would be useful to examine whether shade and shelter in forested areas would result in wind having less influence on precipitation. Future snow surveys where analysis is evenly distributed across the entirety of the watershed are also recommended, as the current research is concentrated in northern areas of the watershed, exclusively around the Coles Lake area. 145 CHAPTER 6: THE WATER BALANCE OF TWO NORTHEASTERN BOREAL WATERSHEDS USING MIKE SHE The physiographic and climatic characteristics of a watershed largely influence its hydrological response and water balance components. Therefore, prior to comparing the water balances of the Coles Lake and Tsea Lake watersheds, it was important to comprehend the range of their physical and climatic characteristics. Thus, a brief description of the physiographic and climatic characteristics of the Coles Lake and Tsea Lake watersheds is described below, prior to summarizing their water balances. 6.1. Characteristics of the Coles Lake and Tsea Lake Watersheds The Coles Lake and Tsea Lake watersheds are located near each other (<40 km) in the moist and cool Boreal White and Black Spruce zone, characterized by wetland complexes of discontinuous permafrost (Johnson, 2012). Their close geographic proximity in the Liard River Basin implies that they are generally under the influence of the same synoptic regimes. However, each watershed has distinct physiographic and climatic characteristics (Das and Saikia, 2012). Thus, even small differences in their characteristics may affect their hydrological responses (Mohamoud, 2004). 6.1.1. Physiographic Characteristics The watershed physiography plays an essential role in the water balance. One of the main surface factors is the type of soil that influences the infiltration rate, the retention capacity, and the runoff coefficient. Northeastern BC is characterized by a thick blanket of clay-rich till deposited by the Laurentide Ice Sheet, overlying the bedrock throughout much of the area (Smith et al., 2005). The prevailing local climate acts upon this geologic setting, determining the quantity and timing of precipitation that falls into the watershed. As a consequence, it influences the vegetation density and patterns. 146 Thus, the type and density of vegetation cover directly determine the quantity of water intercepted and retained by the soil. The Coles Lake watershed consists of 43.3% of open, 51.4% of mixed, and 5.3% of a closed canopy. The Tsea Lake watershed contains 17.0% of open, 56.0% of mixed, and 27.0% of a closed canopy. In addition, the LAI values for the Coles Lake watershed range from 0 to 5.4 while at the Tsea Lake watershed they range from 0 to 9.1 (Table 3-4 and Table 3-5). 6.1.2. Climatic Characteristics The climate dataset is one of the main inputs into the MIKE SHE hydrological model. Among the climate parameters, precipitation and PET are two of its primary inputs. The comparison between precipitation and PET for the Coles Lake and Tsea Lake watersheds is briefly discussed in this section. The historical precipitation data were acquired from the NARR dataset for both the Coles Lake and Tsea Lake watersheds. The precipitation increased for both watersheds over the period of this study (October 1979 to September 2014). The Coles Lake and Tsea Lake watersheds each received 0.0 mm day-1 of minimum daily precipitation, and both watersheds attained 49.0 mm day-1 and 47.6 mm day-1 maximum daily precipitation, respectively. Average daily (1.4 mm day-1) and standard deviation (2.8 mm day-1) in precipitation matched for both watersheds. The overall daily precipitation was similar for both watersheds. The highest annual precipitation (852.4 mm and 849.0 mm) occurred in 2007 (water year), and the lowest annual precipitation (338.0 mm and 342.6 mm) occurred in 1980 (water year), respectively, for the Coles Lake and Tsea Lake watersheds (Table 6-1). The overall mean annual precipitation was also very similar (494.4 mm and 496.6 mm) for both watersheds. 147 Table 6-1.Comparison of precipitation at the Coles Lake and Tsea Lake watersheds (1979_2014). Watershed Maximum (mm) Minimum (mm) Mean (mm) Std. Deviation (mm) Coles Lake Tsea Lake 852.4 849.0 338.0 342.6 494.4 496.6 102.2 99.5 The Coles Lake and Tsea Lake watersheds each had received 0.0 mm day-1 of minimum daily PET, and 6.0 mm day-1 and 6.2 mm day-1 maximum PET, respectively. Average daily PET values (1.6 mm day-1) and standard deviation (1.6 mm day-1) were equal for both watersheds. The overall daily PET was similar for both watersheds. PET is calculated for both watersheds using the FAO Penman-Monteith method. PET is higher at the Tsea Lake watershed compared to the Coles Lake watershed (Table 6-2). The highest PET (640.1 mm and 660.5 mm) calculated in 1998, and the lowest PET (540.4 mm and 553.0 mm) calculated in 2007 for the Coles Lake and Tsea Lake watersheds, respectively. The overall calculated mean PET for both watersheds was quite similar at 583.3 mm and 601.0 mm, respectively for the Coles Lake and Tsea Lake watersheds. Overall, PET and precipitation are correlated; the highest PET occurred in 1998, when precipitation was lowest. The lowest PET occurred in 2007, in the wettest year. Table 6-2. Comparison of PET at the Coles Lake and Tea Lake watersheds (October 1979_ September 2014). Watershed Coles Lake Tsea Lake Maximum (mm) 640.1 660.5 Minimum (mm) 540.4 553.0 Mean (mm) 583.3 601.0 Std. Deviation (mm) 27.3 29.2 6.2. Coles Lake and Tsea Lake Water Balance Simulation The MIKE SHE water balance tool provided detailed information for the Coles Lake and Tsea Lake watersheds. In this study, AET and overland flow to the river form the dominant fluxes of water within the Coles Lake and Tsea Lake watersheds on an annual basis. Annually, the subsurface term was the most significant of the changes in storage, followed 148 by overland storage change, and then canopy and snow storage changes, which were both essentially zero. The yearly water balance components for both watersheds, over the 1979_2014 water year periods, are reported in Appendix E. AET, overland flow, and subsurface storage change patterns increased while overland storage change decreased slightly for both watersheds (Figure 6-1). Figure 6-1. Comparison of ET, overland flow, subsurface storage change and overland storage change between the Coles Lake and Tsea Lake watersheds for 1979 to 2014 water years. Despite the strong and significant correlation of overland flow and subsurface storage change and the moderate and significant correlation between AET and overland storage change between Coles Lake and Tsea Lake (Table 6-3), the simulated AET was higher at the Coles Lake watershed. As well, overland flow was lower at the Coles Lake watershed compared to the Tsea Lake watershed despite the significant correlations. The model simulated regular patterns of overland storage change for the Tsea Lake watershed while it showed an uneven trend for the Coles Lake watershed. 149 Table 6-3. Correlation between daily precipitation, AET, overland flow, subsurface storage change, and overland storage change for the Coles Lake and Tsea Lake watersheds over the 1979 to 2014 water year period. Coles Lake Vs. Tsea Lake AET Overland flow Subsurface storage change Overland storage change Correlation P-Value 0.68 0.00 0.90 0.00 0.98 0.00 0.65 0.00 The maximum AET (532.4 mm) was estimated between 2011_2012 for the Coles Lake watershed, while the maximum AET (441.6 mm) was estimated between 1997_1998 for the Tsea Lake watershed. The minimum AET (413.3 mm) was estimated between 1995_1996 for the Coles Lake watershed, while the minimum AET (368.0 mm) was estimated between 1995_1996 for the Tsea Lake watershed. AET forms a dominant flux of water in boreal watersheds. Overall, the MIKE SHE model reported that the mean monthly AET at the Coles Lake watershed was higher compared to the Tsea Lake watershed (Figure 6-2). Figure 6-2. Comparison of mean monthly AET between the Coles Lake and Tsea Lake watersheds (October 1979_September 2014). Dots outside of the boxes and whiskers are outliers. Maximum and minimum values are shown at the end of each whisker. The upper and lower quartiles are the ends of the box, and the median is the horizontal line in the centre of the box. 150 The AET for the Coles Lake watershed ranged from 0.1 mm month-1 (December) to 124.0 mm month-1 (July). Similar trends were obtained for the Tsea Lake watershed, as the minimum AET (0.5 mm month-1) was estimated in January and maximum AET (96.8 mm month-1) was estimated in July. Overall, the AET rate increased at the beginning of the melting period, reaching its peak in July and reaching its lowest in December and January (Figure 6-2). Mean monthly precipitation ranges were the same for both the Coles Lake and Tsea Lake watersheds (from 21.6 mm month-1 to 78.7 mm month-1). Overall, precipitation increased after the melting period; the highest amount occurred in July, and the lowest amount occurred in February for both watersheds (Figure 6-3). See Appendix F for complete details of mean monthly water balance values for the Coles Lake and Tsea Lake watersheds. Figure 6-3. Comparison of mean monthly precipitation between the Coles Lake and Tsea Lake watersheds (October 1979_September 2014). Dots outside of the boxes and whiskers are outliers. Maximum and minimum values are shown at the end of each whisker. The upper and lower quartiles are the ends of the box, and the median is the horizontal line in the centre of the box. 151 In general, the model suggests the maximum amount of overland flow occurred in May during the melting period while runoff decreased after that for both watersheds (Figure 6-4). The overall amount of overland flow at the Tsea Lake watershed was much higher compared to the Coles Lake watershed. This result is reasonable since the AET was higher at the Coles Lake watershed compared to the Tsea Lake watershed. Although snow storage varies over time, the overall snow storage changes equaled zero for both watersheds during the 1979_2017 time period. The model shows that the maximum amount of snow storage change occurred in May during the melting period (Figure 6-5). In fact, according to the trend in change of snow storage, the snowmelt phase begins in April/May when the maximum amount of overland flow occurs. In addition, the overland storage changes approach zero during the freezing period. The difference of overland storage changes from April to May marks the start of the melting phase. Figure 6-4. Comparison of mean monthly overland flow between the Coles Lake and Tsea Lake watersheds (October 1979_September 2014). Dots outside of the boxes and whiskers are outliers. Maximum and minimum values are shown at the end of each whisker. The upper and lower quartiles are the ends of the box, and the median is the horizontal line in the centre of the box. 152 Figure 6-5. Comparison of mean monthly snow storage between the Coles Lake and Tsea Lake watersheds (October 1979_September 2014). Dots outside of the boxes and whiskers are outliers. Maximum and minimum values are shown at the end of each whisker. The upper and lower quartiles are the ends of the box, and the median is the horizontal line in the centre of the box. The sign convention indicates a negative inflow of precipitation and a positive outflow of melting (MIKE SHE User Guide, 2017). 6.3. Sensitivity Experiments As discussed earlier in Section 5.3, vegetation type plays a significant role in the water balance of the Coles Lake watershed. AET was elevated in areas where there is more vegetation, e.g., dense cover compared to mixed or sparse. Therefore, a higher AET was expected at the Tsea Lake watershed compared to the Coles Lake watershed due to denser vegetation cover (Tables 2-2 and 2-3). Although MIKE SHE simulated similar trends for AET, higher AET was simulated for the Coles Lake watershed compared to the Tsea Lake watershed (Section 6.2). This shows in Section 5.3 that higher AET directly influences the overland and subsurface storage change. As previously discussed in Section 5.3, Devito et al. (2012) indicated differences in vegetation features play a considerable role in the amount of AET. The difference in AET within the boreal watersheds resultsin long-term net moisture surplus for wetlands, and net 153 deficits for forested landscapes (Devito et al., 2012). Further, Devito et al. (2012) found that climate, in addition to vegetation cover, plays a primary role in controlling the water balance in the boreal forest. Therefore, it was essential to examine the influence of climate factors, using the same meteorological forcing dataset for both watersheds. Mohan and Arumugam (1996) identified the influence of RH, temperature and WS variables (in that order) on the AET process completely. Slight changes in these variables may produce a different AET value. Nevertheless, it is essential to note that there are many other variables (e.g., area of the watershed, topography, vegetation classes) in addition to the meteorological inputs that play a significant role in AET. Some of these parameters are discussed in Section 6.3.1. Thus, this sensitivity experiment was applied to extract indicators from the vegetation cover that influence the Tsea Lake water balance by removing the influence of forcing data on the simulation results. To this end, meteorological data (e.g., precipitation, air temperature and PET) from the Coles Lake watershed were applied for the Tsea Lake watershed (October 1979_September 2014) as a sensitivity test while other parameters remained unchanged. After the simulation was completed, significant and strong correlations between the overland flow, subsurface storage change, and overland storage data from the two watersheds were found. Reasonable and significant correlation of the AET data was identified between both watersheds as well (Table 6-4). Table 6-4. Correlation between AET, overland flow, overland storage change, and subsurface storage change for the Coles Lake and Tsea Lake watersheds using the Coles Lake forcing dataset. Coles Lake vs. Tsea Lake AET Overland flow Overland storage change Subsurface storage change Correlation 0.63 0.88 0.84 0.91 P-Value 0.00 0.00 0.00 0.00 154 When the same forcing data for the Tsea Lake watershed are used, simulated AET value increased, while the overland flow, subsurface storage, and overland storage dropped (Figure 6-6). Further, the role of other parameters influencing the water balance components between these watersheds is identified in Section 6.3.1. Figure 6-6. Comparison of AET, overland flow, and subsurface storage change, and overland storage change between the Coles Lake and Tsea Lake watersheds forcing the Coles Lake dataset. 6.3.1. Discussion Northern regions of BC have warmed 1.6°C to 2.0°C from 1900 to 2013, which equals or exceeds twice the global average, and changes mainly occurring in winter months. Thus, the warming effect creates a concern for decision makers regarding whether there will be an intensification of the water cycle in this region and, if so, what the nature of the intensification will be. When the temperature of the atmosphere increases, it may intensify the water cycle and enhance components such as increased evaporation, precipitation and 155 runoff (Jacob and Hagemann, 2007). An intensification of the water cycle may lead to changes in water resource availability (Huntington, 2006), which can have an impact on the oil and gas activities within northeastern BC. Within northeastern BC, annual precipitation and AET are increasing from 1979-2014, while annual net precipitation (P – AET) is decreasing. The simulated annual AET often approaches or even exceeds annual precipitation (P) within both watersheds (see Appendix E for details), and AET and P peak simultaneously during the growing season (see Appendix F for details). These results are consistent with those of Buttle et al. (2009) and Johnson (2010) in similar environments. With AET ≥ P, subsurface storage change was often in deficit with average of -101.9 mm year and -189.5 mm year at the Coles Lake and Tsea -1 -1 Lake watersheds, respectively. Over the 35-year study period, significant inter-annual variability in subsurface storage changes — Coles Lake, from 175.4 to 166.9 mm; Tsea Lake from 170.4 mm to -162.3 mm — was detected in agreement with Buttle et al. (2009), Johnson (2010), and Devito et al. (2012) for a boreal plains ecozone. The majority of the P (60%) falls as rain from May to September (growing season) and returns to the atmosphere through ET. Although the boreal forest remains in a net water deficit overall, a net surplus occurs in the non-growing season (freezing period), and it contributes to lessening the overall deficit (Devito et al., 2012). In contrast to the growing season, during the non-growing season (October-April), a much higher proportion of the precipitation (mostly snowfall) accumulates at the surface. This snowfall accumulates in the seasonal snowpack, and ET is minimized. A higher AET rate compared to P (i.e. declining net precipitation) puts forests and wetlands within northeastern boreal watersheds at risk of significant water stress, due to 156 relatively dry conditions. In addition to water stress, dry conditions in this area can increase wildfire activity. A better understanding of the long-term annual and seasonal climatological trends thus provides essential information to decision-makers in regards to water storage and redistribution (Devito et al., 2012). Though the importance of meteorological forcing has been highlighted, there are other important and influential factors within the hydrological cycle of this region. In the following section, some of these factors are discussed. Vegetation Cover Vegetation cover influences significantly hydrologic patterns. The forest cover type affects evapotranspiration, soil infiltration, and runoff (Ekness and Randhir, 2015). For example, total annual precipitation is highest in the grasslands, compared to deciduous forests, followed by evergreen coniferous forests (Harsch et al., 2009). Consequently, the amount of water recycled to the atmosphere as vapour remains smallest in the grasslands, compared to deciduous forests, while evergreen coniferous forests experience the most evaporation (Harsch et al., 2009). The MIKE SHE simulated higher ET in dense forest cover compared to non-forested areas. The LAI values for the Coles Lake watershed range from 0 to 5.4, while at the Tsea Lake watershed they range from 0 to 9.1 (Table 3-4 and 3-5). Therefore, it was expected to detect higher AET at the Tsea Lake watershed compared to the Coles Lake watershed. However, a higher AET value was simulated at the Coles Lake watershed compared to the Tsea Lake watershed. These results demonstrate that other factors play a significant role within the northeastern boreal watershed in addition to the climatological and vegetation parameters. 157 Water Bodies The variability of land cover and distribution of water bodies in the Coles Lake and Tsea Lake watersheds plays a considerable role in the differences between their water balance components (Table 6-5). A higher AET value at the Coles Lake watershed may arise from the existence of the greater lake surface relative to the Tsea Lake watershed. Overall, AET depends on the availability of water, as more water is evaporated from a lake than from soil. Table 6-5. Water bodies for the Coles Lake and Tsea Lake watersheds (BC Watershed Atlas, 2005). Watershed Lakes (km2) Total area (km2) Coles Lake 3.3 244.1 Tsea Lake 1.2 84.7 Soils The spatial distribution of soil physical properties plays an essential role in hydrological processes (Li et al., 2013). For instance, in the UZ, the type of soil impacts the infiltration/AET, while the soil properties affect the recharge functions at the SZ level (Sandu and Virsta, 2015). However, identical soil features were applied in the simulation of the Coles Lake and Tsea Lake watersheds. Therefore, it can be assumed that the physical properties of soil play a minimal role when comparing the results between the Coles Lake and Tsea Lake watersheds. Watershed Morphometry The watershed morphometry had a dominant influence on the characteristics of its hydrograph such as peak flow, overland flow, and overland timing (Das and Saikia, 2012; Biswas et al., 2014). For instance, the overland flow moved from south to north within both watersheds. However, the overland flow was lower at the Coles Lake watershed (Figure 6-4) 158 as the concentration-time was higher. Hence, for the same rainfall event, it generated a higher flow rate where the watershed has a fan-shape (wider), such as Tsea Lake, since the concentration-time was lower (Musy, 2001). Thus, for the same volume of precipitation, delivery of water volumes occurred earlier and faster leading to a higher peak flow. Section 6.2 confirms this statement, as both watersheds received similar precipitation but the overland and peak flow at the Tsea Lake watershed was higher compared to the Coles Lake watershed. Watershed Orientation In addition, the orientation of a watershed influenced the rate of snowmelt and timing of the spring freshet. For example, the Coles Lake watershed is oriented mainly north to south with good exposure to incoming solar radiation and consequently, slower speed of snowmelt compared to the Tsea Lake watershed, positioned in the east to west direction (Musy, 2001). Slope Slope also controls the distribution of rainfall and its movement within the watershed. The Coles Lake watershed slopes from the south at peak elevations of ~700 m to its lowest of elevations of ~485 m near its outlet in the north (Figure 2-4). The elevation at the Tsea Lake watershed is higher at its southern and eastern section relative to the northern and western part of the watershed (Figure 2-5). The slope of the Coles Lake watershed is about 12% flatter compared to the Tsea Lake watershed. However, there is little significant difference between the two watersheds, as 95% of both are respectively flatter than 2.4° and 2.7°, which may be insufficient to enhance precipitation formation due to orographic effects. Thus, the slope variability between these watersheds can create a different response to the peak flow and timing of discharge. For steeper slopes, flood waters generally drain faster, 159 generating higher peak discharge. The infiltration rate decelerates and promotes surface runoff (Das and Saikia, 2012). In contrast, a gentle slope normally lessens surface runoff, which increases infiltration rate. In general, MIKE SHE was able to simulate reasonably well the hydrology of the Coles Lake and Tsea Lake watersheds. Overall, the water balance tool produced acceptable results for two small sub-watersheds of the Liard River Basin and can be employed to better explore linkages between climate and watershed response in this region. Satisfactory MIKE SHE results were also obtained at a different geographical context by Golmohammadi et al. (2014) for southwestern Ontario’s Grand River Basin (143 km2) and Sandu and Virsta (2015) for the Argesel River catchment (242 km2), Bucharest. In addition to the reasonable regional simulation of the water balance, the MIKE SHE model demonstrated a strong relationship of seasonal water balance components between the Coles Lake and Tsea Lake watersheds. 160 CHAPTER 7: SUMMARY AND CONCLUSIONS 7.1. Summary This chapter summarizes the results of Chapters 1 to 6. In addition, the limitations of this study, recommendations, and possibilities of future research opportunities are discussed in Sections 7.2 and 7.4. Chapter 1 outlined the rationale for this research. According to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC), warming of the global climate system is unequivocal (Barros et al., 2014). It is anticipated that warming in northern BC is amplified relative to the global average. In addition to non-stationarity in hydrological fluxes owing to amplified climate change in northeastern BC, the complex landscape may exacerbate the challenge of the water balance study in this region. This region comprises wetland complexes of discontinuous permafrost, fens, bogs, and swamps on a glaciolacustrine plain that makes studying the northeastern watersheds water balance challenging (Johnson, 2010) but also critical if we want to preserve wetlands while allowing further resource development in the region. Therefore, an improved understanding of the water balance of this region remains a priority for provincial ministries. Accordingly, the current study was conducted to provide both historical and detailed water balance data of two northeastern BC watersheds (the Coles Lake and Tsea Lake watersheds). Chapter 2 discussed data collection processes and fieldwork procedures in detail, which was an integral and critical part of this study. The collected data were used to quantify the detailed water balance of the Coles Lake and Tsea Lake watersheds. The field campaign involved collecting hydrometeorological data such as rainfall and snowfall within three different vegetation canopies, as well as measuring shallow groundwater fluxes and streamflow. 161 Preparation of the datasets was discussed in Chapter 3, as were the complete inputs to the MIKE SHE hydrological model. A complete suite of hydrometeorological and biophysical data (e.g., vegetation cover, elevation, etc.) collected over 35 years (October 1979_September 2014) was used to simulate the Coles Lake and Tsea Lake water balances. This was followed by the model setup and calibration process. Of note, the snowmelt, pressure head, and streamflow data were used for the Coles Lake watershed, while only the snowmelt and streamflow data were available at the Tsea Lake watershed for model calibration. Chapter 4 discussed quantification of the Coles Lake water balance based on the conceptual and numerical model. The measured and simulated water balances results were validated and compared based on available sources. Chapter 5 explored the influence of forest canopy changes (i.e., by altering the LAI value) on the rainfall and snowfall in the Coles Lake watershed. In addition, the role of vegetation cover on the Coles Lake water balance was explored. Chapter 6 described the hydrology of two wetland-dominated northeastern watersheds, with a focus on the simulation of its water balance using a hydrological model. The importance of water balance estimation in the northern region with the support of the hydrological model was highlighted by Johnson (2010), who considered the hydrological model as the most reliable support-tool for water allocation decisions and water resource management. Thus, the MIKE SHE hydrological model was employed to better understand hydrological processes and their interactions in two northeastern BC watersheds. Both the Coles Lake and Tsea Lake watersheds are characterized as wetland dominated watersheds with climate change and large-scale industrial development altering their environments. 162 7.2. Conclusions Based on the comparison between the simulated and measured results, MIKE SHE produced reliable results for the Coles Lake and Tsea Lake watersheds. It can be concluded that the MIKE SHE model simulated similar patterns as the measured data, and detailed information of the fluxes provided by the model can be used to design a sustainable water management plan in the region. However, it was quite challenging to interpret the MIKE SHE results for the lake itself (see Section 4.3 for details). In addition, as a result of vegetation growth, the simulated ET increased significantly, consequently slightly reducing the overland flow and storage change (see Section 5.3 for details). Hence, vegetation type plays a vital role in the water balance of the Coles Lake watershed. This may have implications for the influence of vegetation within all of northeastern BC. The average annual precipitation over 35 years was 499 mm and 501 mm, respectively, for the Coles Lake and Tsea Lake watersheds (see Section 6.1.2 and Appendix E for details). The AET averages for the Coles Lake and Tsea Lake watersheds were estimated at 473 mm and 406 mm, respectively. Overall, inter-annual variations in annual precipitation for both watersheds range from approximately 300 mm to 800 mm. In contrast, PET rarely varies more than 50 mm from the average of about 600 mm a year. As PET was greater than precipitation for most years, these watersheds often experience long-term moisture deficits. Similar results were obtained by Well et al. (2017) for northwestern Alberta, and Devito et al. (2012) for the boreal plains region of the western boreal forest. Fall and winter seasons often account for less than 35% of the annual precipitation, with the most precipitation occurring during the summer when ET rates peak (See Appendix F 163 for details). Similar results were reported by Devito et al. (2012) for the boreal plains region of the western boreal forest. Peak rainfall and maximum ET occurred in July. Thus, excess water was minimal due to the synchronization of peak rainfall with maximum AET (Marshall et al., 1999; Petrone et al., 2007; Devito et al., 2012; Brown et al., 2014). This situation exacerbated seasonal water deficits in wetland dominated watersheds and limits the amount of overland flow throughout the catchment (Devito et al., 2005a; Redding and Devito, 2008). The MIKE SHE model reported that the maximum overland flow occurred in May and then decreased significantly through the remainder of the year. Well et al. (2017) state that, during wet periods (e.g., spring freshet), high water levels can increase the occurrence of surface saturation on upland areas or riparian margins. As a result, it can produce overland or shallow subsurface flow in the area. The hydrologic function of wetlands is constrained by climate, vegetation, and geology (Devito et al., 2005b). The Coles Lake and Tsea Lake watersheds are very similar regarding their climatology and geology. Similar hydrological responses and trends were found by the MIKE SHE model; however, the amount of AET was higher, and overland flow was lower at the Coles Lake watershed compared to the Tsea Lake watershed. Overall, it can be concluded that the Coles Lake and Tsea Lake watersheds demonstrate a strong relationship in terms of water balance fluxes (e.g., precipitation, ET, etc.). 7.3. Limitations 7.3.1. Modelling Framework Limitations In recent years, distributed hydrological models have frequently been used to implement alternative management strategies in the areas of water resources allocation, flood control, land use and climate change impact assessments (Shi et al., 2011). Therefore, distributed 164 hydrological models such as MIKE SHE have been employed extensively. However, some researchers criticize the use of distributed models because many parameters can be adjusted during the calibration process (Abu El-Nasr et al., 2005). Indeed, a critical aspect of distributed models is the problem of over-parameterization (Beven, 1996). In contrast, a flexible parameterization procedure can overcome the problems encountered in the calibration phase (Refsgaard and Storm, 1995). In this study, flexible parameterizations of the MIKE SHE model allowed more accurate and realistic results during the calibration process. 7.3.2. Data Gaps Due to the nature of the distributed model, a significant amount of input data was required, most of which required measurements from field exercises. For example, a major challenge was a lack of river morphology and runoff data. In the initial stages, the runoff was thus underestimated from the overland component to the surface water system during the calibration phase in both watersheds. To correct this, I created additional cross-sections similar to those along the river network and lake to allow more communication between the overland component in MIKE SHE and the stream network in MIKE 11. Runoff is another water balance component, which in this research was impacted by unexpected beaver dams. Several beaver dams were built in the inflow/outflow channels in the summer of 2014. As a consequence, both the inflow and the outflow were almost entirely blocked after July 2014, and no discharge measurement was possible due to the extremely low flow rate. This situation may have further exacerbated the inaccuracy in the calibration process. 165 Further, lack of long-term pressure head and streamflow data could lead to inaccurate results during the calibration process. Specifically, only six months of pressure head data were collected for the Coles Lake watershed, while no measurement was available for the Tsea Lake watershed. Additionally, the hydraulic properties of the SZ and UZ were not available and were estimated from the literature. For this research, the NARR dataset provided the precipitation values. This dataset is believed to underestimate precipitation (Wong et al., 2017), particularly in the fall in the headwaters of the North and South Forks of the Koktuli River and Upper Talarik Creek in Alaska (Wobus et al., 2013). 7.4. Recommendations and Future Opportunities The hydrological model (i.e., MIKE SHE) forms an essential tool for managing and protecting water resources in northeastern BC boreal watersheds. The results of this study can also be utilized for other research with similar objectives in regions with comparable climatic conditions. However, all simulation results are limited to this study site, and the model application would require the appropriate watershed characteristics (e.g., soil, vegetation, and climate data). Overall, the following recommendations are offered for the Coles Lake and Tsea Lake watersheds: § Oil and gas companies need to arrange most of their water withdrawal activities at the beginning of snowmelt (May) when the water level and flow rate are at a maximum and decrease their operations in late summer when the water level is at its lowest. § To assess the influence of water withdrawals on peatland hydrology and vegetation, a water balance needs to be calculated incorporating the treatment values (water withdrawals) and values from control watersheds. Water withdrawal volume and 166 timing can then be placed in the context of the annual hydrograph and vegetation growth (Kabzems et al., 2012). § In situations where oil and gas activities increase during the summer, the following solutions may apply: o Using alternative water resources, like groundwater. o Snowmelt and rainfall can be collected and stored in borrow pits. Pumping water from these dugouts can be an alternative option. o The best practice is to combine all of these solutions to achieve an optimal result. § A major challenge in northeastern BC is a lack of river morphology, runoff data, and long-term pressure head. Therefore, it is recommended to measure the river morphology, runoff data, and groundwater and hydraulic properties within this region. § It is important to note the SWE data were measured at Fort Nelson, which may not entirely be representative of the Coles Lake and Tsea Lake watersheds. Therefore, I recommend the initiation of a long-term snow measurement campaign at both the Coles Lake and Tsea Lake watersheds. § Future snow surveys where analysis is evenly distributed across the entirety of the watershed is also recommended, as the current research is concentrated in northern areas of the watershed, exclusively around the Coles Lake area. § In addition, it could be useful to examine whether shade and shelter in forested areas would result in wind having less influence on precipitation. 167 § This research did not measure the influence of wind activity, which includes wind speed and direction. However, the results demonstrated that the patterns of precipitation in the form of snow were dissimilar to the patterns of precipitation in the observed rainfall. It is possible that wind activity may be present as an additional factor. Therefore, I recommend further investigation into the role of wind in this area, as wind may play a meaningful role in the factors affecting water balance. § I recommend recording a long-term and continuous dataset to examine more robustly differences in precipitation regimes between each canopies. § More measurements with higher temporal resolution are needed to further test the MIKE SHE model for regional applications. § The MIKE SHE model can also be used to identify the impact of other stressors on the northeastern boreal watersheds. For example, the consequences of land use change (e.g., logging) on overland flow dynamics, and the impact of increased water withdrawals 168 BIBLIOGRAPHY Abadzadesahraei, S., Déry, S.J. and Rex, J. (2017): Quantifying the water balance of Coles Lake in northeastern British Columbia using in situ measurements and comparisons with other regional sources of water information. Geoscience BC Summary of Activities 2016, Geoscience BC, Report 2017-1, 69-74. Abbott, M. B., Bathurst, J. C., Cunge, J. A., O'Connell, P. E. and Rasmussen, J. (1986): An introduction to the European hydrological system — systeme hydrologique europeen, “SHE”, 2: Structure of a physically-based, distributed modelling system. Journal of Hydrology, 87(1–2), 61-77. doi: http://dx.doi.org/10.1016/00221694(86)90115-0 Abu El-Nasr, A. Arnold, J.G., Feyen, J. and Berlamont, J. (2005): Modeling the hydrology of a catchment using a distributed and a semi-distributed model. Hydrological Processes, 19(3), 573-587. Allen, R.G., Pereira, L.S., Raes, D. and Smith, M. (1998): Crop EvapotranspirationGuidelines for computing crop water requirements-FAO Irrigation and drainage paper 56. FAO, Rome, 300(9), D05109. Aller, L., Bennett, T., Lehr, J., Petty, R. and Hackett G. (1987): DRASTIC: a standardized system for evaluating ground water pollution potential using hydrogeologic settings. EPA-600/2-87-035, National Water Well Association, Dublin, Ohio/EPA Ada. Oklahoma. Amsel, S. (2005-2017): "Canadian Province - British Columbia" Exploring Nature Educational Resource. URL 169 [November 2017]. Anderson, B., Raabis, T. and Schwarz, C. (2009): Walleye movement and population status in the Petitot River watershed; BC Ministry of Environment, Paper 2004-08. 156 p. Aneljung, M., Sassner, M. and Gustafsson, L.G. (2007): Sensitivity analysis and development of calibration methodology for near-surface hydrogeology model of Laxemar. SKB R-07-52, Svensk Kärnbränslehantering AB. ArcGIS, E. S. R. I. (2008): "9.2 Desktop help." ESRI, Redlands, CA. ArcGIS 9. Barros, V.R., Field, C.B., Dokke, D.J., Mastrandrea, M.D., Mach, K.J., Bilir, T.E., Chatterjee, M., Ebi, K.L., Estrada, Y.O., Genova, R.C. and Girma, B. (2014): Climate change 2014: impacts, adaptation, and vulnerability-Part B: regional aspects-Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. BC Ministry of Environment (2011): Climate Related Monitoring Program (CRMP). CRMP-Data Standards. British Columbia Ministry of Environment. DRAFT. BC Ministry of Forests, Lands and Natural Resource Operations and BC Oil and Gas Commission (2016): BC water portal; BC Ministry of Forests, Lands and Natural Resource Operations and BC Oil and Gas Commission, URL < http://www.bcwatertool.ca/waterportal/#5/55.318/-126.710> [August 2016]. BC Oil and Gas Commission and BC Ministry of Forests, Lands and Natural Resource Operations (2013): NorthEast Water Tool query, Lower Petitot River, British Columbia, Tuesday, November 12, 2013; BC Oil and Gas Commission, query results. 170 BC Oil and Gas Commission and BC Ministry of Forests, Lands and Natural Resource Operations (2016): NorthEast Water Tool; BC Oil and Gas Commission, GIS-based hydrology decision-support tool, URL [December 2016]. Benke, K.K., Lowell, K.E. and Hamilton, A.J. (2008): Parameter uncertainty, sensitivity analysis and prediction error in a water-balance hydrological model. Mathematical and Computer Modelling, 47(11), 134-149. Beven, K.J. (1990): A discussion of distributed hydrological modelling. In Distributed hydrological modelling. Springer Netherlands. 255-278. Biswas, A., Das Majumdar, D. and Banerjee, S. (2014): Morphometry governs the dynamics of a drainage basin: analysis and implications. Geography Journal, 927176. Blessent, D., Therrien, R. and MacQuarrie, K. (2009): Coupling geological and numerical models to simulate groundwater flow and contaminant transport in fractured media. Computers & Geosciences, 35(9), 1897-1906. Bonnaventure, P.P., Lewkowicz, A.G., Kremer, M. and Sawada, M.C. (2012): A permafrost probability model for the southern Yukon and northern British Columbia, Canada. Permafrost and Periglacial Processes, 23(1), 52-68. Bosch, J.M. and Hewlett, J.D. (1982): A review of catchment experiments to determine the effect of vegetation changes on water yield and evapotranspiration. Journal of Hydrology, 55, 3-23. Bowling, L. C. and Lettenmaier, D. P. (2010): Modeling the effects of lakes and wetlands on the water balance of Arctic environments. Journal of Hydrometeorology, 11, 276295. 171 Brazenec, W.A. and Doesken, N.J. (2005): An evaluation of two ultrasonic snow depth sensors for potential use at automated surface weather observing sites. In 13th Symposium on Meteorological Observations and Instrumentation. American Meteorological Society, Savannah, GA. Brazenec, W.A., Doesken, N.J. and Fassnacht, S.R. (2006): January. Ultrasonic Snow Depth Sensors for Measuring Snow in the US. In 10th Symposium on Integrated Observing and Assimilation Systems for the Atmosphere, Oceans, and Land Surface. Brown, S.M., Petrone, R.M., Chasmer, L., Mendoza, C., Lazerjan, M.S., Landhausser, S. M., Silins, U., Leach, J. and Devito, K.J. (2014): Atmospheric and soil moisture controls on evapotranspiration from above and within a Western Boreal Plain aspen forest. Hydrological Processes, 28, 4449-4462. Bullock, A. and Acreman, M. (2003): The role of wetlands in the hydrological cycle. Hydrology and Earth System Sciences, 7 (3), 358-389. Buttle, J.M., Creed, I.F. and Pomeroy, J.W. (2000): Advances in Canadian forest hydrology, 1995–1998. Hydrological Processes, 14(9), 1551-1578. Buttle, J.M., Creed, I.F. and Moore, R.D. (2009): Advances in Canadian forest hydrology, 2003-2007. Canadian Water Resources Journal, 34(2), 113-126. Campbell Scientific. (2013): SR50A Sonic Ranging Sensor Instruction Manual. Canadian Association of Petroleum Producers (2015): Canadian oil and natural gas; Canadian Association of Petroleum Producers website, URL [September 2015]. Canadian Water Network. (2015): 2015 Water and Hydraulic Fracturing, Where knowledge can best support decisions in Canada. Report. CBC News. 172 Center, E.M. (1979): National Centers for Environmental Prediction. National Weather Service/NOAA/US Department of Commerce, NCEP Climate Forecast System Reanalysis (CFSR) Selected Hourly Time-Series Products. Chapman, A., Kerr, B. and Wilford, D. (2012): Hydrological modelling and decision-support tool development for water allocation, northeastern British Columbia. Geoscience BC Report 2012-1, 81-86. Chen, J.M., Chen, X., Ju, W. and Geng, X. (2005): Distributed hydrological model for mapping evapotranspiration using remote sensing inputs. Journal of Hydrology, 305(1), 15-39. Curt, T., Lucot, E. and Bouchaud, M. (2001): Douglas-fir root biomass and rooting profile in relation to soils in a mid-elevation area (Beaujolais Mounts, France). Plant and Soil, 233(1), 109-125. Doi:10.1023/A:1010333308738. Das, M.M. and Saikia, M.D. (2012): Watershed management. PHI Learning Pvt. Ltd. Dekker, P.A. (2013): Evaluating the Potential for Low Impact Development to Mitigate Impacts of Urbanization on Groundwater Dependent Ecosystems using MIKE SHE (Doctoral dissertation). The University of Guelph. DeLong, C., Banner, A., MacKenzie, W.H., Rogers, B.J. and Kaytor, B. (2011): A field guide to ecosystem identification for the Boreal White and Black Spruce Zone of British Columbia. B.C. Min. For. Range, For. Sci. Prog., Victoria, B.C. Land Manag. Handb. No. 65. URL [December 2015]. DeLong, C., Annas, R. M. and Stewart. A. C. (1991): "Boreal white and black spruce zone." Ecosystems of British Columbia. D. Meidinger and J. Pojar (editors). BC 173 Ministry of Forests, Research Branch, Victoria, BC. Special Report Series 6: 237250. Déry, S.J., Taylor, P.A. and Xiao, J. (1998): The thermodynamic effects of sublimating, blowing snow in the atmospheric boundary layer. Boundary-Layer Meteorology, 89(2), 251-283. Déry, S. J. and Yau, M.K. (1999): A climatology of adverse winter-type weather events. Journal of Geophysical Research: Atmospheres, 104(D14), 16657-16672. Déry, S. J. and Yau, M. K. (2001): Simulation of an Arctic ground blizzard using a coupled blowing snow-atmosphere model. Journal of Hydrometeorology, 2, 579-598. Déry, S. J., Clifton, A., MacLeod, S. and Beedle, M. J. (2010): Blowing snow fluxes in the Cariboo Mountains of British Columbia, Canada, Arctic, Antarctic and Alpine Research, 42(2), 188-197. De Schepper, G., Therrien, R., Refsgaard, J.C., He, X., Kjaergaard, C. and Iversen, B.V. (2017): Simulating seasonal variations of tile drainage discharge in an agricultural catchment. Water Resources Research, 53(5), 3896-3920. Devia, G.K., Ganasri, B.P. and Dwarakish, G.S. (2015): A review on hydrological models. Aquatic Procedia, 4, 1001-1007. Devito, K.J., Creed, I.F. and Fraser, C.J.D. (2005a): Controls on runoff from a partially harvested aspen-forested headwater catchment, Boreal Plain, Canada. Hydrological Processes, 19, 3-25. Devito, K., Creed, I., Gan, T., Mendoza, C., Petrone, R., Silins, U. and Smerdon, B. (2005b): A framework for broad-scale classification of hydrologic response units on the 174 Boreal Plain: is topography the last thing consider? Hydrological Processes, 19, 1705-1714. Devito, K., Mendoza, C. and Qualizza, C. (2012): Conceptualizing water movement in the Boreal Plains. Implications for watershed reconstruction. Synthesis report prepared for the Canadian Oil Sands Network for Research and Development, Environmental and Reclamation Research Group, 164 p. DHI Water and Environment (2007): MIKE SHE user manual, volume 2: reference guide; DHI Water & Environment, Horsholm, Denmark, 386 p. Dubé, I. and Rimouski, E.W.S.O. (2003): From mm to cm Study of snow/liquid water ratios in Quebec. Technical Note, SMC Region of Quebec, 1-127. Dunn, S.M. and Mackay, R. (1995): Spatial variation in evapotranspiration and the influence of land use on catchment hydrology. Journal of Hydrology, 171, 49-73. Eckhardt, K., Breuer, L. and Frede, H.G. (2003): Parameter uncertainty and the significance of simulated land use change effects. Journal of Hydrology, 273, 164-176. Edwards, K.A., Classen, G.A. and Schroten, E.H.J. (1983): The water resource in tropical Africa and its exploitation. ILCA Research Report No. 6 International Livestock Centre for Africa. ISBN 92-9053-043-X. Ekness, P. and Randhir, T.O. (2015): Effect of climate and land cover changes on watershed runoff: A multivariate assessment for storm water management. Journal of Geophysical Research: Biogeosciences, 120(9), 1785-1796. Environment and Climate Change Canada (2016a): Canadian climate normals 1981-2010 station data; Environment and Climate Change Canada, URL [August 2016]. Environment and Climate Change Canada (2016b): Wateroffice; Environment and Climate Change Canada, database, URL [January 2016]. Faria, D. A., Pomeroy, J. W. and Essery, R. L. H. (2000): Effect of covariance between ablation and snow water equivalent on depletion of snow-covered area in a forest. Hydrological Processes, 14, 2683-2695. Feiner, K. and Lowry, C.S. (2015): Simulating the effects of a beaver dam on regional groundwater flow through a wetland. Journal of Hydrology: Regional Studies, 4, 675-685. Ferone, J.M. and K.J. Devito. (2004): Shallow groundwater-surface water interactions in pond-peatland complexes along a Boreal Plains topographic gradient. Journal of Hydrology, 292, 75-95. Finch, J. and Calver, A. (2008): Methods for the quantification of evaporation from lakes. Prepared for the World Meteorological Organization’s Commission for Hydrology. Fischer, A. and Durocher, Y. (2007): 6.6 An algorithm deriving snowfall from an ensemble of sonic ranging snow depth sensors. In 14th Symposium on Meteorological Observation and Instrumentation. Fisheries and Environment Canada (1978): Mean annual lake evaporation; in Hydrological Atlas of Canada, Fisheries, and Environment Canada, plate 1 7, URL [August 2016]. 176 Foote, L. and Krogman, N. (2006): Wetlands in Canada’s western boreal forest: Agents of Change. The Forestry Chronicle. 82(6), 825-833. Foster, S. (2014): Characterizing Groundwater-Surface Water Interactions within a Mountain to Ocean Watershed, Lake Cowichan, British Columbia (Master thesis). Simon Fraser University, BC, Canada, 270 p. Fundamental of Environmental Measurements, streamflow measurements. (2016): URL [February 2016]. Gerten, D., Schaphoff, S., Haberlandt, U., Lucht, W. and Sitch, S. (2004): Terrestrial vegetation and water balance—hydrological evaluation of a dynamic global vegetation model. Journal of Hydrology, 286(1), 249-270. Ghosh, A., Tiwari, A.K. and Das, S. (2015): A GIS based DRASTIC model for assessing groundwater vulnerability of Katri Watershed, Dhanbad, India. Modeling Earth Systems and Environment 1(3), 1-14. Gibson, J.J., Edwards, T.W.D., Bursey, G.G. and Prowse, T.D. (1993): Estimating evaporation using stable isotopes: quantitative results and sensitivity analysis for two catchments in northern Canada. Nordic Hydrology, 24(2-3), 79-94. Golder Associates. (2010): Surface Water Study – Horn River Basin. Report submitted to Oil and Gas Division, Ministry of Energy, Mines and Petroleum Resources. URL [October 2013]. 177 Golmohammadi, G., Prasher, S., Madani, A. and Rudra, R. (2014): Evaluating three hydrological distributed watershed models: MIKE-SHE, APEX, SWAT. Hydrology, 1(1), 20-39. Gonsamo, A. and Chen, J.M. (2014): Improved LAI algorithm implementation to MODIS data by incorporating background, topography, and foliage clumping information. IEEE Transactions on Geoscience and Remote Sensing, 52(2), 10761088. Graham, D.N. and Butts, M.B. (2006): Flexible, integrated watershed modelling with MIKE SHE. Watershed models, 245-272. Granger, R. J. and Hedstrom, N. (2011): Modelling hourly rates of evaporation from small lakes. Hydrology and Earth System Sciences, 15(1), 267-277. Gray, D. M. and Male. D. H. (1981): Handbook of Snow: Principles, Processes Management & Use, Pergamon Press, Toronto, 13-31. Greenland, D. (1994): Regional context of the climate of H.J. Andrews Experimental Forest, Oregon. In: Tenth Annual Pacific Climate (PACLIM) Workshop, 4-7 April 1993, Asilomar Conference Center, Pacific Grove, CA, 41-57. Grundling, P., Clulow, A.D., Price, J.S. and Everson, C.S. (2015): Quantifying the water balance of Mfabeni Mire (iSimangaliso Wetland Park, South Africa) to understand its importance, functioning and vulnerability. Mires and Peat, 16 p. Gupta, H.V., Kling, H., Yilmaz, K.K. and Martinez, G.F. (2009): Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling. Journal of Hydrology, 377(1-2). 80-91. 178 Harrelson, C.C., Rawlins, C.L. and Potyondy, J.P. (1994): Stream channel reference sites: an illustrated guide to field technique. General Technical Report. RM-245. Fort Collins, CO: U.S. Department of Agriculture, Forest Service, Rocky Mountain Forest and Range Experiment Station. 61 p. Hartman, C.W. and Carlson, R.F. (1973): Water balance of a small lake in a permafrost region. University of Alaska, Institute of Water Resources. Harsch, N., M. Brandenburg, and Klemm, O. (2009): Large-scale lysimeter site St. Arnold, Germany: Analysis of 40 years of precipitation, leachate and evapotranspiration, Hydrology and Earth System Sciences, 13, 305-317. Hendriks, M. (2010): Introduction to Physical Hydrology. Oxford University Press. Oxford University Press. Holding, S. and Allen, D.M. (2015): Aquifer vulnerability mapping in Northeast BC. Draft Report BC Ministry of Forests, Lands, and Natural Resource Operations. September 2015. Holding, S., Allen, D.M., Notte, C. and Olewiler, N. (2017): Enhancing water security in a rapidly developing shale gas region. Journal of Hydrology: Regional Studies, 11, 266-277. Holland, S.S. (1976): Landforms of British Columbia: a physiographic outline; BC Ministry of Energy, Mines and Petroleum Resources, Bulletin 48, 138 p. Huggett, R. (2007): Fundamentals of Geomorphology. Routledge. eBOOK. PhysicalGeography.net Hughes, J. and David. A. (2015): Clear Look at BC LNG: Energy security, environmental implications, and economic potential. 179 Huntington, T.G. (2006): Evidence for intensification of the global water cycle: review and synthesis. Journal of Hydrology, 319(1), 83-95. Huntley, D., Hickin A. and Chow. W. (2011): Surficial geology, geomorphology, granular resource evaluation and geohazard assessment for the Maxhamish Lake map area (NTS 94O), northeastern British Columbia. Geological Survey of Canada. Open file 6883. Ingram, H.A.P. (1983): Hydrology. In: Gore, A.J.P. (ed.) Mires: Swamp, Bog, Fen and Moor, Ecosystems of the World 4A, Elsevier, New York, 67-158. Interactive Broadcasting Corporation. (2007): The Trees of British Columbia, Canada. Interactive Broadcasting Corporation. Global Marketing Specialists. URL [January 2016]. IPCC. (2007): Summary for Policymakers. In: Climate Change 2007: Impacts, Adaptation and Vulnerability. Contribution of Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, M.L. Parry, O.F. Canziani, J.P. Palutikof, P.J. van der Linden and C.E. Hanson, Eds., Cambridge University Press, Cambridge, UK, 7-22. Jaber, F. H. and Shukla, S. (2012): MIKE SHE: Model use, calibration, and validation. Transactions of the ASABE, 55(4), 1479-1489. Jacob, D. and Hagemann, S. (2007): Intensification of the hydrological cycle: an important signal of climate change. In Global change: Enough water for all? Wissenschaftliche Auswertungen, 170 -173. 180 Johnsen, T., Ferbey T., Levson V. M. and Kerr. B. (2004): Quaternary geology and aggregate potential of the Fort Nelson Airport area. Summary of Activities: 2004-1. Johnson, E.G. and Johnson, L.A. (2012): Hydraulic fracture water usage in northeast British Columbia: Locations, volumes and trends. Geoscience Reports, 41-63. Johnson, E. (2010): Conceptual water model for the Horn River Basin, northeast British Columbia (NTS 094O, parts of 094P, J). Geoscience Report 2010. BC Ministry of Energy, Mines and Petroleum Resources, 99-121. Johnston, G. H. and Brown, R. J. F. (1964): Effect of a lake on distribution of permafrost in the Mackenzie River delta. Arctic, 17(3), 162-175. Kabzems, R., D’Aloia, M., Rex, J., Geertsema, M., Maloney, D. and MacKenzie, W. (2012): Monitoring indicators of wetland dynamics and hydrologic function in Fort Nelson watersheds. FLNRO internal report. Kennedy, M. (2011): BC Oil & Gas Commission _ experiences in hydraulic fracturing; Ministry of Economy, Warsaw, Poland. Khyade, V.B. (2016): Hydraulic Fracturing; Environmental Issue. World Scientific News, 40-58. Kristensen, K. J. and Jensen, S. E. (1975): A model for estimating actual evapotranspiration from potential evapotranspiration. Hydrology Research, 6(3), 170-188. Kwong, Y. T. J. and Gan, T. Y. (1994): Northward migration of permafrost along the Mackenzie Highway and climatic warming, Climatic Change, 26, 399-419. Leij, F. J. (1996): The UNSODA unsaturated soil hydraulic database: User's manual. Cincinnati, Ohio: National Risk Management Research Laboratory, Office of 181 Research and Development, U.S. Environmental Protection Agency. Retrieved from URL [ February 2016]. Letts, M.G., Roulet, N.T., Comer, N.T., Skarupa, M.R. and Verseghy, D.L. (2000): Parametrization of peatland hydraulic properties for the Canadian Land Surface Scheme. Atmosphere-Ocean, 38(1), 141-160. Lewkowicz AG, Etzelmüller B. and Smith SL. (2011): Characteristics of discontinuous permafrost from ground temperature measurements and electrical resistivity tomography, southern Yukon, Canada. Permafrost and Periglacial Processes, 22, 320–342. DOI: 10.1002/ppp.703. Li, X., Zhang, Q. and Ye, X. (2013): Effects of spatial information of soil physical properties on hydrological modeling based on a distributed hydrological model. Chinese Geographical Science, 23(2), 182 p. Liu, Y.B., and De Smedt, F. (2004): WetSpa extension, a GIS-based hydrologic model for flood prediction and watershed management. Vrije Universiteit Brussel, Belgium, 126 p. MacLean, A. (2009): Calibration and analysis of the MESH hydrological model applied to cold regions (Master's thesis, University of Waterloo). 188 p. Madsen, H. (2000): Automatic calibration of a conceptual rainfall-runoff model using multiple objectives. Journal of Hydrology, 235(3-4), 276-288. Marshall, I.B., Schut, P. and Ballard, M. (1999): Canadian Ecodistrict Climate Normal for Canada 1961-1990. A National Ecological Framework for Canada: Attribute Data. Environmental Quality Branch, Ecosystems Science Directorate, Environment Canada and Research Branch, Agriculture and Agri-Food Canada, Ottawa/Hull. 182 Matrix Solutions Inc. (2011): Monitoring Results Dilly and Cordova Leases, April 2012. Calgary Alberta. Matsuzaki, E., Rex, J., Kabzems, R., D’Aloia, M., Maloney, D., and Abadzadesahraei, S. (2013): A hydrological review and water budget estimate for the Tea Lake watershed near Fort Nelson, British Columbia. FLNRO internal report. Matsuzaki, E., Rex, J., Kabzems, R., D’Aloia, M., and Maloney, D. (2012): Water balance estimation in a boreal wetland watershed (Coles Lake) near Fort Nelson, British Columbia. FLNRO internal report. McCuen, R.H. (1973): The role of sensitivity analysis in hydrologic modeling. Journal of Hydrology, 18(1), 37-53. McDowell, J. A. (1996): The ecology of the boreal white and black spruce zone. British Columbia, Ministry of Forests. Mesinger, F., DiMego, G., Kalnay, E., Mitchell, K., Shafran, P.C., Ebisuzaki, W., Jović, D., Woollen, J., Rogers, E., Berbery, E.H., and Ek, M.B. (2006): North American regional reanalysis. Bulletin of the American Meteorological Society, 87(3), 343360. Metcalfe, R.A. and Buttle, J.M. (1998): A statistical model of spatially distributed snowmelt rates in a boreal forest basin. Hydrological Processes, 12, 1701-1722. Miceli, C. M., and Lewkowicz. A. G. (2011): Seasonal cycling in electrical resistivities at thin permafrost sites, southern Yukon and northern British Columbia. In AGU Fall Meeting Abstracts, vol. 1, p. 04. 183 Ministry of Forests and Range. (2007): Tree Book: Learning to Recognize Trees of British Columbia. [February 2016]. Ministry of Sustainable Resource Management. (2002): Vegetation Resource Inventory, the B.C. land cover classification scheme. Terrestrial Information Branch for the Terrestrial Ecosystems Task Force- Vegetation Resources Inventory Committee. Version 1.3. Mohamoud, Y. (2004): Comparison of hydrologic responses at different watershed scales (pp. 1-81). Office of Research and Development, United States Environmental Protection Agency. Mohan, S. and Arumugam, N. (1996): Relative importance of meteorological variables in evapotranspiration: Factor analysis approach. Water Resources Management, 10(1), 1-20. Musy, A. (2001): Watershed Characteristic, www.civil.pdn.ac.lk/ academic%. 05/10/2005. National Aeronautics and Space Administration (NASA), Snow Observations Network: Winter’s Story. Snow Pit Procedures. URL [October 2013]. Naturallywood, (2010): British Columbia’s Boreal Region. URL [November 2017]. Network for Environment and Weather Applications (NEWA). (2016): Accuracy of rain gauge measurements. URL 184 [September 2015]. Oil, B.C. and Gas Commission. (2013): Oil and gas land use in the Fort Nelson LRMP area of British Columbia. Onuchin, A., Burenina, T. and Pavlov, I. (2017): Hydrological consequences of timber harvesting in landscape zones of Siberia. Environments, 4(3), 51 p. Parish, R. and Thomson, S. (1994): Tree book: learning to recognize trees of British Columbia. Pelster, D., Burke J.M., Couling K., Luke S.H., Smith D.W. and Prepas E.E. (2008): Water and nutrient inputs, outputs, and storage in Canadian boreal forest wetlands: a review. Journal of Environmental Engineering and Science, 7(S1), 35-50. Petrone, R.M., Silins, U. and Devito, K.J. (2007): Dynamics of evapotranspiration from a riparian pond complex in the Western Boreal Forests, Alberta, Canada. Hydrological Processes, 21, 1391-1401. Pidwirny, M. (2006): Periglacial Processes and Landforms. Fundamentals of Physical Geography, 2nd Edition. URL [November 2016]. Pike, R.G., Redding, T.E., Moore, R.D., Winker, R.D., and Bladon, K.D. (2010): Compendium of forest hydrology and geomorphology in British Columbia. BC Min. For. Range. In For. Sci. Prog., Victoria, BC and FORREX Forum for Research and Extension in Natural Resources, Kamloops, BC Land Management Handbook (Vol. 66). 185 Pomeroy, J.W. and Goodison, B.E. (1997): Chapter 4, Winter and Snow. In the Surface Climates of Canada, 68-100. Post, V., Kooi, H. and Simmons, C. (2007): Using hydraulic head measurements in variabledensity ground water flow analyses. Groundwater, 45(6), 664-671. Prucha, R.H., Leppi, J., McAfee, S. and Loya, W. (2011): Integrated Hydrologic Effects of Climate Change in the Chuitna Watershed Alaska. 24 p. Raes, D. (2009): The ETo calculator reference manual version 3.2, food and Agriculture Organization of the United Nations. Land and Water Division, Via delle Terme di Caracalla, 153 p. Redding, T.E. and Devito, K.J. (2008): Lateral flow thresholds for aspen forested hillslopes on the western boreal plain, Alberta, Canada. Hydrological Processes, 22, 42874300. Refsgaard, J.C.; Storm, B. (1995): MIKE SHE. In Computer Models of Watershed Hydrology; Singh, V.P., Ed. Water Resources Publications: Highlands Ranch, CO, USA; 809-846. Refsgaard, J.C., Storm, B. and Refsgaard, A. (1995): Recent developments of the Systeme Hydrologique Europeen (SHE) towards the MIKE SHE. IAHS Publications-Series of Proceedings and Reports-International Association Hydrological Sciences, 231, 427 p. Regional Aquatics Monitoring Program. (2017): "The Boreal Forest and Biodiversity" URL [November 2017]. 186 Resources Information Standards Committee. (2009): Manual of British Columbia Hydrometric Standards. British Columbia Ministry of Environment, Province of British Columbia, Version, 1, 204 p. Ridgway, J. W., Higuchi R., Hoffman L., and Pettit R. (2016): Environmental Impacts of Water Withdrawals and Discharges in Six Great Lakes Communities: A Role for Green Infrastructure. Environmental Consulting & Technology Inc. Report, 45 p. Romanovsky V.E., Smith S.L. and Christiansen, H.H. (2010): Permafrost thermal state in the polar northern hemisphere during the International Polar Year 2007–2009: a synthesis. Permafrost and Periglacial Processes, 21, 106-116. DOI: 10.1002/ppp.689 Rosenberry, D.O., LaBaugh, J.W. and Hunt, R.J. (2008): Use of monitoring wells, portable piezometers, and seepage meters to quantify flow between surface water and ground water. Field techniques for estimating water fluxes between surface water and ground water. US Geological Survey Techniques and Methods, 4-D2. Roulet, N. and Woo. M.K. (1986): Wetland and lake evaporation in the low Arctic. Antarctic and Alpine Research, 18 (2), 195-200. Sandu, M.A. and Virsta, A. (2015): Applicability of MIKE SHE to Simulate Hydrology in Argesel River Catchment. Agriculture and Agricultural Science Procedia, 6, 517524. Sanford, S.E., Creed, I.F., Tague, C.L., Beall, F.D. and Buttle, J.M. (2007): Scaledependence of natural variability of flow regimes in a forested landscape. Water Resources Research, 43(8). Schmidlin, T.W. (1989): Climatic summary of snowfall and snow depth in the Ohio snowbelt at Chardon. Ohio Journal of Science, 89, 101-108. 187 Simandl, G.J., Ferbey, T., Levson, V.M., Demchuck, T.E., Hewett, T., Smith, I.R. and Kjarsgaard, I. (2005): Heavy mineral survey and its significance for diamond exploration, Fort Nelson area, British Columbia, Canada; In Hickin, A.S. and Paulen, R.C. (compilers) 2005: Poster presentations of Quaternary studies in northeastern British Columbia and northwestern Alberta; Mineral Exploration Roundup 2005; British Columbia Ministry of Energy and Mines, Petroleum Geology Open-File 2005-1; Alberta Energy and Utilities Board, Alberta Geological Survey Information Series 131; Geological Survey of Canada, Open File 4943; CD-ROM. Shi, P. Chen, C. Srinivasan, R. Zhang, X.S. Cai, T. Fang, X.Q. Qu, S.M. Chen, X. and Li, Q.F. (2011): Evaluating the SWAT model for hydrological modeling in Xixian watershed and comparison with the XAJ model. J. Am. Water Resources Management, 25, 2595-2612. Smith, I.R., Trommelen, M. and Ferbey, T. (2005): Quaternary geology studies in northeastern British Columbia; In Hickin, A.S. and Paulen, R.C. (compilers) 2005: Poster presentations of Quaternary studies in northeastern British Columbia and northwestern Alberta; Mineral Exploration Roundup 2005; British Columbia Ministry of Energy and Mines, Petroleum Geology Open-File 2005-1; Alberta Energy and Utilities Board, Alberta Geological Survey Information Series 131; Geological Survey of Canada, Open File 4943; CD-ROM. Smith MW. and Riseborough, DW. (1996): Ground temperature monitoring and detection of climate change. Permafrost and Periglacial Processes, 7, 301-310. Smith SL, Romanovsky VE, Lewkowicz AG, Burn CR, Allard M, Clow GD, Yoshikawa K. and Throop J. (2010): Thermal state of permafrost in North America: a contribution 188 to the International Polar Year. Permafrost and Periglacial Processes, 21, 117-135. DOI: 10.1002/ppp.690. Sokolov, A.A. and Chapman, T.G. (1974): Methods for water balance computations; an international guide for research and practice-A contribution to the International Hydrological Decade. Spence, C., Guan. X. and Phillips., R. (2011): The hydrological functions of a boreal wetland. Wetlands, 31, 75-85. Spence, C., Rouse, W. R., Worth, D. and Oswald, C. (2003): Energy balance processes of a small northern lake. Journal of Hydrometeorology, 4(4), 694-701. Stats, British Columbia. (2012): Population estimates. URL < http://www. bc stats. gov. bc. ca/Statistics By Subject/Demography/Population Estimates. aspx. > [October 2016]. Stewart, R. E., R. W. Crawford, H. G. Leighton, P. Marsh, G. S. Strong, G. W. K. Moore, H. Ritchie, W. R. Rouse, E. D. Soulis, and B. Kochtubajda. (1998): The Mackenzie GEWEX Study: The water and energy cycles of a major North American river basin. Bulletin of the American Meteorological Society, 79(12), 2665-2683. Stothoff, S. (2012): Literature review of groundwater flow in permafrost. Technical report number 2012:43. ISSN: 2000-0456. Available at www.stralsakerhetsmyndigheten.se Thompson, J. R., Sørenson, H. R., Gavin, H. and Refsgaard, A. (2004): Application of the coupled MIKE SHE/MIKE 11 modelling system to a lowland wet grassland in Southeast England. Journal of Hydrology, 293(1–4), 151-179. Thompson, R.I. (1977): Geology of Beatton River, Fontas River and Petitot River MapAreas, Northeastern British Columbia; Geological Survey of Canada, Paper 75-11, 8 p. 189 United States Department of Agriculture (USDA). Snow Survey Sampling Guide. URL [October 2013]. University of British Columbia. (2016): Centre for Forest Conservation Genetics. URL [July 2016]. Voeckler, H.M., Allen, D.M. and Alila, Y. (2014): Modeling coupled surface water– groundwater processes in a small mountainous headwater catchment. Journal of Hydrology, 517,1089-1106. Wang, T., Spittlehouse, D.L., and Murdock, T.Q. (2012): ClimateWNA-high resolution spatial climate data for western North America. Journal of Applied Meteorology and Climatology, 51, 16-29. Warren, F. J., E., R. Schwartz, J. Andrey, Mills., B. and Riedel. D. (2004): Climate Change Impacts and Adaptation: a Canadian Perspective. URL< http://www.citeulike.org/group/13619/article/6113882.> [August 2015]. Wells, C., Ketcheson, S. and Price, J. (2017): Hydrology of a wetland-dominated headwater basin in the Boreal Plain, Alberta, Canada. Journal of Hydrology, 547, 168-183. Wigmosta, M.S., Vail, L.W. and Lettenmaier, D.P., (1994): A distributed hydrologyvegetation model for complex terrain. Water Resources Research, 30(6), 1665-1679. Winkler, R. D., Spittlehouse, D. L. and Golding D. L. (2005): Measured differences in snow accumulation and melt among clearcut, juvenile, and mature forests in southern British Columbia. Hydrological Processes, 19(1), 51-62. 190 Winkler, R. D. and Moore R. D. (2006): Variability in snow accumulation patterns within forest stands on the interior plateau of British Columbia, Canada. Hydrological Processes, 20(17), 3683-3695. Winkler, R.D., Moore, R.D. Redding, T.E., Spittlehouse, D.L., Carlyle-Moses, D.E. and Smerdon, B.D. (2010a): Hydrologic processes and watershed. Compendium of forest hydrology and geomorphology in British Columbia. BC Min. For. Range, 66. Winkler, R.D., Moore, R.D. Redding, T.E., Spittlehouse, D.L., Smerdon, B.D. and CarlyleMoses, D.E. (2010b): The Effects of Forest Disturbance on Hydrologic Processes and Watershed. Compendium of forest hydrology and geomorphology in British Columbia. BC Min. For. Range, 66. Winsemius, H.C., Schaefli, B., Montanari, A. and Savenije, H.H.G. (2009): On the calibration of hydrological models in ungauged basins: A framework for integrating hard and soft hydrological information. Water Resources Research, 45(12), W12422, doi:10.1029/2009WR007706. Wobus, C.W., Prucha, R.H., Albert, D. and Jones, R. (2013): Climate Change Alters Future Hydrologic Regimes in an Alaskan Salmon Stronghold. In AGU Fall Meeting Abstracts. Wong, J.S., Razavi, S., Bonsal, B.R., Wheater, H.S. and Asong, Z.E. (2017): Intercomparison of daily precipitation products for large-scale hydro-climatic applications over Canada. Hydrology and Earth System Sciences, 21(4), 2163-2185. Woo, M.K. and Waddington, J.M. (1990): Effects of beaver dams on subarctic wetland hydrology. Arctic, 43(3), 223-230. 191 World Meteorological Organization. (2003): Manual on the Global Observing System. Volume I, WMO-No. 544, Geneva. World Wildlife Fund. (2009): HydroSHEDS; World Wildlife Fund, dataset, URL [October 2015]. World Meteorological Organization. (1976): The CIMO International Evaporimeter Comparisons. WMO-No. 449, Geneva. Yan, J. and Zhang, J. (2001): Evaluation of the MIKE-SHE modeling system. Southern Cooperative Series Bulleting, 398. Yang, X. and You, X. (2013): Estimating parameters of van Genuchten model for soil water retention curve by intelligent algorithms. Applied Mathematics & Information Sciences, 7(5), 1977-1983. Zhang, G., Xie, H., Yao, T. and Kang, S. (2013): Water balance estimates of ten greatest lakes in China using ICESat and Landsat data. Chinese Science Bulletin, 58(31), 3815-3829. Zhang, L., Dawes, W.R. and Walker, G.R. (2001): Response of mean annual evapotranspiration to vegetation changes at catchment scale. Water Resources Research, 37(3), 701-708. Zhang, Z., Wang, S., Sun, G., McNulty, S.G., Zhang, H., Li, J., Zhang, M., Klaghofer, E. and Strauss, P., (2008): Evaluation of the MIKE SHE model for application in the Loess Plateau, China. Journal of the American Water Resources Association, 44(5), 1108-1120. Zoltai, S.C. and Pollett, F.C. (1983): Wetlands in Canada, Their Classification, Distribution, and Use. Elsevier Scientific, 246-268. 192 Zotarelli, L., Dukes, M.D., Romero, C.C., Migliaccio, K.W. and Morgan, K.T. (2010): Step by step calculation of the Penman-Monteith Evapotranspiration (FAO-56 Method). Institute of Food and Agricultural Sciences. University of Florida. 193 Ph.D. Dissertation: University of Northern British Columbia APPENDIX A Table. A.1. Meteorological station equipment configuration and accuracy. Hydrological Component Sensor RM Young wind monitor Barometric pressure sensor Temperature and RH probe Model Number 05103-10 Measurements Units Sensitivity Accuracy Wind speed m s-1 0 to 100 m s-1 Wind direction Pressure degrees hPa N/A 0.1 hPa ± 0.3 m s-1 or 1% of reading ± 3° ± 0.5 hPa HMP45C212 Air temperature Relative humidity °C % N/A N/A 61205v Evaporation* Snowfall Sonic ranging sensor SR50 Snow depth cm 0.1 mm Snow depth (Manually) Rainfall Graduated snow probe Tipping bucket rain gauge (T-probe) Snow depth cm N/A TE525 Precipitation mm 0.254 mm per tip Davis tipping bucket 7852 Precipitation mm SWOFFER Model 2100 Velocity m s-1 Odyssey Capacitance Water Level Recorders -----Depth to the water mm Rainfall Streamflow (Manually) Groundwater The rate of rainfall range is 0 to 999 mm hr-1 Range between 0.03 to 7.5 m s-1 N/A ± 0.1°C ± 2% (090% RH) ± 3% (90100% RH) at 20 °C ± 0.4% of distance to target ± 0.5 cm ± 1% precip. rate up to 1 in./hr ±4%, ±1 rainfall count between 0.01" and 2.00" per hour (0.2 mm and 50.0 mm per hour) More than 1% ± 0.8 mm *Only instrumental errors are presented in the above table. The rate of evaporation is defined as the amount of water evaporated from a unit surface area per unit of time. It can be expressed as the mass or volume of liquid water evaporated per area in unit of time, usually as the equivalent depth of liquid water evaporated per unit of time from the whole area. The unit of time is normally a day. The amount of evaporation should be read in millimetres (WMO, 2003). Depending on the type of instrument, the usual measuring accuracy is 0.1 to 0.01 mm. For daily totals, an extreme outer range is 0 to 100 mm, with a resolution of 0.1 mm. The uncertainty, at the 95% confidence level, should be ±0.1 mm for amounts of less than 5 mm, and ±2% for larger amounts. A figure of 1 mm has been proposed as an achievable accuracy. In principle, the usual instruments could meet these accuracy requirements, but difficulties with exposure and practical operation cause much larger errors (WMO, 1976). 194 Ph.D. Dissertation: University of Northern British Columbia APPENDIX B Table B.1. Piezometer installation details and manual water level measurements following installation of the three transects. Piezometer Transect Transect 1 Transect 2 Transect 3 Depth of piezometer (cm bgs)1 Logger SN P1-1 (115 cm) P1-2 (145 cm) P1-3 (175 cm) P2-1 (115 cm) P2-2 (145 cm) P2-3 (175 cm) P3-1 (115 cm) P3-2 (145 cm) P3-3 (175 cm) 43431 43615 43737 43730 43609 43610 43733 43729 43735 GPS coordinates Elevation at Surface (m.a.s.l)2 Latitude Longitude 59°47'11.40 "N 122°36'26.70 "W 472 59°46'56.70 "N 122°37'56.10 "W 474 59°47'0.80" N 122°35'40.00 "W 475 1 Elevation of the midscreen 470.85 470.55 470.25 472.85 472.55 472.25 473.85 473.55 473.25 cm bgs – centimetres below ground surface m.a.s.l – metres above sea level (datum) *land surface slope between nested piezometers and piezometer deployed 5 m 2 195 Depth of Piezometer below ground surface (m) 1.15 1.45 1.75 1.15 1.45 1.75 1.15 1.45 1.75 Land surface slope from nested piezometers (1 & 2) to 3 (%) 6 6 5 *Actual elevation of the screen intake (m) (Adjusted considering surface slope) 470.85 470.55 470.55 472.85 472.55 472.55 473.85 473.55 473.50 Ph.D. Dissertation: University of Northern British Columbia APPENDIX C Once the ice-free period was identified, the Penman-Monteith method was selected to estimate evaporation from the lake surface (Hendriks, 2010). To use the FAO PenmanMonteith method, net radiation at the earth’s surface Rn (MJ m-2 day-1), air temperature and relative humidity RH (%) were all needed, as these determine the atmospheric demand or saturation deficit es – ea (kPa), as well as the resistances r (s m-1) to evaporation. The Penman-Monteith equation is as follows (Hendriks, 2010): 86400 × ρa Cp (es - ea ) 1000 ∆ Rn -G + ra E= × r pδ ∆+ γ 1+ s (Equation C.1) ra Here the soil heat flux density is represented by G (MJ m-2 day-1) and transfers into soil, rock, or water bodies, although it is relatively small and can often be neglected. Thus, if G is set to zero and the value of the constant is inserted into the Penman-Monteith formulation, it can then be written as follows: E = 0.408 × ∆ Rn + 105.028 (es- ea ) (Equation C.2) ra ∆+0.067 1+ rs ra where E is evaporation (mm day-1), ∆ is the gradient of the saturation pressure curve as a function of air temperature (kPa◦C-1), rs is surface resistance (s m-1), and ra is aerodynamic resistance (s m-1). The following steps were applied to estimate the net radiation at the earth’s surface. The following section was prepared using the “Introduction to Physical Hydrology” by Hendriks (2010) and “Procedure for calculating ETo using FAO Penman-Monteith with only minimum and maximum temperature” prepared by the DHI team. Step 1. Convert latitude at the UF station to radians using j ( rad ) = where: lat = latitude of the station in degrees 196 p 180 lat (°) (Equation C.3) Step 2. Calculate solar declination, d (rad) æ 2p ö J - 1.39 ÷ è 365 ø (Equation C.4) d = 0.409 sin ç where: J = Julian day Step 3. Calculate sunset hour angle, ws (rad) ws = arccos [- tan(j ) tan(d )] (Equation C.5) Step 4. Calculate extraterrestrial radiation, Ra (MJ m-2 day-1) Ra = 24 ´ 60 p Gsc d r [ws sin(j ) sin(d ) + cos(j ) cos(d ) sin(ws )] (Equation C.6) where: Gsc is the solar constant = 0.082 (MJ m-2) min-1 and dr (rad) is the inverse relative distance from the Earth to the Sun æ 2p ö d r = 1 + 0.033 cos ç J÷ è 365 ø (Equation C.7) Step 5. Calculate clear sky solar radiation, Rso (MJ m-2 day-1) ( ) (Equation C.8) Rso = 0.75 + 2 ´ 10 -5 z Ra where: z = elevation of the UF station above sea level (m) Step 6. Calculate solar radiation, Rs (MJ m-2 day-1) Use adjustment factor KRs depending on station location, Rs = K Rs (Equation C.9) (Tmax - Tmin ) Ra where: KRs = 0.16 is for interior locations and KRs = 0.19 is for coastal locations Step 7. Calculate net long wave radiation, Rnl (MJ m-2 day-1) Rnl = s (Tmax + 273 .16 ) 4 + (Tmin + 273 . 16 ) 4 2 (0 .34 - 0 .14 e ) æçç1 .35 RR - 0 .35 ö÷÷ a s è so ø (Equation C.10) where: s = 4.903 × 10-9 (MJ K-4 m-2 day-1) 197 Step 8. Calculate net solar radiation, Rns (MJ m-2 day-1) Rns = (1 - a ) Rs (Equation C.11) where: a = 0.15 Step 9. Calculate net radiation, Rn (MJ m-2 day-1) Rn = Rns - Rnl (Equation C.12) where Rns is the net incoming short wave radiation at the earth’s surface and Rnl is the net long wave radiation at the earth’s surface (Hendriks, 2010; Zotarelli et al., 2014). The saturated vapour pressure es can be estimated from T (°C) using equation (Equation C.11) and the actual vapour pressure ea from RH using the following equations: es=0.6108 e&' 17.27 T 237.3+T (Equation C.13) ea =es × RH (Equation C.14) Recall that ∆ represents the gradient of the saturation vapour pressure curve, which is the slope of the tangent line to the curve of the saturation pressure for liquid water, ∆ increases with T, and mathematically, it can be written as follows: ∆= 4098 es (Equation C.15) (237.3+T)2 The aerodynamic resistance ra is the resistance encountered by water vapour as it diffuses into the air from a vegetation surface (interception loss) or water body and is reciprocal to the roughness of the earth’s surface and WS near the surface (Hendriks, 2010). A rougher surface (e.g., a forest has a rougher surface than grasses, which have a rougher surface than open water) combined with stronger winds causes more turbulent mixing of the air and thus smaller resistance to evaporation. The aerodynamic resistance under unstressed conditions commonly varies (Table C.7-1). 198 The surface resistance rs is a physiological resistance imposed by the vegetation stomata on the movement of water vapour by transpiration: rs varies with water availability (soil moisture content). If rs is the minimum surface resistance for a vegetated surface (unstressed conditions), then the computed value is the potential evaporation. Open water has the highest evaporation rate compared to forests and grasslands due to the low albedo value (a low reflection of solar radiation) and the absence of surface resistance rs (Hendriks 2010). Table C.7-1. The possible range of values of aerodynamic resistance ra and surface resistance rs under unstressed conditions (actual evaporation = potential evaporation) for a number of lands uses (Hendriks, 2010). Land Use ra (s m-1) rs (s m-1) Forests 5-10 80-125 Grasses 50-70 40-70 Water 110-125 0 199 APPENDIX D Table D.1. A final comparison of snow survey at the three different canopies – Coles Lake Watershed. Canopy Date Min Snow Depth (cm) Max Snow Depth (cm) Average Snow Depth (cm) Coefficient of Variation (CV) (%) Standard Deviation (SD) (cm) Number of Snow Depth Sampling Number of SWE Sampling Average Snow Water Equivalent (cm) Average Snow Density -3 (g cm ) 8 8 Total Snow Course Samples Weight (g) (Standard) 67.0 87.0 Open Closed 4 Feb 2014 4 Feb 2014 10.5 26.5 52.0 67.0 30.7 46.5 30.0 20.0 9.9 7.7 64 64 8.4 10.9 273.0 235.0 Mixed Open Closed Mixed Open Closed Mixed 4 Feb 2014 3 Mar 2014 3 Mar 2014 3 Mar 2014 10 Mar 2014 10 Mar 2014 10 Mar 2014 25.0 23.0 40.0 25.0 25.0 43.0 25.0 85.0 47.0 70.0 75.0 52.0 74.0 75.0 49.7 34.5 55.2 53.9 36.5 58.3 58.2 20.0 20.0 10.0 20.0 20.0 10.0 20.0 9.5 6.4 6.5 10.8 7.1 7.1 9.1 64 36 36 36 36 36 36 8 6 6 6 6 6 6 99.0 40.0 63.0 89.0 54.0 54.0 55.0 12.4 6.7 10.5 14.8 9.0 9.0 9.2 249.0 194.2 190.2 274.6 246.5 154.4 158.0 Open Closed Mixed Open Closed 17 Mar 2014 17 Mar 2014 17 Mar 2014 24 Mar 2014 24 Mar 2014 9.0 15.0 25.0 7.0 36.0 40.0 70.0 75.0 42.0 66.0 25.8 48.1 52.3 27.6 53.5 30.0 30.0 30.0 30.0 10.0 8.6 13.3 13.1 9.6 7.9 36 36 36 36 36 6 6 6 6 6 54.0 61.0 54.0 79.0 87.0 9.0 10.2 9.0 13.2 14.5 348.8 212.0 172.1 478.3 271.3 Mixed Open Closed Mixed Open Closed Mixed Open Closed Mixed Open Closed Mixed Open 24 Mar 2014 31 Mar 2014 31 Mar 2014 31 Mar 2014 7 April 2014 7 April 2014 7 April 2014 15 April 2014 15 April 2014 15 April 2014 22 April 2014 22 April 2014 22 April 2014 28 April 2014 40.0 5.0 30.0 35.0 0.0 31.0 30.0 0.0 20.0 25.0 -7.0 20.0 -- 76.0 40.0 66.0 68.0 40.0 65.0 70.0 18.0 60.0 65.0 -40.0 70.0 -- 59.4 24.5 50.0 54.4 19.0 48.9 47.6 4.4 39.7 45.2 -25.9 40.5 -- 20.0 40.0 20.0 10.0 60.0 20.0 20.0 30.0 30.0 30.0 -40.0 30.0 -- 10.4 9.1 9.5 7.4 10.4 8.2 10.0 5.7 12.5 11.8 -9.6 12.8 -- 36 36 36 36 36 36 36 36 36 36 -36 36 -- 6 6 6 6 6 6 6 6 6 6 -6 6 -- 81.0 54.0 52.0 71.0 59.0 53.0 55.0 -38.0 72.0 -41.0 42.0 -- 13.5 9.0 8.7 11.8 9.8 8.8 9.2 -6.3 12.0 -6.8 7.0 -- 227.3 367.3 174.0 216.9 515.8 179.9 193.3 -158.7 265.5 -262.5 172.8 -- Closed 28 April 2014 7.0 46.0 26.1 Mixed 28 April 2014 10.0 68.0 46.1 40.0 9.5 36 6 64.0 10.6 406.1 40.0 16.9 36 6 58.0 9.6 208.2 200 APPENDIX E Table E.1. Coles Lake watershed simulated total water balance for each water year and yearly averages. The overland flow to the river is negative as the water leaves the watershed. The model reports small positive and negative values (sums equal to zero) for canopy and snow storage change. Water Year P (mm) AET (mm) Overland flow to river (mm) Canopy storage change (mm) 79-80 80-81 81-82 82-83 83-84 84-85 85-86 86-87 87-88 88-89 89-90 90-91 91-92 92-93 93-94 94-95 95-96 96-97 97-98 98-99 99-00 00-01 01-02 02-03 03-04 04-05 05-06 06-07 07-08 08-09 09-10 10-11 11-12 12-13 13-14 419.0 419.5 481.5 456.8 548.3 504.1 410.7 509.4 601.9 459.7 442.7 473.0 481.1 379.3 431.1 347.8 530.3 589.2 340.7 426.3 421.2 449.7 486.0 455.0 374.2 580.1 555.8 852.4 642.8 651.4 568.1 560.1 491.0 605.7 588.0 -508.2 -486.5 -454.0 -443.8 -458.0 -470.0 -473.0 -480.8 -471.0 -496.9 -481.3 -465.5 -448.8 -466.8 -475.6 -435.7 -413.3 -445.3 -513.4 -437.5 -439.5 -462.2 -419.1 -445.8 -445.0 -466.7 -489.1 -472.4 -500.7 -510.5 -514.4 -514.1 -532.4 -492.4 -521.7 -3.7 -14.3 -5.8 -4.4 -1.8 -23.3 -26.7 -1.3 -26.8 -43.3 -16.9 -0.7 -5.3 0.4 0.4 1.2 -0.5 -2.0 -4.2 0.4 0.8 0.6 -0.4 -0.8 0.6 -0.7 -3.1 -124.7 -200.0 -123.4 -100.7 -44.5 -62.9 -29.3 -53.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 201 Snow Overland storage storage change change (mm) (mm) 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 -0.1 0.0 0.0 0.0 0.0 0.0 0.0 8.8 -3.8 1.5 0.6 8.9 2.5 -9.0 1.2 17.8 -14.2 -4.4 -1.6 2.4 -5.7 -2.2 -2.0 1.3 11.3 -10.1 -1.1 -1.0 -0.3 1.4 0.8 -1.4 2.8 5.9 61.5 -32.9 3.1 -14.2 1.5 -16.3 10.4 -7.6 Subsurface storage change (mm) Total model error (mm) -103.3 -81.6 21.6 7.8 79.6 4.6 -80.3 26.2 83.7 -70.6 -51.0 8.4 22.3 -80.2 -41.9 -84.7 114.9 128.8 -166.9 -9.3 -16.7 -11.7 64.6 6.5 -68.5 109.5 55.5 175.4 -35.5 7.2 -38.6 -2.2 -91.4 67.5 -51.6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Table. E.2. Tsea Lake watershed Simulated total water balance for each water year and yearly averages. The overland flow to the river is negative as the water leaves the watershed. The model reports small positive and negative values (sums equal to zero) for canopy and snow storage change. Water Year P (mm) AET (mm) Overland flow to river (mm) Canopy storage change (mm) 79-80 80-81 81-82 82-83 83-84 84-85 85-86 86-87 87-88 88-89 89-90 90-91 91-92 92-93 93-94 94-95 95-96 96-97 97-98 98-99 99-00 00-01 01-02 02-03 03-04 04-05 05-06 06-07 07-08 08-09 09-10 10-11 11-12 12-13 13-14 421.9 419.9 478.5 461.0 550.4 510.2 408.1 521.7 598.7 466.6 438.2 469.1 476.5 383.6 439.9 354.8 534.3 597.3 342.8 431.7 419.6 450.0 485.7 460.2 388.2 580.6 546.5 849.0 645.8 639.0 555.7 564.9 494.6 614.9 533.6 -436.4 -422.9 -386.1 -384.7 -392.2 -398.9 -401.4 -412.5 -392.8 -407.8 -410.0 -399.4 -388.9 -407.9 -434.0 -407.9 -368.0 -379.7 -441.6 -384.7 -410.6 -423.7 -378.5 -400.9 -414.3 -412.6 -418.3 -385.2 -393.0 -410.0 -412.0 -416.6 -437.3 -398.5 -423.7 -47.7 -56.7 -37.0 -33.3 -35.3 -50.6 -45.2 -27.6 -74.3 -54.8 -44.1 -19.2 -34.3 -12.4 -15.9 -7.6 -17.4 -58.1 -24.7 -15.5 -6.0 -10.0 -23.3 -33.8 -22.1 -27.1 -43.5 -167.8 -153.4 -109.4 -86.0 -67.7 -74.2 -72.2 -77.7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 202 Snow Overland storage storage change change (mm) (mm) 0.5 -0.5 0.0 1.2 -1.2 0.2 -0.2 0.0 0.0 0.1 -0.1 0.0 3.4 -3.4 0.3 -0.3 0.1 -0.1 0.0 0.0 0.0 0.0 0.0 0.1 0.0 -0.1 0.0 0.4 -0.4 0.0 0.0 0.0 0.0 0.0 0.0 5.2 -3.1 0.4 0.5 1.9 0.2 -2.5 1.5 0.8 -1.2 -1.6 1.5 0.2 -1.1 -0.5 -0.9 2.8 1.1 -3.6 0.5 0.3 0.1 1.5 -1.0 -1.6 2.6 0.1 3.7 -2.4 0.7 -0.7 -0.1 -2.7 2.7 -1.8 Subsurface storage change (mm) Total model error (mm) -139.8 -106.3 10.5 1.4 70.4 0.6 -88.4 37.2 54.9 -52.5 -63.1 19.5 7.9 -57.6 -31.4 -77.5 121.8 93.8 -162.3 7.9 -14.0 -4.1 57.1 -13.8 -75.5 108.2 32.2 170.4 -28.8 17.5 -31.0 0.9 -87.8 78.9 -46.7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 APPENDIX F Table F.1. Coles Lake watershed simulated mean monthly water balance results (1979_2014). Month P (mm) AET (mm) Overland flow to river (mm) Canopy storage change (mm) Snow storage change (mm) Overland storage change (mm) Subsurface storage change (mm) Total model error (mm) Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec 28.3 21.6 27.6 27.0 47.4 65.6 77.7 65.2 43.3 36.5 30.7 28.1 -0.2 -0.7 -2.5 -19.6 -83.5 -111.3 -124.0 -83.4 -35.5 -11.5 -0.6 -0.1 -0.6 -0.5 -0.6 -5.6 -12.0 -0.9 -1.9 -0.5 -0.5 -1.4 -1.1 -0.7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -27.2 -20.9 -20.9 103.6 27.9 0.0 0.0 0.0 0.0 -7.6 -27.7 -27.1 0.4 -0.1 -0.1 18.4 -8.3 -4.9 -3.9 -2.6 0.4 0.9 -0.1 0.3 -0.2 -0.5 3.7 84.3 -12.3 -41.1 -44.1 -15.9 7.0 14.8 1.5 -0.1 0.1 0.1 0.0 -1.9 0.3 0.6 0.4 0.2 0.1 -0.1 0.1 0.1 Table. F.2. Tsea Lake watershed simulated mean monthly water balance results (1979_2014). Month P (mm) AET (mm) Overland flow to river (mm) Canopy storage change (mm) Snow storage change (mm) Overland storage change (mm) Subsurface storage change (mm) Total model error (mm) Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec 28.3 21.6 27.9 27.1 48.1 66.4 78.0 65.0 43.4 36.3 30.9 28.1 -0.4 -0.7 -2.8 -29.5 -67.7 -87.7 -96.8 -73.7 -31.7 -12.3 -1.7 -0.5 -1.1 -0.8 -0.8 -11.3 -18.4 -2.8 -4.0 -1.6 -2.3 -2.3 -1.6 -1.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -28.0 -21.5 -24.6 91.7 50.4 0.1 0.0 0.0 -0.2 -11.1 -29.5 -27.3 0.0 0.0 0.0 2.4 -1.4 -0.6 -0.7 -0.2 0.2 0.4 0.0 0.0 -4.2 -3.8 -3.0 67.1 1.2 -26.5 -26.2 -12.8 5.6 6.3 -4.8 -4.3 0.0 0.0 0.0 -0.6 0.4 0.1 0.1 0.0 0.0 -0.1 0.1 0.0 203