ASSESSMENT OF OPTIMUM SNOWMELT MODEL COMPLEXITY, CARIBOO HIGHLANDS, BRITISH COLUMBIA, CANADA. by Kara Francine Przeczek B.Sc, University of Calgary, 2005 THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE IN NATURAL RESOURCES AND ENVIRONMENTAL STUDIES THE UNIVERSITY OF NORTHERN BRITISH COLUMBIA January 2011 ©KaraPrzeczek, 2010 1*1 Library and Archives Canada Bibliotheque et Archives Canada Published Heritage Branch Direction du Patrimoine de I'edition 395 Wellington Street Ottawa ON K1A 0N4 Canada 395, rue Wellington OttawaONK1A0N4 Canada Your file Votre reference ISBN: 978-0-494-75120-6 Our file Notre reference ISBN: 978-0-494-75120-6 NOTICE: AVIS: The author has granted a nonexclusive license allowing Library and Archives Canada to reproduce, publish, archive, preserve, conserve, communicate to the public by telecommunication or on the Internet, loan, distribute and sell theses worldwide, for commercial or noncommercial purposes, in microform, paper, electronic and/or any other formats. L'auteur a accorde une licence non exclusive permettant a la Bibliotheque et Archives Canada de reproduire, publier, archiver, sauvegarder, conserver, transmettre au public par telecommunication ou par I'lnternet, preter, distribuer et vendre des theses partout dans le monde, a des fins commerciales ou autres, sur support microforme, papier, electronique et/ou autres formats. The author retains copyright ownership and moral rights in this thesis. Neither the thesis nor substantial extracts from it may be printed or otherwise reproduced without the author's permission. L'auteur conserve la propriete du droit d'auteur et des droits moraux qui protege cette these. Ni la these ni des extraits substantiels de celle-ci ne doivent etre imprimes ou autrement reproduits sans son autorisation. In compliance with the Canadian Privacy Act some supporting forms may have been removed from this thesis. Conformement a la loi canadienne sur la protection de la vie privee, quelques formulaires secondaires ont ete enleves de cette these. While these forms may be included in the document page count, their removal does not represent any loss of content from the thesis. Bien que ces formulaires aient inclus dans la pagination, il n'y aura aucun contenu manquant. 1+1 Canada Abstract The goal of this study is to assess the performance of melt models of varying complexity to simulate snowmelt under different aspect, forest cover, and input data conditions in the British Columbia interior. Observed snow water equivalent data from the 2008 melt season were used to compare the performance of a basic temperature-index (TI) model, three modified TI models, and a simple energy balance model (EBM). The largest (smallest) values of NSE (RMSE) for the snow course sites were 0.81 (0.0243 m) for the TI models and 0.58 (0.0362 m) for the EBM. At the automatic snow pillow (ASP) the largest (smallest) values of NSE (RMSE) were 0.54 (0.0055 m) for the TI models and 0.65 (0.0048 m) for the EBM. At the snow course sites all TI models performed better than or equivalent to the EBM. At the ASP one EBM version performed better than the TI models. Table of Contents Abstract ii Table of Contents iii List of Tables vi List of Figures viii List of Appendices xii Acknowledgement xiii 1. Introduction 1 2. 3. 1 1 Importance and Characteristics of Snow 1 1 2 Energy Budget in Mountainous Terrain 5 1 3 Snowmelt Modelling 10 1 4 Thesis Objectives and Outline 15 Study Site and Data Collection 18 2 1 Study Area 18 2 2 Data Collection 26 2 2 7 Meteorological Station 26 22 2 AspectPyranometers 27 22 3 Albedo 28 2 2 4 Snow Surveys 28 2 2 5 Snowmelt Lysuneters 35 2 2 6 East and West Aspect Total Season Melt 36 2 27 37 Automatic Snow Pillow 2 3 Data and Analysis Methods 38 2 3 1 Meteorological Data 38 2 3 2 Snow Course Data 38 2 3 3 TidBit Data 41 2 3 4 Automatic Snow Pillow Data 42 2 3 5 Lys imetet Data 42 Field testing and comparison of temperature-index melt models of varying complexity43 3 1 Intioduction 43 3 11 Outline and Objectives 3 2 Methods and Analysis 3 2 1 Snowmelt Models 3 2 2 Radiation Model 3 2 3 Albedo Model 3 2 4 Snowmelt Modelling 3 2 5 Measures of Error 3 3 Results 3 3 1 Meteorological Data 3 3 2 Observed Snowmelt 3 3 3 Model Calibration 3 3 3 1 Basic tempeiature-index Model 3 3 3 2 HBV-EC Model 3 3 3 3 Hock Model 3 3 3 4 Pelhcciotti Model 3 3 4 Model Performance 3 3 4 1 Mt Tom north and south aspect snow 3 3 4 2 BaikerviUe automatic snow pillow 3 3 4 3 TidBit east and west aspect sites 3 3 5 Model Comparison 3 4 Discussion 3 4 1 Model Pai ameter Values 3 4 2 Optimal Melt Model Complexity 3 4 3 Benejit of Including Shortwave Radiation 3 4 4 Sources of Error 3 5 Conclusion Energy balance modelling at Mt Tom 4 1 Introduction 4 11 Outline and Objectives 4 2 Methods and analysis 42 1 All Model 1 Radiation calculations 4 2 12 Latent heat flux 4 2 13 Sensible heat flux 42 2 Running the Model 114 423 Input data 115 42 4 Testing the Model 116 4 3 Results 117 43 1 Model Performance 117 4 3 1 1 Open Sites 117 4 3 1 2 Forest Sites 123 4 3 1 3 Energy Balance 130 4 3 2 Barkerville Automatic Snow Pillow 135 433 136 Energy Balance vs Temperature-Index Models 4 3 3 1 Mt Tom Snow Course Sites 136 4 3 3.2 140 Barkerville Automatic Snow Pillow 4 4 Discussion 142 44 1 Model Performance 142 44 2 Modelled Energy Balance 146 44 3 Energy Balance vs Temperature-Index Models 149 4 5 Conclusion 152 Conclusion 155 5 1 Optimal Melt Model Complexity 155 5 2 Data Availability and Model Performance 157 5 3 Model and Data Limitations 157 5 4 Dncctions for Furthci Reseaich 160 5 5 Summary 162 6. Literature Cited 163 5. V List of Tables Table 2 1 Air temperature and precipitation data from the Environment Canada stations at Quesnel (Q) and Barkerville (B), BC, and from the Mt Tom (MT) meteorological station for the 2008 winter/spring season The 30 year normals (1971-2000) monthly temperature and precipitation for Quesnel and Barkerville stations are italicized Values in bold are at least 2°C cooler than the appropriate station climate normal 20 Table 2 2 Mt Tom study site characteristics Tree height refers to the approximate average height of trees surrounding the opening Tree height was only measured for the main north and south open sites and the meteorological station site Due to their proximity, adjacent north and south aspect forest sites were assumed to have the same average tree heights as measured for the respective open sites 'na' indicates data that were not available 25 Table 2 3 Description of instrumentation at Mt Tom meteorological station 27 Table 2 4 Average tip volume of the tipping bucket for each snowmelt lysimeter 36 Table 2 5 Start and end dates for snow course measurement periods in 2008 for each snow course site S = south aspect, N = north aspect, O = open, F = forest * north aspect forest data were only compared until 30 May, 2008 41 Table 2 6 Start and end dates for melt periods in 2008 for each TidBit site W = west aspect, E = east aspect, O = open, F = forest 42 Table 3 1 Equations used to calculate measures of error for snow melt models tested in this study M P = predicted melt, M 0 = obseived melt, Mbp = melt predicted by the benchmark model, and N = number of samples 60 Table 3 2 Optimized melt factor (Mp) using the basic temperature-index model for each leave-oneout run (coeff set 1 4) and for all four snow course sites together (coeff set 5) 73 Table 3 3 Paiameter values calibrated using the HBV-EC model for each leave-one-out run (coeff set 1 4) and for all four snow couise sites together (coeff set 5) Minimum melt factoi (Cmm), inciease in melt factor (DC), slope-aspect reduction factoi (AM), foiest melt l eduction factoi (MRF) 74 Table 3 4 Paiameter values cahbiated using the Hock model for each leave one-out mn (set 1 4) and all sites pooled (set 5) Sub tables aic for each type of radiation input data used MF = melt factoi, Rr = ladiation factoi, -O = open site, -F = forest site 75 Table 3 5 Parameter values calibrated using the Pellicciotti model for each leave-one-out run (set 1 4) and all sites pooled (set 5) Sub-tables are for each type of radiation input data used MF = melt factoi, Rh = ladiation factoi, O = open site, -F = forest site 78 Tabic 3 6 Calculated mcasuies of enoi foi each model for both leave-one-out and all-sites combined calibration methods Subscript "suiK" indicates the model run with potential surface shortwave radiation Non-subscnpted model runs used measured shortwave radiation Root mean square enor (RJVISE), mean bias error (MBE), relative bias error (RBE), NashSutcliffe efficiency criterion (NSE), and goodness-of-fit indices (G) The basic temperatuie-indcx model (BTI) was used as the benchmark scnes N = 14 foi each loo and all methods 83 Table 3 7 Calculated measuies of error using coefficient set 5 to predict melt at the Baikeiville ASP site Subscript "surK" indicates the model mn with potential surface shortwave radiation Non-subscripted model runs used measuied shortwave ladiation Root mean square eiroi (RMSE), mean bias error (MBE), relative bias error (RBE), Nash-Sutchffe efficiency criterion (NSE), and goodness-of-fit (G) N = 30 87 Table 3 8 Calculated measures of error using coefficient set 5 to predict melt at the TidBit east and west aspect sites Subscript "surK" indicates the model run with potential surface shortwave radiation Non-subscripted model runs used measured shortwave radiation Root mean square error (RMSE), mean bias error (MBE), relative bias error (RBE), NashSutchffe efficiency criterion (NSE), and goodness-of-fit (G) The - N S and -met affixes indicate which measured shortwave radiation was used in modelling as radiation was not measured for the east and west aspects noith and south aspect (NS) or that measured at the meteorological station (met) N = 4 89 Table 3 9 The "best" model for each set of test data used based on computed measures of error Model [resulting test value] Root mean square error (RMSE), mean bias error (MBE), relative bias error (RBE), Nash-Sutchffe efficiency criterion (NSE), and goodness-of-fit (G) The - N S (north and south aspect) and -met (meteorological station) affixes indicate which measured shortwave radiation was used in modelling as radiation was not measured for the east and west aspects 99 Table 3 10 The "worst" model for each set of test data used based on computed measures of error Model [test value] Root mean square error (RMSE), mean bias error (MBE), relative bias error (RBE), Nash-Sutchffe efficiency criterion (NSE), and goodness-of-fit (G) The - N S (north and south aspect) and -met (meteorological station) affixes indicate which measuied shortwave radiation was used in modelling as radiation was not measuied for the east and west aspects 100 Table 4 1 Constant values used in the current version of the Waltei et al (2005) model 112 Table 4 2 Measures of error calculated for each variation of the Walter model and each open site Values in bold indicate the best error measure overall and values in bold italic indicate the best en or measure for the other site SO = south aspect open, NO = north aspect open '1 April' indicates that the model was run staiting 1 April, 2008, and 'Melt' indicates that the model was ran only for the melt period for each site (starting 27 Apnl, 2008 for the south and 9 May, 2008 foi the north) 'RH' indicates that relative humidity was used to calculate the latent heat flux and 'rv' indicates that resistance to heat transfer was used to calculate the tuibulent heat fluxes as in the paper (Eq 12, Walter et al 2005) N = 4 120 Table 4 3 Mcasuics of cnoi calculated foi each vauation of the Walter model and each foiest site using measuied forest covet for each site Values in bold indicate the best en or measure overall and values in bold italic indicate the best enoi measure foi the othei site SF = south aspect forest, NF = north aspect forest 'Melt' indicates that the model was ran only foi the melt period for each site (starting 2 May, 2008 for the south and 9 May, 2008 foi the north) 'RH' indicates that relative humidity was used to calculate the latent heat flux, and 'rv' indicates that the turbulent heat fluxes weie calculated as in the paper (Eq 12, Walter et al 2005) N = 3 126 Table 4 4 Measures of en or calculated for each vauation of the Waltei model and each forest site using only the optimal forest cover for each site and model variation Values in bold indicate the best enoi measuie oveiall and values in bold italic indicate the best enor measure for the other site SF = south aspect foiest, NF - north aspect forest 'Melt' indicates that the model was run only for the melt peuod foi each site (starting 2 May 2008 for the south and 9 May foi the north) 'RH' indicates that lelative humidity was used to calculate the latent heat flux, and 'rv' indicates that the tuibulent heat fluxes were calculated as in the paper (Eq 12, Waltei et al 2005) N = 3 129 Table 4 5 Calculated measures of error for the Walter model run at the Barkerville ASP site 'Melt' indicates model run only for melt period (1 May - 30 May), and "rv" indicates the model run using the original turbulent flux calculation Root mean square error (RMSE), mean bias error (MBE), relative bias error (RBE), Nash-Sutchffe efficiency criterion (NSE), and goodness-of-fit (G) N = 30 136 Table 4 6 Calculated measures of error for each Walter model run using all the snow course sites combined Root mean square error (RMSE), mean bias error (MBE), relative bias error (RBE), Nash-Sutchffe efficiency criterion (NSE), and goodness-of-fit index (G) The basic temperature-index with all-sites-combined coefficient was used as the benchmark series N = 14 137 Table 4 7 Calculated measures of error for each model for both leave-one-out and all-sites-combined calibration methods Subscript "surK" indicates the model run with potential surface shortwave radiation Non-subscripted model runs used measured shortwave radiation Root mean square error (RMSE), mean bias error (MBE), lelative bias error (RBE), NashSutchffe efficiency criterion (NSE), and goodness-of-fit indices (G) The basic temperature-index model (BTI) was used as the benchmark series N = 14 for each loo and all methods 138 Table 4 8 Calculated measures of error using coefficient set 5 to predict melt at the Barkerville ASP site Subscript "surK" indicates the model run with potential suiface shortwave radiation Non-subscripted model runs used measured shoitwave radiation Root mean square erroi (RMSE), mean bias error (MBE), relative bias error (RBE), Nash-Sutchffe efficiency criterion (NSE), and goodness-of-fit (G) N = 30 141 Table 6 1 Mean air temperature (Ta ) and incoming shortwave ladiation (K) as measured at the Mt Tom meteorological station Values were aveiaged across the penod between snow course measmcments End date of each measurement penod shown 184 List of Figures Figuie 2 1 Study area location N = north aspect, E = east aspect, S = south aspect, W - west aspect, O = open Satellite image from ESRI (2009, blip goto aicgi^onlnie com m ips T SRI ftndgo y \\ orld 2D) Map pioduccd by T Mlynowski 19 Figure 2 2 Block map of Mt Tom study area Openmg D, DDD, BB, and J weie used in the snowmelt study The meteorological station was located in opening G Map provided by Patuck Teti, BC Ministiy of Foiests, Mines and Lands 24 Figuie 2 3 Meteoiological station at Mt Tom Photo taken 9 March, 2008 26 Figuie 2 4 Diagram ot south aspect open snow course site The location of snowmelt lysimeteis, pyranometer, and snow survey sample points are shown Tree symbols denote the presence of trees, but do not represent the size oi numbei 30 Figuie 2 5 Diagram of noith aspect open snow course site The location of snowmelt lysimeteis, pyianometer, and snow survey sample points arc shown Tiec symbol denotes the piesence of a single tiees, but does not represent the size 31 Figure 2 6 Relationship between coefficient of variation (CV) and distance fiom start of snow course foi Pat Teti's 2007 snow depth data SWE data were collected on 21 March for the south aspect sites and 31 March for the noith aspect sites Vertical dotted line indicates twelve sample points length Distance fiom stait only indicates the cumulative distance fiom the start of each snow course The distance from the forest edge to the start of each snow course was not necessarily the same for each site, and these data were not available 32 Figure 2 7 Diagram of south aspect forest snow course site The locations of snow survey sample points are shown Tree symbols indicate the edge of the forest on the south side The southeast-most sample point fell on the forest edge The trees were distributed throughout the rest of the snow course, but omitted for clarity in this diagram 33 Figure 2 8 Diagram of north aspect forest snow course site The locations of snow survey sample points are shown Tree symbols indicate the edge of the forest on the north side The trees were distributed throughout the entire snow course, but omitted for clarity in this diagram 34 Figure 2 9 Mean SWE measured at each snow course site throughout the sampling period in 2008 S = south aspect, N = north aspect, O = open, F = forest a) data before correcting for SWE measurement errors at the open sites on 20 April, b) data after correction applied 40 Figure 3 1 Comparison of SR50 measured daily average snow depth at the Mt Tom meteorological station and the Barkerville ASP 62 Figure 3 2 Comparison of the daily Barkerville automatic snow pillow SWE measurements with a) the measuied SWE at each Mt Tom snow course site, and b) the averaged measured SWE for all the Mt Tom snow course sites 63 Figure 3 3 Comparison of the cumulative precipitation at the Environment Canada station and the Barkerville ASP with a) the daily precipitation at the Environment Canada station and the Barkerville ASP, and b) the cumulative SWE at the Barkerville ASP 65 Figuie 3 4 Comparison of the measured daily average air temperatures at the Mt Tom meteorological station, Barkerville ASP, and Barkerville Environment Canada station 66 Figure 3 5 a) daily average incoming shortwave radiation at the south and north aspect and meteorological station pyianometers, and b) the ratio between daily average incoming shortwave radiation at each of the three pyranometeis 67 Figuie 3 6 Change in SWE between each snow course sample date with 95% confidence inteivals Positive values indicate accumulation and negative values indicate ablation of snow S = south aspect, N = noith aspect, O = open, F = foiest 70 Figuie 3 7 a) total pcuod and b) avciage daily melt at each snow course site foi the cahbiation period 27 April, 2008 - 5 June, 2008 70 Figure 3 8 Total melt calculated for each TidBit measurement point E = east aspect, W = west aspect, F forest, O = open 71 Figuie 3 9 a) avciage daily melt rate for each Mt Tom study site as a function of total SWE melted at each site, b) number of days to melt total SWE at each Mt Tom aspect site as a function of total SWE melted at each site 72 Figuie 3 10 Average daily observed melt for each snow course site Melt was averaged acioss the inteivening measurement period, plotted with a) daily average modelled potential direct shoitwave radiation at the suiface foi the same penod, b) measuied shoitwavc incoming radiation at the meteorological station, and c) average daily air temperature at the metcoiological station 77 Figure 3 11 Obscivcd and modelled SWE over time for each of the Mt Tom snow couise sites starting 1 Apnl, 2008 a) south aspect open, b) south aspect foiest, c) noith aspect open, d) noith aspect foiest 80 IX Figure 3 12 Observed and modelled daily SWE over time for the Barkerville automatic snow pillow starting 1 April, 2008 81 Figure 3 13 Predicted vs observed melt calculated for each snow course site using leave-one-out (loo) and all-sites-combined (all) calibration methods Period of comparison was 1 April, 2008 to 5 June, 2008 The straight line indicates a l l relationship a) basic temperature index, b) HBV-EC, c) Hock (measured radiation), d) Hock (modelled surface radiation), e) Pelhcciotti (measured radiation), f) Pelhcciotti (modelled surface radiation) 86 Figure 3 14 Predicted vs observed melt calculated for the Barkerville automatic snow pillow using coefficient set 5 Period of comparison was 1 to 30 May, 2008 The straight line indicates a 1 1 relationship a) basic temperature index, b) HBV-EC, c) Hock (measured radiation), d) Hock (modelled surface radiation), e) Pelhcciotti (measured radiation) Pelhcciotti model results using potential surface shortwave radiation were identical and therefore are not shown here 88 Figure 3 15 Predicted vs obseived melt calculated for the Mt Tom east and west aspect TidBit sites using coefficient set 5 Period of comparison was 1 April, 2008, to 5 June, 2008 The straight line indicates a l l relationship a) basic temperature index, b) HBV-EC, c) Hock (measuied meteoiological station radiation), d) Hock (measuied north and south aspect radiation), e) Hock (modelled surface radiation), f) Pelhcciotti (measured radiation) Pelhcciotti model results using all other radiation data were identical and therefore are not shown here 90 Figure 3 16 Observed and modelled SWE over time for the penod of continuous melt for each of the Mt Tom snow course sites a) south aspect open, b) south aspect forest, c) north aspect open, d) north aspect forest 92 Figure 3 17 Observed and modelled SWE over time for the period of continuous melt for the Barkerville ASP (1-30 May, 2008) 93 Figure 3 18 Observed and modelled SWE over time for each of the TidBit sites a) east aspect open, b) cast aspect forest, c) west aspect open, d) west aspect foiest 94 Figure 4 1 Waltei model tun starting 1 Apnl 2008 Compauson of observed (eucles) veisus modelled (line) SWE ovei time at a) south aspect open site, and b) noith aspect open site 117 Figuie 4 2 Compauson of obseived (eucles) veisus modelled (line) SWE ovei time a) model run liom stait ot melt, 27 April 2008, at the south aspect open site, b) model iun fiom stait of melt, 9 May 2008, at the north aspect open site, c) model iun fiom 27 Apnl, 2008 at the south aspect open site using relative humidity to calculate the latent heat flux, and d) model iun from 9 May, 2008 at the north aspect open site using relative humidity to calculate the latent heat flux 118 Figure 4 3 Compauson of obseived (cnclcs) versus modelled (line) SWE over time The ouginal equation incorporating resistance to heat tiansfer to calculate the turbulent fluxes was used (Eq 12, Waltei et al 2005) a) model run fiom start of melt, 27 April, 2008, at the south aspect open site, b) model run from stait of melt, 9 May, 2008, at the north aspect open site, c) model run fiom 27 Apnl, 2008 at the south aspect open site using relative humidity to calculate the latent heat flux, and d) model run from 9 May, 2008 at the north aspect open site using lclative humidity to calculate the latent heat flux 119 Figuie 4 4 Piedictcd veisus obseived melt for the obseived snow eouise penods Model run starting 1 Apnl 2008 The concsponding penods of compauson aie indicated on the figure Diagonal line indicates a l l relationship a) model run for south aspect open site, only the foui melt penods aie shown, b) model ran for north aspect open site, only the four melt periods are shown, c) model mn for south aspect open site, all periods from 1 April to end of observed snow cover are shown, and d) model run for north aspect open site, all periods from 1 April to end of observed snow cover are shown 121 Figure 4 5 Predicted versus observed melt for the observed snow course periods The corresponding periods of comparison are indicated on the figure Diagonal line indicates a l l relationship a) model run from start of melt, 27 April, 2008, at the south aspect open site, b) model run from start of melt, 9 May, 2008, at the north aspect open site, c) model run from 27 April at the south aspect open site using relative humidity to calculate the latent heat flux, and d) model run from 9 May at the north aspect open site using relative humidity to calculate the latent heat flux 122 Figure 4 6 Predicted versus observed melt for the observed snow course periods The original equation incorporating resistance to heat tiansfer to calculate the turbulent fluxes was used (Eq 12, Walter et al 2005) The conesponding periods of comparison are indicated on the figure Diagonal line indicates a l l relationship a) model run from start of melt, 27 April 2008, at the south aspect open site, b) model run from start of melt, 9 May 2008, at the north aspect open site, c) model run from 27 April at the south aspect open site using relative humidity to calculate the latent heat flux, and d) model run from 9 May at the north aspect open site using relative humidity to calculate the latent heat flux 123 Figure 4 7 Time scries of SWE using the measured forest cover fraction at the south aspect site, for a) model mn from start of melt, 2 May 2008, b) using relative humidity to calculate latent heat flux, c) using the original equation incorporating resistance to heat transfer to calculate the turbulent fluxes, and d) using the resistance to heat transfer and relative humidity to calculate the fluxes 124 Figuie 4 8 Time senes of SWE using the measuied forest covei fraction at the north aspect site, for a) model run from stait of melt, 2 May 2008, b) using lelative humidity to calculate latent heat flux, c) using the original equation incorporating resistance to heat transfei to calculate the turbulent fluxes, and d) using the resistance to heat transfer and relative humidity to calculate the fluxes 125 Figure 4 9 Time scucs of SWE foi each 10% increment in forest cover between 0 and 1 for a) model mn ftom stait of melt, 2 May 2008, at the south aspect foiest site, b) model mn from stait of melt at the south aspect foiest site using lelative humidity to calculate latent heat flux, c) model lun fiom stait ot melt, 9 May 2008, at the north aspect foiest site, and d) model run from stait ot melt at the north aspect foiest site using lelative humidity to calculate latent heat flux 127 Figure 4 10 Time series of SWE foi each 10% increment in forest cover between 0 and 1 The onginal equation mcorpoiating resistance to heat transfei to calculate the turbulent fluxes was used (Eq 12, Walter et al 2005) a) model run from start of melt, 2 May 2008, at the south aspect foiest site, b) model run from start of melt at the south aspect forest site using relative humidity to calculate latent heat flux, c) model run fiom start of melt, 9 May 2008, at the north aspect foiest site, and d) model run fiom start of melt at the north aspect forest site using relative humidity to calculate latent heat flux 128 Figure 4 11 Compauson of the measured and modelled incoming shortwave ladiation to the surface The top panel is for the south aspect open site and the bottom panel is for the north aspect open site The potential incoming shortwave was adjusted foi tiansmissivity using the daily an tempcratuic lange 130 Figure 4 12 Latent, sensible, and net cncigy flux foi each model mn for a) the south aspect open and b) the north aspect open sites 'Melt' indicates model was mn from start of melt, 'RH' indicates that relative humidity was used to calculate the latent heat flux, and 'rv' indicates that the resistance to transfer term was used to calculate the turbulent fluxes 133 Figure 4 13 Latent, sensible, and net energy flux for each model run for a) the south aspect forest and b) the north aspect forest sites 'Melt' indicates model was run from start of melt, 'RH' indicates that relative humidity was used to calculate the latent heat flux, and 'rv' indicates that the lesistance to tiansfei term was used to calculate the turbulent fluxes 134 Figure 4 14 Comparison of observed (circles) versus modelled (line) SWE over time at the Baikerville ASP Model was run only for the melt period 1-May to 30-May, 2008 using a) the new equation to calculate turbulent heat fluxes, and b) the original equation incorporating resistance to heat transfer to calculate the turbulent fluxes (Eq 12, Walter et al 2005) 135 Figure 4 15 Observed and modelled snow course SWE over time using variations of the temperatureindex model a) south aspect open, b) south aspect forest, c) north aspect open, d) north aspect forest 139 Figure 4 16 Obseived and modelled snow course SWE over time using variations of the Walter et al (2005) model a) south aspect open, b) south aspect foiest, c) north aspect open, d) north aspect forest 140 Figure 4 17 Observed and modelled SWE over time for the period of continuous melt for the BarkerviUe ASP a) Walter energy balance model results, and b) temperature-index model results 141 Figure 4 18 The Mt Tom foiest snow couise sites a) south aspect, and b) north aspect 146 Figure 6 1 Cumulative sum of precipitation measuied at the BarkerviUe Environment Canada station and with a tipping bucket rain gauge at Mt Tom, and obtained by converting snow depth (measured with SR50) accumulations to water equivalent assuming a 10 1 density latio 176 Figure 6 2 Comparison of the occurrence of snowfall (denoted by 1) when a snowfall event was extiacted from the SR50 recoid oi when the precipitation data fiom the Envuonmcnt Canada station was used 179 Figuie 6 3 Companson ot daily albedo modelled fot each snow couise site and the albedo measuied at the Mt Tom metcoiological station from 9 May to 30 May 180 Figuie 6 4 Comparison ot the occmience of snowfall (denoted by 1) when the precipitation data from the Baikerville ASP was used oi when the piecipitation data fiom the Envuonmcnt Canada station was used 181 Figure 6 5 Modelled daily albedo foi the Mt Tom snow couise sites and the Barkeiville ASP 182 Figuie 6 6 Modelled daily albedo for the Mt Tom TidBit sites 183 List of Appendices Appendix 1 BarkerviUe Automatic Snow Pillow 175 Appendix 2 Piecipitation at Mt Tom 176 Appendix 3 Temperatmc-index Albedo Model at Mt Tom 178 Appendix 4 Hock and Pelhcciotti Paiameter Comparison 184 Acknowledgement There are many people I would like to thank for their assistance, support, and contribution to the completion of this thesis. First and foremost I would like to thank my supervisors Drs. Stephen Dery and Brian Menounos for their advice, encouragement, and guidance, and for giving me this research opportunity. I would also like to thank my committee member, Dr. Dan Moore from UBC, for all his help from coding models and troubleshooting R to providing valuable references and identifying comma splices. I would like to thank my entire committee for taking the time to answer my many questions and Emails. Finally, I would like to thank my external examiner, Dr. Rita Winkler, for participating in my defense and for her insightful and constructive comments on my thesis. I would like to thank the Northern Hydrometeorology group at UNBC for their assistance in the field and in the office. The members during my time at UNBC included Shane MacLeod, Jinjun Tong, Tullia Upton, Theo Mlynowski, Marco Hernandez, Dr. Andy Clifton, Ian Picketts, Bryndel Fell, and Heidi Knudsvig. I would like to give a special thank you to Pat Teti with the Ministry of Forests, Mines and Lands in Williams Lake for giving me access to his data, the use of his equipment, and his time to answer questions and provide advice. Dr. M.T. Walter from Cornell University kindly sent me his model code in R and provided advice on its use. Drs. Scott Green and Peter Jackson at UNBC provided much of the meteorological equipment. I would also like to thank Scott Jackson with the River Forecast Centre for his assistance and access to data, Brian Hutchinson and Paul Whitfield with Environment Canada for providing data, and the people at Mecca Electronics and Battery Direct in Prince George who went out of their way to help. I also want to recognize all the people who generously volunteered their time to help set up equipment, calibrate sensors, wrestle snowmobiles, and brave winter weather and late nights in the field. Thank you Tullia Upton, Andy Clifton, Melanie Grubb, Theo Mlynowski, Dan Burns, Peter Davies, Jay and Pat from Wells, Lyssa Maurer, Andy's Mom, Alex Swigoda, Kristina Poitras, Audrey Blake, Joanne Lee, Shane MacLeod, Laura Thomson, Nikki Bloomfield, and Denny Straussfogel. Finally I would like to thank my family and friends for their support and encouragement and for enduring the thesis process! This research was supported by the National Sciences and Engineering Research Council of Canada and the British Columbia Pacific Leaders Fellowship. XIII 1. Introduction 1.1 Importance and Characteristics of Snow Snow is a key component of the hydrologic cycle and is an important freshwater reservoir. This storage is especially important in areas where most annual precipitation occurs during the winter. Snowmelt runoff supplies more than 50% of the annual streamflow in regions above 45°N latitude and supplies an estimated 17% of the global population with water (Barnett et al. 2005). Snow cover lasts longest at high elevations and latitudes due to low temperatures. Snow accumulation and distribution patterns are a function of several factors including dominant regional weather, climate conditions during and between snowfalls, frequency of snowfall, physiography, and vegetative cover (McKay and Gray 1981). Land features that affect atmospheric and snow retention processes control the spatial distribution and characteristics of the seasonal snowpack. During the accumulation period, variations in snow depth arise from snow-canopy interactions, snow redistribution by wind, and orographic influences on precipitation. Local variability in these processes results in high spatial heterogeneity of snow cover. Snow possesses a number of unique characteristics that, collectively, make it an important element of the climate system. Fresh snow, for example, has the highest albedo of any natural surface, reflecting up to 95% of incoming solar radiation. Snow is also a good emitter of long wave radiation, which results in high rates of radiative cooling from snow covered areas. Snow is likewise an excellent insulator and limits the movement of heat between the ground and air and keeps snow covered soils warmer than exposed soils and typically much warmer than air temperatures. Finally, the temperature at the snow surface cannot exceed 0°C, which can create large temperature gradients between the snow and the air during the melt period. From a climatic perspective, the most important effect of snow is on the amount of solar energy that is retained by the earth system (Berry 1981). The presence of snow generally reduces energy exchange between the surface and atmosphere. The high albedo and emissivity of snow enhances the energy deficit at the poles which helps drive global circulation, transferring excess energy from the equator and lower latitudes poleward to areas of deficit. Snow cover gives rise to a large difference between summer and winter surface albedo, which significantly affects seasonal surface temperatures. The high albedo and rates of radiative cooling of snow give rise to the cold, high pressure air masses typically seen over the northern continents during winter. These air masses typically move southward, bringing cold air to the mid-latitudes and causing major storms and heavy precipitation in these areas (Berry 1981). From a human perspective, the presence of snow dictates where and when many activities can occur. Snow affects growing season length and agricultural capability and suitability, for example. Flooding from snowmelt runoff can be costly as well as potentially deadly. However, humans also rely heavily on the runoff from snowmelt to supply water for irrigation, municipalities, power generation, and fisheries. Because a smaller proportion of snow is evaporated or transpired compared to rain, snowfall contributes proportionally more to runoff (Dingman 2002). In most of western North America, the spring freshet is the period of greatest discharge during the year with snowmelt supplying 50%-80% of annual river flow volume (Steppuhn 1981, Barnett et al. 2005, Stewart et al. 2005). The timing and magnitude of snowmelt runoff is controlled by meteorological events and the water 2 equivalent of the snowpack. The timing and volume of the freshet is important for reservoir storage, hydroelectric generation, and flood control. The timing of the freshet also affects flow levels during spring and summer with an earlier melt resulting in lower late season flows and a longer low-flow period (Dery et al. 2009; Stewart et al. 2005). British Columbia (BC) has a wide range of climates and topography from semi-desert valley bottoms and arid tundra to coastal rainforests and glacierized mountain ranges. In BC, particularly in the interior of the province, runoff is predominantly snowmelt-driven (Moore and Wondzell 2005; Eaton and Moore 2007); snowmelt contributes more than 50% of the annual river discharge in the interior of BC (Stewart et al. 2004). The pattern of snow accumulation and melt in mountainous catchments is complex as it is modified by atmospheric, biologic, and topographic variations and their interactions (Male and Granger 1981; Olyphant 1986a; Pliiss and Ohmura 1997; Watson et al. 2006; Ellis and Pomeroy 2007). For nival systems, the annual low streamflow typically occurs during winter as the precipitation is stored in solid form, although low streamflows can also occur in late summer or early fall during dry years (Moore and Wondzell 2005; Eaton and Moore 2007). Forest harvesting is an important factor for hydrology in BC due to the vast areas covered by forest and the amount of forest harvesting that occurs in the province. Once trees have been removed the area experiences increases in incoming shortwave radiation and wind speed and decreases in incoming longwave radiation, litter, and canopy interception. These effects typically result in increases in albedo, SWE, and melt rates, with the end effect that forest harvesting tends to increase annual water yield (Moore and Wondzell 2005). In addition, the mountain pine beetle epidemic in BC has raised concerns regarding the potential impact on snow accumulation and melt and watershed hydrology (Rex and Dube 2006; van de Vosse et al. 2008; Boon 2009). Globally averaged air temperatures are predicted to rise 0.4°C in the next two decades (IPCC 2007). Fewer cold events, more frequent warm events, and greater temperature increases at high altitudes and latitudes have been observed and are projected for the future (Bradley et al. 2004, IPCC 2007). Temperatures in BC have been increasing over the twentieth century, with the largest increases seen for minimum temperatures (Bonsai et al. 2001; Mote et al. 2005; Vincent and Mekis 2006; Burford et al. 2009), and temperatures are predicted to increase in the coming decades (Stewart et al. 2004). One of the predicted consequences of this trend is the reduction in the areal extent and persistence of snowpack (Mote et al. 2005). Reductions in snow cover extent and volume have been observed over western North America (Mote et al. 2005; Dery and Brown 2007; Burford et al. 2009). Along with this warming trend, a shift toward an earlier spring freshet for BC rivers has been observed (Loukas et al. 2002; Morrison et al. 2002; Stewart et al. 2005; Burford et al. 2009; Dery et al. 2009) and a resulting decline in summer/autumn flows (Rood et al. 2008). A study on the Fraser River watershed observed earlier occurrence of the spring freshet and predict even earlier freshet timing and declining peak flows with continued warming (Morrison et al. 2002). Unfortunately predictions for changes in precipitation are highly variable (Loukas and Quick 1999; Whitfield and Cannon 2000; Rood et al. 2005; Stewart et al. 2005). In many areas a greater proportion of winter precipitation may fall as rain rather than snow. This decrease in solid precipitation has a two-fold impact on snowpack: there is less snow accumulating and rain increases melt rates above those caused by warmer temperatures alone (Knowles et al. 2006). These changes will have significant consequences for snowpack, snowmelt onset, hydrology, and water resource management. Surface runoff in BC is closely linked to snow accumulation and melt, which are dependent on a number of factors including topography, land cover changes and climate. The variability of accumulation and ablation processes results in variable melt rates and runoff contribution within a basin (Male and Granger 1981; Olyphant 1986a,6; Hock 1999; Fierz et al. 2003). 1.2 Energy Budget in Mountainous Terrain The general consensus in the literature is that net radiation contributes the largest amount of energy for snowmelt (Cline 1997). However, some studies indicate turbulent and radiative exchanges are nearly equal (Aguado 1985), and some studies find that turbulent exchanges are dominant (Prowse and Owens 1982; Moore and Owens 1984; Jackson and Prowse 2009). There is also disagreement among studies as to whether longwave or shortwave inputs are more important for open sites (Aguado 1985; Ohmura 2001). Unfortunately, the presentation of radiation data among different studies confounds the comparison of study results. Some studies present the net values of each energy term whereas other studies present the absolute incoming values. This is an important point because it is only the net energy at the snow surface that is available for melt. Incoming solar radiation is a large term, for example, but once the snow albedo is taken into account, the energy provided by solar radiation is reduced by up to 95%. This flux term is also zero at night. Incoming longwave radiation can also be a large term, especially under certain conditions, but snow has a high emissivity; thus the net longwave contribution is usually negative (Cline 1997). The contribution of the individual energy terms also varies through time and among study sites. Male and Granger (1981) reviewed the contribution of radiative and turbulent fluxes to snowmelt. They observed a high variability in contribution of the individual energy budget components among similar sites and concluded that the dominant flux cannot be predicted for a particular environment. Both topography and large scale circulation patterns affect the relative magnitudes of the fluxes. The radiation balance in mountainous catchments is complex as it is modified by atmospheric and topographic variations and their interaction. The atmospheric factors that reduce incoming radiation are widely variable and random in nature as they are linked to weather conditions and specifically to cloud cover (Cazorzi and Dalla Fontana 1996). The topographic factors include elevation, aspect, slope, and shading, and are much simpler to incorporate. Complex alpine terrain can further complicate the energy balance as it results in topographic reflection, shading, and emission. Albedo is important for determining the shortwave flux and leads to multiple reflections in snow covered terrain (Male and Granger 1981). Olyphant (1986a) modelled radiation inputs to an alpine watershed during the melt season. Aspect, shading, and slope affected the radiation budget at each location. The interaction of topography with radiation resulted in compensating effects such that different aspect-slope combinations received similar overall energy inputs despite variable inputs by the individual radiation components. Comparing direct insolation values over a watershed may be indicative of relative radiation exposure, but actual differences in total incoming radiation are substantially less than suggested by comparisons of direct beam insolation alone. Olyphant (1986a) concluded that terrain heterogeneity was nearly as important as the 6 effects of spectral heterogeneity in determining variations in the surface radiation balance. Pliiss and Ohmura (1997) concluded that longwave radiation from surrounding topography is an important energy term in areas of high relief. Olyphant (19866) compared the distribution of incident and terrain longwave radiation in a mountainous environment. The headwall orientation of the study basins explained part of the variation in total terrain radiation. The differences in incident atmospheric and terrain longwave radiation often compensated resulting in smaller differences in total longwave radiation among basins. The surrounding rock walls increased the radiation balance in cirques by decreasing the net longwave loss by 37-63%. As well, the temperature of terrain features was usually greater than air temperature, which invalidated the common assumption that these terms are equal. Thus estimating incoming longwave radiation without accounting for terrain effects is inaccurate. Forest cover also affects the energy balance at the snow surface (Hardy et al. 1997; Link and Marks 1999; Storck et al. 1999; Koivusalo and Kokkonen 2002; Pomeroy et al. 2009). Trees shelter the surface from wind, intercept and reflect shortwave radiation, and emit longwave radiation. Accumulation of needles, leaves, bark, and other debris on sub-canopy snow can reduce the albedo (Link and Marks 1999). Turbulent fluxes and short wave radiation are lower and longwave radiation is greater in the forest compared to open areas and these effects increase with forest density (Suzuki et al. 1999). The lower energy inputs to the snow surface can sometimes be offset by the increased longwave radiation from the canopy (Hardy et al. 1997; Pomeroy et al. 2009). It is the balance between the incoming and outgoing radiation that determines the relative rate of melt in the forest in comparison to the open (Lopez-Moreno and Stahli 2008). The sub-canopy energy balance thus is affected by forest structure and composition. Boon (2009) observed that shortwave radiation was higher 7 and longwave radiation was lower in beetle killed lodgepole pine stands in comparison to live lodgepole pine stands as a result of the loss of needles, small branches and stems. Shrubs can also affect ablation rates by increasing net radiation to the snow and affecting turbulent fluxes (Pomeroy et al. 2006) Almost every study that incorporated radiation into a snowmelt model discussed cloud cover. Cloud cover has a significant effect on the energy balance by altering the relative contribution of long- and shortwave radiation (Male and Granger 1981; Young 1985). Under clear sky conditions, longwave radiation represents a net loss from the snowpack whereas the fluxes can balance under completely overcast conditions (Cline 1997). Due to the high albedo of snow, the radiative flux is essentially governed by the longwave fluxes under cloudy conditions (Male and Granger 1981). Increasing the cloud cover in a climate change study resulted in greater snowmelt, as the increase in longwave radiation offset the reduction in shortwave radiation (Brubaker et al. 1996). Cloud cover reduces the variability of incoming radiation and thus the influence of aspect (Bloschl et al. 19916). Snow albedo can also be significantly affected by the presence of clouds: a thicker cloud cover usually increases albedo, but it is reduced by greater cloud height (Choudhury and Chang 1981). Thus the energy budget under overcast conditions may be very different from that under clear conditions. Cloud cover is frequently noted as a potential source of error in energy balance or snow melt calculations (e.g. Dunn and Colohan 1999) but it is rarely incorporated in energy balance studies because it is so difficult to quantify (Cazorzi and Fontana 1996). Kustas et al. (1994) measured six-fold differences in global radiation under one-eighth cloud cover, highlighting the variability associated with partly cloudy conditions. Simpson et al. (2004) observed streamflow fluctuations that were due to the passage of transient and 8 scattered clouds over the study basins. They concluded that a seasonal mean could not represent the variability in cloud cover and insolation and the influences of these clouds were not captured by large-scale fluctuations of temperature at weather stations typically used to force the models. However, the inclusion of these clouds via new remote sensing techniques improved runoff simulations (Simpson et al. 2004). The addition of direct beam short- or longwave radiation terms and the effects of shading, slope, and aspect on their receipt (e.g. Hock 1999) are important steps toward improving the simulation of snowmelt within mountainous terrain. However, the effects of cloud cover, reflection, backscatter, vegetation, and terrain emissions will need to be included to completely describe the complex contribution of radiation to snowmelt in mountainous terrain. Air mass characteristics also contribute to the energy balance. The energy and moisture content of an air mass controls the turbulent sensible and latent heat fluxes, respectively. In their review of energy balance studies Male and Granger (1981) noted that net radiation had the lowest contribution during a melt event that was caused by warm air advection, the same period for which latent heat had the greatest contribution. In general, the effect of an air mass on the energy balance is as follows: cold-dry = negative sensible flux and evaporation from snow; warm-dry = positive sensible flux and evaporation; and warm-wet = positive sensible flux and condensation (Male and Granger 1981). The transition from warm to cold air masses leads to high evaporative rates and large latent heat losses relative to the snow surface. Zimmerman (1972) noted that surface air temperature is a poorer index for melt when an air mass is changing significantly, and Male and Granger (1981) therefore suggested using air mass and weather forecast information as input for snowmelt studies. 9 Cline (1997) related the variations in energy input and melt to synoptic weather conditions. Four different synoptic patterns were identified during the study and mean air temperature, humidity and energy fluxes were significantly different between each of the patterns. The different synoptic patterns affected the overall energy available for melt by changing the relative contributions of the individual energy terms. The systems affected the energy budget by altering the air temperature (available energy), the humidity gradients, atmospheric stability, cloud cover, and wind speeds over the snow. Thus the relationship between synoptic patterns and the contribution of energy balance components may be useful for forecasting snowmelt runoff. 1.3 Snowmelt Modelling Snowmelt is determined by the energy balance at the snow surface, and this depends on the relative importance of individual energy components. The energy balance of a snowpack generally takes the form: *Qs = Q * +QH +QE + QG+QP AQs= change in energy storage Q* = net radiation (short- and longwave) QH = sensible heat transfer QE = latent heat transfer QG = heat transfer from the ground QP = heat transfer from precipitation QM = energy used for melt 10 + QM (i-i) 7 -9-1 The above energy terms are commonly in units of either W m" or J m" s" . Individual energy components vary with climate, season, and surface at the large scale and with time of day, terrain, and weather conditions at small scales (Hock 2003). During the melt period, variations in melt processes affect snow distribution through the availability of solar radiation due to slope and aspect, modification of short- and longwave radiation inputs by vegetation and surrounding topography, and local advection. The variability of accumulation and ablation processes results in a mosaic of snow-covered and snow-free areas during the melt period (Liston 2004), and variable melt rates and runoff contribution within a basin (Hock 1999). Snow accumulation and melt are sensitive to climate and land cover changes, and runoff may vary dramatically from year to year. An increasing demand for fresh water has heightened the need for more precise water management strategies. Melt modelling is a critical component of any attempt to model runoff from snow covered or glacierized areas; the success of modelling this runoff depends on the accurate representation of the snowmelt process. Hydrologic models are developed for two primary purposes: 1) to guide water resource management procedures; and 2) as tools for scientific investigation (Dingman 2002). Much hydrologic research is directed toward improving the ability to predict the effects of land use and climate change on aspects of the hydrologic cycle. Modelling is extensively used in hydrology because of the need to address these forecasting issues, because the processes involved are complex and have high spatial and temporal variability, and because of the limited availability of relevant hydrologic data. There are two basic types of snowmelt models. Energy balance models are based on fundamental physical principles and attempt to quantify melt from the energy balance 11 equations (Equation 1.1). Temperature-index (or "degree day") models simulate melt based on an empirical relationship between air temperature and ablation. Energy balance models are more accurate than temperature index models due to their physical basis, and this accuracy has been well established at the site scale (e.g. Anderson 1968; Bloschl et al. 1991a). Energy balance models have been successfully applied for point snowmelt measurements and to small watersheds (e.g. Anderson 1968; Aguado 1985; Bloschl et al. 19916; Cline 1997; Fiertz et al. 2003). However, energy balance models require large amounts of accurate input data that can inhibit their application at the watershed scale or for operational runoff models (Whitaker et al. 2003; Franz et al. 2008). According to Anderson (1973), the minimum data required for an energy balance model are incoming solar radiation, air temperature, vapour pressure, and wind speed. Dingman (2002) adds cloud cover, precipitation, snow surface temperature, and incident and outgoing longand shortwave radiation to the list. These data need to be collected frequently enough to calculate at least daily averages and at high enough density to accurately represent the dominant energy fields of interest. Anderson (1968) commented on the difficulty of obtaining adequate meteorological and snowpack measurements during his study. Even at heavily instrumented sites, the measurement of energy exchanges is difficult, and the use of empirical relationships or modelling is necessitated (Male and Granger 1981). The components of an energy balance model can be calculated using more readily available meteorological measurements (c.f. equations in Dingman 2002, for example), but a number of assumptions and approximations need to be made. Watson et al. (2006) opted for a daily time step in their energy balance study due to the difficulty of estimating meteorological variables at sub-daily time scales. Walter et al. (2005) modelled snowmelt with a physically 12 based energy budget using no more data than required for simple temperature-index models daily maximum and minimum air temperatures and precipitation. The results for the processbased model were generally as good as or better than for the temperature index model. The predicted and observed individual components of the energy budget diverged, however. The simplified energy balance calculations did not accurately represent individual energy fluxes. Acceptable snowmelt predictions resulted due to the overestimation of some parameters compensating for the underestimation of others (Walter et al. 2005). Another energy balance model was developed by Watson and others (2006) that required only daily temperature and precipitation data to run. However, this model still required the use of 30 parameters, 11 of which had to be calibrated for the site. Specific problems related to melt timing were noted under deep snowpack situations, which were hypothesized to be related to excessive accumulation of heat deficit (Watson et al. 2006). However, the timing of snowpack ripening was correctly simulated, which may have been due to the balancing of errors in heat flux simulations to and from the snowpack (Watson et al. 2006). The additional data required by energy balance models combined with uncertainties in input data plus the need for calibration in some models can cause greater uncertainties compared to simpler temperature-index models (Franz et al. 2008). The most basic temperature-index models require only minimum and maximum daily temperature values and usually give snowmelt and runoff results that are comparable to the energy balance approach at the basin scale (WMO 1986; Hamlin et al. 1998). Temperatureindex modelling is used for most operational hydrologic models, as input for ice dynamic models, and to predict the effect of climate change on the cryosphere (Brubaker et al. 1996; Semadeni-Davies 1997; Hock 2003). The temperature index approach enjoys this popularity 13 because of the wide availability of air temperature data, the relatively easy interpolation and forecasting possibilities of air temperature, computational simplicity, and generally good model performance (Hock 2003). One important advantage of using the temperature index approach to model snowmelt runoff is that these models are quick to set up and calibrate for a watershed compared to more rigorous models. Temperature-index models do not account for the physical processes behind snowmelt and are limited by temporal and spatial resolution. The melt factor, which is calibrated and used to convert air temperature into melt, varies in both space and time, yet the value used in modelling usually represents only the average for the basin (Kuusisto 1980; Lang and Braun 1990; Braithwaite 1995; Rango and Martinec 1995; Hock 2003). A temperature index approach is employed in two of the most popular runoff models currently in use: the Snowmelt Runoff Model (SRM) (Martinec and Rango 1986) and the Hydrologiska Byrans Vattenbalansavdelning (HBV) model (Bergstrom and Forsman 1973). Due to its modest data requirements and reliable results it will likely retain its popularity into the future (Rango and Martinec 1995), but it must rise to the demand for spatially distributed snowmelt modelling at shorter time scales. Distributed models give more realistic simulations of snow cover changes over time (Bloschl et al. 19916), and are becoming the norm for both temperature index and energy balance models (Kirnbauer et al. 1994). It has been shown that the incorporation of additional variables to give temperature index models a stronger physical basis often improves model performance. The components of the energy balance that are responsible for snowmelt vary through time as a result of atmospheric and topographic factors. Modifications to temperature index models to address this variability are leading to a gradual transition from the classical temperature index approach to energy- 14 balance-type expressions. Promising areas of advancement incorporate the long- and shortwave components of the energy budget based on digital elevation models (DEMs) and solar geometry. These inputs improve the physical basis of the model without incorporating many or any additional meteorological variables thus maintaining the characteristic low data requirements and usefulness in remote and poorly instrumented basins. 1.4 Thesis Objectives and Outline Snow is an important component of the climate system, and people rely heavily on snowmelt runoff for water supply in much of the Northern Hemisphere. Predicted changes in climate will have significant consequences for the region's snowpack, snowmelt, hydrology, and water resource management. Thus, an understanding of snow hydrology processes and accurate forecasting of runoff is important in the context of global change. Developments in snowmelt modelling have yielded melt models of varying complexity and input data requirements that range from simple temperature-index to complete energy balance models. Micovic and Quick (2009: 873) define "optimal" model complexity as ".. .the minimum watershed model structure required for realistic representation of runoff processes." The optimal complexity necessary for snowmelt models has been debated for some time (Dunn and Colohan 1999; Ferguson 1999; Sivapalan 2003). The intended application of a snowmelt model is also an important consideration in determining "optimal" complexity. When the model is to be applied at the watershed scale for operational purposes the model should be based on sound science but with reasonable data requirements, a level of complexity supported by its performance, and it must be understandable by users (Lindstrom et al. 1997; Ferguson 1999). Additionally, the optimal melt model for a given situation may depend on data availability. Even though a full energy balance model may 15 provide the best results when applied using all of the input data, for example, if the data required to run the model are not available, and most of the terms must be estimated does the energy balance model still provide the best simulation? The trade off between investment in the model and performance improvement should also be considered. Obtaining the data required to run an energy balance model can be costly in terms of manpower, equipment, and computational time. Several studies that compared the performance of energy balance and temperature index models at point and small basin scales concluded that their performance was similar (Kustas et al. 1994, Kane et al. 1997, Pellicciotti et al. 2008, and Singh et al. 2009). The cost of using an energy balance model was deemed too high relative to the minimal improvement in performance obtained in a paper by Zappa et al. (2003). Many snowmelt models were initially developed on glaciers (e.g. Hock 1999, Pellicciotti et al. 2005), which have a different energy balance than the mountainous and forested interior BC. Will these models perform as well for sites with variable slope, aspect, and forest cover conditions? The goal of this study is to assess the performance of melt models of varying complexity to simulate snowmelt under different aspects, forest cover, and input data availability conditions. The specific objectives of this study are to: (1) Determine the optimal melt model complexity to employ for runoff modelling at the plot scale; (2) Determine the effect of data availability on melt model performance; (3) Provide stand level snowmelt and radiation data for runoff model testing and to constrain model parameter ranges. 16 These objectives will be addressed within four chapters of this thesis. Chapter 2 describes the relevant site information and data collection and analysis methods. Chapter 3 compares four different temperature-index models with varying complexity and data requirements. Snow water equivalent (SWE) data collected through snow surveys and an automatic weather station at my research site provide the input, calibration, and validation data. SWE data from an automatic snow pillow is used for further validation. Meteorological, snow depth, SWE, and melt data are discussed with respect to aspect and forest cover. The performance of a simple energy balance model is discussed in Chapter 4. The performance of the model at my research site is assessed and is also compared to the performance of the temperature-index models. The advantages and disadvantages of temperature-index versus energy balance models are briefly discussed. Chapter 5 provides a summary of the findings in this study as well as suggestions for future work. 17 2. Study Site and Data Collection 2.1 Study Area The study area is located within the Willow River Basin in the Southern Interior hydrological zone of British Columbia, Canada (Figure 2.1). The Willow River Basin covers an area of 2870 km2 with 40% old forest, 2 1 % young forest, 34% logged or burned, 3% wetland, and 5% combined open areas including lakes, alpine, and agricultural zones based on Landsat imagery from the early 1990s (Environment Canada, 2002). The basin comprises predominantly sub-boreal spruce (67%) and secondarily Engelmann spruce - subalpine fir (32%) biogeoclimatic zones (Environment Canada, 2002). The basin's mean elevation is 1235 m above sea level (a.s.L), and elevation varies from 800 m to 2200 m. (Environment Canada, 2002). The Willow River flows north and enters the Fraser River east of Prince George. Mean annual temperature at the Environment Canada station near BarkerviUe (1283 m elevation) for the period 1971-2000 is 1.9°C and mean annual precipitation is 1014 mm with almost half (481 mm) falling as snow (Environment Canada, 2009). Between November and March the average monthly temperature varies from -8.8°C in January to 2.9°C in March. From April to June (Table 2.1) average monthly temperatures range from 1.4°C in April to 9.8°C in June (Environment Canada, 2009). 18 122"' W 122" W 122° W 122° W 122° W Figure 2.1 Study area location. N = north aspect, E = east aspect, S = south aspect, W = west aspect, O = open. Satellite image from ESRI (2009, hUp://goto.arcgison1inc.com/maps/ESRi Imagery World 2D). Map produced by T. Mlynowski. 19 Table 2 1 An temperature and precipitation data fiom the Environment Canada stations at Quesnel (Q) and Barkerville (B), BC, and from the Mt Tom (MT) meteorological station for the 2008 winter/spnng season The 30 year normals (1971-2000) monthly temperature and precipitation for Quesnel and Barkeiville stations are italicized Values in bold are at least 2°C cooler than the appropnate station climate normal Year (stn.) 2008(MT) ab Avg AT Avg Max AT Avg Min AT Extr Max AT Extr Min AT Tot Precip [°C] [°C] [°C] [°C] [°C] [mm] Month Jan na na na Na Na na 0 Jan -9 8 -4 8 -14 7 20 -37 0 1150 2008(Q) Jan -8 8 -3 4 -14 1 58 -33 1 29 0 1971-2000(B) Jan -8 8 -4 0 -13 6 10 0 -46 7 99 0 1971-2000(Q) Jan -8 6 -4 3 -12 8 13 9 -46 7 45 5 2008(MT) Feb -3 8 01 -7 2 41 -13 0 na C Feb -5 9 -0 3 -11 6 65 -27 5 108 0 2008 (Q) Feb -4 8 1 7 -11.2 104 -24 5 39 5 1971-2000(B) Feb -5 9 -0 5 -11 3 15 0 -43 3 64 3 1971-2000(Q) Feb -4 2 08 -9 1 15 1 -42 2 24 4 2008(MT) Mar -3 8 03 -7 4 21 -8 3 na 2008(B) 0 Mar -4 8 13 -10.8 70 -18 5 93 5 2008(Q) Mar 04 66 -5 7 11 0 -13 1 27 0 1971-2000(8) Mar -2 9 29 -8 7 172 -37 2 67 6 1971-2000(Q) Mar 11 70 -4 8 22 3 -38 9 28 9 2008(MT) Apr -2 3 20 -6 5 83 -11 2 na 0 Apr -2 9 3.3 -9.0 135 -23 0 25 1 d 2008(Q) Apr 32 9.1 -2 8 20 8 -8 9 150 1971-2000(B) Apr 1 4 74 46 27 8 -26 1 56 0 1971-2000(Q) Apr 63 13 4 -1 0 31 0 20 0 21 9 2008(MT) May 62 107 22 16 1 02 70 1 a 2008(B) 0 May 56 12 0 -0 8 23 0 -8 0 62 3 2008(Q) May 11 7 188 46 28 5 -2 6 64 5 1971-2000(B) May 62 12 6 -0 1 31 5 -75 0 69 9 1971-2000(Q) May 10 9 18 3 35 36 5 -10 0 40 7 2008(MT) Jun 84 136 37 26 8 -0 4 73 7 a 2008(B) 0 Jun 79 148 1.0 29 5 -3 5 62 7* 2008(Q) Jun 137 20 7 67 31 1 20 36 5 1971-2000(B) Jun 98 16 1 35 32 2 -6 7 101 9 1971-2000(Q) Jun 14 3 21 2 73 35 6 -3 3 68 6 2008(B) 2008(B) 2008(B) J Sohd piccipitation data collected at Mt Tom aie not ichable due to icing issues bJanuary an tempeiatuie values were not available due to mstiument failure c2008 Baikcivillc data have not yet undergone official QA check 'Missing one data point in calculating monthly total 20 The study area plots are located within the BC Ministry of Forests and Range Mt Tom Adaptive Management Trial (MTAMT) area. The Mt Tom site is located approximately 30 km north of Wells, BC, and encompasses an area of approximately 3200 ha of old-growth Engelmann spruce (Picea engelmanii) and subalpine fir {Abies lasiocarpa) forest crossing both the Engelmann spruce-subalpine fir Cariboo wet cool variant (ESSFwkl) and Cariboo wet cold variant (ESSFwc3) biogeoclimatic zones and includes a small portion of alpine area. The MTAMT is the third stage of the Quesnel Highland Alternative Silvicultural systems project, the goal of which is to develop sustainable silvicultural systems that preserve mountain caribou habitat while permitting timber harvesting within these high elevation forests (Armleder et al. 2002; Waterhouse 2002). Approximately 1200 ha were cut within the MTAMT area between 2001 and 2006. The method used was group selection harvesting on an 80 year cutting cycle with a maximum of 33% of the area harvested. Openings range from 0.1 - 1.0 ha with a few 3.0 ha clear-cuts and are spread across a range of elevations and aspects within the MTAMT area. The Mt Tom Adaptive Management Trial involves a number of projects including lichen abundance, windthrow, natural and planted stock regeneration, and snow accumulation and melt studies (Armleder et al. 2002). Patrick Teti1 with the BC Ministry of Forests, Mines, and Lands established two snow survey plots in the MTAMT area in 2000 to monitor changes in snow melt and runoff resulting from the proposed harvesting. The two plots are located within CP242-4 on the eastern side of the MTAMT area (Figure 2.2). The elevation of CP242-4 ranges from 1450 m - 1675 m encompassing two hills and straddling the headwaters of Wiley and Two Bit creeks (Waterhouse 2002). The north half of the block was classified as ESSFwkl with ' Emeritus Scientist, British Columbia Ministry of Forests, Mines and Lands, Williams Lake, BC 21 moderate slopes and predominantly south, and warmer aspects. The south half of the block was classified as ESSFwc3 with steep slopes and predominantly north aspects. Engelmann spruce and subalpine fir are the predominant tree species, but a few lodgepole pine {Pinus contortd) are also present (Gabriel 2004). The stand CP242-4 was likely established after a fire approximately 200 years ago and has experienced some partial disturbances throughout its development (Gabriel 2004). Soils are moderately drained, fine-to-medium textured silty loams, and formed on till (Waterhouse 2002). Starting in 2001, pre-harvest snow depth and density data were collected in each plot at 10 m intervals on each aspect along six east-west lines 400 m long and 40 m apart. Measurements were taken approximately every 10 days starting at the peak snow water equivalent (SWE) in March until the melt was complete in June. This area was logged in the winter of 2005/2006 creating openings ranging from 0.2 1.2 ha in size, and replanted that same year. One year of post-harvest data were collected within the same plots in the spring of 2007. Fish eye photographs were also taken at each snow measurement point to characterize canopy closure (Teti 2010, unpublished data, see section 4.2.3). The two main research plots established for the current study are located adjacent to the BC Ministry of Forests and Range's snow survey plots (Figure 2.2). See Table 2.2 for detailed site characteristics. The north plot is in opening D with a south aspect and the south plot is in opening DDD with a north aspect. The meteorological station is located in opening G, on the south side of CP 242-4 east of opening D. Opening D has dimensions of 112 m x 112 m (equating to 5-6 times the surrounding tree height (H)), on an approximate 16% slope (Table 2.2). There is a small patch of short trees at the centre of the plot, but there is little woody debris. Two small streams run through the opening and still had flow at the end of 22 October 2007. Opening DDD is oblong-rectangular with dimensions of 125 m x 75 m (6H x 3H), on an approximate 22% slope (Table 2.2). The terrain is hummocky with complex micro topography and some large shrubs. There are many large slash piles in this opening. Opening G is rectangular with dimensions of 175 m x 75 m (6H x 2.5H), on an approximate 5.1% slope (Table 2.2). It is fairly flat at the top end and is the opening closest to centre between the two research plots. The bottom of the opening drops off steeply and there are a few large slash piles. The meteorological station is located on an approximate 3° slope. The south aspect forest site is located to the east of opening D and the north aspect forest site is located to the west of opening DDD. Based on data from the adjacent Ministry of Forests, Mines and Lands sites the forest cover fraction for the south aspect forest site was 0.60 and the forest cover fraction for the north aspect forest site was 0.46. Two additional openings are used in the study; an east aspect (BB) and a west aspect (J) (Table 2.2). Data from the Ministry of Environment Automatic Snow Pillow (ASP) at BarkerviUe are also used. The snow pillow is located approximately 20 km southwest of the Mt Tom sites in an open area at an elevation of 1483 m (Table 2.2). 23 CP 242-4 Block Map •H CP m-a Onv OKI-fa ^Jptn' c w t U'jm ^Oia. Legend C.iPinc-MA1hn3o« M i f f CfBftf ^GWS S^ta I M M 3>».8) 3nuh tali ^ CP 242-4 lAtfugal ^ 2006 r— CJ V^ Figure 2.2 Block map of Mt Tom study area. Opening D, DDD, BB, and J were used in the snowmelt study. The meteorological station was located in opening G. Map provided by Patrick Teti, BC Ministry of Forests, Mines and Lands. 24 Table 2.2 Mt Tom study site characteristics. Tree height refers to the approximate average height of trees surrounding the opening. Tree height was only measured for the main north and south open sites and the meteorological station site. Due to their proximity, adjacent north and south aspect forest sites were assumed to have the same average tree heights as measured for the respective open sites, 'na' indicates data that were not available. Site Opening Aspect Cover Elev. (m) Slope (deg) Size (ha) Tree ht(m) Latitude (°N) Longitude (°W) SO D south open 1469 9.25 1.2 26 53°11.53' 121°40.18' SF - south forest 1469 9.25 na -26 53°11.53' 121°40.18' NO DDD north open 1520 12.47 1.0 25 53°11.30' 121°39.44' NF - north forest 1520 12.47 na -25 53°11.30' 121°39.44' Met G south open 1487 2.95 1.2 29 53°11.54' 121°39.82' WO BB west open 1427 17.00 0.7 na 53°11.57' 121°40.52' WF - west forest 1427 17.00 na na 53°11.57' 121°40.52' EO J east open 1589 4.80 0.8 na 53°11.75' 121°39.53' EF - east forest 1589 4.80 na na 53°11.75' 121°39.53' ASP - none forest 1483 0.00 -0.0064 na 53°03' 121°29' 25 2.2 Data Collection 2.2.1 Meteorological Station Figure 2 3 Meteoiological station at Mt Tom. Photo taken 9 March, 2008. A meteorological station (Figure 2.3) was set up at a level location in Hepving G approximately 400 m east of the south aspect snow course sites and 625 m north-west of the north aspect sites (Figure 2.2). Air temperature (°C), relative humidity (%), wind speed (m s"1) and direction (°), incoming shortwave radiation (W m" ), and snow depth (cm) were recorded every 10 seconds and averaged every 15 minutes beginning 27 January 2008 (Table 2.3). Daily geometric mean wind speed was calculated to be consistent with methodology in Walter et al. (2005). Solid and liquid precipitation (mm) were totaled daily and precipitation data were also available throughout the study at the Environment Canada station at 26 Barkerville (1283 m) and the River Forecast Centre automatic snow pillow (1483 m) both within ca. 30 km of the Mt Tom site. Table 2.3 Description of instrumentation at Mt Tom meteorological station. Climate Variable Instrument air temperature Vaisala/Campbell Scientific HMP45C212 probe relative humidity Vaisala/Campbell Scientific HMP45C212 Mounting Height (cm) 209 Accuracy ±0.1° C at 20°C: 209 ±2% (0-90%) +3% (90-100%) incoming shortwave radiation Kipp & Zonen CMP-3 thermopile pyranometer wind speed and direction RM Young 05103-10A wind monitor snow depth Campbell Scientific SR50A sonic ranging sensor 230 Texas Electronics TE525M* with CS705 precipitation adapter ca. 215 (top of adapter) precipitation ca. 210 290 ±5%(-10°C-40°C) +10% (daily sums) speed: ±0.3 m s"1 or 1% of reading direction: ±3° depth: max. of: ±1cm or 0.4% of distance to target temp.:<±0.75°C ±1% for rainfall up to 10 mm h"1 *Piecipitation gauge was unshielded and mounted to tower mast 2.2.2 Aspect Pyranometers Incoming shortwave radiation was also measured at both the south and north aspect open sites with Kipp & Zonen CMP-3 thermopile pyranometers. A single pyranometer was set up near the centre of each opening and was mounted parallel to the respective average site slope to measure radiation incident to the surface. Average incoming shortwave radiation (W m"2) was measured every minute, and 15 minute averages were logged on Campbell Scientific data loggers. Data were collected between 2 February 2008 and 5 June 2008. 27 2.2.3 Albedo Two LiCor LI200S silicon pyranometers were set up to measure albedo adjacent to the meteorological station beginning 9 May 2008. Albedo data were collected until all the snow had melted at the station. Incoming and reflected shortwave radiation was recorded every 10 minutes with a CR10 data logger. Measurement of albedo using silicon pyranometers is not recommended without a comparison correction to thermopile pyranometers (Henneman and Stefan, 1998). This correction is required since silicon pyranometers only sub-sample the shortwave radiation spectrum (350-1000 nm), and are calibrated to predict all of the solar radiation (280 to 2800 nm). Due to equipment restrictions it was neither possible to use thermopile pyranometers to measure albedo nor to compare the albedo measurements made with the silicon pyranometers with the thermopile pyranometers. It is expected that albedo measured with silicon pyranometers will be overestimated (Henneman and Stefan, 1998). 2.2.4 Snow Surveys Snowmelt was measured using four methods: weekly snow surveys of depth and SWE throughout the melt season, continuous snowmelt lysimeter measurements of runoff from the snowpack, total melt over the season using temperature sensors at the ground surface, and hourly automatic snow pillow measurements from a nearby station. Snow survey SWE and depth data were collected ten times between 1 April 2008 and 5 June 2008 at four sites: north aspect open, north aspect forest, south aspect open, and south aspect forest (Figure 2.2). A snow pit was dug on the south aspect on 9 March 2008 and on the north aspect on 18 March 2008 to obtain temperature and density profiles. 28 At each site, a snow course transect was made up of four measurement points along an east-west line. There were three transects per site, giving a total of 12 sampling points per site. The sampling points were evenly distributed across the open area, resulting in different spacing at each aspect (Figure 2.4 and Figure 2.5). To reduce edge effects, transects started at least 16 m away from the forest edge for all sites. Each snow course was at least 100 m long. Snow course length was based on the relationship between the coefficient of variation (CV) and transect length for Ministry of Forest and Range's 2007 snow depth data (Figure 2.6). For the south aspect open site, measurements were taken starting at 16.5 m away from the forest edge on both the east and west side of the opening at an east to west spacing of 20 m and a north to south spacing of 10 m (Figure 2.4). For the north aspect open site, measurements were taken starting at 16 m away from the forest edge on both the east and west side of the opening at an east to west spacing of 10 m and a north to south spacing of 29 m (Figure 2.5). For the south aspect forest site, measurements were taken starting at 17 m inside the forest on the east side of opening D at an east to west spacing of 20 m and a north to south spacing of 7 m (Figure 2.7). For the north aspect forest site, measurements were taken starting at 20 m inside the forest on the west side of opening DDD at an east to west spacing of 10 m and a north to south spacing of 29 m (Figure 2.8). Measurements were taken approximately 80 cm from the marker. At each measurement point, a snow core was taken at arm's length in the predetermined direction from the marker. Each measurement within a site was made in the same direction for a given sampling date. The direction changed for each sampling date to prevent sampling in a previous core hole. If an obstruction was encountered another sample was taken adjacent to the first until successful. If a forest course sample landed in a tree well it was still taken, but a note was made. 29 South Aspect Open @ *£ — oo Core point Lysimeter Road MFR survey boundary 7m N 15m A W-«- -•E 14 m ^JIV t & t 7m 3 m •* 38 m • Forest Course X 3 5 nt 20 m Iff ;imete Sip** Lyiimete ^B 1 3 5 m 48 m 20 m 112m Figure 2.4 Diagram of south aspect open snow course site. The location of snowmelt lysimeters, pyranometcr, and snow survey sample points are shown. Tree symbols denote the presence of trees, but do not represent the size or number. 30 North Aspect Open @ Core point $• Lysimeter ---• Road oo MFR survey boundary 7m S • E-4- *!m* ->W 3 m Forest Course 29 m 18m Ly simetef B Ly >imete(-A 21 m 2m -> •*— 8.5 m 22 m 10m 10m 67 m Figure 2.5 Diagram of north aspect open snow course site. The location of snowmelt lysimeters, pyranometer, and snow survey sample points are shown. Tree symbols denote the presence of trees, but do not represent the size or number. 31 o TO & Ol A. A' O O O ^—i A -^. o south aspect open -°south aspect forest - A o o north aspect open north aspect forest ~I 1 I 50 100 150 Distance from start (rn) Figure 2.6 Relationship between coefficient of vaiiation (CV) and distance fiom start of snow couise for Pat Tcti's 2007 snow depth data. SWE data were collected on 21 Maich for the south aspect sites and 31 Maich for the north aspect sites. Vcitical dotted line indicates twelve sample points length. Distance from stail only indicates the cumulative distance fiom the stait of each snow course. The distance from the forest edge to the start of each snow course was not necessarily the same for each site, and these data were not available. 32 Figure 2.7 Diagram of south aspect forest snow course site. The locations of snow survey sample points are shown. Tree symbols indicate the edge of the forest on the south side. The southeast-most sample point fell on the forest edge. The trees were distributed throughout the rest of the snow course, but omitted for clarity in this diagram. 33 North Aspect Forest © Core point - - - Road o o MFR survey boundary S • 29 m <4i 20m H E<«- _&2 .BLOCK DDD 10 m Figure 2.8 Diagram of north aspect forest snow course site. The locations of snow survey sample points are shown. Tree symbols indicate the edge of the forest on the north side. The trees were distributed throughout the entire snow course, but omitted for clarity in this diagram. The empty weight of the Federal snow sampler was recorded at the start and end of each snow course. If snow was freezing within the tube the weight was recorded after each sample. Snow depth was measured against the outside of the tube while it was inserted. The length of the soil plug was subtracted from this measurement to obtain snow depth. If there was no soil plug the core was retaken. After removing the soil plug the snow core was weighed. Soil plugs were dumped at the marker or outside of the snow course boundary to avoid accelerated melting around the debris. Snow survey data were used to generate a time series of SWE for each site throughout the melt season and provided snow melt data for comparison with model simulations. 34 -•W 2.2.5 Snowmelt Lysimeters Two snowmelt lysimeters were installed in each open aspect site (Figure 2.4 and Figure 2.5). Snowmelt lysimeters were not installed in the forested sites because the high spatial variability of snow accumulation and melt within forested areas could not be captured without the installation of a prohibitively large number of lysimeters (Kattelman 2000; Storck et al. 1999). Each lysimeter had a 2.97 m2 collection area. In a review of lysimeter studies, Kattelman (2000) concluded that when the collection area of the lysimeter is less than 2 m , more than three snowmelt lysimeters should be deployed at each site to accurately sample runoff. The lysimeters used in this study needed to represent snowmelt runoff variability over an area of approximately 112 m x 112 m and 125 m x 75 m. Thus the two lysimeters together represented slightly less than 0.05% of the opening area. Since the snow courses cover a smaller area within the centre of each opening (60 m x 28 m for the south aspect and 30 m x 58 m for the north aspect), the two lysimeters represented approximately 0.35% of the snow course area for each site. The purpose of two lysimeters at each site is to provide replication and to provide backup if one of the lysimeters fails. The lysimeters were constructed on advice from Patrick Teti. Each lysimeter collection box is a 1.22 m x 2.44 m piece of plywood with a support base and 15.2 cm tall sides. The plywood base is flush to the ground surface and anchored in place with rebar. The box is lined with a grey tarp and water drains from the bottom centre of the box through a length of buried 3.8 cm ABS pipe into a tipping bucket rain gauge located down slope of the box. The tipping bucket rain gauges used were designed by Patrick Teti and custom built to handle the large volume of water runoff during snowmelt. Each bucket has a tip volume of-440 mL (Table 2.4). Calibration was performed for each side of each tipping bucket at two different flow rates to 35 simulate flow volumes likely experienced during runoff. Each box was surrounded by 3' tall wire fence erected at a distance of 3 ' from all sides of the box to keep animals away from the plywood. The tipping buckets were fitted into wooden frames to protect them from damage from the snow and were surrounded with hardware cloth to protect them from damage by animals. In mid-March, the tipping buckets were dug out, the ice cubes were removed from the buckets, the data loggers were tested, and the holes were filled back in to keep the buckets from freezing. Data were downloaded from the lysimeters in late May once the snow had melted from the lysimeter pan. Table 2.4 Average tip volume of the tipping bucket for each snowmelt lysimeter. Lysimeter Location Tip volume (mL) SA west D 451.8 SB east D 425.8 NA west DDD 421.2 NB east DDD 465.6 2.2.6 East and West Aspect Total Season Melt In addition to the north and south aspect snow surveys, the total melt and average melt rate between 27 April 2008 and the end of the melt were estimated at an additional four sites at Mt Tom: east aspect open, east aspect forest, west aspect open, and west aspect forest. HOBO TidBit v2 data loggers (Onset Computer Corporation, 2010) were used to estimate snow free date at the east and west aspect sites. This was then used to calculate the total and average daily melt. Melt was estimated at two locations within each forest or open aspect site for a total of eight estimates. Opening J and BB (Figure 2.2) were divided into four quadrants. Sensors were deployed at the centre of the bottom left quadrant and the 36 centre of the top right quadrant. The forest sites were located to the east side of each opening. Sensors were deployed at two points in the forest randomly located at least 12 m beyond the forest edge. At each location a SWE measurement was taken with the Federal sampler. The TidBit, a small, waterproof temperature sensor, was then placed at the bottom of the core hole and the core of snow was returned to the hole to cover the logger. Once the snow had melted in the spring the TidBits were collected and the temperature data extracted. As soon as the snow melted, exposing the TidBit, the temperature record showed clear diurnal cycles. This date was used to determine the amount of time to melt the total initial SWE for each site. 2.2.7 Automatic Snow Pillow The automatic snow pillow at Barkerville (1A03P) is located ca. 19.8 km southeast from the Mt Tom snow course sites in a small, level opening, and has been maintained at this location since 1968. The opening is approximately 8 m in diameter and surrounded by mature forest (Scott Jackson pers. comm.). More detailed information on tree height and stand characteristics were not available. Quality checked hourly meteorological data are available through the River Forecast Centre up to the end of 2006 (River Forecast Centre, 2009). Hourly air temperature (°C), accumulated precipitation (mm), snow depth (cm), and SWE (mm) were obtained for the period 1 January - 31 May 2008. These data have undergone preliminary quality assurance, but were suitable for the purpose of this study. See Appendix 1 for further details on the snow pillow station and measurements. 2 Technician, River Forecast Centre, British Columbia Ministry of Environment, Victoria, BC 37 2.3 Data and Analysis Methods 2.3.1 Meteorological Data The 15 minute average Mt Tom meteorological data underwent routine quality control analysis. Spikes in SR50 snow depth data and missing or en-oneous air temperature, relative humidity and shortwave radiation readings were replaced with the average of the measurement that preceded and followed the errant reading. The data were mostly error-free other than noise in the snow depth data. Daily averages of shortwave radiation and air temperature and daily sums of precipitation were calculated from the 10 minute meteorological data. Shortwave radiation averages incorporated the full 24 hour period. The precipitation gauge at the Mt Tom meteorological station was neither correctly mounted nor shielded due to time and funding constraints. The precipitation adapter was kept topped up with antifreeze throughout the winter. Issues with undercatch, evaporation, freezing, and recording lag (MacDonald and Pomeroy, 2007) prohibited the use of any winter precipitation data from this station (Appendix 2). Cumulative precipitation is also available at the Barkerville ASP, but the gauge is unshielded and there is no wind speed monitor on-site to enable calculation of under-catch corrections. Precipitation data from the Barkerville Environment Canada station located at 53°4.200' N, 121°30.600' W and 1283 m elevation were used for this study. 2.3.2 Snow Course Data Snow density (kg m"3) was calculated from the depth and SWE measurements as: swe Ps=—r-*P„ a.. 38 (2-1) -1 Where ps is snow density [kg m" ], swe is snow water equivalent [m], ds is snow depth [m], andpw is density of water [kg m"3]. Average SWE was calculated for each site and each date. As snow melted from the sampling points the point SWE was set at zero. On subsequent sample dates these points were still recorded as zero SWE. To account for differences in snow disappearance dates between survey stations I calculated the average site melt by first calculating the change in SWE at each sample point and then averaging the resulting values for each site and sample period. On 20 April, 2008, the measured SWE values in the openings were anomalously low (Figure 2.9a). The day was characterized by cold (-15°C), sunny conditions. Due to the freezing air temperatures and high solar radiation the snow was freezing inside the Federal snow tube between samples, an effect that was worse at the open sites. The snow corer was zeroed after each sample, but anomalous readings still resulted. To replace the anomalous measured values a linear model was created for each of the four sites using the snow density values from the two sample dates before and the two dates following 20 April, 2008. These models were used to interpolate snow density for 20 April, 2008 for each site and the SWE was back calculated for each site using the interpolated snow density and the measured snow depth (freezing did not affect depth measurements). Resulting corrected SWE values were more reasonable and were used for modelling purposes (Figure 2.9b). 39 — I Apr 01 I 1 I 1 I 1 April Apr21 May01 May 11 May21 May31 — I j , \ Apr01 1 1 1 1 I 1 April Apr21 May01 May 11 May21 May31 Figure 2.9 Mean SWE measured at each snow course site throughout the sampling period in 2008. S = south aspect, N = north aspect, O = open, F = forest, a) data before correcting for SWE measurement errors at the open sites on 20 April; b) data after correction applied. Due to the difficulty in determining precipitation at the Mt Tom site (Appendix 2), the snow course data used for melt modelling were chosen so that the SWE data used from each site were only for the period of sustained melt. It was assumed that once melt was occurring the snowpack was isothermal at 0°C and any liquid precipitation would run off rather than be retained in the snowpack. Based on the precipitation data taken from the Barkerville Environment Canada station data no snowfall occurred after 27 April, 2008 (Appendix 2 and 4). The contribution of snowfall to SWE should be accounted for in the calculation of melt for studies with significant accumulation during the simulation period. The melt period for each site was chosen based on the time series of SWE (Figure 2.9b) such that SWE only decreased at each site between sample periods. Change in SWE was calculated as the difference between the average site SWE on subsequent sampling dates. Start and end dates for melt periods varied among the snow course sites as a result of variable melt rates (Table 2.5). 40 Table 2.5 Start and end dates for snow course measurement periods in 2008 for each snow course site. S = south aspect, N = north aspect, O = open, F = forest. * north aspect forest data were only compared until 30 May, 2008. SF SO NO/NF Period Start Date End Date Start Date End Date Start Date End Date 1 27-Apr 02-May 02-May 09-May 09-May 17-May 2 02-May 09-May 09-May 17-May 17-May 23-May 3 09-May 17-May 17-May 23-May 23-May 30-May 4 17-May 23-May 23-May 30-May 30-May* 05-Jun 2.3.3 TidBit Data At the Mt Tom TidBit sites the total SWE recorded at each site on 27 April, 2008, was assumed to equal the total melt between that date and the time of snow disappearance. The SWE was averaged for the two TidBits at each site. The end of the melt period for each TidBit pair was determined from the temperature record and varied for each site (Table 2.6). The snow melt was considered complete once a diurnal cycle appeared in the TidBit temperature record. The end date was the same for each of the two TidBits at a site except for the east aspect open site. One TidBit melted out on 25 May, 2008 and the other TidBit melted out on 29 May, 2008. The average of the two end dates was used (27 May, 2008). 41 Table 2.6 Start and end dates for melt periods in 2008 for each TidBit site. W = west aspect, E = east aspect, O = open, F = forest. Start Date End Date Site Initial SWE [mm] EO 475 27-Apr 27-May EF 325 27-Apr 26-May WO 305 27-Apr 17-May WF 125 27-Apr-2008 16-May-2008 2.3.4 Automatic Snow Pillow Data To calculate daily precipitation and change in SWE at the Barkerville automatic snow pillow (ASP), the value at midnight was subtracted on subsequent dates. This method was recommended by the River Forecast Centre to reduce the effect of noise in the data. Daily average air temperature was calculated from the hourly data measured at the ASP. Due to the small opening size and surrounding forest the Barkerville ASP was treated as a forested site for modelling. 2.3.5 Lysimeter Data Unfortunately, the lysimeter data collected at the end of the 2008 melt season did not match with observed melt at the sites. The lysimeters did not record runoff until only a few days before the snow had completely melted from each site and the total volume recorded equated to less than half the SWE that was measured on site. It is suspected that the drain covers froze over preventing runoff from entering the tipping bucket and being measured. The melt water likely overflowed the pan during the time the drain was frozen. 42 3. Field testing and comparison of temperature-index melt models of varying complexity 3.1 Introduction Modelling is a powerful tool to predict how changes in land use will affect surface runoff. The HBV model is a conceptual precipitation-runoff model that is currently used around the world for flood forecasting, hydroelectric generation, and climate change studies (Bergstrom 1976). The HBV model, like many hydrologic models, uses the empirical relation between air temperature and snowmelt (the degree day approach). The advantage of a degree-day approach to model snowmelt runoff is that these models typically require only air temperature as input data. As a result, the degree day approach is ideal for operational applications or for studies which examine hydrologic scenarios for many basins. The physical basis of temperature index models is the relationship between air temperature and the energy exchanges driving snowmelt. Braithwaite (1995) observed a fairly strong correlation between temperature and ablation (r=0.78), although there was considerable scatter in the relationship. Zuzel and Cox (1975) showed that air temperature was the best index for melt. Air temperature is successful as the sole index for melt because of the correlation between temperature and several other energy balance components. Air temperature is correlated with incoming shortwave radiation; a good correlation between melt and solar radiation (i-0.75) was observed by Anderson (1968). A review by Ohmura (2001) of energy balance studies at multiple glacier sites concluded that incoming longwave radiation made up the greatest percentage of energy input and together incoming longwave and sensible heat accounted for approximately 75% of energy input. Both of these fluxes are 43 strongly influenced by near-surface air temperatures hence providing an explanation for the close relationship between air temperature and melt (Ohmura 2001). To further investigate this relationship, atmospheric longwave radiation was modelled using profiles from radiosonde measurements. About 66% and 88% of the total atmospheric longwave radiation was found to originate from the first 100 m and 1 km of the atmosphere, respectively. Air temperature measured at standard screen height is representative of the first 100 m of the atmosphere and the conditions in the next 100 m are highly correlated with the first. Thus, air temperature is a good indicator of melt because the first 100 m of the atmosphere contribute more than half of the incident longwave radiation and incoming longwave radiation makes up the largest contribution to the overall energy budget (Ohmura 2001). In practice the amount of melt over a given time period is related to air temperature by a factor of proportionality. Melt can be related to average air temperature or to positive temperature deviations (degree days). Either method yields accurate results; the success is dependent on the appropriateness of the melt factor or degree-day factor chosen. A constant melt or degree-day factor (Mp, DDF) was originally used for most temperature index models. This factor was assumed to account for average basin melt. However, there is a large body of evidence that documents the spatial and temporal variability of Mp and provides grounds for the use of a varying Mp (e.g. Kuusisto 1980; Lang and Braun 1990; Braithwaite 1995). The melt factor is known to vary with vegetation cover, precipitation (humidity), snow density, albedo, and daily and seasonal cycles (Kuusisto 1980). Rango and Martinec (1995) and Hock (2003) give good reviews of melt factor variability and the influencing factors. 44 A daily time step is most commonly used for calculating melt. However, the use of daily average air temperature may introduce error when the temperature fluctuates around 0°C. There may have been periods during the day when temperatures were conducive to melt even though the daily average temperature was below 0°C. The threshold temperature for melt is often assumed to be 0°C (Hamlin et al. 1998; Sloan et al. 2004). However, this is not necessarily true. Kuhn (1987) showed theoretically that melt could occur at temperatures as low as -10°C and may not occur at temperatures up to +10°C. Melt was observed at temperatures at and below 0°C by Braithwaite (1995). Since air temperature is only one of several factors affecting snowmelt neither the melt factor nor the threshold temperature can be expected to remain fixed (Kuusisto 1980). These issues have been addressed to varying degrees through research on snowmelt and energy balance processes and with modification to the classical temperature index model. The relation between temperature and melt was recognized as early as 1887 and was used to simulate melt for models beginning in the 1930s (Rango and Martinec 1995; Hock 2005). The temperature index approach is still the most widely used approach to melt modelling today and dominates in operational runoff models (Bergstrom and Forsman 1973; Quick and Pipes 1977; Martinec and Rango 1986; Micovic and Quick 1999; Sorman et al. 2009). The strength of the temperature index model is in simulating melt on a watershed scale over longer time periods in areas with little meteorological data. For early models it was observed that melt was sufficiently simulated at a daily or longer time step, but these models were not recommended for simulations at a shorter time step (WMO 1986; Rango and Martinec 1995). 45 There has been a number of developments that improve the ability of the temperature index model to represent the snowmelt process. Almost all models employ a seasonally variable melt factor to capture the changes in insolation and albedo throughout the melt season. The UBC Watershed model uses a monthly variable melt factor (Quick and Pipes 1977) and the HBV model uses a sinusoidal melt factor function that increases from a minimum in December to a maximum in June (Lindstrom et al. 1997). Attempts have also been made to improve the physical basis of temperature index models by incorporating additional variables to approximate additional energy terms important for snowmelt. Zuzel and Cox (1975), for example, concluded from a statistical analysis that the additions of wind speed, vapour pressure, and net radiation would improve melt estimates. Dunn and Colohan (1999) added a number of factors including wind redistribution, slope and aspect classes, albedo, and snow density. Anderson (1973) introduced a straightforward modification by using the temperature index approach for nonrainy conditions and a simplified energy balance during rainy periods. Cazorzi and Dalla Fontana (1996) simply used a different melt factor during precipitation events. A radiation or energy balance component has been incorporated into a number of temperature index models and has generally been seen to improve model performance (Kustas et al. 1994; Cazorzi and Dalla Fontana 1996; Brubaker et al. 1996; Hamlin et al. 1998; Hock 1999; Verbunta et al. 2003; Kling et al. 2006). However, even in recent years the temperature index method is still employed in its simplest form (e.g. Daly et al. 2000; Woo and Thome 2006). Hock (1999) identified two major drawbacks of classical temperature index models. First, they have a restricted temporal resolution. The daily time step employed by most temperature index models is insufficient to capture the diurnal melt cycles. Second, they 46 have limited spatial variability in melt rates. The melt factor is usually assumed to be invariant in space, or perhaps different for discrete surface, vegetation, elevation, or aspect classes. The melt factor thus represents the average melt rates of the basin, but cannot take the spatial or temporal dynamics of the melt process into account. In complex topography melt rates can be highly variable (Male and Granger 1981; Olyphant I986a,b; Fierz 2003). Cline (1997) showed that melt would stop overnight as net energy fluxes became negative and would not begin again until the energy deficit in the snowpack had been compensated for the following day. The same observation was made by Singh et al. (2000) despite air temperatures being above 5°C throughout the study. A simple temperature-index method would likely overestimate melt in this situation, since with this method melt occurs as long as temperatures remain above the critical temperature (usually 0°C). In this situation a lower melt factor would be needed to correctly simulate the overall melt. This demonstrates one of the reasons why simple temperature-index models do not perform well for periods less than one day and are unable to capture diurnal cycles in melt (Rango and Martinec 1995; Hock 1999). Thus there has been a push for more distributed models that can describe melt at finer temporal and spatial scales, and for more physically based models that can belter account for the variability in the contribution of the individual energy balance components. Since the melt factor accounts implicitly for all the energy balance processes affecting snowmelt it can be highly variable in space and time. This variability has been cited as a concern (Kustas et al. 1994) as it is difficult to accurately represent the melt factor throughout the melt season. The melt factor is only accurate as long as the relative contribution of the energy fluxes remains the same. Current studies have been aimed at 47 extracting these energy balance processes from the melt factor to better represent the spatial and temporal variability for the basin. One of the current deficits of many common runoff models is that radiation is not used in calculating snowmelt runoff (e.g. CHC 2007). Shortwave radiation is an important component of the energy budget available to melt snow, particularly in mountainous terrain (Cline 1997). Under overcast conditions shortwave radiation is not as strongly related to air temperature and has a fairly poor relationship with melt. Further, the turbulent exchange processes are difficult to reliably quantify and assess (Fierz 2003). It therefore follows that radiation was the first to be incorporated as an index. Cazorzi and Dalla Fontana (1996) used monthly potential clear sky solar radiation maps constructed from digital elevation models (DEMs) and solar geometry to determine an energy index. This index could vary as a function of local topography. The modified model simulated the spatial and temporal distribution of SWE accurately across their basin of study whereas the classical temperature-index model was unable to reproduce the spatial dynamics of snow cover. Kustas et al. (1994) and Brubaker et al. (1996) added a net radiation term, which improved the performance of the new model over the original version. However the net radiation term needed to either be measured or modelled, for which measurements of additional meteorological variables and components of the radiation budget were necessary to obtain accurate estimates. Hamlin et al. (1998) also included a net radiation term, but found that it made little improvement over the original model. The lack of improvement was attributed to the lack of measured radiation data or good data for estimating the radiation terms. Hock (1999) developed a grid-based temperature-index model that accounted for the spatial and daily variability of melt rates by incorporating potential clear-sky direct solar 48 radiation. This information was estimated from DEMs and solar geometry and thus was a simplification over the net radiation indices added in previous studies as it did not require additional meteorological input. Also the hourly calculations gave better estimates of radiation than the monthly averages used by Cazorzi and Dalla Fontana (1996). The addition of the radiation term produced diurnal cycles, improved discharge values, and simulated a more realistic pattern of melt rates (i.e. increased melt on south-facing slopes) whereas the classical approach did not (Hock 1999). The addition of a radiation term to the temperature index model has been seen to considerably reduce the variation in the restricted melt factor used in these models (Kustas et al. 1994; Brubaker et al. 1996; Hock 1999), providing a significant improvement in the models. 3.1.1 Outline and Objectives The optimal melt model for a given situation may depend on data availability. For example, even though the addition of other measured meteorological variables improves melt forecasting, if the data required to run the model are not available and the terms must be estimated does the modified model still perform as well or does the simple temperature index perform better? As well, how does the performance of modified temperature index models compare? Although many studies compare the performance of a new modified temperature index model against observations or to the basic temperature index model, modified temperature index models are not as often compared against each other. Finally, the goal is to apply these melt models in operational runoff models at the watershed scale. Therefore the performance of these models under variable terrain conditions is an important component in their selection for incorporation into watershed scale runoff models. Many snowmelt models were originally developed and tested on glaciers (Hock 49 1999; Pellicciotti et al. 2005). The presence and structure of forest canopy affects the energy balance at the snow surface (Male and Granger 1981; Storck et al. 1999; Koivusalo and Kokkonen 2002; Pomeroy et al. 2008) which results in different rates and timing of snowmelt amongst stands (Storck et al. 1999). Slope angle and aspect affect the amount of incident shortwave radiation received at a point resulting in variability of shortwave radiation distribution across a basin as a function of time of day and time of year (Olyphant 1986a; Cazorzi and DallaFontana 1996) and differing melt rates (Bloschl et al. 19916). The majority of snowmelt modelling studies use river discharge data to calibrate and test model performance. However, running snowmelt models within a rainfall-runoff model smoothes the basin response and reduces the effect of an individual model component on the overall output (Kustas et al. 1994; Kane et al. 1997). Comparison of the snowmelt models at the small plot scale avoids the need to extrapolate data across a large area, which could introduce uncertainties not directly related to the performance of the model. Comparison of the snowmelt models independently from a runoff model removes the interference of parameters related to other watershed processes as well as the dampening effect of basin discharge. This study seeks to assess the accuracy of temperature-index melt models of varying complexity to simulate snowmelt under different aspects, forest cover, and input data availability conditions. The objectives of this study are to: (1) Determine the optimal temperature-index melt model complexity to simulate melt at the plot scale for mixed open and forest conditions (2) Assess whether the inclusion of shortwave radiation is beneficial for snowmelt prediction at a daily simulation time step 50 (3) Provide stand level snowmelt and radiation data for runoff model testing and to constrain model parameter ranges. 3.2 Methods and Analysis See Chapter 2 for detailed description of the study area and data collection methods. 3.2.1 Snowmelt Models Basic Temperature Index The basic temperature index (TI) model requires only daily air temperature and a calibrated melt factor (Mp) to predict melt. The basic TI approach usually takes the form: M = MF(Ta-T0) (3-1) Where M= the amount of melt (mm d"1) MF = the melt factor (mm °C"' d"1) Ta = daily average screen level air temperature (°C) To = threshold temperature below which melt does not occur (°C) M=0whenTfl _ \\ o \ o - * - Mt Tom north aspect open -°- Mt Tom north aspect forest -A_ MtTom south aspect open - & Mt Tom south aspect forest — Barkerville automatic snow pillow tD \ " Mt Torn snow course average i i i i Apr 01 April Apr 21 May 01 May 11 \° Barkerville snow pillow CD i - \ y_ — l May 21 May 31 Apr01 April Apr21 May01 May 11 May 21 Figure 3.2 Comparison of the daily Barkerville automatic snow pillow SWE measurements with: a) the measured SWE at each Mt Tom snow course site; and b) the averaged measured SWE for all the Mt Tom snow course sites Precipitation data are a large source of uncertainty in this study. Accurate precipitation data were necessary for several modelling steps. Measurement of solid precipitation in cold and remote locations is particularly challenging (Goodison et al. 1998). Precipitation measured at Mt Tom was unusable due to issues with undercatch, evaporation, freezing, and recording lag (i.e. MacDonald and Pomeroy, 2007) and extraction of snowfall amounts from the SR50 at Mt Tom was subjective and unsatisfactory (Appendix 2). In addition, precipitation data varied inconsistently between the gauge at the Barkerville ASP and measurement at the nearby Environment Canada station, which made the representativeness of these data to the study site questionable (Appendix 2). The precipitation data recorded at ASP and Barkerville Environment Canada station (EC) do not match well with observed SWE changes and snowfall events at Mt Tom (Appendix 2). It is particularly noticeable for the north aspect open site that the measured precipitation cannot account for the observed increase in SWE. This observed accumulation 63 May31 o is odd given that this increase is not seen at any other snow course site. It may have resulted from drifting snow as there were signs of drifting around some of the slash piles, but this was not monitored. Wind speeds at MT were greater than 1.5 ms" 1 for the period 1 2 - 1 7 April, 2008, and 26 April - 2 May, 2008, with the maximum daily geometric mean wind speed of 2.71 m s" measured on 12 April, 2008. Wind speeds were relatively low during the period of accumulation seen at the north aspect open site during 2 - 1 2 May, 2008. Also, fresh snow falling in the open on a north aspect might adhere to the surface. The recorded precipitation was also compared between ASP and EC (Figure 3.3a). The records differed although the cumulative precipitation was similar. Agreement between the records was better in April than in May. There were several occasions when precipitation occurred at one station, but was not recorded at the other or the precipitation differed by a considerable amount (i.e greater than double). Since the ASP precipitation gauge is unshielded it is expected that the ASP precipitation would be less than that at EC due to undercatch. This was not always the case. When the precipitation records were compared to the cumulative increase in SWE measured at the snow pillow, the SWE increased more rapidly than either of the precipitation records (Figure 3.3b). Since the snow pillow measured the net change in SWE, this record was not valid if any melt occurred. 64 Environment Canada Station Automatic Snow Pillow Environment Canada Station Precipitation Automatic Snow Pillow Precipitation Automatic Snow Pillow SWE X Apr 01 April Apr 21 May May 11 May 21 May 31 Apr 01 April Apr 21 May 01 May 11 May 21 a) Figure 3.3 Comparison of the cumulative precipitation at the Environment Canada station and the Barkerville ASP with a) the daily precipitation at the Environment Canada station and the Barkerville ASP, and b) the cumulative SWE at the Barkeiville ASP. The ASP and MT daily average air temperature data were similar to each other but vary from the EC temperatures on a number of occasions (Figure 3.4). The MT air temperatures are slightly warmer than at ASP on most days. There are numerous cases of temperature inversions with air temperatures at the higher elevation ASP and MT stations exceeding those measured at the EC station. The temperatures in May at the EC station are consistently colder than at either of the higher elevation stations. The EC station may be subject to cold air pooling or shading due to its narrow valley bottom location. Daily average air temperatures rise well above the freezing point on 12 and 13 April, 2008, slightly on 16 April, 2008, and warm considerably again on 27, 28, and 29 April, 2008, and then remain above freezing from 2 May, 2008, onward. The EC station gets above the freezing point on 27 April, 2008, and remains positive from that point on. These observations provide support for using on-site air temperature to determine precipitation type. 65 May 31 b) u TO m a. E ID TO < TO Q •o^- Apr 01 April Apr 21 May 01 Mt Torn (1487 m) Barkerville ASP (1483 m) Barkerville E.C. (1283 m) l l l May 11 May 21 May 31 Figure 3.4 Comparison of the measured daily average air temperatures at the Mt Tom meteorological station, Barkerville ASP, and Barkerville Environment Canada station. During the study period hourly solar radiation values up to 1129 W m" were measured at the south aspect pyranometer (SPyr) and up to 1118 W m" were measured at the meteorological station pyranometer (MET). The high values occurred during midday and were thus assumed acceptable. The north and south aspect pyranometer data were compared against those measured at the meteorological station. SPyr and MET radiation values were similar and the north aspect pyranometer (NPyr) radiation was less than either SPyr or MET radiation on most days (Figure 3.5a). On overcast days all three measures of incoming shortwave radiation approached the same value. A comparison of the ratio of each radiation measurement showed that SPyr and MET received similar amounts of radiation and NPyr 66 received less radiation than either of the other two stations (Figure 3.5b). The difference in radiation receipt between sites decreased over time as an increasing solar angle reduced the effect of aspect. Anomalous ratios were observed on three dates in April. On 15 April, 2008, NPyr received relatively less radiation than the other sensors and on 18 April, 2008, this pyranometer received relatively more radiation than the other sensors (Figure 3.5b). Both of these days were overcast (based on low overall radiation) and some snow may have accumulated on 15 April, 2008. On 30 April, 2008, MET received relatively more radiation than the other two sensors. Snowfall occurred on this date and snow may have accumulated or ice may have built up on the two aspect pyranometers, which are located in sheltered locations compared to the meteorological station, while blowing clear from MET. SAspect MetStn NAspect MetStn SAspect NAspect j^V» b) n Apr 01 1 April 1 Apr 21 r May 01 May 11 May 21 May 31 Apr 01 April Apr 21 May 01 May 11 May 21 May 31 Figure 3.5 a) daily average incoming shortwave radiation at the south and north aspect and meteorological station pyranometers, and b) the ratio between daily average incoming shortwave radiation at each of the three pyranometers. 67 3.3.2 Observed Snowmelt Snow surveys done between 1 April, 2008, and 5 June, 2008, showed the expected effect of forest cover and aspect on snow accumulation and melt (Figure 2.9b). The north aspect sites had more snow and melted later in the year compared to the south aspect sites. The open sites accumulated more snow and the melt proceeded more rapidly than at the forest sites, as can be seen by the steeper slope in the plot of SWE versus time (Figure 2.9b). A lower accumulation of snow (Winkler et al. 2005) and smaller melt rate in forested compared to open areas is common (Boon 2009). However, the forest sites in this study melted out sooner than their open counterparts. The earlier melt may have resulted from the much lower initial SWE at the forest sites combined with the increase in longwave radiation to the snowpack from the surrounding canopy, trunk, branches, and shrubs (Pomeroy et al. 2006; Pomeroy et al. 2008). During my snow surveys 1 observed large wells around the trunks, especially on the north aspect forest site. Lopez-Moreno and Stahli (2008) observed both an extension and a shortening of the snow season in forested areas in their study, and theorized that the balance between the reduction in shortwave and increase in longwave radiation beneath a forest canopy largely determines the relative timing of melt-out. At my study sites the forest on the north and south aspect sites appeared to differ in stand structure, which could affect melt, but these characteristics were not monitored in this study. Based on data collected by Patrick Teti at his sites (see section 4.2.3 for more information) the south aspect forest had a denser canopy closure (60%) compared to the north aspect forest (46%). The change in SWE was calculated between subsequent sample dates for each snow course site (Figure 3.6). Negative values indicate net loss and positive values indicate net accumulation of SWE during the intervening period. At the north aspect sites accumulation 68 dominated for the first few sample periods. At the south aspect sites the change in SWE was nearly neutral. After 9 May, 2008, all sites experienced melt. The 95% confidence intervals are plotted for the change in SWE for each period (Figure 3.6). There are only three instances where the change in SWE differs significantly between sites. For the first measurement period the south aspect forest site experienced melt whereas the north aspect forest site experienced accumulation. The south aspect open site experienced significantly greater melt in comparison to the south aspect forest site for the period 17 to 23 May, 2008. The north aspect open site experienced significantly greater melt than any of the other sites for the period 23 to 30 May, 2008. However, this is partly due to the fact that there was very little snow left to melt for the south aspect sites by this sample period. Total period and daily average melt are shown in Figure 3.7. These data were not used in calibrating or testing melt models, but are useful as they show a standardized comparison of melt among different sample periods. This is particularly important given the different timing of melt for each site and the different lengths of each sample period. Melt started earlier and ended earlier at the south aspect sites, but the maximum daily melt was similar among aspect sites with the same cover (Figure 3.7b). 69 o o J E E, Ilr tu -t: CO c .1.1 J i en w i_ a> o o A" ZP 3 i A so o o A SF • NO • NF 100 i 120 110 130 i ~~i i r 140 150 160 Julian Date Figure 3 6 Change in SWE between each snow course sample date with 95% confidence intervals Positive values indicate accumulation and negative values indicate ablation of snow S = south aspect, N = north aspect, 0 = open, F = forest o ID CM • O O A 0 - A E E o m - £ • D m 6 o o - as A n A a m o A A & I May 11 A A SO * SF • NO a NF o I May 01 g A I I May 21 May 31 m A SO A SF • NO a NF A A g) I May 01 I May 11 I I May 21 May 31 Figure 3 7 a) total peuod and b) aveiagc daily melt at each snow couise site foi the cahbiation penod 27 Apt ll 2008 - 5 June, 2008 70 b) The two TidBit sites had similar initial SWE and time-to-melt values (Figure 3.8). Snow disappeared on the same date for each site pair, except for the east aspect open site. One location became snow free on 25 May, 2008, and the other location became snow free on 29 May, 2008 (Figure 3.8). o in o 0) o CO o o CM 1 May 17 r May 21 Figure 3.8 Total melt calculated for each TidBit measurement point. E = east aspect, W = west aspect, F forest, O = open. Average daily melt rates for the TidBit sites (calculated from total melt and length of melt period) were greatest for the east open site (ca. 15 mm d" ) and lowest for the west forest site (ca. 6 mm d"1). The average daily melt at the east forest site was slightly lower than that at the west open site (ca. 11 mm d"' and 14.5 mm d"1 respectively). Although the melt finished earlier on the south compared to the north aspect sites (Figure 2.9b) and on the 71 west compared to the east aspect sites (Figure 3.8), the average melt rates were greater for the north and east aspect sites compared to the other sites (Figure 3.9a). The north aspect sites started melting later in the season yet a greater amount of snow was melted in approximately the same amount of time as for the south aspect sites. A south aspect open A south aspect forest a north aspect open a north aspect forest £ £ ° east aspect open e east aspect forest O west aspect open • west aspect forest 200 n~ ~r 300 400 200 500 Total Melt (mm) 300 400 500 Total Melt (mm) Figure 3.9 a) average daily melt rate for each Mt Tom study site as a function of total SWE melted at each site, b) number of days to melt total SWE at each Mt Tom aspect site as a function of total SWE melted at each site. 3.3.3 Model Calibration 3.3.3.1 Basic temperature-index Model The basic temperature-index model (bTl) had similar melt factors among all the calibration runs (Table 3.2). The melt factor was lower when two forest sites were used in calibration and higher when two open sites were used in calibration. 72 Table 3.2 Optimized melt factor (MF) using the basic temperature-index model for each leave-one-out run (coeff. set 1:4) and for all four snow course sites together (coeff. set 5). MF coeff. set withheld 1 SO 2.33 35.61 2 SF 2.64 38.83 3 NO 2.30 34.80 4 NF 2.78 35.13 5 none 2.52 37.26 (mm °C 1 d"1) RMSE (mm) 3.3.3.2 HBV-EC Model The calibrated parameter values for the HBV-EC model (Table 3.3) show how the parameters depend on which sites are used for calibration. A non-zero value for the melt factor at winter solstice (Cmm) is only calibrated when the south aspect open (SO) site is withheld. In this situation the increase in the melt factor (DC) is much lower than for the other four coefficient sets. When the SO snow course data are included, Cmm does not improve model performance. However, the sum of Cmm and DC, which is the melt factor at summer solstice, is consistent among all coefficient sets. The effect of aspect and forest cover is reflected in the slope-aspect (AM) and the forest melt (MRF) reduction factors. When there is only one open site data set in the calibration data (coefficient sets 1 and 3) the effect of slope is reduced and the AM is lower. Interestingly the MRF appears to be responding more to the effect of aspect than to the presence of forest cover. The MRF value is the same for coefficient set 1 and 2, despite the first set having two forested and one open snow course in the calibration data and the second set having one forested and two open snow courses in the calibration data. Coefficient sets 3 and 4 show a similar pattern, but with a slightly lower MRF value. The MRF is supposed to represent the ratio between the MF in the forest to the 73 MF in the open. An MRF of 1.0 indicates that melt rates in the forest are equal to those in the open. When the model was calibrated with two north and one south aspect sites there was a slightly higher forest melt rate compared to the open melt rate than when the model was calibrated with one north and two south aspect sites. Table 3.3 Parameter values calibrated using the HBV-EC model for each leave-one-out run (coeff. set 1:4) and for all four snow course sites together (coeff. set 5). Minimum melt factor (C,„,„); increase in melt factor (DC); slope-aspect reduction factor (AM); forest melt reduction factor (MRF) DC coeff. set withheld 1 _1 (mm°C- d ) (mm°C- 1 d" 1 ) AM MRF RMSE (mm) 1 SO 1.54 1.58 0.17 0.65 26.19 2 SF 0 3.72 0.69 0.65 25.83 3 NO 0 3.98 0.23 0.55 15.94 4 NF 0 3.72 0.69 0.55 25.54 5 none 0 3.68 0.58 0.61 24.26 3.3.3.3 Hock Model The Hock model melt factors for forest sites were fairly consistent among loo runs and type of radiation input (Table 3.4). When potential surface radiation was used, all but one M/.-0 value calibrated to zero. The Rf-F was zero for all cases. This indicates shortwave radiation was not important in describing either the spatial or temporal variation in melt at forested sites. The RF-0 was greater than zero for all cases except one (to be discussed later), which indicates that shortwave radiation had at least some effect on the spatial pattern of melt at open sites in this study. The value of RF-0 was dependent on both the withheld site and the type of radiation input used in the model run. 74 Table 3.4 Parameter values calibrated using the Hock model for each leave-one-out run (set 1:4) and all sites pooled (set 5). Sub-tables are for each type of radiation input data used. Mf = melt factor, Rf = radiation factor, -O = open site, -F = forest site. Measured RF-0 RF-F (mm°C" 1 d" 1 ) (mm (W m"2)"1 °C"1 d"1) (mm (W rrf2)"1 c C -i J-1J RMSE (mm) MrF MF-0 (mm °C_1 d_1) set withheld 1 SO 1.86 0.80 0.0105 0 25.65 2 SF 1.83 1.89 0.0057 0 28.22 3 NO 1.86 3.47 0 0 18.39 4 NF 1.91 1.89 0.0056 0 28.32 5 none 1.86 1.89 0.0056 0 26.20 Modelled Mp-O set withheld 1 RF-0 RP-F (mm°C" 1 d~ 1 ) (mm°C- 1 d- 1 ) (mm (W m"2)"1 °C"1 d_1) (mm (W m"2)"1 °C' 1 d"1) RMSE (mm) SO 1.86 2.28 0.0017 0 26.32 2 SF 1.83 0 0.0088 0 27.95 3 NO 1.86 0 0.0099 0 17.26 4 NF 2.08 0 0.0088 0 28.27 5 none 1.86 0 0.0088 0 25.98 Rf-0 was zero for coefficient set 3 when the measured radiation input values were used (Table 3.4). Thus, when considering the south aspect open site the simple temperature index worked as well as or better than a model including radiation. For all other cases Rp-0 was greater than zero indicating that radiation was needed to explain the spatial variation in melt for cases where both open site aspects were used for calibration (sets 2,4, and 5) and was needed to explain the temporal variation in melt when only one site was used in calibration (sets 1 and 3). 75 The open and forest site melt factors were similar when measured radiation was used. When potential surface radiation was used, the Mf-0 equalled zero for all cases except when the south open site was withheld (Table 3.4). In these situations, shortwave radiation explained all the variation in melt and accounted for melt variation better than air temperature for these conditions. When the south aspect open site data were withheld (set 1), both air temperature and shortwave radiation were necessary to explain the observed melt at the north aspect open site. The melt at the south aspect open site follows a similar pattern to the daily average potential incoming shortwave radiation: a relatively continuous increase over time (Figure 3.10a, open triangles). Since it is potential incoming radiation, the only variation in the daily average radiation used as a model input is a steady increase as the season progresses. This good fit between the modelled potential radiation and melt at the south aspect site may explain the calibration result, but the relationship may be coincidental. Measured daily radiation (Figure 3.10b) and daily air temperature (Figure 3.10c) follow a more variable increase over time, as weather conditions change, which does not fit as cleanly with the observed melt. Flowever, air temperature does appear to show a stronger increasing temporal trend compared to measured shortwave radiation. Melt rates at the north aspect snow course sites tend to fluctuate, whereas at the south aspect sites melt rates increase with time. This may help explain the low or absent radiation factors calibrated for the models using the modelled radiation values (Figure 3.10a, Table 3.4); modelled potential radiation did not help explain variation in melt among sites. In these situations the basic temperature index approach was able to predict melt as well or better than a model including potential incoming shortwave radiation. 76 May 01 May 11 May 21 May 31 May 01 g) May 01 May 11 May 21 May 11 May 21 May 31 May 31 Figure 3.10 Average daily observed melt for each snow course site. Melt was averaged across the intervening measurement period, plotted with: a) daily average modelled potential direct shortwave radiation at the surface for the same period, b) measured shortwave incoming radiation at the meteorological station, and c) average daily air tcmperatuie at the meteorological station. 3.3.3.4 Pellicciotti Model The Pellicciotti model (Pellicciotti et al. 2005) includes a measured or modelled albedo term. Albedo was not measured on site until the last week of melt and attempts to model the albedo using the methods in Pellicciotti et al. (2005) could not be confirmed (section 3.2.3 Albedo Model). However, the calibrated RF values using these albedo values were zero for almost every case (and zero for all coefficient set 5 values). Therefore, neither the albedo value nor the type of radiation used mattered for calculating melt. The Pellicciotti model 77 essentially became a basic temperature-index model for this study The averages of the open and forest melt factors for the Pelhcciotti model were similar to the melt factors calibrated for the basic temperature-index model (Table 3 2) except when a radiation factor was calibrated for the Pelhcciotti model (1 e for set 1, Table 3 5) Table 3 5 Parameter values calibrated using the Pelhcciotti model for each leave-one-out mn (set 1 4) and all sites pooled (set 5) Sub-tables are for each type of radiation input data used MF = melt factor, RF = radiation factor, -O = open site, -F = forest site Measured RF-0 MF-0 set withheld 1 RF-F RMSE (mm) (mm°C1d1) (mm°C1d1) (mm (W m ) d1) (mm (W m 2 ) 1 d1) SO 1 86 0 90 0 1533 0 25 25 2 SF 1 83 3 03 0 0 28 70 3 NO 1 86 3 47 0 0 18 39 4 NF 1 91 3 03 0 0 28 79 5 none 1 86 3 03 0 0 26 61 2 1 Modelled set withheld 1 RF-F RrO MF-F 2 1 RMSE (mm) ( m m ° C 1 d1) (mm°C1d1) (mm (W m ) d1) (mm(Wm2)1 d1) SO 1 86 2 68 0 0073 0 26 30 2 SF 1 83 3 03 0 0 28 70 3 NO 1 86 3 47 0 0 18 39 4 NF 1 91 3 03 0 0 28 79 5 none 1 86 3 03 0 0 26 61 The cahbiated Pelhcciotti model coefficients were identical between the two radiation data sets used, except when a radiation factoi was cahbiated (set 1, Table 3 5) The forest radiation factor (Rp-F) was zero for every case and the open ladiation factor (RF-O) was 78 greater than zero only when calibration was done with just the north aspect open site (set 1, south aspect open site withheld). When measured shortwave radiation was used the set 1 open site melt factor was smaller and the forest radiation factor was larger than when modelled shortwave radiation was used. Shortwave radiation had only a small effect on melt for the Pellicciotti model. It is expected that radiation would more likely be required to explain melt when two open sites are used in calibration (i.e. set 2, 4, and 5). However, The RF-0 was significant only when the north aspect open site was used in calibration (set 1) indicating that shortwave radiation was necessary to explain the temporal variation in melt only for the north aspect. When only the south aspect or when both open sites were used in calibration, shortwave radiation no longer improved model performance. Further details on the Hock and Pellicciotti model parameters can be found in Appendix 5. 3.3.4 Model Performance None of the temperature index models compared in this study account for the heat content of the snowpack. Melt simulations started before the date of peak SWE overestimated melt for most sites (Figure 3.11 and Figure 3.12). All models were therefore run from the start of the period of continuous melt and comparison of model performance was based on these results. 79 •St o a- '—V> ° ° V y ^ . . . l _ I J ••• I'A O ° ^^\N \ ° Observed Basic temperature-index HBV-EC Hock measured K Hock modeled K - - - Pellicciotti \\° \t "i Apr 01 Apr 21 i Apr 01 r Apr 21 May 11 i i May 11 a) \ v»°V \ h \ 1 8 r May 31 Apr 01 r May 31 Apr 01 Apr 21 "i i Apr 21 r i May 11 i ] i May 11 r May 31 r May 31 Figure 3.11 Observed and modelled SWE over lime for each of the Mt Tom snow course sites starting 1 April, 2008. a) south aspect open; b) south aspect forest; c) north aspect open; d) north aspect forest. 80 m OOO33COO0OCO0OCQ;) CO tfCocP o co _ m CO 're Q CD CM - CD CN — m cu 3 CD in o _ o o Observed Basic temperature-index HBV-EC Hock measured K Hock modeled K Pellicciotti —i i i Apr01 April Apr 21 1 1 1 r May 01 May 11 May 21 May 31 Figure 3.12 Observed and modelled daily SWE over time for the Barkerville automatic snow pillow starting 1 April, 2008. 3.3.4.1 Mt Tom north and south aspect snow course sites When the leave-one-out (loo) calibration method was used, the Hock model run with modelled potential surface shortwave radiation had the best performance based on the NSE, RMSE and G values (Table 3.6). Using all the sites, the HBV-EC model had the best performance based on these same values (Table 3.6). Compared to the other three models the basic TI model had higher RMSE and smaller NSE values. The measures of bias (MBE and RBE) indicate that all the models moderately overpredict the melt for the snow course sites. Using all the sites in model calibration produced a better fit with lower bias compared to the leave-one-out approach. The relative bias error (RBE) indicates the percentage by which predicted values of the mean differ from observed. For example, the Pellicciotti model run with measured radiation had an RBE of 14.74% 81 using loo calibration coefficients (Table 3.6). Thus, this model will over-predict melt by ca. 15% of the mean observed melt (average melt was 104.87 mm in the case of the snow course sites, therefore an average over prediction of 0.1474 x 104.87 = 15.46 mm for each measurement period). According to McCuen et al. (2006) an RBE greater than 5% indicates noteworthy bias in the model. Using this criterion, all the models except for the basic temperature-index would be considered biased using the loo calibration coefficients. When the all calibration coefficients were used the RBE did not exceed 2.25% for any model. 82 Table 3 6 Calculated measures of en or for each model for both leave-one-out and all-sites-combincd calibration methods Subscript "surfC" indicates the model lun with potential surface shoitwave radiation Nonsubscripted model luns used measured shoitwave radiation Root mean squaie enor (RMSE), mean bias eiror (MBE), relative bias eiror (RBE), Nash-Sutchffe efficiency criterion (NSE), and goodness-of-fit indices (G) The basic temperatuie-index model (BTI) was used as the benchmark senes N = 14 for each loo and all methods Leave-one-out calibration method MODEL RMSE (mm) MBE (mm) RBE (%) NSE G basicTI 44 30 1 10 1 05 0 37 0 HBV-EC 31 95 5 96 5 69 0 67 0 48 Hock 32 31 8 74 8 34 0 66 0 47 HocksurK 31 85 7 67 7 32 0 68 0 48 Pellicciotti 38 34 1546 14 74 0 53 0 25 Pellicciotti surK 32 39 7 97 7 60 0 66 0 46 All-sites-combined calibration method MODEL RMSE (mm) MBE (mm) RBE (%) NSE G basicTI 37 26 2 25 2 15 0 56 0 HBV-EC 24 26 2 35 2 24 0 81 0 58 Hock 26 20 1 90 1 82 0 78 0 50 HocksurK 25 98 1 74 1 66 0 78 0 51 Pellicciotti* 26 61 2 04 1 95 0 77 0 49 PelllCCIOttlsur|<* 26 61 2 04 1 95 0 77 0 49 * The calculated measures of enor aic the same between the Pellicciotti model lun with measured shoitwave ladiation and that lun with modelled potential surface shoitwave ladiation because with all sites used in calibration shoitwave ladiation was not used in delcimining melt, thus, the type of ladiation used is inconsequential to model perfoimance 83 Even though the Pellicciotti model using all coefficients is essentially a basic temperature-index model (no radiation factor), the error values for the Pellicciotti model using all calibration coefficients are more similar to the Hock model and are better than the error measures for the basic temperature-index model (Table 3.6). The only difference between the Pellicciotti and basic temperature-index models is that the effect of forest cover is considered by calibrating separate melt factors for open and forest sites for the Pellicciotti model. This appears to improve model performance without requiring additional meteorological data. Based on the goodness-of-fit indices, all models performed better than the basic TI for both calibration methods (Table 3.6). Overall, the HBV-EC model gave the greatest improvement over the basic temperature-index model (G = 0.58). In addition to the improved measures of error for the all-sites-combined calibration method, a graphical comparison of predicted and observed melt for each model shows that the all-sites-combined calibration method resulted in better prediction of melt than the leaveone-out method (Figure 3.13). When the leave-one-out method was used, melt was generally best predicted for the forest sites. The two or three smallest melt amounts for the south aspect open site are consistently predicted accurately by all models using both calibration methods, but the one or two highest melt amounts for that site were usually under predicted (Figure 3.13). Melt at the north aspect open site was the least accurately predicted (hence the lowest RJV1SE for coefficient calibration occurred for all four models whenever the north aspect open site was withheld). The smaller melt amounts were better predicted than the larger amounts. All the models show less variation between predicted and observed melt in comparison to the basic temperature-index model (Figure 26). The basic temperature-index 84 model poorly predicts melt when greater than ca. 100 mm occurred; melt amounts were mostly over-predicted for the forest sites and under-predicted for the open sites (Figure 3.13a). The high bias observed for the Pellicciotti model run with loo coefficients and measured shortwave radiation resulted primarily from the over-prediction of melt at almost all the open sites (Figure 3.13f). 85 alt loo all loo South Open • South Forest North Open North Forest r •a Z< a. i d) O E i S E e) 50 100 200 Observed melt tmmi 300 f) 0 50 100 200 300 Observed melt (rrifn) 0 50 100 200 300 Observed melt (mm) 50 100 Observed melt (mm) Figure 3.13 Predicted vs. observed melt calculated for each snow course site using leave-one-out (loo) and allsites-combincd (all) calibiation methods. Period of comparison was 1 April, 2008 to 5 June, 2008. The straight line indicates a 1:1 relationship, a) basic temperature index, b) HBV-EC; c) Hock (measured radiation); d) Hock (modelled surface radiation); e) Pellicciotti (measuied radiation); f) Pellicciotti (modelled surface ladintion) 3.3.4.2 Barkerville automatic snow pillow Each model was tested at the Barkerville automatic snow pillow (ASP) site, using the ASP station record of air temperature and SWE. Daily air temperature data were highly correlated between the Barkerville ASP and Mt Tom meteorological station (r = 0.99) (Figure 3.4). Model performance was similar among the modified temperature-index models at the Barkerville ASP, and all models performed better than the basic TI as indicated by positive G indices (Table 3.7). The Hock model run with modelled potential and measured 86 200 300 shortwave radiation and the Pelhcciotti model had identical measures of error (Table 3 7) This occurred because all of these models calibrated a zero radiation factor for forested conditions Shortwave radiation did not improve model performance over using air temperature by itself The HBV-EC model had marginally smaller RMSE, MBE and RBE and larger NSE values compared to the other temperature-index models As seen in the plots of observed against predicted daily melt (Figure 3 14) the basic temperature-index model overestimated snowmelt whereas the modified temperatuie-index models underestimated snowmelt for the BarkerviUe ASP Table 3 7 Calculated measures of cnor using coefficient set 5 to piedict melt at the BarkerviUe ASP site Subsciipt "surK." indicates the model run with potential suiface shoitwavc ladiation Non-subscripted model mns used measured shortwave radiation Root mean square error (RMSE), mean bias error (MBE), relative bias enor (RBE), Nash-Sutchffe efficiency cntenon (NSE), and goodness of fit (G) N = 30 MODEL RMSE (mm) MBE (mm) RBE (%) NSE G basicTI 6 84 2 57 22 11 0 29 0 HBV-EC 5 49 -0 05 -0 44 0 54 0 35 Hock 5 58 -1 15 -9 90 0 53 0 33 HOCKsur|< 5 58 -1 15 -9 90 0 53 0 33 Pelhcciotti* 5 58 -1 15 -9 90 0 53 0 33 * The Pelhcciotti calculated measmes of enoi using modelled potential surface shoitwavc radiation weic identical and thercfoic aie not shown heie 87 o J~ a) 10 20 b) 30 10 Observed melt (mm) E 20 30 Observed melt (mm) •=> -, c) 10 20 SQ Observed melt (mm} d) 10 i i 20 30 Observed melt (mm) 10 20 30 Observed melt (mm) Figure 3.14 Predicted vs. observed melt calculated for the Barkerville automatic snow pillow using coefficient set 5. Period of comparison was 1 to 30 May, 2008. The straight line indicates a 1:1 relationship, a) basic temperature index; b) HBV-EC; c) Hock (measured radiation); d) Hock (modelled surface radiation); e) Pellicciotti (measured radiation). Pellicciotti model results using potential surface shortwave radiation were identical and therefore aie not shown heie. 3.3.4.3 TidBit east and west aspect sites All models performed better for the TidBit sites than for the Barkerville ASP data, with NSE values mostly in the 0.90 range (Table 3.8). However, only four data points were used for model testing at the TidBits sites. Each data point represented the total melt at each TidBit site with no reflection of daily or weekly melt variation. All models exhibited low bias, with no RBE values exceeding 3% (Table 3.8). The Hock model run with measured north and south aspect radiation (Hock-NS) or with modelled radiation (HockSLU|<) under- predicted melt whereas the Hock model run with radiation measured at the meteorological station (Hock-met) over-predicted melt All models provide improved melt estimation in comparison to the basic temperatuie-index model, with G ranging from 0 74 to 0 84 HockNS had the smallest RMSE and the largest NSE and G (Table 3 8) Since season total melt rates were used at the TidBit sites the absolute errors in prediction were large compared to RMSE values at the Barkerville ASP (Table 3 7 and Table 3 8) Table 3 8 Calculated measures of eiror using coefficient set 5 to predict melt at the TidBit east and west aspect sites Subscript "surK" indicates the model mn with potential suiface shortwave radiation Non-subscripted model runs used measuied shoitwave ladiation Root mean square error (RMSE), mean bias en or (MBE), lelative bias error (RBE), Nash-Sutchffe efficiency cnterion (NSE), and goodness of-fit (G) The -NS and met affixes indicate which measuied shortwave ladiation was used in modelling as ladiation was not measured for the east and west aspects noith and south aspect (NS) oi that measuied at the meteorological station (met) N-4 MODEL RMSE (mm) MBE (mm) RBE (%) NSE G basicTI 65 61 6 49 2 11 0 72 0 00 HBV-EC 30 71 8 71 2 83 0 94 0 78 Hock-NS 25 99 -3 93 -1 28 0 96 0 84 Hock-met 33 68 9 01 2 93 0 93 0 74 HocksurK 29 29 -4 28 -1 39 0 94 0 80 Pellicciotti* 28 72 1 44 0 47 0 95 0 81 * The calculated measuies of ciroi aie the same among all the Pellicciotti model luns because all the ladiation factois in coefficient set 5 are zeio, thus, the type of ladiation used is inconsequential to model peifoimance 89 v East Open ^ East Forest - West Open / West Forest "> East Open '- East Forest - West Open • West Forest irt ~* '> East Open — East Forest _ ~ West Open West Forest a) 100 200 300 400 500 m "" - East Open • *"* East Forest - West Open S -, - West forest £ o c) b) 0 100 200 300 400 500 ' East Open • East Forest West Open West Forest tQO 200 300 400 500 '" East Open <- East Forest - West Open " West Forest , e) d) 100 20O 300 Observed melt [mmj 400 600 0 100 200 300 Observed ineft (mm) 400 5D0 f) 0 1O0 200 300 400 Observed mett (mm) Figure 3.15 Predicted vs. observed melt calculated for the Mt Tom east and west aspect TidBit sites using coefficient set 5. Period of comparison was 1 April, 2008, to 5 June, 2008. The straight line indicates a 1:1 relationship, a) basic temperature index; b) HBV-EC; c) Hock (measured meteorological station radiation); d) Hock (measured north and south aspect radiation); e) Hock (modelled surface radiation); f) Pellicciotti (measured radiation). Pellicciotti model results using all other radiation data were identical and therefore are not shown here. The basic temperature-index model overestimated melt for both the east and west aspect forest sites and underestimated melt for both the open sites because the basic temperature-index model does not account for differences in forest cover (Figure 3.15a). Since melt rates are typically lower in the forest when compared to open sites the averaged melt factor will simulate too much melt for the forest sites and too little melt for the open sites. Thus these results are as expected from this model. Including the effect of radiation (through slope, aspect, and forest cover) resulted in better prediction of melt at all TidBit sites (Figure 3.15b-f). The modified temperature-index models had results similar to each other. Melt at the east aspect open and west aspect forest 90 500 sites was overestimated and melt at the east aspect forest and west aspect open sites was underestimated (Figure 3.15b-f). The east open and west forest sites experienced the greatest (475 mm) and least (125 mm) amount of melt respectively, whereas the east forest and west open sites experienced similar, mid-range amount of melt (325 mm and 305 mm respectively). The melt at the Mt Tom east and west aspect sites is a direct reflection of the amount of snow at the sites as of 27 April, 2008 - more snow equates to more melt volume. This does not necessarily mean that daily melt rates followed the same pattern, as the length of the melt period varied for each site. In this case melt rates do appear to follow this same pattern (Figure 3.9). Average daily melt rates for the TidBit sites (calculated from total melt and length of melt period) were greatest for the east open site (ca. 15 mm d"1) and lowest for the west forest site (ca. 6 mm d"1). The average daily melt at the east forest site was slightly lower than that at the west open site (ca. 11 mm d" and 14.5 mm d"1 respectively). 3.3.5 Model Comparison To calculate SWE over time the predicted daily melt for each site and each model was subtracted from the initial measured SWE and daily snow fall was added to the daily SWE. The basic temperature-index model predicts SWE reasonably well during the melt period for the open sites (Figure 3.16a,c), but under-predicts SWE for the forest sites (Figure 3.16b,d). The same under-prediction of SWE is seen for the Barkerville ASP, which is a small gap in the forest and was modelled as a forested site (Figure 3.17). The melt factor appears to be more appropriate for the open snow course sites. At the TidBit sites, however, the predicted SWE is too high for the open sites and matches the observed forest site SWE well (Figure 3.18). The other three models (HBV-EC, Hock, and Pellicciotti) performed similarly for all the sites being compared. Only the Hock model used shortwave radiation to help predict melt, but the HBV-EC accounted for the effects of slope, aspect, and forest cover with melt reduction parameters and the Pellicciotti model calibrated separate melt factors for open and forest sites. Based on the similarity in SWE predictions among these three models despite their differences in structure, it would appear that the effect of cover (being either open or forested) had a strong influence on the melt factor. a) south open E o Observed Basic temperature-indexY HBV-EC \ \ Hock measured K Hock modeled K Pellicciotti T I May 01 May 11 May 21 May 01 May 31 north open May 11 May 21 May 31 c) \ E o r^o\o \\ May 11 I I I i I May 16 May 21 May 26 May 31 Jun 05 May 11 ~\ i 1 I r May 16 May 21 May 26 May 31 Jun 05 Figuie 3 16 Obscived and modelled SWE ovei time for the period of continuous melt foi each of the Mt Tom snow couise sites a) south aspect open, b) south aspect foiesl, c) north aspect open, d) noith aspect foicst 92 o — — --•• — — 1 Observed Basic temperature-index HBV-EC Hock measured K Hock modeled K Pellicciotti I I 1 1 1 May 01 May 06 May 11 May 16 May 21 May 26 May 31 Figure 3.17 Observed and modelled SWE over time for the period of continuous melt for the Barkerville ASP (1-30 May, 2008). 93 <*. east o p e n o Observed V \^ Basic temperature-index N ^ \ HBV-EC N^ \ HockNS \ \ Hock met N» \ --- Hock modeled K V. — Pellicciotti \ ,\ May 01 May 11 May 21 3 ) May 31 May 01 May 11 May 21 May 31 May 01 May 11 May 21 May 31 west open Hfc^>^ ^ \ X. o^A I I I May 01 May 11 May 21 e« c) i May 31 Figure 3 18 Observed and modelled SWE over time for each of the TidBit sites, a) east aspect open; b) cast aspect fotcsl; c) west aspect open; d) west aspect forest. The ability of the models to correctly predict the snow free date can be examined for the Barkerville ASP and TidBit data sets. Unfortunately, the snow-free date for the Mt Tom snow course sites was not known precisely. The site was snow free on the last sample date, but it may have gone snow-free on any day between the previous measurement and the last. Thus the accuracy of the predicted snow-free date cannot be assessed for the Mt Tom snow course sites. However, simulation was not good for the north aspect forest site since all models, except the basic temperature-index, predicted the snow free date when snow was still observed on site (Figure 3.16d). 94 At the Barkerville ASP the basic TI model predicted the snow free date approximately three days before observed, the HBV-EC predicted the snow free date accurately, and the Hock and Pellicciotti models predicted approximately 4 cm of SWE on the observed snow free date (Figure 3.17). Prediction of snow free date was most accurate at the TidBit sites. For the east aspect open and the west aspect forest sites all models, except the basic temperature-index, correctly predicted the snow free date (Figure 3.18a,d). For the east aspect forest and west aspect open sites all models, except for the basic temperature-index, predicted the snow free date one to two days later than observed (Figure 3.18b,c). 3.4 Discussion 3.4.1 Model Parameter Values The melt factor values in this study are smaller than calibrated for the basic temperature-index models in Hock (1999), 6.3 mm ° C ' d"1; and Pellicciotti et al. (2005), 7.68 mm °C~' d"1, but within the range calculated for a forested mountain basin in southwestern Idaho, USA of 0.2 to 3.4 mm "C"1 d"1 (Franz et al. 2008). Hock (2003, Table 1) presented a number of degree-day factors calculated for snow both on and off glaciers. On glacier factors ranged from 2.7 to 11.6 mm °C~' d"1, while off-glacier factors ranged from 2.5 to 5.5 mm °C~' d"1. Hamlin et al. (1998) found melt factors ranging from 1.8 to 3.36 mm °C~' d"1 for various classes of mixed, transitional, and coniferous forests. Kuuisisto (1980) found melt factors ranging from 1.75 - 3.36 mm °C~1 d"1 for forested sites, and ranging from 2.82 4.94 mm °C~ d" for open sites located across Finland. Kane et al. (1997) found melt factors for a small Arctic watershed ranged between 1.2-3.9 mm °CA d"1 for different years. The 95 melt factor calibrated for the combined open and forest sites at Mt Tom fits within the range of melt factors found in other off-glacier studies incorporating forest cover. In the HBV-EC model (CHC 2007), the default values are: Cmin = 2.0 mm °Cl d"1, DC = 2 mm °C"1 d"1, and MRF = 0.7 (typical range 0.6 - 1.0). Hamilton et al. (2000) calibrated the following parameter values for a forested drainage basin in the Yukon: Cmin = 1.65 mm°C"' d"\ DC= 2.55 m m ' C ' 1 d"1, and MRF = 0.706. A Cmi„ value of zero was calibrated in this study, and the DC value was larger than in other studies. The difference in values might be due to calibrating the model only to the period of continuous melt when the Cmm is intended to represent the melt at winter solstice and DC reflects how the melt factor changes throughout the entire snow season. In Hock's study (Hock 1999), the calibrated parameters were: MF = 6.3 mm °C_1 d"1, and RF = 0.0144 mm (W m"2)"1 °C"' d"1 for the model using potential direct clear sky shortwave radiation, and MF = 2.1 mm °C"' d"1 and RF = 0.0168 mm (W m"2)"1 °C"1 d"1 for the model using measure shortwave global radiation. Pellicciotti et al. (2005) found the parameters for the Hock model using potential direct clear sky shortwave radiation to be MF = 1.97 mm °C~' d"1 andRF= 0.0125 mm (W m"2)"1 °C"' d"1. Schuler et al. (2002) found the parameters for the Hock model on a valley glacier in the Swiss Alps to be MF = 0.45 mm °C"' d"1 and RF =0.012 mm (W m"2)'1 °C"' d"1. The radiation factors calibrated in this study are an order of magnitude lower than observed in other studies indicating a lower importance of shortwave radiation to melt variability. The comparison studies found were for glacier sites. The melt factors calibrated for glacier sites were larger than off-glacier sites in basic temperature-index studies (Kuuisisto 1980, Kane et al. 1997, Hamlin et al. 1998, Hock 1999, Hock 2003, Pellicciotti et al 2005, Franz et al. 2008). Radiation and 96 temperature melt factors might also be higher on glacier sites. In this study the melt factors (radiation factors) were smaller (larger) when modelled potential shortwave radiation was used compared to measured shortwave radiation which is consistent with Hock (1999). The RF was zero for forest sites. This indicates shortwave radiation was not important in describing either the spatial or temporal variation in melt at forested sites. This is expected as the canopy excludes a portion of shortwave radiation, dependent on the canopy structure and density (Male and Granger 1981; Storck et al. 1999; Koivusalo and Kokkonen 2002; Pomeroy et al. 2008). The canopy cover fraction at my forest sites ranged from 0.46 to 0.60. In Pellicciotti et al. (2005) the calibrated parameters were: Mp = 1.2 mm °C"1 d"1, and Rp = 0.226 mm (W m 2 ) 1 d"1 for the model using both modelled and measured shortwave radiation and albedo. For a small Arctic watershed Kane et al. (1997) found the Mp ranged between 1.68 - 3.12 mm "C"1 d"1 and the Rp ranged between 0 - 0.072 mm (W m"2)"1 d~' for different years. A radiation factor of zero was calibrated for one year of Kane et al.'s (1997) study. The zero radiation factor calibrated in this study may have resulted from specific conditions in 2008. It would be interesting to calibrate the Pellicciotti model to other years of data to test if a positive radiation value will be calibrated for this site. In this study the Mp values were larger than Pellicciotti et al. (2005), but fit within the range observed by Kane et al. (1997) over a four year period. 3.4.2 Optimal Melt Model Complexity No melt model was conclusively the "best". Performance depended on the measure of error used, the test data, and the model variation. For each test data set there was a different model that outperformed the others based on the measures of error. The Flock model run with modelled radiation had the lowest error values for the snow course loo 97 calibration data. The HBV-EC had the lowest error values for the snow course all calibration data. The HBV-EC model had the lowest error values for the Barkerville ASP data set. The Hock model run with measured radiation had the lowest error values for the TidBit data set. In all cases the modified temperature-index models performed better than the basic temperature-index model. In a recent intercomparison study of various snow models, a "best" model could not be identified either due to inconsistencies in model performance between different sites or years of simulation (Essery et al. 2009). The Hock and HBV-EC models received the most top scores in this study. The HBVEC had eight and the two Hock model variations combined had seven top scores out of 16 (Table 3.9). The Hock model run with modelled potential shortwave surface radiation had the smallest RMSE and the largest NSE for the loo calibration data and also the smallest MBE/RBE for the all calibration data. Having a low bias error score, however, does not indicate that the model accurately predicted melt only that it was not biased. It is important to note that the "best" score is determined relative to the performance of the other models at a particular site and therefore does not necessarily indicate good performance for a given site. The basic temperature-index model had the "worst" scores for every site and the basic temperature-index model was outperformed overall by all the modified models based on the G index. In addition, the differences in error scores among the models were often small. For example, with the exception of the basic temperature-index model, all models using the TidBit data had a NSE score within the range 0.93 to 0.96 and had RMSE values within 5 mm of each other over the entire season (Table 3.8). 98 The Barkerville ASP is located in a small opening in a mature forest and the models were run using the same coefficients as for the forested snow course and TidBit sites. Model performance was reasonable for the modified temperature-index models. The basic temperature-index model, which did not account for forest cover, over-predicted melt for the site. The Pellicciotti model used in this study was essentially a basic temperature-index model as a zero radiation factor was calibrated. However, the performance of the Pellicciotti model at all study sites was more similar to the Hock model than the basic temperature-index model. The Pellicciotti model was at times better (based on RMSE and NSE values) than the Hock model. The only difference between the structure of the basic temperature-index and the Pellicciotti-based model in this study is that the effect of forest cover is being taken into account by calibrating separate melt factors for open and forest sites for the Pellicciotti model. Model performance was improved without the need for additional meteorological data however, more calibration data were required. Table 3.9 The "best" model for each set of test data used based on computed measures of error. Model [resulting test value]. Root mean square error (RMSE), mean bias error (MBE), relative bias error (RBE), Nash-Sutcliffe efficiency criterion (NSE), and goodness-of-fit (G). The - N S (north and south aspect) and-met (meteorological station) affixes indicate which measured shortwave radiation was used in modelling as radiation was not measured for the east and west aspects. test data RMSE (mm) MBE (mm)/RBE (%) NSE G Calibration - loo HocksurK [31.85] basicTI [1.10/1.05] HocksurK [0.68] HBV-EC/Hock surK [0.48] Calibration - all HBV-EC [24.26] HocksurK [1.74/1.66] HBV-EC [0.81] HBV-EC [0.58] Barkerville ASP HBV-EC [5.49] HBV-EC [-0.05/-0.44] HBV-EC [0.54] HBV-EC [0.35] TidBit E/W aspect Hock-NS [25.99] Pellicciotti [1.44/0.47] Hock-NS [0.96] Hock-NS [0.84] 99 Table 3.10 The "worst" model for each set of lest data used based on computed measures of error. Model [resulting test value]. Root mean square error (RJVISE), mean bias error (MBE), relative bias error (RBE), Nash-Sutcliffe efficiency criterion (NSE), and goodness-of-fit (G). The -NS (north and south aspect) and -met (meteorological station) affixes indicate which measured shortwave radiation was used in modelling as radiation was not measured for the east and west aspects. test data RMSE (mm) MBE (mm)/RBE (%) NSE G Calibration - loo basicTI [44.30] Pellicciotti [15.46/14.74] basicTI [0.37] basicTI [0.00] Calibration - all basicTI [37.26] HBV-EC [2.35/2.24] basicTI [0.56] basicTI [0.00] Barkerville ASP basicTI [6.84] basicTI [2.57/22.11] basicTI [0.29] basicTI [0.00] TidBit E/W aspect basicTI [65.61] Hock-met [9.01/2.93] basicTI [0.72] basicTI [0.00] It is difficult to interpret the measures of error among the different test data sets. Model performance can only be assessed relative to other models within each test data set. The NSE value is affected by a number of factors, and high values can result even for a poor fit, for example if the sample variance is large (McCuen et al. 2006). When considering the Nash-Sutcliffe criterion, the larger the amount of variance in the observed data, the larger the denominator and the smaller the ratio. Therefore, all other things being equal, greater variance in the observed data results in a better NSE score as long as the modelled-observed variance is smaller than the observed variance. Small differences between the observed and predicted melt values also results in a higher NSE score. If there is small sample variance and/or a large difference between the observed and predicted values the NSE will be smaller. Thus, interpretation of good or bad fit depends on sample size. The comparison of errors should only be used within each test data set due to variation in the time scales used. The Barkerville ASP had daily data, the Mt Tom snow courses were measured approximately weekly, and the TidBit data were only available as a season total value. The RMSE is strongly dependent on the number of observations in the data set. Therefore it is difficult to compare the RMSE between the Barkerville ASP results 100 (which used 30 data points) and the TidBit results (which used 4 data points). For example, the RJvlSE ranged from ca. 26 - 65 mm for the season total TidBit data whereas it was less than 7 mm for the Barkerville ASP daily measurements, yet the daily equivalent for the TidBit sites was ca. 0.9 - 3 . 2 mm. Model performance at these study sites appeared to improve as the time scale of comparison increased. Micovic and Quick (2009) observed that a simple runoff model performed well at a daily simulation time step but a more complex model structure was necessary for good performance at an hourly time step. They concluded that optimal model complexity depended on the simulation period and computational time step. At the Barkerville ASP site melt is predicted and compared at the daily scale, the Mt Tom snow course sites were compared at the weekly time scale and the Mt Tom TidBit sites were compared at the season total scale. The best model performance was seen at the TidBit sites and the poorest was seen at the Barkerville sites. The Barkerville ASP site was also the furthest from the snow course sites where the models were calibrated, and it was located in a small opening. The sloped snow course sites were used to calibrate all the models thus the model parameters should be the most accurate for those sites and should give better performance. In addition, the models were calibrated to approximately weekly melt totals from the snow course sites, yet melt was modelled at the daily scale (and summed up to appropriate totals as required). 3.4.3 Benefit of Including Shortwave Radiation In this study, the inclusion of shortwave radiation did not appear to improve model performance. Only the Hock model used shortwave radiation to predict melt, as a zero value 101 radiation factor was calibrated for the Pellicciotti model. Neither the measured nor modelled radiation version of the Hock model consistently out-performed the other models. In Hock (1999) and Zappa et al. (2003), including potential shortwave radiation improved model performance. The three modified temperature-index models (HBV-EC, Hock, and Pellicciotti) had similar performance at each site, and, based on visual inspection and the measures of error, performed better than the basic temperature-index model. The Hock model run using modelled shortwave radiation performed better than the run using measured radiation. This may be due to the steady increase in daily modelled radiation over time, which mirrors the increase in melt over time (see section 3.3.3.3 Hock Model). Measured radiation had more variation due to changes in cloud cover and lacked a strong pattern. Since the model is calibrated to shortwave radiation the relationship does affect the parameters however, the close, increasing relationship between melt and modelled radiation may just be coincidental. Due to the low frequency of snow course measurements in this study, modelling was done at the daily scale. The Hock (1999) and Pellicciotti et al. (2005) studies looked at the benefit of shortwave radiation on hourly melt modelling. Diurnal cycles in melt are controlled by shortwave radiation (Martinec 1989; Munro 1990; Bales et al. 1993). It may be that shortwave radiation has more influence on melt variability at an hourly scale compared to a daily scale. An hourly scale incorporates the effect of night and day on shortwave energy input. In addition, in glacier studies such as Hock (1999) and Pellicciotti et al. (2005) the melt period is longer as it incorporates snow and then ice melt as well as the progression of melt up the glacier. The melt period at Mt Tom (< 30 days) may not have been long enough to be clearly influenced by the increasing trend in daily average incoming radiation. 102 The model run using measured shortwave radiation measured at the snow course sites and the meteorological station is more representative of the Mt Tom sites than for the Barkerville ASP site due to distance and differences in cloud cover. Differences in slope, aspect, and shading for the east and west aspect Tidbit sites could also have affected the representativeness of the measured radiation. On the other hand, modelled radiation was simulated directly for each site. In addition, the Pellicciotti model includes net shortwave radiation whereas the Hock model only includes direct incoming radiation. The inclusion of albedo in the Pellicciotti model may have reduced the energy of the incoming shortwave such that it was no longer important for melt modelling and thus no radiation factor was calibrated for that model. 3.4.4 Sources of Error Precipitation data were a large source of uncertainty in this study. Accurate precipitation data were necessary for several modelling steps, yet precipitation measured at Mt Tom was unusable and the representativeness of the Environment Canada station data to the study site was questionable (see section 3.3.1 and Appendix 2). Daily observations and manual measurement of snowfall on-site (Jackson and Prowse 2009) would have been beneficial. However, this would not be practical particularly for application to larger areas. Albedo was another source of uncertainty in this study. Albedo varies over time with snowfall, grain size, impurity concentration, and exposure of vegetation and ground features, and generally decreases as the melt progresses resulting in greater absorbance of shortwave radiation. Including only the incoming shortwave radiation misses this important modifier to the actual energy received at a given site. Pellicciotti et al. (2005) concluded that using albedo in their temperature-radiation index model eliminated the need for a seasonally 103 varying melt factor. Unfortunately the albedo measurements and simulations available for this study could not be verified, and error in albedo would have affected the Pellicciotti model results. Measured on-site albedo would have been ideal, and albedo should have been measured in several locations in the open and forest sites to be representative. However, this would not have been feasible given time and financial constraints and would not be feasible for regional melt modelling. Alternatively, a study to test the applicability of the albedo model for my study site would be beneficial. The late season accumulation of SWE at the north aspect open site could not be explained by any of the precipitation records. In other studies the effect of radiation differences due to aspect was only seen during the melt period and did not affect the accumulation of SWE (Murray and Buttle 2003; Watson et al. 2006). However, Jost et al. (2007) found that aspect did affect snow accumulation in a partly forested watershed. The late season accumulation at my site may have resulted from drifting snow. Drifting snow can significantly affect SWE distribution and this should be considered in future studies even for areas that may not appear to be prone to blowing snow (Dery et al. 2010). Data collection at remote sites during winter is difficult. Thus testing and calibration of models on detailed data sets is a daunting task for large areas or areas with dangerous terrain or difficult access. It can be difficult and costly to collect the amount of data required to calibrate or test models using snow surveys, unless the surveys are already established for an area. The majority of snowmelt modelling studies use river discharge data to test the model's performance, often supported by point snow course or snow stake measurements of change in SWE. River discharge data enable snowmelt runoff testing at the sub-daily time step and provide a continuous measure of runoff. The snow melt lysimeters in this study were 104 supposed to fill this role, supplying a continuous measure of snow melt runoff at the plot scale. Unfortunately due to equipment failure these data were not available. River discharge incorporates several dampening effects for a watershed and can mask differences in melt model performance (Rango and Martinec 1995). The snow melts as well as evaporates in situ, water then runs through and beneath the snowpack. Runoff can be absorbed by the underlying surface or evaporated when exposed to the air at any point on route to the river gauge. In addition the time lag between the release of water until it reaches the measuring gauge can vary. Many hydrologic processes are not well understood so that if discharge records are used for calibration, parameterization of the real physics of snowmelt and runoff are not captured. Even if considerable errors are introduced in the melt simulation the basin discharge can still be sufficiently accurate (Rango and Martinec 1995). Kustas et al. (1994) observed a reduction in the RMSE values of nearly 50% when models were used to predict runoff in comparison to lysimeter outflow predictions. Testing models at the plot scale is different than testing at the watershed scale using river discharge. Kane et al. (1997) found that different coefficient values were obtained when the model was calibrated to snow course measurements versus discharge measurements. There is greater variability in the melt at the plot scale which is not uniformly due to specific slope, aspect or cover conditions (Sivapalan 2003). The slope varies depending on the scale being examined and micro-topography has a considerable impact at the plot scale. Shading from surrounding trees is important when dealing with small clearings, and shrubs, stumps, and slash piles all contribute to melt variability within each plot. Watson et al. (2006) found a burned forest site to have the greatest SWE variation due to large woody debris. For the forested sites in this study forest cover varied between the two 105 aspects, with larger spacing and a more mature forest on the north aspect compared to the south. The influence of edge effects was not quantified in this study. A number of studies have observed the effect of forest edges on snow accumulation and melt both in the clearing and in the adjacent forest (e.g. Golding and Swanson 1986; Stather et al. 2001; Spittlehouse et al. 2004). Shading from the surrounding trees affects the amount of shortwave radiation received within forest openings, and this effect is greatest for the north- and south-facing edges of the opening (Spittlehouse et al. 2004). Golding and Swanson (1986) observed significantly lower SWE on the north edge of all openings larger than 3/4H in diameter. In my study average tree heights were only recorded for the main north and south aspect sites and the meteorological station site. The north aspect and meteorological openings were 3H in the east-west direction and 6H in the north-south direction. The south aspect opening was 5-6H in both directions. For both south and north aspects the furthest east and west snow course measurements were taken at a distance less than one tree height from the forest edge (Figure 2.4 and Figure 2.5). The furthest south and north snow course measurements were taken at a distance of greater than one tree height from the forest edge (data not presented). The nearest south aspect forest snow course sample points were taken less than one tree height distance in to the west-facing edge. The nearest north aspect forest snow course sample points were taken just short of one tree height distance into the eastfacing edge. In a survey of several studies comparing forested and open areas Spittlehouse et al. (2004) concluded that most of the influence on the microclimate occurs within one tree height (H) distance on either side of the forest edge. Thus it is likely that a portion of the snow course measurements for both open and forest sites on both aspects were affected by 106 the forest edge. However, Golding and Swanson (1986) observed that SWE in the east and west sides of an opening were nearly identical and Spittlehouse et al. (2004) observed that the east and west edges of openings received similar energy for melt. Since the sample points within one tree height distance were on the east and west edges of the opening and forest, the edge effects may be less. However, the influence of edge effects on snow accumulation and melt should be assessed for this study area in the future. The ASP was located in an 8 m diameter opening in the forest. Tree height information was not available for the station, but if tree heights are assumed similar to the nearby Mt Tom area, the opening is approximately 1/3H in diameter. I treated the ASP site as a forested site for modelling purposes. Spittlehouse et al. (2004) and Golding and Swanson (1986) found that small clearings (~1H and smaller) had a microclimate similar to the forest and no difference in ablation among the different edges of the opening. This information supports the assumption that snow melt would proceed as for a forest site for the ASP. However, the influence of edge effects on snow accumulation and melt should be assessed for this study area in the future. Finally, the number of observations of SWE may not have been sufficient. Random effects within the plot scale may have overwhelmed the variation due to aspect and vegetation, even though the coefficient of variation data indicated that the variability could be represented by the number of samples used. The 95% confidence intervals indicate that there were few periods where the change in SWE differed significantly among the snow course sites (Figure 3.6). Watson et al. (2006) found that random effects were greatest at scales less than 100 m and were larger than the effect of radiation and vegetation. They concluded that intense sampling of SWE (up to four times as many samples as were taken in 107 this study), would be necessary to be able to isolate the effects of radiation and vegetation on SWE at their study site. In addition it would have been beneficial to know the date of snow disappearance for each site as this introduced uncertainty in the calculation of melt rates for the final observation period. The models did not perform well for the accumulation period. This probably resulted from a combination of the poor precipitation data and the lack of snowpack temperature or ripeness being considered in these models. These temperature-index models are restricted by their inability to assess snowpack ripeness. When only applied to the continuous melt period they performed better compared to when they were applied to the entire period from 1 April, 2008 onward as seen in comparison of Figure 3.11 with Figure 3.16 and Figure 3.12 with Figure 3.17. During early April no change in SWE was observed on site, yet the models predicted melt whenever air temperatures exceeded the threshold temperature. Although surface melt may have been occurring it refroze in the pack without resulting in loss of SWE. This resulted in over-prediction of melt and the end of the melt predicted earlier than it was observed. 3.5 Conclusion The optimal melt model could not be defined for this study. However, the HBV-EC and Hock models consistently gave good performance. Considering the effect of slope, aspect and/or using separate melt factors for open and forested areas appeared to improve model performance without requiring any additional meteorological data. It is common in glacier studies to use separate melt factors for snow and ice due to their different melt properties (e.g. Table 1, Hock 2003), and other studies have used separate land-cover based melt factors (e.g. Hamlin et al. 1998). 108 The idea behind including shortwave radiation is attractive: shortwave radiation is a large energy term; inclusion of shortwave radiation reduces melt factor variability, and may bypass the need for a seasonally varying melt factor. The inclusion of shortwave radiation in temperature-index models has been seen to improve melt simulations in other studies (Cazorzi and Dalla Fontana 1996; Kane et al. 1997; Dunn and Colohan 1999; Hock 1999; Pelhcciotti et al. 2005). It may need to be applied at the sub-daily scale to be informative for short melt seasons such as in this study. The strength of other non-radiative factors in this study could be causing the lack of an apparent effect of shortwave radiation. The inclusion of shortwave radiation is complicated by interactions with vegetation and the balance between long and shortwave radiation. There is a considerable body of research on canopy radiation models and the inclusion of these effects in snowmelt models (e.g. Ellis and Pomeroy 2007; Pomeroy et al. 2008; Esseiy et al. 2009; Pomeroy et al. 2009; Rutter et al. 2009). A number of input data uncertainties, which were primarily due to lack of sufficient measurements or instrumentation, may have confounded the results. However, one of the purposes of this study was to compare these melt models for use in operational runoff models. It is not feasible to make exhaustive measures of these variables for entire watersheds. In addition, site SWE measurements may have been insufficient to evaluate the influence of aspect, forest cover, and radiation on melt in this study. Future research should incorporate these snow melt algorithms into a runoff model and test them against watershed discharge data. 109 4. Energy balance modelling at Mt Tom 4.1 Introduction There are two basic types of snowmelt models. Temperature-index (or "degree day") models simulate melt based on an empirical relation between air temperature and ablation. Energy balance models are based on fundamental physical principles and calculate melt from the energy balance equation (Equation 1.1). Energy balance models are more accurate than temperature-index models due to their physical basis; however they typically require large amounts of input data to run (See Chapter 1, section 1.3). 4.1.1 Outline and Objectives The additional data required by energy balance models combined with uncertainties in input data can cause greater uncertainties compared to simpler temperature-index models (Franz et al. 2008). However, if the model can be run without calibration it is an advantage. Areas of interest for snowmelt modelling are often remote and may not have adequate river discharge or snow water equivalent (SWE) data to calibrate melt factors. A simple energy balance model that does not require input data beyond that of a basic temperature-index model and that does not require calibration is appealing. This model would be easy to run and easy to transfer among sites without concern about the transferability of melt factors from site to site or year to year. The goal of this study is to assess the performance of a simple energy balance model that does not require calibration (Walter et al. 2005) at a partly forested site with varying slope and aspect conditions. The specific objectives of this study are to: (1) Assess the ability of the simple energy balance model to predict melt at a partly forested site with an undulating surface; and 110 (2) Compare the performance of a simple energy balance model against various temperature-index models. 4.2 Methods and analysis For site and data collection information please see Chapter 2 (sections 2.1 to 2.3). 4.2.1 Model The simple energy balance model created by Walter et al. (2005) requires only daily maximum and minimum air temperatures and precipitation as input data. Wind speed is required, but a constant value can be substituted if daily measurements are not available (Walter et al. 2005). Relative humidity can be used if it is available, but the model can be run without it. The model assumes the snowpack behaves as a single slab of uniform temperature. All the terms necessary to run an energy balance model are approximated using these variables. Dr. M.T. Walter3 shared the new code for his simple energy balance model. This code is based on his publication (Walter et al. 2005), but differs on several accounts. This new code has not been extensively tested, although comparison with data from Danville, Vermont indicates it is performing reasonably well (Walterpers. comm.) with assumed theoretical and empirically-derived constants (Table 4.1). Assistant Professor, Department of Biological and Environmental Engineering, Cornell University, Ithaca, NY 111 Table 4.1 Constant values used in the current version of the Walter et al. (2005) model. Value Units kJ m"2 d"1 Solar constant 117500 Snow emissivity 0.94 - Latent heat of vapourization 2500 kJ kg"1 Latent heat of freezing 333.3 kJ kg"1 Heat capacity of snow 2.1 kJ kg"1 °C"1 Heat capacity of air 1.25 kJ m"3 °C"1 Heat capacity of water 4.2 kJ kg"1 °C"1 Density of water 1000 kgm" 3 Maximum albedo 0.98 - Bare ground albedo 0.25 - Stefan-Boltzmann constant 4.89 x 10'6 kJ m"2 K"4 d"1 Thermodynamic vapour constant 0.4615 kJ kg"1 K"1 For the original calculation of resistance to heat and vapour transfer Height of zero-plane displacement 0.0 m Momentum roughness parameter 0.001 m Heat and vapour roughness parameter 0.0002 m Von Karman's constant 0.41 - 4.2.1.1 Radiation calculations: Shortwave radiation (K [kJ m"2]) was calculated as in Walter et al. (2005) using modelled albedo (A), calculated atmospheric tiansmissivity (T,), and daily potential extraterrestrial solar radiation (S0 [kJ m"2]). K = {l-A)T,So 112 (4-1) Longwave radiation (L [kJ m" ]) was calculated with the Stefan-Boltzmann equation as in Walter et al. (2005). L = soTAK (4-2) Emissivity (e) was assumed to equal 0.94 for snow and was calculated for the atmosphere using cloud cover fraction. Cloud cover fraction was back-calculated from the calculated atmospheric transmissivity (Walter et al. 2005). The Stefan-Boltzmann constant (a [kJ m~2 K" d"1]) is listed in Table 4.1. Temperature (TK [°C]) was assumed to be equal to daily average air temperature to calculate incoming longwave radiation and was assumed to be equal to daily average snow surface temperature to calculate outgoing longwave radiation. The remaining energy balance equations were also those in Walter et al. (2005) except for changes to the turbulent flux calculations discussed below. 4.2.1.2 Latent heat flux: The saturation vapour density of the surface and air are calculated according to Walter et al. (2005). The resistance to vapour exchange and the latent heat of vapourization were replaced with a wind function (Walter pen', comm.) such that, E = 86400 x (ps - pa)x (5.3(1 + W)) (4-3) Where E is the latent heat flux [kJ m" d" ], ffis the average daily wind speed [m s" ], andp s andpo are the vapour densities at the surface and air respectively [kg m" ]. If relative humidity (RH) is available, then the vapour density of air can be calculated using the saturation vapour density (p0) calculated with average daily air temperature (Ta [°C]) as: 113 pa=RHxp0[Ta] (4-4) Otherwise it is assumed that the minimum air temperature (T„ [°C]) approximates the dew point temperature and the vapour density of air is calculated as the saturation vapour density calculated with the minimum daily air temperature: Pa=P0[Tn] (4-5) 4.2.1.3 Sensible heat flux: Surface and air temperature were used as in the paper. The resistance to heat exchange term was replaced with a wind function (M.T. Walterpers. comm.) H = 86400xC a x{Ts - T a ) x ^ \ + W ^ Where H is the latent heat flux [kJ m"2 d"1], W is the average daily wind speed [m s"1], Ca is the heat capacity of air [kJ m"3 °C"'] from Table 4.1, Xv is the latent heat of vapourization [kJ kg" ] from Table 4.1, and Ts and Ta are the temperatures at the surface and air respectively [°C]. 4.2.2 Running the Model The model was run both with and without using relative humidity in the calculation of latent heat flux, and with and without the new turbulent heat flux calculations. The relative humidity option was not used for the Barkerville ASP since relative humidity was not measured at the ASP station. The model was initially run starting 1 April, 2008, for each site, then was subsequently run from the start of active melt for each site (Table 2.5). The start of the melt was defined as the date of observed peak SWE. 114 (4-6) 4.2.3 Input data The Mt Tom meteorological station provided daily average air temperature, maximum and minimum air temperatures, wind speed, and relative humidity data for the snow course sites. The model runs for the Barkerville ASP used daily temperature from the ASP station. Daily precipitation totals from the Barkerville Environment Canada station were used for all study sites. Snow survey and automatic snow pillow SWE were used to compare against model output. Analysis and preparation of observed meteorological and SWE data were executed as outlined in section 2.3. Initial SWE was taken as the Mt Tom snow course or Barkerville automatic snow pillow SWE [m] measured on the start date for the model run. Initial snow surface temperature was -5°C, the same value used in Walter's code (M.T. Walterpers. comm.), but an initial albedo of 0.7 was used instead of 0.5 from Walter's code. A higher initial albedo was chosen based on the high albedo measured later in the season at the site and because it was a deeper snowpack than in Walter's study. Predicted SWE was adjusted through time within the model based on simulated melt and measured precipitation. I used daily geometric wind speed for the Mt Tom snow course sites since the data were available at my station and recommended by Walter et al. (2005). Station records showed periods where wind speed was equal to zero. Since zero values cannot be used in calculating geometric wind speed I used 0.0001 m s" in place of the 0 m s~' measurements. Since wind speed was not measured at the Barkerville ASP I used season total geometric mean wind speed measured at my station as discussed in Walter et al. (2005). In the Walter model forest cover is incorporated using the fraction of forest cover between 0 and 1. Forest cover fraction was obtained from the adjacent Ministry of Forest 115 and Range study area (Teti 2010, unpublished data). Forest cover was interpreted using Gap Light Analyzer from fisheye photographs taken at 1.1 m above the ground (Teti 2008). The averages of the pre-harvest percent cover within a 30 degree zenith angle for each south and north aspect were used as model input (Teti 2003). The forest cover fraction used for the south aspect site was 0.60 and for the north aspect forest site it was 0.46. No forest cover measurement was available for the Barkerville ASP. The site was described as 'mature forest' (Scott Jackson RFC,pers. comm.). I decided to use the average measured forest cover values from Mt Tom and used a value of 0.53 for the forest cover fraction at the Barkerville ASP. The effect of forest cover is only taken into consideration in the calculation of shortwave radiation in the model. Forest cover is applied much like albedo as (1-forest) multiplied by the potential incoming shortwave radiation. Neither wind speed nor longwave radiation is altered by forest cover and thus is identical for all sites. 4.2.4 Testing the Model Observed SWE and melt at the Mt Tom snow course sites and the Barkerville ASP site were compared to the modelled SWE and melt for each of the model variations. Snow melt was calculated as the difference between SWE at subsequent snow survey dates. Evaluation of model performance was based on graphical representations of modelled and observed SWE data as well as the calculation of root mean square error (RMSE), mean bias error (MBE), relative bias error (RBE), and the Nash-Sutcliffe efficiency criterion (NSE) and goodness-of-fit (G) (Table 3.1). The G index is used to evaluate model performance against that of a benchmark model. In this study the basic temperature-index model was used as the benchmark (see section 3.2.5). 116 4.3 Results 4.3.1 Model Performance 4.3.1.1 Open Sites Model performance from 1 April, 2008, was poor (Figure 4.1 a and b). Observed precipitation was not sufficient to balance the simulated melt early in the season. south open » north open 9 * ® 9 9 •a ® © © « « « • Apr 01 Predicted swe Observed swe April Apr 21 * May 01 May 11 May 21 a) * b) Predicted swe Observed swe ® May 31 Apr 01 April Apr 21 May 01 May 11 May 21 May 31 Figure 4.1 Walter model run starting 1 April 2008. Comparison of observed (circles) versus modelled (line) SWE over time at a) south aspect open site, and b) north aspect open site. When the model was run from the start of the continuous melt for each site (Figure 4.2 a and b), model performance improved. When relative humidity was used to calculate the latent heat flux model performance improved for both the north and south aspect (Figure 4.2 c and d). The melt proceeded more slowly which resulted in a better estimation of SWE over time. 117 north open south open V E, \ o £ a * a) Predicted swe Observed swe May 01 May 11 May 21 —* May 31 May 11 May16 May21 May26 May31 Jun 05 north open south open V b) Predicted swe Observed swe * \. « N. * Predicted swe Observed swe May 01 * \ May 11 c) « May 21 d) Predicted swe Observed swe May 31 Mjy11 May 16 May 21 May 26 May 31 Figure 4.2 Comparison of observed (circles) versus modelled (line) SWE over time, a) model run from start of melt, 27 April 2008, at the south aspect open site, b) model run from stall of melt, 9 May 2008, at the north aspect open site, c) model run from 27 April, 2008 at the south aspect open site using relative humidity to calculate the latent heat flux, and d) model run from 9 May, 2008 at the north aspect open site using relative humidity to calculate the latent heat flux. Figure 4.3 a and b show the effect of using the original latent heat flux equations. Again, melt is slowed and SWE prediction appears somewhat more representative than in the original model (Figure 4.2 a and b). The additional use of relative humidity with the original turbulent flux calculations minimizes melt even more providing the most realistic SWE for the south aspect site (Figure 4.3 c), although on the north aspect it appears that melt is too low and results in an over-prediction of SWE. 118 Jun 05 north open south open Predicted swe Observed swe May 11 May 01 Predicted swe Observed swe a) May 21 May 11 May 31 May 16 May 21 May 26 May 31 Jun 05 north open south open \ b) 9 X^ 9 \ « —• Predicted swe Observed swe Ma/01 Predicted swe Observed swe c) M'iy 11 Ma/21 May 11 May 31 May 16 d) May 21 Ma/26 May 31 Figure 4.3 Comparison of observed (circles) versus modelled (line) SWE over time. The original equation incorporating resistance to heat transfer to calculate the turbulent fluxes was used (Eq. 12; Walter et al. 2005). a) model run from start of melt, 27 April, 2008, at the south aspect open site, b) model run from start of melt, 9 May, 2008, at the north aspect open site, c) model run from 27 April, 2008 at the south aspect open site using relative humidity to calculate the latent heat flux, and d) model run from 9 May, 2008 at the north aspect open site using relative humidity to calculate the latent heat flux. The measures of error for each open site were calculated using only four melt values. Thus comparison of the error values among the south and north aspect in Table 4.2 is reasonable, but comparison with other sites that have a larger data set may be inaccurate. Running the model from the start of the continuous melt period for each site resulted in a clear improvement in model performance (Table 4.2). The RMSE is reduced from 0.10 m to 0.08 m for the south aspect and is reduced from 0.14 m to 0.06 m for the north aspect. The 119 NSE values improve somewhat, but still remain negative for both sites indicating a poor fit between the model and observations (Table 4.2). Negative values of NSE can result if the predictions of a linear model are biased (McCuen et al. 2006). Table 4.2 Measures of error calculated for each variation of the Walter model and each open site. Values in bold indicate the best error measure overall and values in bold italic indicate the best error measure for the other site. SO = south aspect open; NO = north aspect open. '1 April' indicates that the model was run starting 1 April, 2008, and 'Melt' indicates that the model was run only for the melt period for each site (starting 27 April, 2008 for the south and 9 May, 2008 for the north). 'RH' indicates that relative humidity was used to calculate the latent heat flux and 'rv' indicates that resistance to heat transfer was used to calculate the turbulent heat fluxes as in the paper (Eq. 12; Walter et al. 2005). N = 4. Model Run RMSE [m] MBE [m] RBE [%] NSE [-] SO, lApril, 4 melt periods 0.1015 -0.0324 -32.22 -1.39 NO, lApril, 4 melt periods 0.1398 -0.1218 -83.01 -7.02 SO, lApril, all periods 0.0780 -0.0014 -2.58 -0.36 NO, lApril, all periods 0.1000 -0.0154 -25.75 -0.72 SO, melt 0.0814 0.0040 3.92 -0.54 NO, melt 0.0614 -0.0096 -6.57 -0.55 SO, melt, RH 0.0488 0.0040 3.92 0.45 NO, melt, RH 0.0440 -0.0096 -6.57 0.20 SO, melt, rv 0.0516 0.0040 3.92 0.38 NO, melt, rv 0.0460 -0.0207 -14.10 0.13 SO, melt, rv, RH 0.0393 0.0012 1.23 0.64 NO, melt, rv, RH 0.0500 -0.029 -19.79 -0.02 The observed and predicted melt are plotted against each other clearly showing when melt was over- or under-predicted for each site and model variation (Figures 4.4 - 4.6). Based on the RMSE and NSE calculated measures of error (Table 4.2) the 'best' model performance is for the south aspect site using the model run from the start of the melt using the original turbulent flux equations and using relative humidity to calculate the latent heat flux. The RMSE is 0.04 m and the NSE is 0.64. The 'best' model performance for the north 120 aspect site occurred when the model was run from the start of the melt using the new turbulent flux calculation and using relative humidity to calculate the latent heat flux The RMSE is also 0 04 m although the NSE is only 0 20 0 00 0 05 0 10 0 15 0 20 0 00 obs^r%pd melt(m) 0 05 0 10 0 15 0 20 obsened melt(m) Figuic 4 4 Picdicted veisus obseived melt foi the observed snow couise peuods Model lun staiting 1 April 2008 The conesponding peuods of companson are indicated on the figure Diagonal lines indicate a 1 1 relationship a) model lun foi south aspect open site, only the foui melt periods aie shown, b) model urn foi noith aspect open site, only the foui melt peuods aie shown, c) model lun for south aspect open site, all periods fiom 1 Apnl to end of obseived snow covei aie shown, and d) model urn for north aspect open site, all peuods from 1 April to end of obseived snow cover are shown 121 0 25 north open / o May 23-May 30 May 17 - May 23 / o May^ May 17 / May 30 - Jun 5 b) 000 005 0 00 0 05 010 015 020 0 00 0 05 0 10 0 1? 0 20 0 00 0 05 observed melt (m) 0 10 015 0 20 0 25 0 10 0 15 0 20 0 25 observed melt (m) Figure 4.5 Predicted versus observed melt for the observed snow course periods. The corresponding periods of comparison are indicated on the figure. Diagonal lines indicate a 1:1 relationship, a) model run from start of melt, 27 April, 2008, at the south aspect open site, b) model run from stait of melt, 9 May, 2008, at the north aspect open site, c) model run from 27 April at the south aspect open site using relative humidity to calculate the latent heat flux, and d) model run from 9 May at the north aspect open site using relative humidity to calculate the latent heat flux. 122 ^ 0 00 1 1 1 1— 0 05 0 10 0 15 0 20 ooseivpij melt (m) ^ 1 1 1 1 r~ 0 00 0 05 0 10 0 15 0 20 0 25 observed melt (m) Figure 4.6 Predicted versus observed melt for the observed snow course periods. The original equation incorporating resistance to heat transfer to calculate the turbulent fluxes was used (Eq. 12; Walter et al. 2005). The coiresponding periods of comparison are indicated on the figure. Diagonal lines indicate a 1:1 relationship, a) model run from start of melt, 27 April 2008, at the south aspect open site, b) model run from start of melt, 9 May 2008, at the north aspect open site, c) model run from 27 April at the south aspect open site using relative humidity to calculate the latent heat flux, and d) model run from 9 May at the north aspect open site using relative humidity to calculate the latent heat flux. 4.3.1.2 Forest Sites Model performance at the south aspect forest site was improved by using measured relative humidity to calculate the latent heat flux (Figure 4.7b). Use of the original turbulent flux calculation also improved melt performance (Figure 4.7c) relative to the original model 123 (Figure 4.7a). However, using both the relative humidity and the original turbulent flux calculation resulted in over-prediction of SWE later in the season (Figure 4.7d). CO _ o CM _ o o o © o o o \s in CD _ O o $ May 01 Predicted swe Observed swe May 11 a) « May 21 May 01 May 11 0 May 31 May 01 c) Predicted swe Observed swe May 21 b) Predicted swe Observed swe May 11 May 21 May 31 d) Predicted swe Observed swe May 31 May 01 May 11 May 21 May 31 Figuic 4 7 Time series of SWE using the measured foiest covei fi action at the south aspect site, foi a) model lun fiom stall of melt, 2 May 2008, b) using relative humidity to calculate latent heat flux, c) using the oiigmal equation mcoipoiating tesistance to heat liansfer to calculate the turbulent fluxes, and d) using the resistance to heat tiansfci and lclalivc humidity to calculate the fluxes Model performance at the north aspect forest site was also improved by using measured relative humidity to calculate the latent heat flux (Figure 4.8b) However, use of the original turbulent flux calculation or both the relative humidity and the original turbulent flux calculation resulted in over-prediction of SWE later in the season (Figure 4.8c,d) 124 CO _ CO o o in CN _ m o o \ * o \ * o CN o o m m o o o _ \« o ® o o o » Predicted swe Observed swe 9 i May 11 i May 16 in a) o __ o _ o _ o o i May 21 i May 26 i May 31 \ Predicted swe Observed swe « b) * ® m r Jun 05 i May 11 Predicted swe Observed swe i May 16 i May 21 i May 26 i May 31 r Jun 05 Predicted swe Observed swe 1 1 1 1 1 r 1 i i i i r May 11 May 16 May 21 May 26 May 31 Jun 05 May 11 May 16 May 21 May 26 May 31 Jun 05 Figuie 4 8 Tunc seaes of SWE using the mcasuied toicsl covci fiaction at the noith aspect site, loi a) model lun fiom stait of melt, 2 May 2008, b) using relative humidity to calculate latent heat flux, c) using the onginal equation mcoipoidting lesistance to heat tiansfei to calculate the turbulent fluxes, and d) using the lesistance to heat tiansfei and lclative humidity to calculate the fluxes The measures of error foi each forest site were calculated using only three melt values Thus companson of the error values among the south and north aspect in Table 4.4 is reasonable, but comparison with other sites that have a larger data set may be inaccurate The model run with lelative humidity gave the best measuies of eiror for both aspects, with smallei RMSE and laigei NSE values relative to the othei model runs (Table 4 3) The model mn using the original tuibulent flux calculation still gave modeiate performance for 125 the south aspect site with an NSE of 0.54. The 'best' model performance for the north aspect site only gave an NSE of 0.28, but the RMSE was lower than for the south aspect site. Table 4.3 Measures of error calculated for each variation of the Walter model and each forest site using measured forest cover for each site. Values in bold indicate the best error measure overall and values in bold italic indicate the best error measure for the other site. SF = south aspect forest; NF = north aspect forest. 'Melt' indicates that the model was run only for the melt period for each site (starting 2 May, 2008 for the south and 9 May, 2008 for the north). 'RH' indicates that relative humidity was used to calculate the latent heat flux, and 'rv' indicates that the turbulent heat fluxes were calculated as in the paper (Eq. 12; Walter et al. 2005). N = 3. Model Run Forest Cover [-] RMSE [m] MBE [m] RBE [%] NSE [-] SF, melt 0.60 0.0310 -0.0079 -11.83 0.27 NF, melt 0.46 0.0288 0.0008 0.90 -1.75 SF, melt, RH 0.60 0.0151 -0.0079 -11.83 0.83 NF, melt, RH 0.46 0.0147 0.0008 0.90 0.28 SF, melt, rv 0.60 0.0248 -0.0163 -24.41 0.54 NF, melt, rv 0.46 0.0229 -0.0194 -20.87 -0.75 SF, melt, rv, RH 0.60 0.0305 -0.0248 -37.15 0.29 NF, melt, rv, RH 0.46 0.0310 -0.0277 -29.90 -2.20 Using zero forest cover resulted in a large over-prediction of melt and subsequent under-prediction of SWE over time at the forest sites. To examine the sensitivity of the model to forest cover fraction 1 ran the model with forest cover set at increments of 0.1.1 also determined the optimal forest cover fraction based on RMSE in this manner. As forest cover was increased the melt decreased, which improved the model prediction to a point, after which melt was too small and over-prediction of SWE resulted (Figure 4.9). The optimal forest cover changed depending on the model variation used. Using relative humidity to calculate the latent heat flux resulted in a lower optimal forest cover for each site (Figure 4.9 b and d), and for larger forest cover values the SWE would not melt by the end 126 of the season. Using both the resistance function and relative humidity resulted in extremely low melt such that much of the SWE was still on the ground at the end of the melt period (Figure 4.10). south forest • Observed swe V \ May 11 May 12 Ma/26 May 11 M3y16 May 31 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% * \\ \ May 01 - * May 21 Ma/21 b) May 31 M iy 2b Figuic 4 9 Time seues of SWE foi each 10% incicmcnt in forest covei between 0 and 1 for a) model lun fiom stait of melt, 2 May 2008, at the south aspect toicst site, b) model uin from stait of melt at the south aspect foiest site using lelative humidity to calculate latent heat flux, c) model uin fiom stait of melt, 9 May 2008, at the noiih aspect forest site, and d) model um fiom stait of melt at the noith aspect foicst site using lelative humidity to calculate latent heat flux 127 * Observed swe south forest - V 0% 10% 20% 30% 40% 50% b0% 70% 80% 90% 100% south forest ° Observed swe V- — * \ 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% \ \ X W \ May 01 May 11 v * May 21 a) X X. May 31 May 11 y 01 north forest b) . \ May 21 May 31 * Observed swe north forest Nr-—_ xX \ E „ — 0% - c) May 11 May 10 Ma/ 21 Ma/31 May11 \ X \ x^ * \ 20% 30% 40% 50% 60% 70% 80% 90% 100% X X \ \ X? May16 May21 \ V May26 Ma/31 L Figuie 4 10 Time scnes of SWE for each 10% incicment in foiest covci between 0 and 1 The onginal equation incorpoiatmg icsistancc to heat tiansfei to calculate the tuibulcnt fluxes was used (Eq 12, Waltci et al 2005) a) model inn fiom stait of melt, 2 May 2008, at the south aspect foiest site, b) model uin fiom stait ot melt at the south aspect foiest site using lelative humidity to calculate latent heat flux, e) model lun fiom stait ot melt, 9 iViay 2008, at the i.oith aspect foie>t sAc, and ci) model uin l.om stait of melt at the noith aspeet lotcst -.ite using lelative humidity to calculate latent heat flux For each model variation the optimal forest cover was lower for the north forest compared to the south (Table 20). For each progressive model variation the optimal forest cover decreased (Table 4.4). The 'best' north aspect forest site simulation was obtained usin^ the basic model with a forest cover of 0.7. This also lesulted in the smallest overall RMSE (0.011 m). The 'best' south aspect forest site simulation was obtained using the basic model and relative humidity to calculate latent heat flux with a forest cover of 0.6. This also 128 d) resulted in the largest overall NSE (0 83) In comparison to the measured forest cover fraction values, the model run using relative humidity gave the closest optimized forest cover values - measured south aspect forest was 0 6 and optimized was 0 6, measured north aspect forest was 0 46 and optimized was 0 5 (Table 4 4 SF and NF "melt, RH") Note that the optimization was only done for increments of 0 1 Using the optimized forest cover of 0 7, the basic Walter model run gave the lowest RMSE and highest NSE for the north aspect forest site (NSE = 0 61)- better values than using the measured forest cover of 0 46 For the south aspect forest site the model run with measured relative humidity gave the best results using both the measured and optimized forest cover Table 4 4 Measures of cnor calculated foi each variation of the Waltei model and each forest site using only the optimal foiest cover for each site and model vaiiation Values in bold indicate the best enoi measure overall and values in bold italic indicate the best eiroi measure for the othci site SF = south aspect foiest, NF = north aspect foiest 'Melt' indicates that the model was tun only for the melt period for each site (starting 2 May 2008 for the south and 9 May for the noith) 'RH' indicates that relative humidity was used to calculate the latent heat flux, and 'rv' indicates that the turbulent heat fluxes weie calculated as in the paper (Eq 12, Walter et al 2005) N = 3 Model Run Forest Cover [-] RMSE [m] MBE [m] RBE [%] NSE [-] SF, melt 08 0 0179 -0 0092 -13 82 0 76 NF, melt 0.7 0.0108 -0 0048 -5 15 0.61 SF, melt, RH 0.6 0.0151 -0 0079 -11 83 0.83 NF, melt, RH 05 0 0128 0 0008 0 89 0 46 SF, melt, rv 05 0 0222 -0 0079 -11 83 0 62 NF, melt, rv 03 0 0129 -0 0012 -1 25 0 44 SF, melt, rv, RH 05 0 0223 -0 0138 -20 70 0 62 NF, melt, rv, RH 02 00151 0 0008 0 89 0 24 129 4.3.1.3 Energy Balance Each of the energy balance terms was output from the model. These are compared to available measured data and among sites to give some additional msight on the performance of the Walter model. The modelled incoming shortwave radiation to the surface tracks the measurements made parallel to each open aspect site (Figure 4.11), however, modelled shortwave for the south aspect is under-predicted throughout the melt season, whereas for the north aspect it is over-predicted early in the season. o o o South aspect LO o o o LO ON o o o o LO o o o o 10 CO Measured--*- Walter o o — LO l l l l l l l 4pr01 April Apr 21 May 01 May 11 May 21 May 31 North aspect ^~ o _ i . LO Ol -r> v ~' < ~ l cz o to TJ ro o o o LO LT (11 CO e o sz CO o o \$fm LO Measured Walter I l l 1 I I l Apr 01 April Apr 21 May 01 May 1 1 May 21 May 31 Figure 4.11 Comparison of the measured and modelled incoming shortwave radiation to the surface. The top panel is for the south aspect open site and the bottom panel is for the north aspect open site. The potential incoming shortwave was adjusted for transmissivity using the daily air temperature range. 130 Net energy flux to the snow surface is positive from the start of the simulation for almost every site and model variation (bottom panels, Figure 4.12, Figure 4.13). One exception is the south aspect forest site for the original model variation run (Figure 4.13a, bottom panel). Thus melt is not predicted for two periods near the beginning of May. Otherwise melt is continuous for all sites. The largest energy input occurs around the middle of May with high values of sensible and latent heat. Equally high energy input occurs at the south aspect sites near the end of May. Modelled SWE had completely melted at these sites by this time, which resulted in low albedo values and high shortwave radiation. After 10 May, 2008, the latent and sensible heat fluxes were similar among the four sites for a given model variation since the minimum daily air temperature remains above 0°C after this point. The calculation of latent heat flux is dependent on both relative humidity and the resistance term thus the values for latent heat vary for each model run (top panels Figure 4.12 and Figure 4.13). Using relative humidity has a noticeable effect on the latent heat flux. At times the latent heat flux changes from a large positive value to a negative value (top panels Figure 4.12 and Figure 4.13). For most days using measured relative humidity resulted in a smaller calculated latent heat flux compared to that predicted using minimum daily air temperature. Using relative humidity resulted in smaller values and less variability in the air vapour density (used to calculate the latent heat flux); the correlation between the two methods of calculating air vapour density was around 0.4. Use of the resistance term also decreased the latent heat flux. The calculation of the sensible heat flux is only dependent on the resistance term thus only two variations are seen (middle panels Figure 4.12 and Figure 4.13). Use of the resistance term also tends to decrease the sensible heat flux. The patterns seen in the latent 131 heat flux for the model variations that use relative humidity are strongly influenced by the minimum daily air temperature. The patterns are nearly identical when latent heat flux is greater than zero, but when the latent heat flux is less than zero it tends to mirror the daily minimum temperature instead (results not shown). For the models that do not use relative humidity the pattern in latent heat flux is not tied to air temperature. As expected, the sensible heat flux also follows air temperature, but the net energy is also strongly related to air temperature patterns. 132 V?7 V M~t,rv Melt Melt.RH MB./V,RH ,•21 May 11 A ^ V tr XWATW May 31 fvA a) ^ A^ 1 May 31 May 16 May 26 May 31 T Jun 05 Figure 4.12 Latent, sensible, and net energy flux for each model run for a) the south aspect open and b) the north aspect open sites. 'Melt' indicates model was run from start of melt, 'RK' indicates that relative humidity was used to calculate the latent heat flux, and 'rv' indicates that the resistance to transfer term was used to calculate the turbulent fluxes. 133 T May 11 May 01 Ma/11 May 31 Ma/16 1 May 21 May 26 Melt Melt.RH Melt,rv Melt,rv,RH May 26 May 11 May 16 May 21 I 1 May 11 May 16 May 21 May 26 May 31 JunOS May 31 May 31 Jun05 Figure 4.13 Latent, sensible, and net energy flux /or each model run for a) the south aspect forest and b) the north aspect forest sites. 'Melt' indicates model was run from start of melt, 'RH' indicates that relative humidity was used to calculate the latent heat flux, and 'rv' indicates that the resistance to transfer term was used to calculate the turbulent fluxes. 134 4.3.2 Barkerville Automatic Snow Pillow The Walter model was also run for the Barkerville automatic snow pillow (ASP). As seen with the Mt Tom snow course sites, running the model for the period of active melt improved model performance. Model performance was good using the new turbulent flux calculation (hereafter called 'new Walter') (Figure 4.14a), but melt was underestimated using the original turbulent calculations (hereafter called 'Walter rv') resulting in an overestimation of SWE over time (Figure 4.14b). b) a) g o Predicted swe Observed swe Predicted swe Observed swe "i 1 i i i i r May 01 May 06 May 11 May 16 May 21 May 26 May n i i i i i Figure 4.14 Comparison of observed (circles) versus modelled (line) SWE over time at the Barkerville ASP. Model was run only for the melt period 1-May to 30-May, 2008 using a) the new equation to calculate turbulent heat fluxes, and b) the original equation incorporating resistance to heat transfer to calculate the turbulent fluxes (Eq. 12; Walter etal. 2005) The goodness-of-fit indices indicate that the new Walter model performed better than the basic temperature-index model, but the Walter rv model performed more poorly than the basic temperature-index model (Table 4.5). I also ran the model for each 10% increment in forest cover to examine the effect of forest cover on model performance. The forest cover increment that resulted in the 'best' model results for the new Walter model was 50%, close to the estimated value of 53% used to run the model. However, the 'best' model results for 135 r May 01 May 06 May 11 May 16 May 21 May 26 May 31 the Walter rv model run were obtained using a 30% forest cover as opposed to the estimated value of 53% used to run the model. Table 4.5 Calculated measures of error for the Walter model run at the Barkerville ASP site. 'Melt' indicates model run only for melt period (1 May - 30 May), and "rv" indicates the model run using the original turbulent flux calculation. Root mean square error (RMSE), mean bias error (MBE), relative bias error (RBE), NashSutcliffe efficiency criterion (NSE), and goodness-of-fit (G). N = 30. Model run RMSE [m] MBE [m] RBE [%] NSE G Melt 0.0048 -0.0008 -7.24 0.65 0.52 Melt, rv 0.0076 -0.0049 42.17 0.13 -0.23 4.3.3 Energy Balance vs. Temperature-Index Models 4.3.3.1 Mt Tom Snow Course Sites The Mt Tom snow course sites were also used to calibrate and test a series of temperature-index models (see Chapter 3). The measures of error were calculated for both the Walter model runs and the temperature-index models using the results from all four snow course sites pooled together (Table 4.6 and Table 4.7). The range in RMSE values tended to be slightly higher (by 0.01 to 0.02 m) for the Walter model compared to the temperatureindex models. In most cases the NSE values for the Walter model were smaller compared to the temperature-index models. A negative MBE was calculated for the Walter model whereas a positive value was calculated for all the temperature-index models. Melt was overpredicted at the open sites and under-predicted at the forest sites with the Walter model, but resulted in an overall under-prediction of melt. Finally, the G indices indicated that most versions of the Walter model performed similarly to the basic temperature-index model with values close to zero (Table 4.6). The temperature-index models all performed better than the basic TI at the snow course sites with a minimum value of G of 0.25. Based on these 136 measures of error it appears that the modified temperature-index models outperformed the simple energy balance model in predicting melt for the Mt Tom snow course sites. Table 4 6 Calculated measuies of eiror for each Walter model run using all the snow course sites combined Root mean squaie error (RMSE), mean bias eiror (MBE), relative bias enor (RBE), Nash-Sutcliffe efficiency criterion (NSE), and goodness-of-fit index (G). The basic temperature-index with all-sites-combined coefficient was used as the benchmark series N = 14. Model run RMSE [m] MBE [m] RBE [%] NSE G Melt 0.0554 -0.0046 -4.42 0.02 -1.21 Melt, RH 0.0363 -0.0031 -3.00 0.58 0.05 Melt, rv 0.0388 -0 0067 -6.41 0.52 -0.08 Melt, rv, RH 0.0362 -0.0107 -10 23 0.58 0.06 137 Table 4 7 Calculated measures of error for each model for both leave-one-out and all-sites-combined calibration methods Subscnpt "surK" indicates the model run with potential suiface shoitwave ladiation Nonsubscnptcd model runs used measuied shoitwave radiation Root mean square eiror (RMSE), mean bias error (MBE), lelative bias eiror (RBE), Nash-Sutchffe efficiency criterion (NSE), and goodness-of-fit indices (G) The basic temperature-index model (BTI) was used as the benchmark series N = 14 for each loo and all methods Leave-one-out calibration method MODEL RMSE [m] MBE [m] RBE [%] NSE G basicTI 0 0443 0 0011 1 05 0 37 0 HBV-EC 0 0320 0 0060 5 69 0 67 0 48 Hock 0 0323 0 0087 8 34 0 66 0 47 Hocksur|< 0 0318 0 0077 7 32 0 68 0 48 Pellicciotti 0 0383 0 0155 14 74 0 53 0 25 Pellicciotti surK 0 0324 0 0080 7 60 0 66 0 46 All-sites-combined calibration method MODEL RMSE [m] MBE [m] RBE [%] NSE G basicTI 0 0373 0 0022 2 15 0 56 0 HBV-EC 0 0243 0 0024 2 24 0 81 0 58 Hock 0 0262 0 0019 1 82 0 78 0 50 HocksurK 0 0260 0 0017 1 66 0 78 0 51 Pellicciotti* 0 0266 0 0020 1 95 0 77 0 49 Pellicciotti surK * 0 0266 0 0020 1 95 0 77 0 49 * The calculated mcasuics of enor aie the same between the Pellicciotti model um with measuied shoitwave ladiation and that lun with modelled potential suiface shortwave ladiation because with all sites used in calibration shoitwave radiation was not used in determining melt, thus, the type of radiation used is inconsequential to model peifoimance The tempeiatme-index models appealed to follow variations in the rate of change in SWE bettei than the eneigy balance model Foi example, the pattern of observed SWE at the noith open aspect site showed an acceleration followed by slowing in melt rate (Figure 4 15c versus Figure 4 16c) The eneigy balance model appealed to pi edict the SWE bettei foi the foiest sites than seen from the temperatuie-index models and the tempeiature-index models 138 appeared to perform better for the open sites. For example, the energy balance model predicted a rapid initial decline in SWE for the south aspect open site that was not observed (Figure 4.16a), whereas this was not predicted by any of the temperature-index models (Figure 4.15a). south open ° — — a) Observed Basic temperature index HBV-EC Hock measured K Hock modeled K Pellicciotti May 01 May 11 M3y21 May 31 May 01 May 11 May 21 May 31 north forest d) "X \ May 11 May 16 May 21 May 26 May 31 May 11 May 16 ^ v, \ "^^N. May 21 May 26 May 31 Figure 4.15 Observed and modelled snow course SWE over time using variations of the temperaturc-index model, a) south aspect open; b) south aspect forest; c) north aspect open; d) north aspect forest. 139 south open V south forest a) b) o • o • E o ~ 5 > o >a.f • Obseived Melt Melt RH Melt iv Melt iv RH May 01 w^.. ^X o \ V. t May 11 May 21 • Xr— o o May 31 May 11 May 01 May 11 May 11 May 10 May 21 May 21 May 26 May 31 May 31 Figuie 4 16 Obscivcd and modelled snow couise SWE ovei time using vaiiations of the Waltei et al (2005) model a) south aspect open, b) south aspect toiest, c) noith aspect open, d) noith aspect foicst 4 3 3 2 Barkerville Automatic Snow Pillow All models, except for the Waltei rv model, perfoimed bettei than the basic tempeiatme-index model (Table 4 5 and Table 4 8) The measmes of error weie similai among all models The RMSE values weie slightly smaller foi the tcmpeiatuie-mdex models (mean 0 0059 m) compaied to the Waltei model mns (mean 0 0062 m) The NSE and goodness-of-fit values were slightly higher for the Walter model mn with the new tuibulent flux calculations compaied to the tempeiatuie-index models 140 Ju 05 Table 4 8 Calculated measures of error using coefficient set 5 to predict melt at the BarkerviUe ASP site Subscupt "surK" indicates the model run with potential suiface shortwave radiation Non-subscripted model runs used measuied shortwave radiation Root mean square error (RMSE), mean bias error (MBE), relative bias en or (RBE), Nash-Sutcliffe efficiency criterion (NSE), and goodness-of-fit (G) N = 30 MODEL RMSE [m] MBE [m] RBE [%] NSE G basicTI 0 0068 0 0026 22 11 0 29 0 HBV-EC 0 0055 -0 0001 -0 44 0 54 0 35 Hock 0 0058 -0 0012 -9 90 0 53 0 33 HocksurK 0 0058 -0 0012 -9 90 0 53 0 33 Pellicciotti* 0 0058 -0 0012 -9 90 0 53 0 33 * The Pellicciotti calculated measures of error using modelled potential surface shortwave radiation were identical and therefore are not shown here The basic temperature-index model over-predicted melt at the BarkerviUe ASP resulting in a predicted snow free date occurring eaiher than it was observed, the HBV-EC and new Walter models predicted the snow-free date accuiately, and the remaining models under estimated the snow-free date (Figure 4 17) The new Walter model had the closest fit to the observed SWE (Figure 4 17a) a-P_y •SK, ^ ^ ^ c a) o 0 — Observed Melt Melt rv May01 May06 May11 May16 May21 May 26 May 31 Observed Basic temperature-index HBV-EC Hock measured K Hock modeled K Pellicciotti \ T~ T~ T i r May 01 May 06 May11 May16 May 21 May 26 May 31 Figure 4 17 Obseivcd and modelled SWE ovei time toi the period of continuous melt foi the Barkeivillc ASP a) Waltei eneigy balance model lesults, and b) tempeiatuie index model lcsults 141 4.4 Discussion 4.4.1 Model Performance The Walter model tended to overestimate melt for certain sites and certain model variations. Melt was over-estimated in particular for the south aspect open site resulting in too early of a predicted snow-free date. However, the overall mean bias error was negative for several model runs where the decline in simulated SWE proceeded more rapidly than observed. This was primarily due to the over-estimation of melt resulting in an early snowfree date and thus a period of under-prediction of melt at the end of the season when simulated melt was zero yet observed melt was positive. Over-prediction of melt was also observed in several other energy balance studies (Whitaker et al. 2003; Pellicciotti et al. 2008; Zeinivand and De Smedt 2009). On the other hand, melt was often under-predicted in general for the north aspect open site and the forest sites. The model run using the original turbulent heat flux calculations was prone to under-predicting melt for these sites. This under-prediction is obvious for the Barkerville ASP. This simple energy balance model essentially assumes an isothermal snowpack (M.T. Walter, pers. comm.). Therefore as soon as the snow surface temperature reached 0°C, melt commenced. This early season melt resulted in large model errors at all sites for the initial model run starting 1 April since simulated energy balance conditions conducive to surface melt occurred long before the snowpack had ripened. Pellicciotti et al. (2008) observed the same overestimation in melt with a single-layer energy balance model that assumed an isothermal snowpack. In reality, before the snowpack is ripe the surface layer may melt during the day, but the meltwater percolates into the deeper cold snow and refreezes. This releases latent energy into the snowpack that raises the snow temperature (Dingman 2002). 142 Alternatively, once melt has commenced sub-zero overnight temperatures cool the surface layer and it could take several hours for the snow to warm and resume melting the following day (Tseng et al. 1994). This latter occurrence is a concern mostly for sub-daily time steps. In this study model performance was thus drastically improved by running the model from the start of the melt period for each site. The need to know when the snowpack is isothermal complicates the running of the model. If this involves periodic snow surveys, the application of a model to predict snow evolution is somewhat redundant as it is being measured in the field. Remotely sensed microwave data as brightness temperature or backscatter have been used to identify snowmelt onset (Wang et al. 2008; Takala et al. 2009). A potentially easier approach is the addition of a second snowpack layer to the model structure (Marks et al. 2002; Whitaker et al. 2003; Garen and Marks 2005; Watson et al. 2006; Pellicciotti et al. 2008). The top layer changes temperature rapidly in response to diurnal changes to the energy balance while the internal layer adjusts more slowly (Garen and Marks 2005). Once the entire snowpack reaches a temperature of 0°C melt commences. This produces a more realistic snowpack response and hopefully removes the need to measure when the snowpack becomes isothermal. There are methods to track the energy content of the snowpack with a single layer approach as well (e.g. Coughlan and Running 1997; Watson et al. 2006; Franz et al. 2008;Debele^a/. 2009). Using measured relative humidity improved model performance for all snow course sites in this study. This improvement makes physical sense since using relative humidity provides a more accurate measure of the vapour pressure of the air than assuming that minimum daily air temperature is equal to the dew point temperature. In the paper by Walter 143 et al. (2005), the atmospheric vapour density was consistently over-estimated by the model, and this was likely due to the minimum daily air temperature actually being higher than the dew point temperature. Using relative humidity at Mt Tom resulted in a smaller modelled latent heat flux to the snowpack and eliminated the extreme latent heat flux inputs to the pack occurring on days with high minimum air temperatures. The resulting simulated daily melt is lower and is a better reflection of the observed melt. Relative humidity was not measured at the Barkerville ASP station thus minimum daily temperature had to be assumed to represent the saturation temperature. Using the original equation for the turbulent fluxes, which included the resistance to heat transfer term (Eq. 12; Walter et al. 2005), suppressed the daily melt even more than using relative humidity alone. The original equation for the turbulent fluxes also resulted in smaller fluxes than using the new function. In Walter et al. (2005) they observed that their simple energy balance model tended to over-predict the turbulent fluxes. Given this tendency for over-prediction, the original function may be a better choice. Using the original resistance function improved model performance for the south aspect open site, but inhibited melt too much for the other three snow course sites and the Barkerville ASP. Thus, SWE was overestimated using the resistance to heat transfer for four of the five test sites. The model appears to perform better for the forest sites in comparison to the open sites. The lowest RMSE and highest NSE values for the forest sites were lower and higher than for the open sites. The model also appears to perform better for the south aspect sites in comparison to the north aspect. Better RMSE and NSE values are seen for the south aspect. Melt was often over-predicted for the south aspect and well- or under-predicted for the north site. Thus, model variations (such as using the original turbulent flux calculation) that 144 lowered modelled melt improved the performance for the south aspect but made the overprediction of melt for the north aspect worse. Hence the 'best' model variation for the north aspect was not the same as that for the south aspect. Using the original turbulent flux calculations and relative humidity decreased the modelled daily melt for the forest sites as seen for the open sites. The measured forest cover fraction was smaller for the north aspect than for the south. This is in agreement with field observations of the two sites. The forest on the south aspect is a younger, denser stand (Figure 4.18a) whereas on the north aspect it is an older spruce and fir forest with large gaps between trees (Figure 4.18b). When forest cover was optimized the optimal forest cover was different for each model variation. For model variations where changes in the code decreased the melt, the optimized forest cover decreased. Forest cover also lowered melt and therefore less cover was required as the model output decreased. The optimal forest cover was lower for the north forest compared to the south forest for each model variation. The optimized forest cover fraction was the same as the measured forest cover fraction for the south aspect site, but not for the north aspect. A higher forest cover fraction in conjunction with the original model variation, which output larger melt values, gave a better prediction of melt for the north aspect site than using the measured forest cover fraction. Much as for temperatureindex models, using calibration parameters can improve energy balance model performance by tuning the model output to specific site conditions (Rutter et al. 2009); thus it is expected that better simulations were obtained from optimized values. 145 Figure 4 18 The Mt Tom foiest snow course sites a) south aspect and b) noith aspect A A 9 *« >H^n ^rt >- W ' I l l " ^ When the Walter model was urn fiom 1 April, 2008, the individual eneigy balance outputs weie sometimes difficult to inteipret on a day-to-day basis froi example, on some days when the average air temperatuie was below zero the modelled snow surface tempeiature was still equal to 0°C, SWE deci eased, the net energy balance was positive, and the sensible heat flux was positive This helps to explain laige oveiestimation of melt when the model was urn from 1 Apnl, 2008, which is still within the accumulation penod 146 For a given model variation the net energy flux varies among the four sites only as a result of changes in the shortwave radiation flux. The shortwave radiation is affected by forest cover, slope, aspect, and albedo. All the other energy inputs are identical among the sites. Since the modelled snow surface temperature remains above 0°C during the continuous melt there is no feedback among the energy fluxes which is why they are identical among all the sites despite changes in the shortwave radiation. The simulated melt thus varies as a result of shortwave radiation, the initial SWE, and the date on which the model is started for each site. In this model forest cover is used in the calculation of net shortwave radiation only. Therefore the observed modification to model output for the forest sites is due solely to changes in the net shortwave radiation. In the model equations, forest cover is applied much like albedo as (1-forest) multiplied by the potential incoming shortwave radiation. Thus a 10% increase in forest cover translates into a 10% reduction in shortwave radiation to the surface. This is likely an invalid assumption as shortwave radiation is not reduced linearly by the canopy density. This is valid only when the sun is directly above the canopy, and then only if there is no reflection by the canopy. During the winter, when solar angles are typically low and snow is often retained in the canopy, the reduction in shortwave radiation by the canopy is greater than the percent canopy cover due to additional shading and increased canopy albedo. There is no accounting for the effect of forest cover on the turbulent or longwave fluxes. Longwave radiation to the snow surface is increased due to contribution from the surrounding canopy, trunk, branches, and shrubs (Pomeroy et al. 2006; Pomeroy et al. 2008), and the increase in net longwave can even outweigh the reduction in shortwave due to the canopy (Sicart et al. 2008). Wind speed in the forest is lower than in 147 the open resulting in lower turbulent fluxes (Hardy et al. 1997, Storck et al. 1999, Suzuki et al 1999, Koivusalo and Kokkonen 2002). Even if the effect of canopy on longwave is considered, the canopy and trees are not necessarily at air temperature as commonly assumed (Pomeroy et al. 2009). Insolation intercepted by the canopy increases the temperature above air temperature resulting in greater longwave radiation input. This effect is enhanced with increasing shortwave radiation and should be considered for canopy models (Pomeroy et al. 2009). Thus consideration of the effect of forest cover on the net longwave radiation and the sensible heat flux are of particular importance for snowmelt modelling. In their study Walter et al. (2005) were able to compare the modelled energy balance terms to on-site measurements. Unfortunately only incoming shortwave radiation in the open was measured at Mt Tom and therefore there is no way to validate individual energy components. It would be beneficial to compare the individual outputs to measured values in a future study. Comparison between the simulated and measured incoming shortwave radiation at Mt Tom indicated good agreement. Walter et al. (2005) also observed that shortwave radiation was well simulated, although with a slight tendency for under-prediction by the model. The atmospheric longwave radiation in the Walter et al. (2005) study was strongly underestimated while the terrestrial longwave was overestimated, but to a lesser degree. Atmospheric stability was not considered in the Walter et al. (2005) equations for sensible and latent heat exchanges (equations 11 and 13). Stable atmospheric conditions, where the temperature gradient near the surface is less steep than the adiabatic lapse rate, suppress turbulent transfer (Dingman 2002). These conditions typically occur over snowpacks because warm air is overlaying a cold surface, a condition that strongly represses 148 the turbulent fluxes. The Richardson number (Oke 1987) is one method to incorporate atmospheric stability. Measurements of wind speed and temperature at two heights are compared to determine the atmospheric stability then the appropriate correction is applied (Oke 1987; Dingman 2002). Use of atmospheric stability in this simple energy balance model would likely decrease the simulated turbulent fluxes bringing those more in line with observations. 4.4.3 Energy Balance vs. Temperature-Index Models Since the Mt Tom snow course sites were used for both calibration and validation of the temperature-index models the optimal performance from the temperature-index models was expected for these Mt Tom sites. Walter et al. (2005) found during their study that this simple energy balance model outperformed the basic temperature-index model based on correlation and standard error values. In my study the results were mixed. The temperatureindex models using all-site-combined calibration (all) had smaller RJV1SE values compared to the Walter model. When one of the snow course sites was left out during calibration of the temperature-index models and used as an 'independent' data set for testing, the RMSE values were larger compared to the all calibration method and closer to those calculated for the energy balance model. Based on NSE, RMSE, and G indices, the new Walter model results were better than all of the temperature-index models for the Barkerville ASP. However, the Walter rv model results were poorer than all the models at the Barkerville ASP. At the snow course sites all temperature-index models were better than or equivalent to the Walter model performance. Studies by Kustas et al. (1994) (open snowfield), Walter et al. (2005) (open), and Singh et al. (2009) (70% forested) compared energy balance models of varying complexity 149 with basic temperature-index models and concluded that the energy balance models performed better. Ichii et al. (2008) (Columbia River Basin) found that a simple energy balance model based on Walter et al. (2005) performed considerably better than a net radiation-temperature-index model. Studies by Kustas et al. (1994) (open snowfield), Kane et al. (1997) (arctic), Pellicciotti et al. (2008) (glacier), and Singh et al. (2009) (70% forested) concluded that energy balance model results were approximately the same as those for temperature-index models. It is notable that all but one of the temperature-index models considered in these studies had some sort of modification from the basic approach. Zeinivand and De Smedt (2009) compared a slightly modified basic temperature-index model with a simple energy balance model using only minimum, maximum, and average daily temperature, wind speed, and precipitation and found both provided similar results with model efficiencies greater than 80%. Zappa et al. (2003) (open basin) compared an energy balance model with three variations of the basic temperature-index model and found that all four models provided good discharge simulations but the shortwave radiation based temperature-index model had the highest model efficiencies. The energy balance model simulated SWE slightly better than the shortwave radiation temperature-index model, but they concluded that the increase in computational time outweighed the small increase in model performance (Zappa et al. 2003). Only one study that I found concluded that the energy balance model was inferior to the temperature-index model (Franz et al. 2008) (forested). However, model perfonnance depended on the year of simulation. For some years the temperature-index model performed better than the energy balance model, on other years the energy balance model performed better than the temperature-index model, and for other 150 years the performance was nearly the same (Franz et al. 2008). The single year of comparison done in this study is not sufficient to confidently conclude the 'best' model type. The energy balance model performed better for the forest sites than the temperatureindex models. The temperature-index models generally overestimated melt at the forest sites. The energy balance model considers the impact of forest cover on incoming shortwave radiation; melt can still respond to aspect effects and daily variation in radiation due to cloud cover. Thus the melt simulated for the forest sites was appropriately reduced on a day-to-day basis using the energy balance model. The temperature-index models, however, use a reduced melt factor for both forest sites that does not depend on differences in aspect or daily changes in insolation. The temperature-index model appeared to simulate melt better for the open sites than the energy balance model. This result could be due to the use of melt factors calibrated specifically for the Mt Tom snow course sites. In a comparison study of 33 snowpack models model performance was more consistent for open sites compared to forest sites and simulation of forest site SWE without calibration was worse than for the open sites (Rutter et al. 2009). The Barkerville ASP site was independent from the temperature-index calibration data and also had less input data for the energy balance model. The ASP site provided information on how the two model types performed when being transferred to another site away from the measurement site. Snow melt at the Barkerville ASP site was reasonably simulated by almost all models considered. All models performed better than the basic temperature-index model except for the Walter rv model, which under-estimated melt. The new Walter model simulated melt better than the modified temperature-index models at the 151 Barkerville ASP, but the HBV-EC model gave close results. The NSE values at the Barkerville ASP were not as large as calculated for the snow course sites. Snow melt at the Barkerville ASP is being compared at a daily time step whereas for the snow course sites it is being compared to approximately weekly melt totals. Perhaps the temperature-index models are too simple to be applied as accurately at a daily time step. In addition the Barkerville ASP site was furthest away from the snow course sites where the models were calibrated, and it was located in a flat, forested area. The temperature-index models, however, were calibrated to sloped sites with south and north aspects and were calibrated with approximately weekly melt totals. Perhaps the difference in calibration site conditions resulted in the poorer performance for the temperature-index models when transferred to the Bakerville ASP. This is typical of temperature-index models; they only work well for the location where they were calibrated. The energy balance model does not require calibration and should in theory perform as well for either the snow course or Barkerville ASP sites. As discussed in section 3.4.2, comparison of model performance among the different test sites is not straightforward since different data sets and sizes were used for each. However, based on visual inspection of the figures of SWE over time the Walter model appears to perfonn better at the Barkerville ASP site in comparison to its performance at the Mt Tom snow course sites. The Barkerville ASP site also had less input data than the Mt Tom sites - no measured RH or windspeed. 4.5 Conclusion The Walter model overestimated melt for the south aspect open site; thus model adjustments that decreased the melt tended to improved model performance for that site. However this resulted in under-prediction of melt for the other sites. Model performance at 152 Mt Tom was clearly improved for all sites by incorporating observed relative humidity to calculate the latent heat flux. Relative humidity is commonly measured along with air temperature at remote stations and should be used in the model when available. The Walter model was limited by the assumption of an isothermal snowpack. It was necessary to know the date at which the pack at each location became ripe to obtain reasonable melt simulations. Use of a two-layer model or the inclusion of corrections to account for the energy content of the entire snowpack would be beneficial (Watson et al. 2006). The individual energy balance components could not be compared to measured terms. However, since the relative relationship among the energy balance terms was the same as in Walter et al. (2005) it is reasonable to worry that the energy balance tenns in this study were not well simulated. As discussed in both Walter et al. (2005) and Watson et al. (2006) the inaccuracies within the modelled energy fluxes are cancelling each other out resulting in good model output for the wrong reasons. If this is the case the argument for using an energy balance model to more realistically represent the physical processes is no longer accurate. Overall this simple energy balance model reasonably simulated the SWE for open and forest sites on north and south aspects and for a flat site. Some model variations performed better than others for certain sites. Although melt simulations were not as good, periods of over- and under-prediction of melt cancelled each other out resulting in fairly good overall SWE prediction for this study. Other studies using similar simple energy balance models also reported satisfactory simulation of ablation and SWE (Ichii et al. 2008, Zeinivand and De Smedt 2009). 153 Based on the measures of error used in my study, the Walter model performed similarly to the basic temperature-index model, and the modified temperature-index models had better performance in general than the Walter model. In the study by Walter et al. (2005) the simple energy balance model outperformed the basic temperature-index model. Since the modified temperature-index models performed as well or better than the energy balance model and did not require any more information to run than the basic temperature-index model, it is thus debatable whether it is better to use a physically based model where a number of assumptions had to be made and there were substantial discrepancies in the estimates of the individual energy components (Walter et al. 2005, pg 73), or an empirical index-based model which relies on site-specific calibration. Further research should attempt to quantify the errors associated with the simple energy balance approach and with the various temperature-index models. How much error is accumulated through the assumptions in the energy balance approach? Which model type has the smallest uncertainty associated with the final melt simulation? 154 5. Conclusion 5.1 Optimal Melt Model Complexity The optimal temperature-index model complexity could not be determined for this study as individual model performance varied for each test data set. No single model performed better in all situations, but some models worked well for specific cases. Incorporating the effect of slope aspect and/or forest cover in some manner improved model performance in comparison to the basic temperature-index model without necessarily requiring any additional meteorological data. The simple energy balance model tended to over-predict melt for the open snow course sites, in particular the south aspect open site, and tended to slightly under-estimate melt for the forest sites. Model results varied for the different model runs and model variations provided improved performance for some sites. The use of measured relative humidity, as opposed to the assumption of dewpoint temperature, gave better results for almost every model run. Since relative humidity is commonly measured along with air temperature at remote stations it is recommended that this be used in the model whenever possible. Melt simulations were not as accurate as the basic temperature-index model. Yet, periods of overand under-prediction of melt cancelled each other out resulting in fairly good overall simulation of SWE decline over time for this study. The simple energy balance model may be obtaining an acceptable result for the wrong reasons. The optimal model type could not be defined for this study. Based on the measures of error used in my study, the simple energy balance model performed similarly to the basic temperature-index model and the modified temperature-index models had better performance in general than the basic temperature-index and energy balance model. In the 155 study by Walter et al. (2005) the simple energy balance model outperformed the basic temperature-index model. Since the modified temperature-index models performed as well or better than the energy balance model and did not require any more meteorological data to run than the basic temperature-index model, it is thus debatable whether it is better to use this Walter et al. (2005) physically based model where a number of assumptions had to be made and there were substantial discrepancies in the estimates of the individual energy components (Walter et al. 2005, pg 73), or an empirical index-based model which relies on site-specific calibration. The assumptions that were used for the study sites in Walter et al. (2005) may not be applicable to the Mt Tom area. These assumptions should be tested for other locations. The temperature-index models actually required more data than the simple energy balance model because of the need for calibration data. The amount of calibration data used for this study was minimal and more data over multiple years are usually used. It was much easier to apply the simple energy balance model in this study as only basic meteorological data were required, and different melt factors were not calibrated for each site type within the study area. Not requiring calibration is an advantage, making a simple energy balance model quick and easy to set up and run and easily transferable among sites. In addition there was no need for subjective modification of the calibration parameters, which can occur with temperature-index models. However, calibration of some parameters, such as the relationship between incoming shortwave radiation beneath the canopy and canopy density, might improve model performance. Calibration is frequently used for energy balance modelling. 156 The energy balance model was expected to perform better than the calibrationdependent temperature-index models when transferred to an independent site outside the study area. One version of the simple energy balance model performed better than the temperature-index models when used to predict melt for the Barkerville ASP site. However, the other version of the simple energy balance model performed more poorly than all other models considered and under-estimated the observed melt. 5.2 Data Availability and Model Performance In this study the inclusion of shortwave radiation did not appear to improve model performance. Of the three modified temperature-index models tested, only one incorporated shortwave radiation directly, yet all three models had similar performance. The simple energy balance model also performed similarly to models that did not include shortwave radiation. Other non-radiative factors may have overwhelmed the effect of shortwave radiation at the site. The use of either measured or modelled radiation data therefore did not obviously affect model performance. Albedo was only considered in the simple energy balance model due to temperature-index model calibration. Including only the incoming shortwave radiation misses this important modifier to the actual energy received at a given site. Better measurement or simulation of albedo may have changed the calibration outcome and perhaps a difference between measured and modelled radiation would have been clear if albedo were considered. 5.3 Model and Data Limitations Both model types were limited by the inability to track snowpack ripeness. Knowledge of the date on which the snowpack became ripe for each site was therefore necessary to run 157 both of these model types. There are various ways to track snowpack temperature or ripeness that could be utilized in these models either through inclusion of a two-layer snowpack structure or some variation of heat content accounting (Coughlan and Running 1997; Kane et al. 1997; Whitaker et ah 2003; Zappa et al. 2003; Garen and Marks 2005; Watson et al. 2006; Franz et ah 2008; Zeinivand and De Smedt 2009). It is also possible to determine snowpack ripeness through remote sensing methods (Wang et ah 2008; Takala et ah 2009). Input data quality is one of the primary problems for snowmelt modelling (Ferguson 1999). It is important for any modelling exercise to have high quality input data, but energy balance models are more sensitive to input data and thus more affected by data bias than temperature-index models (Zappa et ah 2003, Franz et ah 2008). In addition, the more data required by a model the greater the opportunity for input uncertainties to be propagated through the model and manifested in SWE and melt output (Franz et al. 2008). This is an important consideration when applying these two types of models to an area larger and more complex than the small site used for this study. It was even problematic at this study scale. Precipitation data are a large source of uncertainty in this study. Precipitation measurements are difficult to obtain in the field, particularly in winter, and it is not always appropriate to use the measurements from a nearby weather station. Precipitation is an important input for the albedo model as a light snowfall can increase the albedo considerably and have an important influence on the shortwave energy flux. Daily observations and manual measurement of snowfall on-site (Jackson and Prowse 2009) would have been beneficial. However, this would not be practical particularly for application to larger areas. Albedo was another source of uncertainty in this study. Good albedo measurements were not available and, thus, the representativeness of the modelled albedo values could not 158 be tested. The albedo models used were simple, however, more advanced albedo models can be prohibitively complex and data intensive (e.g. Gardner and Sharp 2010) for application to remote areas. Data collection at remote sites during winter is difficult. Thus testing and calibration of models on detailed data sets is a daunting task for large areas or areas with dangerous terrain or difficult access. It can be difficult and costly to collect the amount of data required to calibrate or test models using snow surveys, unless the surveys are already established for an area. The majority of snowmelt modelling studies use river discharge data to test the models performance, often supported by point snow course or snow stake measurements of change in SWE. The use of river discharge data enables snowmelt runoff testing at the sub-daily time step and provides a continuous measure of runoff. The snow melt lysimeters in this study were supposed to fill this role, supplying a continuous measure of snow melt runoff at the plot scale. Unfortunately due to equipment failure these data were not available. Testing models at the plot scale is different than testing at the watershed scale using river discharge. Kane et al. (1997) found that different coefficient values were obtained when the model was calibrated to snow course measurements versus discharge measurements. There are additional issues, such as the dampening effect, that need to be addressed at the watershed scale. The number of observations of SWE used for this study may not have been sufficient to identify a shortwave radiation effect or to conclusively differentiate among the tested models. Random effects within the plot scale may have overwhelmed the variation due to aspect and vegetation, even though the coefficient of variation data indicated that the variability could be represented by the number of samples used. Watson et al. (2006) found 159 that random effects were greatest at scales less than 100 m and were larger than the effect of radiation and vegetation. They concluded that intense sampling of SWE (up to four times as many samples as were taken in this study), would be necessary to be able to isolate the effects of radiation and vegetation on SWE at their study site. 5.4 Directions for Further Research The model was limited by the assumption of an isothermal snowpack. It was necessary to know the date at which the pack at each location became ripe to obtain reasonable melt simulations. Use of a two-layer model or the inclusion of corrections to account for the energy content of the entire snowpack would be beneficial for this study area. In addition it would be beneficial to know the date of snow disappearance for each site as this introduced uncertainty in the calculation of melt rates for the final observation period. A number of input data uncertainties, which were primarily due to lack of sufficient measurements or instrumentation, may have confounded the results. However, one of the purposes of this study was to compare these melt models for use in operational runoff models. It is not feasible to make exhaustive measures of these variables for entire watersheds. In addition, site SWE measurements may have been insufficient to evaluate the influence of aspect, forest cover, and radiation on melt in this study. Future research should incorporate these snow melt algorithms into a runoff model and test them against watershed discharge data. The input data were easy to apply for this study since all sites were located proximal to the meteorological station. For watershed studies meteorological data have to be extrapolated throughout the watershed. The primary reason quoted for use of temperature-index models in operational forecasting is that air temperature is easy to extrapolate. On the other hand, wind 160 speed and relative humidity, which were required for this simplest of energy balance models, are notoriously difficult to extrapolate (Ferguson 1999). Garen and Marks (2005) give a detailed description of spatial interpolation methods for some common energy balance input variables in their paper. It would therefore be beneficial to compare these models at the watershed scale, taking into consideration the effects of input data interpolation or extrapolation on the performance and error associated with the different model structures. Forest cover clearly influences the energy and mass balance of the snowpack. In this study the effect of forest was considered simplistically. Extensive research has been done on canopy radiation models and the inclusion of these effects in snowmelt models (e.g. Ellis and Pomeroy 2007; Pomeroy et al. 2008; Essery et al. 2009; Pomeroy et al. 2009; Rutter et al. 2009). British Columbia is a mosaic of forest and natural and manmade clearings. A model that incorporates the influence of forest cover should be tested in place of the simple approach taken with these models. Is the extra computation and data required balanced by improved performance? As previously described the models compared in this study were only tested using one season worth of melt data from a small selection of proximal sites. Model performance varies with time, location, and time scale (Kane et al. 1997; Essery et al. 2009; Micovic and Quick 2009; Rutter et al. 2009) thus models should be tested at more sites and with more than one season of data to obtain a better comparison among the different models. This study only gives insight on the relative performance of these models for the 2008 season and for this particular study site. 161 5.5 Summary In conclusion, performance was similar among the basic temperature-index and the simple energy balance models compared in these two studies. The modified temperatureindex models generally performed better than the simple energy balance model for the Mt Tom snow course sites, but performed more poorly in comparison to one version of the simple energy balance model for the independent Barkerville automatic snow pillow data set. There are benefits to the use of each type of model. From a scientific standpoint it is better to move toward a physically-based model for accuracy. In addition, a model that does not require calibration can make application easier especially when considering a number of different sites. From an operational standpoint the cost of obtaining data, maintaining instrumentation and computational efficiency must be considered. Once these models are used within a runoff model, calibration must be performed regardless of the snowmelt model type. 162 6. Literature Cited Aguado, E. 1985. Snowmelt energy budgets in southern and east-central Wisconsin. Annals of the Association of American Geographers 75:203-211. Anderson, E.A. 1968. Development and testing of snowpack energy balance equations. Water Resources Research 4:19-38. Anderson, E.A., 1973. National weather service river forecast system/snow accumulation and ablation model, NOAA Technical Memorandum NWS HYDRO-17, US Department of Commerce, Silver Spring, MD, 217 pp. Armleder, H., M. Waterhouse, T. Newsome, et al. 2002. Mount Tom Adaptive Management project overview: Quesnel Highland Alternative Silvicultural systems. Cariboo Forest Region Research Section Extension Note #37. httpVwww.for.gov.bc.ca/cariboo/research/research.htm. Bales, R.C., R.E. Davis, and M.W. Williams. 1993. Tracer release in melting snow: diurnal and seasonal patterns. Hydrological Processes 7:389-401. Barnett, T., J. Adam, and D. Lettenmaier. 2005. Potential impacts of a warming climate on water availability in snow-dominated regions. Nature 438:303-309. Bergstrom, S. 1976. Development and application of a conceptual runoff model for Scandinavian catchments. Department of Water Resources Engineering. Lund Institute of Technology, University of Lund. Bulletin Series A No. 52. Bergstrom, S., and A. Forsman. 1973. Development of a conceptual deterministic rainfallrunoff model. Nordic Hydrology 4:147-170. Berry, M. 1981. Snow and climate. In: Handbook of Snow: Principles, Processes, Management and Use. D. Gray and D. Male, (Eds.). The Blackburn Press, New Jersey, USA, pp 32-59. Bloschl, G., B. Kimbauer, and D. Gutknecht. 1991a: Distributed snowmelt simulations in an Alpine catchment. 1. Model evaluation on the basis of snow cover patterns. Water Resources Research 27:3171-79. Bloschl, G., Kimbauer, B. and Gutknecht, D. 19916: Distributed snowmelt simuations in an Alpine catchment. 2. Parameter study and model predictions. Water Resources Research 27:3181-3188. Bonsai, B.R., X. Zhang, L.A. Vincent, and W.D. Hogg. 2001. Characteristics of daily and extreme temperatures over Canada. Journal of Climate 14:1959-1976. 163 Boon, S. 2009. Snow ablation energy balance in a dead forest stand. Hydrological Processes 23: 2600-2610. Bradley, R., F. Keimig, and H. Diaz. 2004. Projected temperature changes along the American Cordillera and the planned GCOS network. Geophysical Research Letters 31:1-4. Braithwaite, R.J. 1995. Positive degree-day factors for ablation on the Greenland ice sheet studied by energy-balance modelling. Journal of Glaciology 41:153-160. Brock, B., I. Willis, M. Sharp. 2000. Measurement and parameterization of albedo variations at Haut Glacier d'Arolla, Switzerland. Journal of Glaciology 46:675-688. Brubaker, K., A. Rango, and W. Kustas. 1996. Incorporating radiation inputs into the snowmelt runoff model. Hydrological Processes 10:1329-43. Burford, J.E., S.J. Dery, and R.D. Holmes. 2009. Some aspects of the hydroclimatology of the Quesnel River Basin, British Columbia, Canada. Hydrological Processes 23:12591536. Cazorzi, F., and G. Dalla Fontana. 1996. Snowmelt modelling by combining air temperature and a distributed radiation index. Journal of Hydrology 181:169-187. CHC. 2007. EnSim Hydrologic Reference Manual. Canadian Hydraulics Centre, National Research Council. Ottawa, Ontario, Canada. htlp:/'Www.nrccnrc.gc.ca/obj/chc/doc/kenuc/GreenKenucManual.pdf. Choudhury B., and A. Chang. 1981. The albedo of snow for partially cloudy skies. Boundary Layer Meteorology 20:371-389. Cline, D. 1997. Snow surface energy exchanges and snowmelt at a continental, midlatitude Alpine site. Water Resources Research 33:689-701. Coughlan, J.C., and S.W. Running. 1997. Regional ecosystem simulation: a general model for simulating snow accumulation and melt in mountainous terrain. Landscape Ecology 12:119-136. Daly, S., R. Davis, E. Ochs, and T. Pangbum. 2000. An approach to spatially distributed snow modelling of the Sacramento and San Joaquin basins, California. Hydrological Processes 14:3257-3271. Debele, B., R. Srinivasan, and A.K. Gosain. 2009. Comparison of process-based and temperature-index snowmelt modeling in SWAT. Water Resources Management DOI 10.1007/sl 1269-009-9486-2. 164 Dery, S. J., A. Clifton, S. MacLeod, and M.J. Beedle, 2010: Blowing snow fluxes in the Cariboo Mountains of British Columbia, Canada. Arctic, Antarctic and Alpine Research, 42: 188-197. Dery, S.J., and R.D. Brown. 2007: Recent Northern Hemisphere snow cover extent trends and implications for the snow-albedo feedback, Geophysical Research Letters 34: L22504, doi: 10.1029/2007GL031474. Dery, S.J., K. Stahl, R.D. Moore, P.H. Whitfield, B. Menounos, and J.E. Burford. 2009. Detection of runoff timing changes in pluvial, nival, and glacial rivers of western Canada. Water Resources Research 45, W04426, doi:10.1029/2008WR006975. Dingman, S. 2002. Physical Hydrology 2 nd ed. Prentice-Hall, Inc., New Jersey, USA. 646 PgsDunn, S., and R. Colohan. 1999. Developing the snow component of a distributed hydrological model: a step-wise approach based on multi-objective analysis. Journal of Hydrology 223:1 -16. Eaton, B., and R.D. Moore. 2007 [draft]. Chapter 4 - Regional hydrology. In R.G. Pike et al. (eds) Compendium of forest hydrology and geomorphology in British Columbia. B.C. Ministry of Forests and Range, Victoria, BC and FORREX, Kamloops, BC. Land Management Handbook. Available at: http;//www. forrex.org/prograiii/wa_ler/coirjp_enduiin.asp_. Ellis, C.R., and J.W. Pomeroy. 2007. Estimating sub-canopy irradiance to melting snow on forested slopes. Hydrological Processes 21:2581-2593. Environment Canada. 2002. British Columbia and Yukon Environmental Monitoring Networks Station Information Centre. hltp.//scitcch.pyr.cc.gc.ca.pyr.ec.gc.ca/climhydio/ Accessed June 18,2009. Environment Canada. 2009. Canadian Climate Normals 1971-2000. http://\\\y w.climatc.wcatherofiice.ec.gc.ca/climate normals Accessed June 18, 2009. Essery, R., and others. 2009. An evaluation of forest snow process simulations. Bulletin of the American Meteorological Society 90:1120-1135. Ferguson, R.I. 1999. Snowmelt runoff models. Progress in Physical Geography 23:205-227. Fierz, C , P. Riber, E. Adams, A. Curran, P. Fohn, M. Lehning, C. Plilss. 2003. Evaluation of snow-surface energy balance models in alpine terrain. Journal of Hydrology 282:76-94. 165 Franz, K.J., T.S. Hogue, and S. Sorooshian. 2008. Operational snow modelling: Addressing the challenges of an energy balance model for the National Weather Service forecasts. Journal of Hydrology 360: 48-66. Gabriel, C. 2004. Structure of three high elevation Engelmann spruce-subalpine fir stands on Mt. Tom, British Columbia. Unpublished report. On file Ministry of Forests and Range, Forest Science Section, Williams Lake, B.C. Gardner, A.S., and M.J. Sharp. 2010. A review of snow and ice albedo and the development of a new physically based broadband albedo parameterization. Journal of Geophysical Research-Earth Surface 115:doi:10.1029/2009JF001444. Garen, D.C., and D. Marks. 2005. Spatially distributed energy balance snowmelt modelling in a mountainous river basin: estimation of meteorological inputs and verification of model results. Journal of Hydrology 315:126-153. Golding, D.L., and R.H. Swanson. 1986. Snow distribution patterns in clearing and adjacent forest. Water Resources Research 22: 1931-1940. Goodison, B.E., P.Y.T. Louie, and D. Yang. 1998. WMO Solid Precipitation Measurement Intercomparison: Final Report. WMO Instruments and Observing Methods Report No. 67, WMO/TD-No.872, 300 pp. Hamilton, A.S., D.G. Hutchinson, and R.D. Moore. 2000. Estimating winter streamflow using a conceptual streamflow model. Journal of Cold Regions Engineering 14: 158175. Hamlet, A., P. Mote, M. Clark, and D. Lettenmaier. 2005. Effects of temperature and precipitation variability on snowpack trends in the western United States. Journal of Climate 18: 4545-4561. Hamlin, L., A. Pietroniro, T. Prowse, R. Soulis, and N. Kouwen. 1998. Application of indexed snowmelt algorithms in a northern wetland regime. Hydrological Processes 12:1641-1657. Hardy, J.P., R.E. Davis, R. Jordan, X. Li, C. Woodcock, W. Ni, and J.C. McKenzie. 1997. Snow ablation modelling at the stand scale in a boreal jack pine forest. Journal of Geophysical Research 102: 29,397-29,405. Henneman, H.E., and H.G. Stefan. 1998. Snow and ice albedo measured with two types of pyranometers. Journal of the American Water Resources Association, 34: 1487-1494. Hock, R. 1999. A distributed temperature-index ice- and snowmelt model including potential direct solar radiation. Journal of Glaciology 45:101-111. 166 Hock, R. 2003. Temperature index melt modelling in mountain areas. Journal of Hydrology 282:104-115. Hock, R. 2005. Glacier melt: a review of processes and their modelling. Progress in Physical Geography 29:362-391. Ichii, K., M.A. White, P. Votava, A. Michaelis, and R.R. Nemani. 2008. Evaluation of snow models in terrestrial biosphere models using ground observation and satellite data: impact on terrestrial ecosystem processes. Hydrological Processes 22: 347-355. DOI:10.1002/hyp.6616. 1PCC. 2007. Climate change 2007: the physical science basis: summary for policy makers. http://www.ipcc.ch Jackson, S.I., and T.D. Prowse. 2009. Spatial variation of snowmelt and sublimation in a high-elevation semi-desert basin of western Canada. Hydrological Processes 23:26112627. Jost, G., M. Weiler, D.R. Gluns, and Y. Alila. 2007. The influence of forest and topography on snow accumulation and melt at the watershed-scale. Journal of Hydrology 347: 101115. Judson, A., and N. Doesken. 2000. Density of freshly fallen snow in the Central Rocky Mountains. Bulletin of the American Meteorological Society 81: 1577-1587. Kane, D.L., R.E. Gieck, and L.D. Hinzman. 1997. Snowmelt modelling at small Alaskan arctic watershed. Journal of Hydrologic Engineering 2:204-210. Kattelman, R. 2000. Snowmelt lysimeters in the evaluation of snowmelt models. Annals of Glaciology 31: 406-410. Kirnbauer, R., G. Bloschl, and D. Gutknecht. 1994: Entering the era of distributed snow models. Nordic Hydrology 25:1-24. Kling, H., J. Fiirst, and H. Nachtnebel. 2006. Seasonal, spatially distributed modelling of accumulation and melting of snow for computing runoff in a long-term, large-basin water balance model. Hydrological Processes 20:2141-2156. Knowles, N., M. Dettinger, and D. Cayan. 2006. Trends in snowfall versus rainfall in the western United States. Journal of Climate 19: 4545-4559. Koivusalo, H., and T. Kokkonen. 2002. Snow processes in a forest clearing and in a coniferous forest. Journal of Hydrology 262: 145-164. 167 Krause, P., D.P. Boyle, and F. Base. 2005. Comparison of different efficiency criteria for hydrological model assessment. Advances in Geoscience 5: 89-97. Kuhn, M. 1987. Micro-meteorological conditions for snow melt. Journal of Glaciology 33:263-72. Kustas, W.P., A. Rango, and R. Uijlenhoet. 1994. A simple energy budget algorithm for the snowmelt runoff model. Water Resources Research 30:1515-1527. Kuusisto, E., 1980. On the values and variability of degree-day melting factors in Finland. Nordic Hydrology 11:235-242. Lang, H., and L. Braun. 1990. On the information content of air temperature in the context of snow melt estimation. In: Hydrology of Mountainous Areas, Proceedings of the Strbske Pleso Symposium 1990. Molnar, L., (Ed.). IAHS Publication no. 190, pp. 347354. Leach, J.A., and R.D. Moore. 2010. Above-stream microclimate and stream surface energy exchanges in a wildfire-disturbed riparian zone. Hydrological Processes DOI: 10.1002/hyp.7639. Lindstrom, G., B. Johansson, M. Persson, M. Gardelin, and S. Bergstrom. 1997. Development and test of the distributed HBV-96 hydrological model. Journal of Hydrology 201:272-288. Link, T.E., and D. Marks. 1999. Point simulation of seasonal snowcover dynamics beneath boreal forest canopies. Journal of Geophysical Research 104: 27,841-27,857. Liston, G. 2004. Representing sub grid snow cover heterogeneities in regional and global models. Journal of Climate 17: 1381-1397. Lopez-Moreno, J.I., and M. Stahli. 2008. Statistical analysis of the snow cover variability in a subalpine watershed: assessing the role of topography and forest interactions. Journal of Hydrology 348: 379-394. Loukas A., and M. Quick. 1999. The effect of climate change on floods in British Columbia. Nordic Hydrology 30: 231-256. Loukas, A., L. Vasiliades, and N.R. Dalezios. 2002. Climatic impacts on the runoff generation processes in British Columbia, Canada. Hydrology and Earth System Sciences 6:211-227. MacDonald, J.P., and J.W. Pomeroy. 2007. Gauge undercatch of two common snowfall gauges in a prairie environment. Proceedings of the 64" Eastern Snow Conference. 1 1 9 126. 168 Male, D.H., and R.J. Granger. 1981: Snow surface and energy exchange. Water Resources Research 17:609-627. Marks, D., A. Winstral, and M. Seyfried. 2002. Simulation of terrain and forest shelter effects on patterns of snow deposition, snowmelt and runoff over a semi-arid mountain catchment. HydrologicalProcesses 16:3605-3626. Martinec, J. 1989. Hour-to-hour snowmelt rates and lysimeter outflow during an entire ablation period. Snow Cover and Glacier Variations (Proceedings of the Baltimore Symposium). IAHS Publications No. 193, pp. 19-28. Martinec, J., and A. Rango. 1986. Parameter values for snowmelt runoff modelling. Journal of Hydrology 84:197-219. McCuen, R.H., Z. Knight, A.G. Cutter. 2006. Evaluation of the Nash-Sutcliffe efficiency index. Journal of Hydrologic Engineering, 11:597-602. McKay, G., and D. Gray. 1981. The distribution of snowcover. In: Handbook of Snow: Principles, Processes, Management and Use. D. Gray and D. Male, (Eds.). The Blackburn Press, New Jersey, USA, pp 153-190. Micovic, Z., and M.C. Quick. 1999. A rainfall and snowmelt runoff modelling approach to flow estimation at ungauged sites in British Columbia. Journal of Hydrology 226:101120. Moore, R., and I. Owens. 1984. Controls on advective snowmelt in a maritime alpine basin. Journal of Climatology and Applied Meteorology 23:135-142. Moore, R.D., and S.M. Wondzell. 2005. Physical hydrology and the effects of forest harvesting in the Pacific Northwest: A review. Journal of the American Water Resources Association 41:763-784. Morrison, J., M.C. Quick, and M.G.G. Foreman. 2002. Climate change in the Fraser River watershed: flow and temperature projections. Journal of Hydrology 263:230-244. Mote, P. 2006. Climate-driven variability and trends in mountain snowpack in western North America. Journal of Climate 19: 6209-6220. Mote, P., A. Hamlet, M. Clark, and D. Lettenmaier. 2005. Declining mountain snowpack in western North America. Bulletin of the American Meteorological Society 86: 39-49. Munro, D.S. 1990. Comparison of melt energy computations and ablatometer measurements on melting ice and snow. Arctic and Alpine Research 22:153-162. 169 Murray, C D . , and J.M. Buttle. 2003. Impacts of clearcut harvesting on snow accumulation and melt in a northern hardwood forest. Journal of Hydrology 271:197-212. Nash, J.E., and J.V. Sutcliffe. 1970. River flow forecasting through conceptual models. Part I: A discussion of principles. Journal of Hydrology 10: 282-291. Ohmura, A. 2001. Physical basis for the temperature-based melt index method. Journal of Applied Meteorology 40:753-761. Oke, T. R. 1987. Boundary Layer Climates, 2nd ed., Methuen, London. Olyphant, G. 1986a. The components of incoming radiation within a mid-latitude alpine watershed during the snowmelt season. Arctic and Alpine Research 18:163-169. Olyphant, G. 19866. Longwave radiation in mountainous areas and its influence on the energy balance of alpine snowfields. Water Resources Research 22:62-66. Onset Computer Corporation. 2010. TidBit v2 Data Logger product page. htlpr/'/wwvv.on^etcomp.com'products/data-loggers/iilbi-OOl Accessed April 20, 2010. Pellicciotti, F., B. Brock, U. Strasser, P. Burlando, M. Funk, and J. Corripio. 2005. An enhanced temperature-index glacier melt model including the shortwave radiation balance: development and testing for Haut Glacier d'Arolla, Switzerland. Journal of Glaciology 51:573-587. Pellicciotti, F., J. Helbing, A. Rivera, V. Favier, J. Corripio, J. Araos, J. Sicart, and M. Carenzo. 2008. A study of the energy balance and melt regime on Juncal Norte Glacier, semi-arid Andes of central Chile, using melt models of different complexity. Hydrological Processes 22:3980-3997. Pliiss, C , and A. Ohmura. 1997. Longwave radiation on snow-covered mountainous surfaces. Journal oj Applied Meteorology 36:818-824. Pomeroy, J., A. Rowlands, J. Hardy, T. Link, D. Marks, R. Essery, J.E. Sicart, and C. Ellis. 2008. Spatial variability of shortwave irradiance for snowmelt in forests. Journal of Hydrometeorology 9:1482-1490. Pomeroy, J.W., D. Marks, T. Link, C. Ellis, J. Hardy, A. Rowlands, and R. Granger. 2009. The impact of coniferous forest temperature on incoming longwave radiation to melting snow. Hydrological Processes 23:2513-2525. Pomeroy, J.W., D.S. Bewley, R.L.H. Essery, N.R. Hedstrom, T. Link, R.J. Granger, J.E. Sicart, CR. Ellis, and J.R. Janowicz. 2006. Shrub tundra snowmelt. Hydrological Processes 20:923-941. 170 Prowse, T., and I. Owens. 1982. Energy balance over melting snow, Craigieburn Range, New Zealand. J Hydrol. N. Z. 21:133-147. Quick, M., and A. Pipes. 1977. U.B.C Watershed Model. Hydrological Sciences Bulletin 22:153-161. Rango, A., and J. Martinec. 1995. Revisiting the degree-day method for snowmelt computations. Water Resources Bulletin 31:657-669. Rex, J., and S. Dube. 2006. Predicting the risk of wet ground areas in the Vanderhoof Forest District: Project description and progress report. BC Journal of Ecosystems and Management 7:57-71. River Forecast Centre. 2009. Automated Snow Pillow Data Archive. http://www.cnv.gov.bc.ca/rfc/data/asp/archive.htm Accessed June 6, 2009. Rood, S., G. Samuelson, J. Weber, and K. Wywrot. 2005. Twentieth-century decline in streamflows from the hydrographic apex of North America. Journal of Hydrology 306:215-233. Rood, S., J. Pan, K. Gill, C. Franks, G. Samuelson, and A. Shepherd. 2008. Declining summer flows of Rocky Mountain rivers: changing seasonal hydrology and probable impacts on floodplain forests. Journal of Hydrology 349: 397-410. Rutter, N., and others. 2009. Evaluation of forest snow processes models (SnowMIP2). Journal of Geophysical Research-Atmospheres 114:D06111. Schuler, T., U.H. Fischer, R. Sterr, R. Hock, and G.H. Gudmundsson. 2002. Comparison of modelled water input and measured discharge prior to a release event: Unteraargletscher, Bernese Alps, Switzerland. Nordic Hydrology 33: 27-46. Seibert, J. 2001. On the need for benchmarks in hydrological modelling. Hydrological Processes 15: 1063-1064. Semadeni-Davies, A. 1997. Monthly snowmelt modelling for large-scale climate change studies using the degree day approach. Ecological Modelling 101:303-323. Sicart, J.E., R. Hock, and D. Six. 2008. Glacier melt, air temperature, and energy balance in different climates: the Bolivian tropics, the French Alps, and northern Sweden. Journal of Geophysical Research 113:D24113. Simpson, J., M. Dettinger, F. Gehrke, T. Mclntire, and G. Hufford. 2004. Hydrologic scales, cloud variability, remote sensing, and models: implications for forecasting snowmelt and streamflow. Weather and Forecasting 19:251-276. 171 Singh, P. R., T. Y. Gan, and A. K. Gobena. 2009. Evaluating a hierarchy of snowmelt models at a watershed in the Canadian Prairies. Journal of Geophysical Research 114:doi:10.1029/2008JD010597. Singh, P., N. Kumar, and M. Arora. 2000. Degree-day factors for snow and ice for Dokriani Glacier, Garhwal Himalayas. Journal of Hydrology 235:1-11. Sivapalan, M. 2003. Process complexity at hillslope scale, process simplicity at the watershed scale: is there a connection? Hydrological Processes 17:1037-1041. Sloan, W., C. Kilsby, and R. Lunn. 2004. Incorporating topographic variability into a simple regional snowmelt model. Hydrological Processes 18:3371-3390. Sorman, A.A., A. Sensoy, A.E. Tekeli, A.U. Sorman, and Z. Akyurek. 2009. Modelling and forecasting snowmelt runoff process using the HB V model in the eastern part of Turkey. Hydrological Processes 7:1031 -1040. Spittlehouse, D.L., R.S. Adams, and R.D. Winkler. 2004. Forest, edge and opening microclimate at Sicamous Creek, B.C. Ministry of Forests Research Branch, Victoria, BC. Res. Rep. 24. http://www.for.gov.bc.ca/hfd/pubs/Docs/Rr/Rr24.htin Stathers, R.J., T. Newsome, M.J. Waterhouse, and C. Sutherland. 2001. Microclimate studies on a group selection silvicultural system in a high-elevation ESSFwc3 forest in the Cariboo Forest Region. B.C. Ministry of Forests,Victoria, B.C. Working Paper 58. Steppuhn, H. 1981. Snow and agriculture. In Handbook of Snow: Principles, Processes, Management and Use. D. Gray and D. Male, Eds. The Blackburn Press, New Jersey, USA. pp 60-125. Stewart, I., D. Cayan, and M. Dettinger. 2005. Changes toward earlier streamflow timing across western North America. Journal of Climate 18: 1136-1155. Stewart, I.T., D.R. Cayan, and M.D. Dettinger. 2004. Changes in snowmelt runoff timing in western North America under a 'Business as Usual' climate change scenario. Climatic Change 62:217-232. Storck, P., T. Kern, and S. Bolton. 1999. Measurement of differences in snow accumulation, melt, and micrometeorology due to forest harvesting. Northwest Science 73:87-101. Suzuki, K , T. Ohta, A. Kojima, and T. Hashimoto. 1999. Variations in snowmelt energy and energy balance characteristics with larch forest density on Mt Iwate, Japan: observations and energy balance analysis. Hydrological Processes 13:2675-2688. 172 Takala, M., J. PuUiainen, S.J. Metsamaki, and J.T. Koskinen. 2009. Detection of snowmelt using spaceborne microwave radiometer data in Eurasia from 1970 to 2007. IEEE Transactions on Geoscience and Remote Sensing 47: 2996-3007. Teti, P. 2003. Relations between peak snow accumulation and canopy density. The Forestry Chronicle 79: 307-312. Teti, P. 2008. Effects of overstory mortality on snow accumulation and ablation Mountain Pine Beetle Working Paper 2008-13. MPBI Project #8.39. B.C. Minstry of Forests and Range, Williams Lake, B.C., Canada. 37 p. Tseng, P.H., T.H. Illangasekare, and M.F. Meier. 1994. Modelling of snow melting and uniform wetting front migration in a layered subfreezing snowpack. Water Resources Research 30:2363-2376. van de Vosse, H. et al. 2008. Mountain pine beetle infestation: hydrological impacts. http://www.llbc.leg.bc.ca/public/pubdocs/bcdocs/445817/finishdownloaddocument.pdf. Verbunta, M., J. Gurtz, K. Jasper, H. Lang, P. Warmerdam, and M. Zappa. 2003. The hydrological role of snow and glaciers in alpine river basins and their distributed modelling. Journal of Hydrology 282:36-55. Vincent, L.A., and E. Mekis. 2006. Changes in daily and extreme temperature and precipitation indices for Canada over the twentieth century. Atmosphere-Ocean 44:177193. Walter, M.T., E.S. Brooks, D.K. McCool, L.G. King, M. Molnau, and J. Boll. 2005. Processbased snowmelt modelling: does it require more input data than temperature-index modelling? Journal of Hydrology 300:65-75. Wang, L., C. Derksen, and R. Brown. 2008. Detection of pan-Arctic terrestrial snowmelt from QuikSCAT, 2000-2005. Remote Sensing of Environment 112: 3794-3805. Waterhouse, M. 2002. Group selection silvicultural systems for high elevation forests (ESSFwc3) to maintain caribou habitat in the Cariboo Forest Region. EP1104.02. Mt. Tom Adaptive Management Trial Working Plan. B.C. Min. For. Williams Lake, B.C. Watson, F.G.R., W.B. Newman, J.C. Coughlan, and R.A. Garrott. 2006. Testing a distributed snowpack simulation model against spatial observations. Journal of Hydrology 328:453-466. Whitaker, A., Y. Alila, J. Beckers, and D. Toews. 2003. Application of the Distributed Hydrology Soil Vegetation Model to Redfish Creek, British Columbia: model evaluation using internal catchment data. Hydrological Processes 17:199-224. 173 Whitfield, P., and A. Cannon. 2000. Recent variations in climate and hydrology in Canada. Canadian Water Resources Journal 25:19-65. Wilmott, C.J. 1982. Some comments on the evaluation of model performance. Bulletin American Meteorological. Society 63: 1309-1313. Winkler, R.D., D.L. Spittlehouse, and D.L. Golding. 2005. Measured differences in snow accumulation and melt among clearcut, juvenile, and mature forests in southern British Columbia. Hydrological Processes 19:51-62. WMO, 1986. Intercomparison of models for snowmelt runoff. Operational Hydrology Report 23 (WMO No. 646). Woo, M., and R. Thorne. 2006. Snowmelt contribution to discharge from a large mountainous catchment in subarctic Canada. Hydrological Processes 20:2129-2139. Young, G. 1985. Overview. In: Techniques for Prediction of Runoff from Glacierized Areas. Young, G., (Ed.). IAHS Publication no. 149, pp. 3-23. Zappa, M., F. Pos, U. Strasser, P. Warmerdam, and J. Gurtz. 2003. Seasonal water balance of an alpine catchment as evaluated by different methods for spatially distributed snowmelt modelling. Nordic Hydrology 34: 179-202. Zeinivand, H., and F. De Smedt. 2009. Hydrological modeling of snow accumulation and melting on river basin scale. Water Resources Management 23: 2271-2287. Zimmerman, A. 1972. Air temperature observations and forecasts-their relationship to the prediction of spring snowmelt in the Eagel River Basin, Colorado. In: Proceedings of the 40th Annual Meeting of the Western Snow Conference, 1972, Phoenix, Arizona. Colorado State University, Colorado, USA, pp 30-36. Zuzel J., and L. Cox. 1975. Relative importance of meteorological variables in snowmelt. Water Resources Research 11:174-176. 174 Appendix 1 Barkerville Automatic Snow Pillow At the Barkerville Automatic Snow Pillow (ASP) SWE is measured over a 3 m diameter bladder containing antifreeze solution As snow accumulates on the pillow, the weight of the snow pushes an equal weight of the antifreeze solution from the pillow up a standpipe in the instrument house Accumulated precipitation is measuied in an unshielded PVC standpipe with a 1500 mm capacity including antifreeze and mineial oil to prevent evaporation The measurements are provided by a 2 psi pressure transducer mounted in the base A battery powered pump circulates the liquid in the gauge once an hour at 20 minutes past the hour for 20 seconds to prevent fieezc-up The readings are taken on the hour - the offset is to ensure that ruibulence from the pump does not influence the pressure transducer Snow depth is measured using a Campbell Scientific SR50 mounted on a 6 m tower Generally, the error margin for the automatic snow pillows is ±10% compared to the manual measurements (Scott Jackson, River Forecast Centre pers comm, 2009) The precipitation gauges aie accurate given the challenges associated with these measuiements, and are usually within 2-5 mm of the total accumulated precipitation value The snow depth measurements arc moie susceptible to enor - the sonic pulse used to make the measuiement is affected by wind, heavy precipitation and solai radiation These measurements are usually within a few millimeters of the manually measured snow depth directly beneath them 175 Appendix 2: Precipitation at Mt Tom. I suspected that delays in precipitation recording occurred due to slush build up and possible icing in the precipitation gauge. As a result, the timing of precipitation measurement at the gauge likely differed from the occurrence of precipitation. To get an idea of the accuracy of the tipping bucket gauge measurements I compared the precipitation measured by the Mt Tom met station gauge to the accumulation of snow as measured by the SR50.1 assumed a snow density of 10:1 for conversion of snow depth to water equivalent. There are obvious issues with using the density assumption for the SR50, for example fresh snow density in the Rocky Mountains can range from 10 to 257 kg m"3 (Judson and Doesken 2000). An hourly increase in snow depth of 1.0 cm was the minimum to be counted as a precipitation event. Sensitivity of the SR50 is ±1.0 cm. The gauge did not agree well with the SR50 measurements (Figure 6.1). Both the amount and timing of the precipitation were offset. I also compared the SR50 measurements from the Barkerville automatic snow pillow and found these were similar to that measured at Mt Tom (Figure 6.4). — BkSR50 --- MtTomSRSO MtTom TB - - EC precip o o E E a. Apr M,-w Figure 6 1 Cumulative sum of piecipitation measured at the BarkciviUc Environment Canada station and with a lipping bucket tain gauge at Mt Tom, and obtained by converting snow depth (mcasuied with SR50) accumulations to water equivalent assuming a 10.1 density nUio The problems with the tipping bucket agree with findings by MacDonald and Pomeroy (2007) and are largely inherent in the design of the precipitation adapter. Unfortunately I could not correct for the wind related gauge undercatch because of the event timing issue. I also compared cumulative SWE and precipitation data at Barkerville ASP with the precipitation data as measured at the Barkerville Environment Canada station (Figure 3.3b). Cumulative precipitation measured at ASP and the Environment Canada station were in fairly good agreement 176 indicating that precipitation measured at the Environment Canada station should be fairly representative of the Mt Tom area. I concluded that the precipitation data as measured at Mt Tom should not be used for this study, and I opted to use the precipitation data from the BarkerviUe Environment Canada station instead. Since the BarkerviUe precipitation gauge was not shielded, I expect there to be undercatch issues. The data agree with this as SWE was greater than precipitation during the accumulation period. As there were no wind speed measurements at the ASP site it would be difficult to correct for this. Finally, since the objective is to use these models for operational purposes it would be best to use data that are commonly available. 177 Appendix 3 Temperature-index Albedo Model at Mt Tom The albedo model developed by Brock et al (2000) used different equations for conditions of deep (>5 cm w e ) and shallow snowpack (<0 5 cm w e ) For a deep snowpack 0 ^ = 0 7 1 3 - 0 1121og 1 0 r (6 Where adi is the albedo for a deep snowpack (>5 cm w e ) and T is the accumulated daily maximum air temperature >0 °C since the last snowfall The values 0 713 and 0 112 are the calibrated empirical coefficients For a shallow snowpack a S i = « „ + 0 442e(-00W) (6 Where ais is the albedo for a shallow snowpack (<0 5 cm w c ) and a„ is the albedo of the underlying surface The lowest albedo measured at Mt Tom, 0 17, was used for the giound siuface albedo To deal with the discontinuity that arises when changing from the deep to the shallow snow model Block et al (2000) recommended calculating the albedo as a weighted function of decieasing SWE ( V _£ l-edt -d\ «* + e * / \ ) 0 f d Where d is the depth of SWE and d* is a scaling length for d, which Block et al (2000) found to be 2 4 cm w e for their glacier site Precipitation was required to ran the a model and it was necessary to know piccipitation foim The occuncnce of snow fall reset the a to its maximum value Since a reliable mcasuic of piccipitation was not available on-site at Mt Tom, and to maintain consistency among all the model runs, the daily piccipitation amount mcasuied at the Envuonmcnt Canada station at BaikeiviUe was used The Envuonmcnt Canada station was located appioximately 200 m lower in elevation and 31 km southwest of the study sites (Ftgutc 2 1) Piccipitation that fell as lain at the Envuonmcnt Canada station may have fallen as snow at the highei elevation Mt Tom study sites Therefore, the air tempciatuie measured at each site was used with a cut off tempeiature of 0°C to sepaiatc lain and snow Daily SWE was lequircd to dcteimine the scaling factoi in equation 6-3 and to switch from the deep to the shallow snow models As SWE was not measuicd daily on site, a simple linear function was used to interpolate the daily SWE values between snow course sample dates The SWE measurements at each of the four snow couise sites were used to run the model The fust pioblem encountcicd was in deteimimng 'snowfall events', which were used to leset the albedo to its maximum value Only daily snow depth data, measuied with the SR50, which has a 1 cm lcsolution, were available Since using the SR50 data to determine snowfall events lequued a subjective analysis, the precipitation lecord at the BaikeiviUe Environment Canada station was used, a snowfall event occuned on any day when piccipitation fell at an on-site daily average an tempciature less than 0°C The numbei of snowfall events vaued depending on the method used with fcwei events being pinpointed using the Envuonmcnt Canada measured piecipitation compaied to the analysis ot the SR50 data collected on site (Figuie 6 2) 178 TV T':> «x? »« Apt 11 .y> Apr 21 NW/01 May 1 I May 21 May 31 Figure 6 4 Compaiison of the occuiience of snowfall (denoted by 1) when the piccipitation data fiom the Barkervillc ASP was used oi when the piccipitation data fiom the Envnonmcnt Canada station was used There were no diffeienees in the number oi timing of snowtall events between the Mt Tom snow coutsc sites and the Baikerville ASP The modelled albedo values vancd slightly between the two sites, paiticulaily on 2 Apul (Figuie 6 5) On this date the difference was due to the maximum daily temperatuic at each site At Mt Tom the maximum daily an lempciahire was above 0°C, thcicioie the albedo decayed ovei that day However, at the BarkeiviUe ASP site the maximum daily air tcmpcratuie was below 0°C, thcrefoie the albedo lemained at the maximum value The small diffeienees later in the simulation pcnod aie due to variation in maximum daily temperatures, and at the end of the pcnod diffeienees aie due to the date on which SWE dropped below the cntical value of 5 cm w c 181 o — — — — — south open south forest north open north forest BkASP 1 1 1 1 I I 1 Apr 01 April Apr 21 May 01 May 11 May 21 May 31 Figure 6.5 Modelled daily albedo for the Mt Tom snow course sites and the Barkerville ASP. It was assumed that the TidBit sites experienced the same precipitation regime as the snow course sites due to their proximal location. The daily SWE had to be approximated for the tidbit sites as only an initial SWE measurement was taken. The SWE is important for determining the contribution of the deep and shallow albedo models to the overall albedo. The date on which the threshold SWE of 5 cm w.e. is crossed marks a strong change in the albedo pattern for the site. Since only an initial SWE value was measured for the tidbit sites the daily SWE was approximated using a simple linear decay function, the same as that applied to the snow course sites. The appropriateness of the albedo value at the tidbit sites is thus dependent on this assumption. The varying end dates of the melt at the tidbit sites resulted in variation in the end albedo values (Figure 6.6). 182 oo O o CD O o o oo o — --— --- CM O Apr 01 west open west forest east open eastfoiest i April 1 Apr 21 1 May 01 Figure 6 6 Modelled daily albedo for the Mt Tom TidBit sites 183 1 May 11 1 May 21 r May. Appendix 4 Hock and Pelhcciotti Parameter Comparison Modelled variation in melt between sites on different aspects should be due to a factor other than air temperature as the same air temperature data were used in modelling melt for all snow course sites This variation is explained by the radiation factor in these models Therefore, where there is significant difference in melt between sites on opposite aspects, radiation factors greater than zero aie expected to be optimized for the model Unfortunately only two measurement periods overlap for snow course measurements on opposing aspects 9-17 May and 17 - 23 May (Figure 3 10a) For the open sites, the first overlapping measurement period melt rates were nearly the same between north and south aspect sites For the second period melt at the south aspect site was much larger than that observed at the north aspect (Figure 3 10a) The radiation factor was greater than zero for all open site calibrations using the Hock model except for a single case where measured radiation was used and the NO site was withheld (Table 3 4) However, the radiation factors were zero for all open site calibration runs using the Pelhcciotti model except where the south aspect open site was withheld (Table 3 5) Since radiation is a factor in determining the spatial variability in melt at the open sites, short wave radiation is expected to have been low (cloudy conditions) for the period pieceding the 17 May measurement to explain the similar melt on both aspects For the second period radiation is expected to have been quite high to explain the high melt on the south aspect compared to the north Examining the daily and penod averaged values for ladiation and temperatuie do not clearly indicate this (Figure 3 10, Tabic 6 1) Air temperatuie was slightly lower and shortwave radiation was gieater than in the second period (Table 6 1) This is the opposite of expected The three days of unusually high an temperatures (11-16°C) that immediately picceded the snow course measurement on 17 May might have affected the expected relationship by causing high melt mtes compared to the pievious five days which experienced daily average temperatures less than 5°C However, the north aspect site had more SWE than the south aspect site for the same measurement date as well as colder average snowpack temperature (SO = -0 9 °C on 9 March 2008 and NO = -2 5 °C on 18 Maich 2008) The difference in melt may be reflecting the difference in ripeness of each snowpack, and thus the capability foi melt in response to cneigy inputs, on those particular dates Table 6 1 Mean an tcmpeialuic (Ta ) and incoming shoilwave ladiation (K) as measuicd at the Mt Tom mcleoiological station Values wcic avciaged acioss the pcnod between snow couisc measurements End dale of each measuiemcnt penod shown MEAN K Date fC) (Wm2) 5/2/2008 23 167 5/9/2008 31 191 5/17/2008 67 198 5/23/2008 75 182 5/30/2008 87 217 6/5/2008 89 217 184 The MF F values were nearly identical between the Hock and Pelhcciotti models (Table 3 4, Table 3 5) For both models the lowest RMSE value occuired when the north aspect open site was withheld For the forest sites, daily average melt at the south aspect site was slightly higher than at the north aspect for the first period and slightly lower during the second period (Figure 3 10a) No radiation factor was calibrated for any foiested site using either the Hock or Pelhcciotti model (Table 3 4, Table 3 51 (Table 3 4, Table 3 5) Either the diffeience in melt was too small to be significant oi the reversal in aspect experiencing greater melt cancelled out the effect Of consequence for both forest and open sites, there were only two data points to compare melt among sites making it difficult to attribute an effect or even significant difference in melt The melt at the north aspect foiest site is higher than the noith open and south foiest sites for the period 17 - 23 May The south aspect forest site was at the tail end of the melt during this period, few measurement points still had snow From 9-17 May, four out of the 12 sites completely melted, and from 17-23 May, seven out of the sites points went to 0 SWE (theie was only one point remaining with snow on 23 May) Unfortunately, as snow couise measuiements were only taken weekly it is impossible to know the exact date on which snow melted out for each sampling point Therefore for any point that melted out before 23 May, the daily average was still calculated over the entire period and would give a lowci melt rate than actually occurred at a given point This could be significant if enough sites melted out before the sampling on 23 May, since only one point still had snow this could be the reason for the relatively low melt rate for the south aspect forest site In contrast the north aspect forest site had only had thice 0 measuiements on 23 May, and full snow for the previous two measurement dates 185