Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Hydrometeorology and streamflow response during rain-on-snow events in a coastal mountain region Trubilowicz, Joel William 2016

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

Item Metadata


24-ubc_2016_september_trubilowicz_joel.pdf [ 14.6MB ]
JSON: 24-1.0303146.json
JSON-LD: 24-1.0303146-ld.json
RDF/XML (Pretty): 24-1.0303146-rdf.xml
RDF/JSON: 24-1.0303146-rdf.json
Turtle: 24-1.0303146-turtle.txt
N-Triples: 24-1.0303146-rdf-ntriples.txt
Original Record: 24-1.0303146-source.json
Full Text

Full Text

Hydrometeorology and streamflow response duringrain-on-snow events in a coastal mountain regionbyJoel William TrubilowiczB.Sc. Environmental Engineering, Michigan Technological University, 2004M.A.Sc., Forest Hydrology, The University of British Columbia, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Geography)The University of British Columbia(Vancouver)May 2016c© Joel William Trubilowicz, 2016AbstractRain-on-snow, in which rainfall occurs upon a previously existing snowpack, com-plicates runoff response to rain events. In some situations the snowpack absorbsrainfall; in others, runoff is enhanced by significant snowmelt. Rain-on-snow hasgenerated major floods around the world, particularly in coastal, mountainous re-gions such as southwestern British Columbia, where the rugged topography cancause rapid and varying transitions between rainfall and snowfall within the samewatershed, and warm, subtropical storms known as atmospheric rivers can deliverlarge quantities of precipitation. This thesis sought to further the understandingof rain-on-snow at the regional scale, including its role in runoff response to awide spectrum of rain-on-snow events. Tools were developed and assessed to helpachieve this goal and for use by others.First, the hydrologic utility of output from a regional weather reanalysis modelwas tested. Results showed air temperature and vapour pressure were likely tobe useful, whereas other variables were not accurate enough to be of use. Airtemperature, in particular, showed potential ability for more accurate specificationof temperature gradients for hydrologic forecasting of rain-on-snow runoff.An analysis of rain-on-snow events across five automated snow pillow sitesover 10 years illustrated the importance of understanding the amount of rainfalloccurring at high elevations during rain-on-snow, and the relatively consistent en-hancement of water available for runoff (WAR) by 25-30% due to snowmelt duringlarge rain-on-snow events. For smaller events, a range of antecedent and meteoro-logical factors influenced WAR generation, particularly the antecedent liquid watercontent of the snowpack.A probabilistic method for infilling cloud obscured pixels in optical remotelyiisensed snow cover imagery showed strong performance compared to more stan-dard methods, and illustrated spatial changes in snow cover during the largest floodevent within the analysis period. This method was developed to maximize infor-mation gain from satellite snowcover imagery while minimizing the transfer ofdisinformation.Finally, high elevation rainfall during atmospheric river events was found to bethe dominant predictor of runoff response in six study catchments in the region.Antecedent snowcover provided only minimal increases in the ability to predictrunoff during these events compared to rainfall alone.iiiPrefaceThis thesis was completed under the supervision of Dr. R.D. Moore, who providedguidance on data analysis and editing for all chapters of this work. R.D. Moorealso helped in the conceptualization of methods for downscaling reanalysis outputto surface observations in Chapter 3, and the mass balance estimation in Chapter 4.I am the lead author of all chapters of this thesis, and led the conceptualization,data analysis, and writing of all chapters.A version of Chapter 3, entitled ’Suitability of North American Regional Re-analysis (NARR) output for hydrologic modelling and analysis in mountainousterrain’ is currently available as an early view article for the journal HydrologicalProcesses (DOI:10.1002/hyp.10795). Co-authors are J.M. Shea, of the Interna-tional Centre for Integrated Mountain Development, G. Jost of BC Hydro, andR.D. Moore.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . 41.1.1 High elevation observations . . . . . . . . . . . . . . . . 41.1.2 Snowpack processes during rain-on-snow . . . . . . . . . 61.1.3 Remote sensing of snow . . . . . . . . . . . . . . . . . . 91.1.4 Rain-on-snow at the watershed scale . . . . . . . . . . . . 111.1.5 Atmospheric rivers . . . . . . . . . . . . . . . . . . . . . 131.2 Research outline . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 Regional description and data sources . . . . . . . . . . . . . . . . . 172.1 Regional description . . . . . . . . . . . . . . . . . . . . . . . . 172.2 Data sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19v2.2.1 Point data . . . . . . . . . . . . . . . . . . . . . . . . . . 192.2.2 Spatial data . . . . . . . . . . . . . . . . . . . . . . . . . 232.2.3 Stream discharge . . . . . . . . . . . . . . . . . . . . . . 302.2.4 Atmospheric river events . . . . . . . . . . . . . . . . . . 333 Suitability of North American Regional Reanalysis (NARR) outputfor hydrologic analysis and modelling . . . . . . . . . . . . . . . . . 353.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.2.1 Weather observations . . . . . . . . . . . . . . . . . . . . 373.2.2 NARR products . . . . . . . . . . . . . . . . . . . . . . . 373.2.3 Alternative methods for estimating incoming radiation . . 393.2.4 Performance measures . . . . . . . . . . . . . . . . . . . 393.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.3.1 Example time series . . . . . . . . . . . . . . . . . . . . 403.3.2 Air temperature . . . . . . . . . . . . . . . . . . . . . . . 423.3.3 Vapour pressure . . . . . . . . . . . . . . . . . . . . . . . 453.3.4 Wind speed . . . . . . . . . . . . . . . . . . . . . . . . . 473.3.5 Incoming shortwave radiation . . . . . . . . . . . . . . . 503.3.6 Incoming longwave radiation . . . . . . . . . . . . . . . . 523.4 Case study: use of NARR products to supplement air temperatureobservations for modelling rain-on-snow events in mountain catch-ments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 533.4.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . 533.4.2 Study catchment and modelling system . . . . . . . . . . 543.4.3 Using NARR to estimate vertical air temperature profiles . 553.4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 563.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 603.5.1 Potential for use of NARR output for hydrologic analysis . 603.5.2 Use of NARR in hydrologic modelling . . . . . . . . . . 623.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64vi4 Point-scale analysis of the role of the snowpack in generating wateravailable for runoff during rain-on-snow events . . . . . . . . . . . . 654.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 654.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 684.2.1 Data sources . . . . . . . . . . . . . . . . . . . . . . . . 684.2.2 Discrimination between rain and snow . . . . . . . . . . . 684.2.3 Identification of rain-on-snow events . . . . . . . . . . . . 704.2.4 Calculation of the mass balance at a snow pillow . . . . . 704.2.5 Modelling of the cold content and liquid water content ofthe snow pack . . . . . . . . . . . . . . . . . . . . . . . . 724.2.6 Calculation of the influence of the energy content of rain-fall on melt . . . . . . . . . . . . . . . . . . . . . . . . . 744.2.7 Analysis of rain-on-snow events . . . . . . . . . . . . . . 754.2.8 Recursive partitioning to predict WAR production . . . . . 764.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 774.3.1 Discrimination between rain and snow . . . . . . . . . . . 774.3.2 Point modelling at snow pillow sites . . . . . . . . . . . . 774.3.3 Event scale analysis . . . . . . . . . . . . . . . . . . . . 784.3.4 Snowmelt rate and advective energy exchange . . . . . . . 814.3.5 Regression tree analysis . . . . . . . . . . . . . . . . . . 814.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 854.4.1 WAR generation during large events . . . . . . . . . . . . 854.4.2 Seasonal variation of WAR composition . . . . . . . . . . 874.4.3 Effects of antecedent conditions . . . . . . . . . . . . . . 904.4.4 Controls on melt generation . . . . . . . . . . . . . . . . 904.4.5 Sources of uncertainty . . . . . . . . . . . . . . . . . . . 914.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 935 Probabilistic estimation of snow cover from MODIS imagery . . . . 955.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 955.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 985.2.1 Data sources . . . . . . . . . . . . . . . . . . . . . . . . 985.2.2 Uncertainty model for unobscured pixels . . . . . . . . . 100vii5.2.3 Uncertainty model for cloud-obscured pixels . . . . . . . 1025.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1055.3.1 Probabilistic classification of observed pixels . . . . . . . 1055.3.2 Infilling of cloud-obscured pixels . . . . . . . . . . . . . 1075.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1105.4.1 Regional accuracy and adjustment . . . . . . . . . . . . . 1105.4.2 Uncertainty in observed pixels . . . . . . . . . . . . . . . 1145.4.3 Infilling of cloud-obscured pixels . . . . . . . . . . . . . 1155.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1176 The influence of antecedent snow and high elevation rainfall on runoffgeneration in atmospheric rivers . . . . . . . . . . . . . . . . . . . . 1186.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1186.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1206.2.1 Data sources . . . . . . . . . . . . . . . . . . . . . . . . 1206.2.2 Flood frequency analysis . . . . . . . . . . . . . . . . . . 1226.2.3 River discharge and rain-on-snow events . . . . . . . . . 1226.2.4 Changes in snow covered area (SCA) during atmosphericrivers . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1236.2.5 Regional response to atmospheric rivers . . . . . . . . . . 1246.2.6 Influence of rainfall and initial snow cover on discharge . 1246.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1256.3.1 Flood frequency analysis and atmospheric rivers . . . . . 1256.3.2 High elevation rainfall during atmospheric rivers . . . . . 1266.3.3 Snow covered area during atmospheric rivers . . . . . . . 1296.3.4 Individual event investigation . . . . . . . . . . . . . . . 1296.3.5 Regional response to atmospheric rivers . . . . . . . . . . 1306.3.6 Influence of rainfall and initial snow cover on discharge . 1326.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1346.4.1 Flood frequency analysis . . . . . . . . . . . . . . . . . . 1346.4.2 Floods and high elevation runoff generation . . . . . . . . 1366.4.3 Snow covered area . . . . . . . . . . . . . . . . . . . . . 1376.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139viii7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1407.1 Summary of key findings . . . . . . . . . . . . . . . . . . . . . . 1407.2 Significance and implications of the findings . . . . . . . . . . . . 1437.3 Recommendations for future research . . . . . . . . . . . . . . . 144References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148A Globcover 2009 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164B Infilled time series at ASP Sites . . . . . . . . . . . . . . . . . . . . . 166ixList of TablesTable 2.1 Hourly point data . . . . . . . . . . . . . . . . . . . . . . . . 21Table 2.2 Basin Elevation and Area . . . . . . . . . . . . . . . . . . . . 24Table 2.3 Study basin land cover . . . . . . . . . . . . . . . . . . . . . . 24Table 2.4 Biogeoclimatic zones of study basins . . . . . . . . . . . . . . 25Table 2.5 Ground based photo georeferencing dates . . . . . . . . . . . . 28Table 2.6 Viewpoints used for photo georeferencing . . . . . . . . . . . 28Table 2.7 Stream gauge locations . . . . . . . . . . . . . . . . . . . . . 32Table 2.8 Atmospheric river events . . . . . . . . . . . . . . . . . . . . 34Table 3.1 Air temperature performance measures . . . . . . . . . . . . . 44Table 3.2 Vapour pressure performance measures . . . . . . . . . . . . . 45Table 3.3 Wind speed performance measures . . . . . . . . . . . . . . . 48Table 3.4 K↓ performance measures . . . . . . . . . . . . . . . . . . . . 51Table 3.5 Potential K↓ performance measures . . . . . . . . . . . . . . . 52Table 3.6 L↓ performance measures . . . . . . . . . . . . . . . . . . . . 53Table 4.1 Summary Variables . . . . . . . . . . . . . . . . . . . . . . . 76Table 4.2 Error statistics for SWE simulation . . . . . . . . . . . . . . . 80Table 4.3 Rain-on-snow event summary . . . . . . . . . . . . . . . . . . 81Table 5.1 Terra naive Bayes probabilities . . . . . . . . . . . . . . . . . 106Table 5.2 Aqua naive Bayes probabilities . . . . . . . . . . . . . . . . . 106Table 5.3 Combined naive Bayes probabilities . . . . . . . . . . . . . . 106Table 5.4 Cloud filter test days . . . . . . . . . . . . . . . . . . . . . . . 107Table 5.5 Test Cloud Masks . . . . . . . . . . . . . . . . . . . . . . . . 110xTable 6.1 ASP and watershed association . . . . . . . . . . . . . . . . . 121Table 6.2 Categorization of P50 results . . . . . . . . . . . . . . . . . . . 123Table 6.3 Percent of annual floods associated with ARs . . . . . . . . . . 127Table 6.4 Correlation in WAR during atmospheric rivers . . . . . . . . . 132Table 6.5 Correlation in discharge during atmospheric rivers . . . . . . . 132Table 6.6 Multiple regression results for discharge prediction . . . . . . . 133Table A.1 Globcover 2009 categories . . . . . . . . . . . . . . . . . . . 165xiList of FiguresFigure 2.1 South coast region . . . . . . . . . . . . . . . . . . . . . . . 19Figure 2.2 Example atmospheric river . . . . . . . . . . . . . . . . . . . 20Figure 2.3 Wolf River automated snow pillow (ASP) site . . . . . . . . . 22Figure 2.4 Whistler peak viewshed . . . . . . . . . . . . . . . . . . . . 29Figure 2.5 DEM projection . . . . . . . . . . . . . . . . . . . . . . . . . 30Figure 2.6 Cheakamus snow cover 6 August, 2012 . . . . . . . . . . . . 31Figure 2.7 Annual hydrographs for study watersheds . . . . . . . . . . . 33Figure 3.1 Timeseries at Bridge . . . . . . . . . . . . . . . . . . . . . . 41Figure 3.2 Timeseries at Place . . . . . . . . . . . . . . . . . . . . . . . 42Figure 3.3 Air temperature correlation . . . . . . . . . . . . . . . . . . . 43Figure 3.4 Air temperature seasonality . . . . . . . . . . . . . . . . . . 44Figure 3.5 Vapour pressure correlation . . . . . . . . . . . . . . . . . . 46Figure 3.6 Vapour pressure seasonality . . . . . . . . . . . . . . . . . . 47Figure 3.7 Wind speed correlation . . . . . . . . . . . . . . . . . . . . . 49Figure 3.8 Wind speed seasonality . . . . . . . . . . . . . . . . . . . . . 50Figure 3.9 K↓ correlation . . . . . . . . . . . . . . . . . . . . . . . . . 51Figure 3.10 L↓ correlation . . . . . . . . . . . . . . . . . . . . . . . . . . 52Figure 3.11 NARR high elevation temperatures . . . . . . . . . . . . . . 58Figure 3.12 Modelled and NARR temperature profiles . . . . . . . . . . . 59Figure 4.1 Rain-snow model fit . . . . . . . . . . . . . . . . . . . . . . 78Figure 4.2 Observed and simulated SWE . . . . . . . . . . . . . . . . . 79Figure 4.3 Rain-snow event summary . . . . . . . . . . . . . . . . . . . 82xiiFigure 4.4 WAR generation . . . . . . . . . . . . . . . . . . . . . . . . 83Figure 4.5 Site specific WAR generation . . . . . . . . . . . . . . . . . 84Figure 4.6 Distribution of melt/WAR ratio . . . . . . . . . . . . . . . . . 85Figure 4.7 Example event timeseries . . . . . . . . . . . . . . . . . . . . 86Figure 4.8 Advective heat from rainfall . . . . . . . . . . . . . . . . . . 87Figure 4.9 Large event regression tree . . . . . . . . . . . . . . . . . . . 88Figure 4.10 Small event regression tree . . . . . . . . . . . . . . . . . . . 89Figure 5.1 Probabilistic model flow chart . . . . . . . . . . . . . . . . . 99Figure 5.2 Probabilistic test days . . . . . . . . . . . . . . . . . . . . . . 108Figure 5.3 Test cloud masks . . . . . . . . . . . . . . . . . . . . . . . . 109Figure 5.4 Probabilistic infilling results . . . . . . . . . . . . . . . . . . 111Figure 5.5 Model accuracy comparison . . . . . . . . . . . . . . . . . . 112Figure 5.6 High confidence model predictions . . . . . . . . . . . . . . . 113Figure 6.1 Flood frequency analysis . . . . . . . . . . . . . . . . . . . . 126Figure 6.2 Precipitation during atmospheric rivers . . . . . . . . . . . . 127Figure 6.3 High elevation rainfall during atmospheric rivers . . . . . . . 128Figure 6.4 SCA at the onset of ARs . . . . . . . . . . . . . . . . . . . . 130Figure 6.5 October 2003 change in P50 by watershed . . . . . . . . . . . 131Figure B.1 Chilliwack River Ta correlation . . . . . . . . . . . . . . . . . 167Figure B.2 Chilliwack River Ta correlation by month . . . . . . . . . . . 168Figure B.3 Chilliwack River Ta filled time series . . . . . . . . . . . . . 169Figure B.4 Chilliwack River P correlation . . . . . . . . . . . . . . . . . 169Figure B.5 Chilliwack River P filled time series . . . . . . . . . . . . . . 170Figure B.6 Disappointment Lake Ta correlation . . . . . . . . . . . . . . 171Figure B.7 Disappointment Lake Ta correlation by month . . . . . . . . . 172Figure B.8 Disappointment Lake Ta filled time series . . . . . . . . . . . 173Figure B.9 Disappointment Lake P correlation . . . . . . . . . . . . . . 173Figure B.10 Disappointment Lake P filled time series . . . . . . . . . . . 174Figure B.11 Jump Creek Ta correlation . . . . . . . . . . . . . . . . . . . 175Figure B.12 Jump Creek Ta correlation by month . . . . . . . . . . . . . . 176Figure B.13 Jump Creek Ta filled time series . . . . . . . . . . . . . . . . 177xiiiFigure B.14 Jump Creek P correlation . . . . . . . . . . . . . . . . . . . . 177Figure B.15 Jump Creek P filled time series . . . . . . . . . . . . . . . . . 178Figure B.16 Upper Squamish River Ta correlation . . . . . . . . . . . . . 179Figure B.17 Upper Squamish River Ta correlation by month . . . . . . . . 180Figure B.18 Upper Squamish River Ta filled time series . . . . . . . . . . 181Figure B.19 Upper Squamish River P correlation . . . . . . . . . . . . . . 181Figure B.20 Upper Squamish River P filled time series . . . . . . . . . . . 182Figure B.21 Wolf River Ta correlation . . . . . . . . . . . . . . . . . . . . 183Figure B.22 Wolf River Ta correlation by month . . . . . . . . . . . . . . 184Figure B.23 Wolf River Ta filled time series . . . . . . . . . . . . . . . . . 185Figure B.24 Wolf River P correlation . . . . . . . . . . . . . . . . . . . . 185Figure B.25 Wolf River P filled time series . . . . . . . . . . . . . . . . . 186xivGlossaryAMSU Advanced Microwave Sounds UnitANOVA analysis of varianceASP automated snow pillowAWS automated weather stationBC British ColumbiaCMS Cheakamus synthetic climate stationDSLR digital single lens reflexDCP data collection platformDEM digital elevation modelEm Nash-Sutcliffe model efficiencyENSO El Nin˜o Southern OscillationFTP file transfer protocolGIS geographic information systemsGOES Geostationary Operational Environmental SatelliteLiDAR light detection and rangingMBE mean bias errorxvMODIS MODerate resolution Imaging SpectroradiometerNARR North American Regional ReanalysisNCEP National Center for Environmental PredictionNDSI Normalized Difference Snow IndexNDVI Normalized Difference Vegetation IndexNSIDC National Snow and Ice Data CenterNOAA National Oceanic and Atmospheric AdministrationNOHRSC National Operational Hydrologic Remote Sensing CenterOLI Operational Land ImagerRMSE root mean squared errorSAMM Sensitive Area Mapping ModuleSCA snow covered areaSSMI Special Sensor Microwave ImagerSWE snow water equivalentTRIM Terrain Resource Information ManagementUBCWM UBC Watershed ModelWAR water available for runoffWSC Water Survey of CanadaxviAcknowledgmentsI have many people to thank for their assistance in the completion of this thesis.First, this thesis would not have been possible without the ideas, assistance, andencouragement of Dr. R. Dan Moore, my Ph.D. supervisor. Along with Dan,I would like to thank my committee members, Dr. Sean Fleming and Dr. IanMcKendry who provided input on the scope and direction of this thesis and myentire education process. Sean also organized the initial agreement and fundingfrom BC Hydro and got me on board as a PhD student. Dr. Nicholas Coopsprovided input in the evolutionary stages and during my comprehensive exams.Georg Jost, Stephanie Smith, Frank Weber, Adam Gobena, Scott Weston, and TimAshman of BC Hydro all provided valuable help, computer use, access to data,and excellent perspective on the scope and goals of this work. In particular, Georgprovided model code and ongoing support in the access of data and model outputfrom BC Hydro. David Campbell of the BC River Forecast Centre, and Karl Jonesand Jessica Byers of the BC Ministry of Environment, provided access to data andinformation and insights for the analysis of the snow pillow data.Joseph Shea of the International Centre for Integrated Mountain Developmentprovided high elevation meteorological observations, and he, along with BrianMenounous of UNBC, provided guidance for access of MODIS snowcover data.Faron Anslow of the Pacific Climate Impacts Consortium provided guidance foraccess of North American Regional Reanalysis model output. Dave Hutchinsonand Ruping Mo of Environment Canada provided perspectives and records of at-mospheric rivers and access to hydrometric data from the Water Survey of Canada.Bill Floyd of the BC Ministry of Forests Lands and Natural Resource Operationsprovided perspectives on understanding rain-on-snow in the field. Justin Knudsenxviiand Spencer Pasieka provided field assistance for photography used in the groundbased photo reprojection, and Vincent Kujala of UBC Geography provided com-putational assistance. My office mates for a majority of this work, Derek van derKamp and Jason Leach, were both great sounding boards and great friends through-out this work.Along with funding from BC Hydro, I have been fortunate to receive fundingfrom the Natural Sciences and Engineering Research Council (NSERC), throughan Industrial Postgraduate Scholarship, and Mitacs, through an Accelerate intern-ship, through this journey. Additionally, Jill Young, Kira Cailes, and Anton Hor-vath from Whistler-Blackcomb helped me with mountain access and insights intosnow safety and forecasting issues in the Coast mountains of BC.I’d like to thank my parents, Bill and Marty Trubilowicz, for all of the encour-agement and access to a great education from a young age. Without them, noneof this would have ever been possible. Finally, I’d like to thank my partner AmyThompson for her support, positivity and flexibility as I completed this thesis. Iknow that the late nights and weekends that I spent working on this thesis haverequired a lot of understanding as we look forward to our lives ahead.xviiiChapter 1IntroductionThe influences of moist maritime air and rugged topography in mountainous coastalregions complicate the more typical snow-dominant hydrograph of much of west-ern North America, as the effect of the rainfall-runoff response is present at low ele-vations, while snowmelt is dominant at high elevations. Depending on weather andclimate conditions, the elevation of the transition between rain and snow can varybetween weather systems and between years. This transition results in some catch-ments in coastal regions being characterized as nivo-pluvial, or ”hybrid” catch-ments, which experience annual peaks in the winter due to rainfall at lower ele-vations, and in the spring due to snowmelt (Wade et al., 2001). Because the hy-drographs of hybrid catchments are highly sensitive to the transition between rainand snow, they are anticipated to be one of the most sensitive hydrologic regimesto climate change as both the winter and summer portions of the hydrographs willbe affected (Schnorbus et al., 2014; Whitfield et al., 2002). Additionally, interan-nual climate variability such as that associated with the El Nin˜o Southern Oscilla-tion (ENSO) can cause these catchments to experience relatively greater peaks inwinter due to rain in warm-phase years (a pluvial-dominant hybrid regime) and rel-atively greater peaks in spring due to snowmelt (a nival-dominant hybrid regime)in cool-phase years (Fleming et al., 2007).Hybrid catchments and maritime snow dominant catchments commonly expe-rience peak flows due to heavy rainfall onto pre-existing snow cover, known as rain-on-snow. Though most common in maritime regions, such as the Pacific North-1west of the United States (McCabe et al., 2007), British Columbia (BC) (Eatonet al., 2010; Loukas et al., 2000), Japan (Whitaker and Sugiyama, 2005), and NewZealand (Fitzharris et al., 1999), rain-on-snow has been observed at nearly everyclimate station that receives snow in the western United States (McCabe et al.,2007), in the eastern United States (e.g. Kroczynski, 2004; Leathers et al., 1998;Pradhanang et al., 2013), and in Europe (e.g. Garvelmann et al., 2015; Merz andBlo¨schl, 2003; Ro¨ssler et al., 2014; Singh et al., 1997; Sui and Koehler, 2001).While peak flow events due to rain-on-snow are the most noteworthy, rain-on-snow can also result in reduced runoff generation due to storage of rainfallwithin the snowpack. This ability of the snowpack to either damp runoff produc-tion through water storage or refreezing (Kroczynski, 2004), or enhance productionthrough rapid snow melt (Marks et al., 1998) and rapid lateral movement of waterthrough the snowpack (Singh et al., 1997), presents a difficult situation for hydro-logic forecasters (e.g Kroczynski, 2004; Ro¨ssler et al., 2014). Forecasting dur-ing rain-on-snow is especially challenging given the typically limited informationavailable for driving operational hydrological models (usually only air temperatureand precipitation). Compounding this difficulty in predicting runoff response dueto rain-on-snow in mountainous environments is a lack of high elevation meteo-rological observations (Moore et al., 2013), which means that hydrologists oftendo not know where or how much rain is falling within a watershed (Ro¨ssler et al.,2014).The high levels of spatial and temporal variability of snow distribution fur-ther complicate the process of understanding rain-on-snow events at the watershedscale. The spatial distribution of runoff generation within a catchment during arain-on-snow event is still not well understood (McCabe et al., 2007), and in re-gions with transient snow cover, hydrologists may not even know where snow ispresent within a watershed. Research has pointed both to runoff generation fromupper elevations with deeper snowpacks as the most important source of runoff(Garvelmann et al., 2015), and large scale synchronized melt from all areas ofa watershed as the dominant control on generation of extreme floods (Jones andPerkins, 2010). Point observations of snow provide little detail on the spatial vari-ability of snow coverage within a watershed and must be extrapolated with caution(Perkins and Jones, 2008), limiting the understanding of the amount of snow stored2in a catchment that is available to melt. Remote sensing offers promise in fillingthis gap in observation, and significant work has already been done in assessing itsusefulness in hydrologic modelling. However, the utility of optical remote sensingof snow in situations with extensive cloud cover (which is common during rain-on-snow) has been somewhat limited (e.g. Andreadis and Lettenmaier, 2006; Lavalle´eet al., 2006).Most research into rain-on-snow has focused on understanding processes thatoccur during individual (or a small number) of rain-on-snow events (e.g. Garvel-mann et al., 2015; Kroczynski, 2004; Leathers et al., 1998; Marks et al., 1998,2001; Ro¨ssler et al., 2014; van Heeswijk et al., 1996), but it is useful to gain an un-derstanding of what differentiates large events from smaller, less significant events.Studies that include many rain-on-snow events of a range of sizes are rare, due tothe extensive long term data collection effort that is required. One of the few sitesin which data records have allowed longer term analysis which includes a range ofevent sizes is the H.J. Andrews Forest in the Oregon Cascades (e.g Jennings andJones, 2015; Jones and Perkins, 2010; Mazurkiewicz et al., 2008), which is primar-ily located in the transient snow zone. The topography of watersheds in southwestBritish Columbia which experience rain-on-snow differs from the H.J. Andrewsforest due to more persistent seasonal snowpacks at upper elevations, wider eleva-tion ranges, and significant occurrence of glacierized and alpine terrain.The broad objective of this research is to improve our understanding of rain-on-snow processes at the regional scale in the mountainous maritime environmentof southwest BC. Because major regional rain-on-snow events are significant tothe hydrology of mountainous coastal areas, but their timing is unpredictable, theresearch objectives have been addressed through a combination of analysis of his-torical data collected by operational hydrologic forecasting agencies and other re-searchers, and incorporation of weather reanalysis model output and snow coverderived from satellite remote sensing.Operational flood forecasting in southwest BC is an underlying motivationfor this thesis. In BC, floods can be generated due to a range of different pro-cesses, including winter rain, winter rain-on-snow, spring/summer rain-on-snow,and spring/summer snowmelt (Loukas et al., 2000). For logistical reasons, in-cluding tight time deadlines and the necessity of producing forecasts for many3watersheds at nearly the same time, operational hydrologic forecasters need to pro-duce accurate forecasts in the face of relatively minimal data collection using con-ceptually based hydrologic models which cannot fully account for the physicalprocesses taking place within a watershed (Weber et al., 2012). Thus, forecast-ing across a range of potential flood generating processes, all which have varyingdominant driving physical processes at work, presents a challenge. For many rea-sons that will be explored throughout this thesis, operational forecasting duringrain-on-snow is particularly difficult. The potentially a-typical snow energy bal-ance and meteorological conditions within many rain-on-snow events, along withan incomplete knowledge of the antecedent snow conditions and snowcover com-bine to make hydrologic forecasting during these events very difficult. Often, theprofessional experience of the hydrologic forecasters makes up for the lack of skillof the model system itself. However, there are cases of rain-on-snow events thathave been mis-forecast (e.g. Ro¨ssler et al., 2014), and practicing hydrologic fore-casters agree that decreasing subjectivity of forecasts is important for increasingconsistency of forecasts (Weber et al., 2012).1.1 Literature review1.1.1 High elevation observationsDue to the logistical difficulties and hazards of installation and maintenance, aswell as the volatile weather that can affect the ability to collect accurate observa-tions, environmental observation networks are often missing high-elevation data(Moore et al., 2013). In situations in which high-elevation data are lacking, obser-vations from low-elevation stations are typically extrapolated to higher elevations,with orographic effects addressed by application of vertical gradients and other to-pographic adjustments. Routines and software packages have been developed toestimate weather variables at individual points (e.g Running et al., 1987) and ingridded spatial fields (e.g. Liston and Elder, 2006; Thornton et al., 1997).Numerous studies have focused on the interpolation/extrapolation of air tem-perature (e.g Bolstad et al., 1998; Garen and Marks, 2005; Hasenauer et al., 2003;Jarvis and Stuart, 2001). In British Columbia, for example, Stahl et al. (2006)4found that knowledge of the vertical temperature gradient was the dominant factorin estimating high-elevation temperatures. A major challenge in estimating high-elevation air temperature is that the temperature gradient can vary substantially inspace and time, and even within and between individual storm events. At the mostextreme, temperature inversions can occur, with significant impacts on the hydro-logic response of a watershed; these situations are not captured well by typicalextrapolation methods.In comparison with air temperature, less research has addressed the extrap-olation of other variables. Water vapour content is often assumed to be relativelyspatially homogeneous, and closely related to the minimum daily air temperature ata site (Running et al., 1987). Kimball et al. (1997) estimated water vapour contentthrough an empirical method based on the daily air temperature, mean annual pre-cipitation, and potential evapotranspiration. Garen and Marks (2005) distributedvapour pressure by elevation through a combination of upper air soundings andlimited surface measurements. Wind speed and direction are highly variable inspace and time, and closely linked to local land cover and topography (Barry,2008). Most work on interpolation/extrapolation of wind observations in moun-tainous areas has involved modelling the effects of topography on wind speed forthe purposes of predicting snow deposition, redistribution, and melt (e.g. Mott andLehning, 2010; Winstral and Marks, 2002). Similar to their approach for vapourpressure, Garen and Marks (2005) used upper air soundings to create three-hourwind speed fields for different elevations. Methods for estimating incoming shortand longwave radiation typically require air temperature observations at the highelevation site. Daily incoming shortwave radiation is often estimated using an ap-proach based on that of Bristow and Campbell (1984), which uses daily air temper-ature range to estimate atmospheric transmissivity. Bristow and Campbell (1984)noted that the accuracy of this method is dependent on site-specific parameter cali-bration. Clear sky longwave radiation is often estimated using air temperature andvapour pressure at a site using empirical formulae such as that presented by Prata(1996).An alternative to interpolation/extrapolation from observational networks is todownscale output from atmospheric models, either to infill missing data (Way andBonnaventure, 2015) or to drive hydrologic models (e.g. Bastola and Misra, 2014;5Wilby et al., 2000), but many atmospheric models are too coarse spatially to repre-sent meteorological variation within mountainous environments at the sub-regionalscale. Assimilation of high resolution weather models offers one potential solu-tion. Deng and Stull (2007) assimilated output from a numerical weather predic-tion model and valley bottom surface temperature observations to improve surfacepotential temperature estimates in mountainous terrain. Customized high resolu-tion point-scale weather forecasts, from mesoscale numerical weather predictionmodels are also available (e.g. those available through the University of BritishColumbia Department of Earth and Ocean Sciences1) as another potential optionfor dynamic tracking of high elevation weather conditions.Regional reanalysis models also offer potential in mountainous environmentsand at smaller scales. Reanalysis models assimilate surface observations withmeteorological model output to produce multivariate state snapshots of the atmo-spheric conditions of the region to which they are applied, using a consistent assim-ilation method over the entire analysis period (as opposed to operational forecastmodels, which may change assimilation methods) (Dee et al., 2011). Fiddes andGruber (2014) developed the TopoSCALE application to downscale output fromthe ERA-Interim 0.75◦ resolution reanalysis model (Dee et al., 2011) and notedthat the same techniques could likely be applied to other regional reanalysis prod-ucts. The North American Regional Reanalysis (NARR) model (Mesinger et al.,2006) is the likeliest candidate for providing high elevation meteorological esti-mates in mountainous areas of North America. However, primarily due to a lack ofground-based observations to test against, the output from this model has not beendirectly tested for hydrologic applications in this manner.1.1.2 Snowpack processes during rain-on-snowThe modern understanding of rain-on-snow processes in North America primarilyoriginates from studies undertaken by the U.S. Army Corps of Engineers (USACE,1956). This report focused on empirical observations of snow processes at multiplesites in the Pacific Northwest and Sierra Nevada regions of the United States thatwere conducted in the 1940s and 50s. Key observations from these original inves-1 were that the majority of water available for runoff (WAR) during rain-on-snow is due to the percolation of rainwater, the turbulent heat fluxes dominatethe energy available for snowmelt, and that advective heat transfer from rainfall isa secondary source of energy. Building on the work of the USACE, Harr (1981)estimated that snowmelt supplemented runoff generation by an average of 17%over rainfall alone for 16 events in the H.J. Andrews watershed; however, theseestimates were calculated from the empirical indices from the U.S. Army Corps ofEngineers (1956), not direct observation.A focus of research into rain-on-snow has been the impact of forest cover onrunoff production during rain-on-snow. In the Sierra Nevada mountains, Kattel-mann (1987) observed 14 rain-on-snow events; in most cases, open sites producedmore water outflow from the snowpack than forested sites, and none of the ob-served sites showed a tendency to delay snowpack outflow. At the watershedscale, Harr (1986) reviewed and reanalyzed streamflow records in the H.J. An-drews Experimental Forest and concluded that forest harvesting increased peakflows over pre-logging conditions. Energy balance modelling of individual rain-on-snow events indicates that the dominant factor in the contrasts in snowmeltbetween forested and open sites is wind speed differences, which drive turbulentheat exchange (e.g. Berris and Harr, 1987; Marks et al., 1998).The snowpack energy balance has been used as an heuristic tool to further theprocess based understanding of snowmelt during rain-on-snow. Most of the energyinvolved in snowmelt during rain-on-snow is thought to be due to the latent andsensible heat fluxes (Harr, 1986; Kormos et al., 2014; Marks et al., 1998, 2001; U.S.Army Corps of Engineers, 1956; van Heeswijk et al., 1996) which are dependenton temperature, humidity, and wind speed. In an analysis of major individual rain-on-snow events using a multi-layer energy and mass balance snow model, Markset al. (1998) found that 60-90% of snowmelt energy came from latent and sensibleheat exchanges due to high temperatures, wind speeds, and humidity.The amount of precipitation falling during rain-on-snow has been found to havea relatively small effect on the amount of snowmelt generated during a rain-on-snow event. van Heeswijk et al. (1996) found that doubling precipitation rates hadthe effect of increasing snowmelt in an energy balance snow model by between2% and 27%. Kormos et al. (2014) found that advective energy due to rainfall7accounted for more than 17% of total energy fluxes during two mid-winter rain-on-snow events in a small watershed in Idaho, U.S.A., and Jennings and Jones(2015) found that advective transport of heat due to rainfall comprised between 29and 44% of the energy balance during events with persistent slow melt.Mazurkiewicz et al. (2008) is one of the few studies to analyze lower mag-nitude, higher frequency rain-on-snow events. When these smaller rain-on-snowevents were included, the traditional view of turbulent energy exchange as thedominant portion of the energy balance during rain-on-snow becomes less clear.Mazurkiewicz et al. (2008) found that radiation (long and shortwave), not the tur-bulent energy exchange, dominated the energy balance and was the largest con-tributor to melt during rain-on-snow over the course of eight years at three sitesin the H.J. Andrews Experimental Forest in Oregon, U.S.A. In an analysis of 20rain-on-snow events that occurred over six years, Berg et al. (1991) found that thetotal WAR produced for rain-on-snow events correlated significantly with precip-itation amount, duration and rate, antecedent snow depth, and snowmelt potential(an index based on air temperature) in the Sierra Nevada Mountains, U.S.A.The complex interactions of the energy balance of the snowpack and internalsnowpack processes makes the prediction of the timing and magnitude of WARgeneration (through percolation of both rainwater and snowmelt) difficult. Alongwith the importance of the energy balance, the ability of a snowpack to store wa-ter is dependent on the depth, density, and internal structure of the snow, and therainfall temperature and intensity. The effects of these processes have been inves-tigated primarily either through plot scale studies (e.g. Kattelmann, 1986; Singhet al., 1997), which do not account for other external energy sources beyond rain-fall, or through combined energy and mass balance modelling of the snowpack (e.g.Marks et al., 1998; Mazurkiewicz et al., 2008; van Heeswijk et al., 1996), whichdepend on our ability to create realistic process based models of these complexsystems.The ability of a snowpack to store rainwater is often considered similar to thatof a coarse sand layer, in which rainfall can percolate and be stored within the layer(Gerdel, 1954). However, this process is complicated by the freezing of rainwaterinto a sub-freezing snowpack and the ongoing metamorphism of the snowpack atnear-freezing temperatures, especially in the presence of liquid water. Understand-8ing of the ability of a snowpack to store rainfall is also complicated by the stratifiednature of a snowpack, where multiple snowfalls have accumulated and formed thecurrent snow layer (Singh et al., 1997).Rainfall percolation through the snowpack during rain-on-snow can be viewedas a threshold process. It is thought that once a snowpack is conditioned, or ’ripe’for release of rainwater, high bulk permeability, augmented by development ofpreferential flow paths, contributes to a rapid contribution of precipitation (as rain-fall) directly to water available for runoff (Singh et al., 1997); however, some fieldresults indicate that the snowpack does not need to reach a ripe stage before out-flow of rainwater occurs (Kattelmann, 1986). Additionally, Conway and Benedict(1994) found that release of latent heat through the refreezing of liquid water wasthe primary energy source for ripening of the snowpack during a rain-on-snowevent, but that this refreezing created ice layers which temporarily impeded theflow of liquid water.It is difficult to determine how much rainwater a snowpack can store beforeit releases WAR based on measurable snowpack parameters. Singh et al. (1997)found that test snow plots in Austria had an average liquid water holding capacityof 6.8% (by volume), but the capacity increased to over 14% when ice layers werepresent.1.1.3 Remote sensing of snowRemote sensing is highly useful for monitoring the spatial dynamics of snow overwide areas, but current shortcomings limit satellite remote sensing for measure-ment of snow covered area (SCA), particularly in hybrid rain-snow catchments.The MODerate resolution Imaging Spectroradiometer (MODIS) instruments, onboard NASA’s Aqua and Terra satellites, provide an appealing combination of spa-tial resolution and short return periods (twice daily) to account for significant cloudobscuration and capture rapid snow dynamics (Dietz et al., 2012).A major limitation of satellite measurements of SCA is blocking due to cloudcover (Dietz et al., 2012; Parajka and Blo¨schl, 2006). In Austria, Parajka andBlo¨schl (2006) found that 63% of days when snow was measured at 754 climatestations were obscured by cloud cover, and this obscuration may be greater in9coastal regions. This limitation has been circumvented by various spatial and tem-poral filters. Eight-day composite images of snow cover are available from theNational Snow and Ice Data Center (NSIDC), which can limit this cloud problemby showing the maximum daily coverage within the eight day period. However,composite images can misrepresent SCA changes during accumulation and ab-lation periods by damping changes that occur under heavy cloud cover. HybridAqua/Terra images are another option. Parajka and Blo¨schl (2008) found that hy-brid images from the two different MODIS instruments (which orbit at differenttimes of day) provided a decrease in cloud obscuration, with only a slight decreasein overall accuracy when compared to ground measurements. Combining the spa-tial filters and hybrid Aqua/Terra images decreased cloud coverage to 46% from63% and decreased accuracy from 95.5% down to 94.2%.Infilling of cloud obscuration often follows a tiered approach, in which in-creasing amounts of cloud removal are traded for decreasing levels of accuracy(e.g. Gafurov and Ba´rdossy, 2009; Parajka and Blo¨schl, 2008). Temporal informa-tion has been incorporated to infill cloud obscuration as well; Dozier et al. (2008)used a space-time infilling method to produce smoothed and infilled ’data cubes’ offractional SCA observations that were produced using the spectral unmixing meth-ods outlined in Painter et al. (2009). Parajka et al. (2010) introduced the SNOWLmethod, in which cloud-obscured pixels are infilled based on the mean snow lineelevation, as determined by non-obscured pixels in the region. This method wasfound to be robust even when cloud coverage was as great as 90% for a region. TheSNOWL algorithm has the advantage of being easily calculable and based on phys-ical reasons for pixel similarity (i.e. the same elevation). Additionally, SNOWLimplicitly includes a middle, ’grey area’, where a pixel is classified as ’partiallysnow-covered’ when the evidence is unclear.Though the SNOWL method implies different accuracy levels in cloud infilledregions, most cloud infilling methods for MODIS snow cover do not provide es-timates of the reliability of these data. Infilled pixels are presented the same asobserved pixels, which may over-represent the confidence in this data and leadto the potential transfer of misinformation about hydrologic processes (Beven andWesterberg, 2011). One obstacle to our understanding of the contribution of snowcover to runoff generation is our inability to directly observe the changes in the spa-10tial distribution of snow cover before, during, and after rain-on-snow events haveoccurred. Maximizing the amount of useful information that can be obtained fromremotely sensed snow cover in commonly cloud-obscured regions is necessary togain a greater understanding of the spatial snow dynamics during these events.Remote sensing of the spatial variability of snow depth and snow water equiv-alent (SWE) is an ultimate goal of snow remote sensing. However, outside of cus-tom data collection campaigns, such as aerial light detection and ranging (LiDAR)(e.g. Hopkinson et al., 2012), which can incur considerable costs, this goal has notyet been fully attained. Passive microwave sensors likely have too coarse a spa-tial resolution for use in mountainous regions (e.g. Tedesco and Narvekar, 2010),and active microwave sensors for remote sensing of snow are still in development(Nolin, 2010). Hence, optical remote sensing of snow via instruments such asMODIS still has considerable importance, and has potential for improved infor-mation gain (as will be explored further in this thesis). Snow covered area fromthe MODIS instruments provide an important complement to combined observa-tional and model based distributed SWE products such as the National OperationalHydrologic Remote Sensing Center (NOHRSC) national snow analysis products2, and spatial snow observations from MODIS have been combined with surfacesnow accumulation and melt models to produce custom high resolution spatiallydistributed snow datasets (Kuchment et al., 2010).1.1.4 Rain-on-snow at the watershed scaleRain-on-snow investigations at the watershed scale have focused on the effectsof forest harvesting (Harr, 1981; Jones and Perkins, 2010), the influence on thehydrograph shape (Perkins and Jones, 2008), and operational forecasting of rain-on-snow events (e.g. Kroczynski, 2004; Ro¨ssler et al., 2014).The importance of rain-on-snow on flood generation for the Willamette Riverin Oregon was illustrated by Harr (1981) who used a range of historical reportsto determine that the majority of annual floods for the river from 1814-1977 weredue to rain-on-snow events. Harr (1981) also found that approximately half of theannual peak floods in a small watershed in H.J. Andrews Experimental Forest in2 Oregon Cascades were due to rain-on-snow.The hydrograph response during rain-on-snow events is generally thought tobe steeper and of a greater magnitude than either purely rainfall or purely snowfallbased runoff (Harr, 1981, 1986; Pradhanang et al., 2013). This has been hypothe-sized to be due to precipitation and air temperature conditions at higher elevationsduring major rain-on-snow events (Ro¨ssler et al., 2014), the rapid lateral move-ment of water through saturated snow packs that can cause a rapid runoff response(Singh et al., 1997), and synchronized snowmelt over an entire catchment duringmajor events (Jones and Perkins, 2010).The physiography of a watershed has been shown to have a profound influ-ence on the shape of the hydrograph during smaller rain-on-snow events. In smallwatersheds (< 1 km2) in Oregon, USA, rain-on-snow floods were analyzed overa 30 year period (Perkins and Jones, 2008). Perkins and Jones (2008) found thatin steeper watersheds, smaller rain-on-snow events had the effect of advancingthe hydrograph response as compared to purely rain events. In flatter watersheds,smaller rain-on-snow events had a delayed response as compared to pure rainfallevents. The catchment topography did not affect the response to large, regionalrain-on-snow events.Field measurements and modelling exercises have demonstrated that open sitesgenerally exhibit higher melt rates than forested sites, primarily due to higher windspeeds, which affect the turbulent exchanges (e.g Harr, 1981; Marks et al., 1998).These results indicate that the spatial contribution of melt during rain-on-snowat the watershed scale may depend significantly on landcover. Areas with forestcover should be expected to contribute less snowmelt than clear areas (either alpineor clearcut), and indeed Jones and Perkins (2010) found that snowmelt increasedduring rain-on-snow following forest harvesting.At the watershed scale, a focus has been on determining why rain-on-snowflooding is often poorly predicted by operational hydrologic models (e.g. Kroczyn-ski, 2004; Ro¨ssler et al., 2014). Numerous potential reasons for these poor pre-dictions are possible. Ro¨ssler et al. (2014) attributed under-forecasting of a majorrain-on-snow event in Switzerland to a lack of observation and process representa-tion of high elevation meteorological conditions, an inability of operational hydro-logic models to represent rapid lateral transport of water within snowpacks, and a12lack of representation of turbulent heat fluxes that drive snowmelt during rain-on-snow. Kroczynski (2004) attributed wide ranging response during rain-on-snow inthe eastern United States to differences in the antecedent cold content and liquidwater content within the snowpack.1.1.5 Atmospheric riversA significant cause of the major rain-on-snow events in southwest British Columbiaare landfalling atmospheric rivers. Atmospheric rivers are concentrated bands ofwater vapour that represent the moisture rich component of the warm conveyor beltof the mid-latitude cyclone, which was first conceptualized by Carlson (1980). Zhuand Newell (1998) first used the term atmospheric river to describe this phenomenadue to the long, narrow (4-500 km wide) band formation. Earlier, using the term’tropospheric rivers’, Newell et al. (1992) found that water vapour in a single riverwas similar to the water flux of the Amazon River (≈ 165 · 106 kg s−1). Zhu andNewell (1998) found that the 4 to 5 atmospheric rivers present on the globe atany given time account for more than 90% of the total net poleward transfer ofmoisture in the atmosphere. The high winds of the pre-cold front low level jet,along with the large amounts of moisture within the atmospheric river and the moistneutral static stability associated with these conditions, can lead to heavy amountsof orographically enhanced precipitation when an atmospheric river makes landfallduring winter (Ralph et al., 2005).Neiman et al. (2008b) catalogued landfalling atmospheric rivers in North Amer-ica, identifying atmospheric rivers through analysis of integrated water vapouridentified by the Special Sensor Microwave Imager (SSMI) instrument (Hollingeret al., 1990) on four polar orbiting United States defense satellites froom 1997 to2005. Neiman et al. (2008b) classified atmospheric rivers as areas with integratedwater vapour depths greater than 2 cm, with widths of less than 1000 km, andlengths of greater than 2000 km. Neiman et al. (2008b) identified 301 atmosphericriver events that made landfall over the Pacific Northwest coast from 1997 to 2005.However, the majority of these occurred during summer, and were not associatedwith significant precipitation due to the ability of warmer air to remain unsatu-rated. Winter atmospheric river events were associated with an offshore trough and13a ridge over the intermountain west, and brought a large onshore vapour flux dueto the association with high winds and the air being near saturation. Atmosphericrivers with greater amounts of integrated water vapour were associated with moreintense storms.Atmospheric river weather patterns have been associated with some of the mostextreme floods on the west coast of North America. The California ’Great Flood’of 1861-1862 is thought to be due to continual atmospheric river conditions thatresulted in constant rainfall for 43 days (Dettinger and Ingram, 2013). This rainfallflooded the entire central valley of California, killing hundreds of thousands of cat-tle and causing millions of dollars worth of damage. More recently, atmosphericriver conditions have resulted in major flood damages in the Pacific Northwest andCalifornia. A major event in November 2006 in western Washington and Oregoncaused extensive flood and debris flows, with approximately 50 million dollars indamage. Neiman et al. (2008a) found that this particularly strong event entrainedmoisture directly from tropical regions of the Pacific ocean and had a maximumintegrated water vapour content of 5 cm. Another major event in January 2009resulted in record flood levels on the Snoqualmie River (1562 km2 drainage area),and disaster conditions to be declared in eight counties in western Washington(Mastin et al., 2010). This event caused record high reservoir levels in Washing-ton’s Howard Hanson Dam and highlighted potential seepage issues in the reservoir(Neiman et al., 2011). Neiman et al. (2011) found that 46 of the 48 annual peakfloods in four basins in western Washington from 1998 to 2009 were associatedwith atmospheric rivers, and Neiman et al. (2008a) found that all of the majorflood events in the Russian River of California from 1997 to 2005 were associatedwith atmospheric river conditions.The well known (along the west coast of North America) ’Pineapple Express’storms represent a subset of particularly strong atmospheric river events that bringin water vapour from subtropical areas near Hawaii (Dettinger et al., 2011). Det-tinger (2004) identified historic pineapple express storms (from 1948 to 1999)across the west coast of North America by identifying moisture fluxes from nearHawaii to western North America in the National Center for Environmental Pre-diction (NCEP) reanalysis model output and found 206 days that this pineappleexpress circulation occurred in this time period. Spry et al. (2014) found that the14precipitation resulting from the intense ’Pineapple Express’ storms, were signifi-cantly larger than other rain storms in southwest BC, and resulted in 11 to 17% ofthe total annual precipitation and water yield at the Capilano River watershed, nearVancouver, BC.Along with the heavy amounts of rainfall that they bring, the warm tempera-tures, neutral stability and high winds of an atmospheric river are conditions thatcan promote rapid snowmelt due to enhanced levels of latent and turbulent heatexchange (see section 1.1.2). Atmospheric river conditions have been associatedwith major rain-on-snow events in the Pacific Northwest and British Columbia (e.g.Marks et al., 1998; Spry et al., 2014), and in Europe (e.g. Ro¨ssler et al., 2014).While major rain-on-snow events have been influenced by atmospheric river con-ditions, the role of snow in these events is somewhat unclear. Neiman et al. (2011)found that atmospheric rivers were a driver for major floods even when the water-shed did not contain significant snow cover, and Perkins and Jones (2008) foundthat the hydrograph shape during major, regional winter rain events was not influ-enced by the presence or absence of snow cover.1.2 Research outlineThe preceding section has provided a brief review of rain-on-snow hydrology, witha particular focus on maritime mountain regions, and has identified a number ofknowledge gaps and challenges in understanding and predicting hydrologic re-sponse to rain-on-snow, especially at the catchment scale. The broad goal of thisthesis is to advance our understanding of rain-on-snow hydrology at a range ofspatial scales to enhance our ability to predict runoff from rain-on-snow events andto improve our understanding of how these processes may be affected by land usechange, climate change, and forest disturbance. Analysis focuses on the southernCoast Mountains of British Columbia, but improving the understanding of rain-on-snow processes and the ability to predict melt will be beneficial to the wider regionof the Pacific Northwest of North America and other mountainous maritime cli-mates which experience rain-on-snow, such as the Southern Alps of New Zealand(Moore and Prowse, 1988).This thesis is organized into six chapters, including four analysis chapters.15Chapter 2 provides an outline of the research area and the data sources used withinthis thesis. Following this, four analysis chapters are included.In the first analysis chapter (Chapter 3), the potential utility of output from theNorth American Regional Reanalysis (NARR) model for hydrologic analysis andmodelling is assessed through direct comparison with high elevation observationsin the Coast Mountains of BC. Following a direct comparison of the performanceof NARR output compared to surface observations and commonly used estimationtechniques, an analysis of the potential for NARR to assist in operational hydro-logic forecasting during an atmospheric river event is undertaken.The second analysis chapter (Chapter 4) focuses on observations taken at auto-mated snow pillow sites in southwest British Columbia during rain-on-snow eventsfrom 2003-2013. This analysis seeks to determine the primary factors governingthe generation of WAR at the point scale during rain-on-snow over a wide spectrumof sizes of rain-on-snow events. Point observations are supplemented with outputfrom a simple snow accumulation and melt model.The third analysis chapter (Chapter 5) focuses on satellite observations of snowcover in the south coast region of British Columbia from the MODIS instruments.A probabilistic snow cover classification method is developed and tested in south-west British Columbia to address the limitations of satellite remote sensing in re-gions with common cloud obscuration. This chapter addresses a common lack ofinclusion of the confidence in estimates of snow cover (both from observed andinfilled pixels) in most remote sensing products.In the final analysis chapter (Chapter 6), knowledge resulting from addressingresearch gaps in the previous chapters is then used at the watershed scale to improveour understanding of the relationship between flood generation, rain-on-snow, andatmospheric river events, which are a major driver of mid-winter rain-on-snowevents and floods in British Columbia and many maritime regions. The role ofatmospheric rivers in the flood generation of selected watersheds in the south coastof British Columbia is assessed, and the potential influences of both high elevationrunoff generation and loss of SCA on runoff generation is assessed during largeatmospheric river events.In the final chapter (Chapter 7), the key conclusions from the four analysischapters are summarized, and potential avenues for future research are proposed.16Chapter 2Regional description and datasourcesThis chapter outlines the study area and the data sources used. Research in this the-sis focused on events at a regional scale and over multiple years; hence, historicaldata from multiple sources were necessary to compile a dataset suitable to answerthe specific research questions.2.1 Regional descriptionThe study was conducted in the southern Coast Mountains of British Columbia(BC), Canada, a rugged, mountainous region with elevations ranging from sea levelto greater than 2500 m (Figure 2.1). The general area shown in Figure 2.1 willhereafter be referred to as the ”south coast region.” Topography in this region isdominated by U-shaped main valleys fed by steep tributary streams, which aretypically headed by cirque basins. The region contains substantial ice cover inthe form of cirque glaciers, as well as high-elevation icefields drained by valleyglaciers.The climate is generally moist maritime, with mild temperatures, wet wintersand relatively dry summers. Low-elevation areas near sea level generally receivemost precipitation in the form of rain, whereas seasonally continuous snow coverwith upwards of 2000 mm of SWE commonly occurs above 1000 m. There is a17marked west-east climatic gradient, with strong maritime conditions on the upwindslopes and the coastal divide, grading to drier conditions east of the divide. Asdescribed in Section 1.1.5, the south coast region can experience atmospheric riverstorms, which are associated with heavy rainfall, high winds, and warm, humid air(Neiman et al., 2008b). A satellite image of precipitable water from an atmosphericriver event on December 10, 2014 is shown in Figure 2.2Elevations near sea level are dominated by the Coastal Western Hemlock (CWH)biogeoclimatic zone, shifting first to the Mountain Hemlock (MH) zone and thenalpine vegetation at higher elevations (Klinka et al., 1991). In the drier, eastern por-tions of the Coast Mountains, the CWH and MH zones are replaced by the InteriorDouglas-Fir and Montane Spruce zones, respectively.In the Campbell River watershed of Vancouver Island, BC, Loukas et al. (2000)found that floods were primarily generated due to either autumn rain or winterrain-on-snow events. However, due to the maritime climate and mountainous to-pography, the south coast region experiences a wide range of complex hydrologicresponses, sometimes all within the same catchment. Watersheds in this regioncan experience high flows in winter months due to rainfall, in the spring and earlysummer due to snowmelt, and in late summer due glacial melt (Eaton et al., 2010).In delineating the hydroclimatic regimes of the southwest British Columbia, Wadeet al. (2001) found that the elevation and distance inland were primary factors gov-erning whether a region was more dominated by rainfall or snowfall in winter, withhigher elevations and areas further inland more likely to be snow dominated. Flem-ing et al. (2007) found that the dominant runoff response could change betweenyears for an individual basin. For example, one catchment was rainfall-dominatedduring the ENSO warm phase, but exhibited a minor spring snowmelt pulse duringthe ENSO cool phase. Catchments in this region typically exhibit a flashy responseto rainfall and/or snowmelt due to the generally steep hillslope gradients and shal-low soils (Winkler et al., 2010).18Figure 2.1: The BC south coast region, including portions of Vancouver Is-land and the mainland of British Columbia. Filled circles indicate thelocations of automated weather stations (AWS), data collection plat-forms (DCP), and hybrid auto/manual weather stations. Filled trianglesindicate locations of automated snow pillows (ASP). Coloured polygonsindicate catchments used in the analysis. Note that the ’CheakamusDaisy’ watershed also includes the ’Cheakamus Millar’ watershed.2.2 Data sources2.2.1 Point dataIn order to characterize rain-on-snow in the south coast of BC, a long term datarecord was compiled through a combination of historic datasets from multiplesources and agencies. The point data sources used in this research project aresummarized in Table 2.1, where Ta is air temperature (◦C), RH is relative humidity(%), w is wind speed (m/s), θ is wind direction (◦), K ↓ is incoming shortwaveradiation (Wm−2), L ↓ is incoming longwave radiation (Wm−2), P is precipitation19Figure 2.2: Example atmospheric river weather pattern, as indicated by pre-cipitable water in the atmosphere. Imagery from the NOAA AMSU in-strument, obtained 12 December, 2014 from, pa is air pressure (kPa), Wxobs are manual weather observations, and SWEis SWE depth (mm).BC Hydro Point data collected at the Daisy Lake data collection platform (DCP),and at the Disappointment Lake and Wolf River automated snow pillow (ASP)sites (see Table 2.1) were supplied by BC Hydro. At present, BC Hydro onlyperforms quality control on daily data; hence the hourly data used in this analysishad to be quality controlled manually. In order to account for small (≤ 1 mm)fluctuations in gauge readings, hourly precipitation was standardized to fit the dailytotal precipitation by proportionally scaling hourly values so that the daily sumequalled the quality controlled daily values. Though data from Disappointment20Table 2.1: Hourly point data used in the south coast region. See text for defi-nitions of variable symbols. RFC = River forecast centre, EC = Environ-ment Canada, Metro Van = Metro VancouverSite Name Source Elev. (m) Date range Variables usedBridge Glacier Shea (2010) 1640 2006-08 Ta, RH, w, θ , K ↓Helm Glacier Shea (2010) 2192 2006-08 Ta, RH, w, θ , K ↓Place Glacier Shea (2010) 2075 summer 2007-08 Ta, RH, w, θ , L ↓Weart Glacier Shea (2010) 2220 2006-08 Ta, RH, w, θ , K ↓Whistler EC 659 Oct. 2003 - Dec. 13 Ta, WxobsUpper Squamish Riv. BC RFC 1387 Oct. 2003 - Dec. 13 Ta, P, SWEChilliwack Riv. BC RFC 1600 Oct. 2003 - Dec. 13 Ta, P,SWEJump Creek BC RFC 1160 Oct. 2003 - Dec. 13 Ta, P,SWEDisappointment Lk. Metro Van 1040 Oct. 2003 - Sep. 09 Ta, P,SWEWolf Riv. BC Hydro 1492 Oct. 2003 - Dec. 13 Ta, P,SWEDaisy Lk. BC Hydro 390 Oct. 2003 - Dec. 13 Ta, PLake were supplied by BC Hydro, the ASP is operated by Metro Vancouver. ASPsites are typically flat and located in small clearings of subalpine treed areas. Asan example of a typical site, a photo from the Wolf River ASP site is shown inFigure 2.3.BC River Forecast Centre The BC River Forecast Centre provided data fromthree ASP sites (Upper Squamish River, Jump Creek, and Chilliwack River). Aswith data supplied by BC Hydro, the River Forecast Centre does not perform qual-ity control on hourly data, so manual quality control was performed. As with pre-cipitation from the BC Hydro sites, hourly precipitation was standardized to matchthe daily total precipitation; however, since daily quality controlled data was notavailable for the full period of record from the River Forecast Centre sites, dailyvalues were calculated as the difference between the first and last cumulative pre-cipitation values for each day.Environment Canada Environment Canada’s Whistler station, officially knownas Whistler - Nesters, was used primarily due to hourly manual precipitation ob-servations that were collected at the site during working hours. The hourly datawere quality controlled by Environment Canada.21Figure 2.3: The meteorological station at the Wolf River ASP site, operatedby BC Hydro. Photo courtesy of Frank Weber, HydrometeorologicalField Programs Scientist, BC Hydro.22Ridgetop automated weather stations Data collected from ridgetop automatedweather stations (AWS) were collected as part of a doctoral research project byJ.M. Shea (Shea, 2010). Site locations were chosen in order to obtain high eleva-tion weather readings and to minimize boundary-layer effects, particularly thoseassociated with katabatic flow (Shea and Moore, 2010). One-second scan intervalswere averaged to produce 10-minute timestep output. The 10-minute-interval datawere integrated to mean hourly values for use in this research. Data were recordedyear-round at the Bridge, Helm and Weart glacier sites, but only during summerat the Place Glacier site. Further details on the instrumentation and collection ofsurface observations can be found in Shea and Moore (2010).2.2.2 Spatial dataWatersheds Six watersheds in the study region were chosen for study (Figure 2.1).The watershed boundary for ’Cheakamus Daisy’ was supplied by BC Hydro. Allother boundaries were supplied by the Water Survey of Canada (WSC). Basin char-acteristics are shown in Table 2.2. Watershed boundaries represent the contributingarea to the outlet of the Daisy Lake Reservoir, and for each corresponding WSCgauge location, accordingly.Digital elevation model The Sensitive Area Mapping Module (SAMM) digitalelevation model (DEM) from Rosin (2010) was used as the base DEM for all anal-ysis. The SAMM-DEM is a 25-m resolution DEM based on the DEM from the BCTerrain Resource Information Management (TRIM) data set that has been hydro-logically corrected to eliminate edge effects between map sheets that may causeproblems in DEM-based flow path modelling. The DEM was resampled to 100-mresolution in order to calculate elevation statistics for the study watersheds.Land cover The Globcover 2009 data sets (Arino et al., 2009) and the BC biogeo-climatic zones of British Columbia (Klinka et al., 1991) were used to classify landcover within the study watersheds. The Globcover 2009 product was further sim-plified into six categories for use in other sections of this research. The resulting23Table 2.2: Summary elevation and area statistics for study basins. Elevationstatistics were calculated from a resampled 100-m resolution version ofthe SAMM DEM.Basin Area (km2) Min. Elev. (m) Median Elev. (m) Max. Elev (m)Cheakamus Daisy 721 352 1401 2629Cheakamus Millar 285 606 1656 2629Chemainus 355 11 621 1493Coquitlam 55 281 1054 1758Cruickshank 214 146 988 2007Elaho 1250 259 1547 2758classes are (1) Alpine/Open, (2) Moderate Canopy, (3) Closed Canopy, (4) Perma-nent Snow/Ice, (5) Open Water, and (6) Urban. The simplified categories assignedto the full Globcover 2009 land cover classification are shown in Table A.1 (seeAppendix A). Land cover statistics for the study basins are shown in Table 2.3.Table 2.3: Simplified land cover (by %) for the study basins by the categoriesdescribed in Table A.1 from the Globcover 2009 dataset.Basin Open Moderate Closed Perm. Snow/Ice Water UrbanCheakamus Daisy 26 37 27 8 1 0Cheakamus Millar 31 25 26 17 2 0Chemainus 1 98 0 0 0 0Coquitlam 14 45 41 0 1 0Cruickshank 13 69 17 0 1 0Elaho 28 29 15 28 0 0Biogeoclimatic zones (Klinka et al., 1991), which represent a combination ofphysiographic and ecological characteristics, are summarized in Table 2.4. Tru-bilowicz et al. (2013) found that the fraction of the catchment covered by CoastalWestern Hemlock and Coastal Douglas Fir was a strong indicator of the contribu-tion of winter rainfall to streamflow. Catchments dominated by Engelmann Spruce- Subalpine Fir, Mountain Hemlock, Montane Spruce, and Alpine had a strongernival influence.MODIS satellite instruments Binary snow cover and individual band data wereobtained from the MODIS instrument on board NASA’s Aqua and Terra satellites.The MODIS instrument has 36 spectral bands, ranging from 0.405 to 14.1 µm,24Table 2.4: Percent coverage of biogeoclimatic zones for each study basin.ESSF = Engelmann Spruce - Subalpine Fir, MH = Mountain Hem-lock, MS = Montane Spruce, CWH = Coastal Western Hemlock, CDF= Coastal Douglas FirBasin ESSF MH MS CWH CDF AlpineCheakamus Daisy 2 23 0 37 0 38Cheakamus Millar 4 25 0 19 0 52Chemainus 0 11 0 87 2 0Coquitlam 0 36 0 37 0 27Cruickshank 0 37 0 44 0 20Elaho 0 19 0 26 0 55with spatial resolution ranging from 250 m to 1000 m.Snow cover delineated via the MODIS satellites is operationally useful sourcefor satellite SCA due to the combination of high temporal resolution (twice daily)and relatively high spatial resolution of 500 m. The NSIDC uses a classificationmethod based on the Normalized Difference Snow Index (NDSI) to detect snowcover:NDSI =b4−b6b4+b6(2.1)where b4 is the 0.55 µm (green) band and b6 is the 1.25 µm (shortwave infrared)band. Additional classification rules are supplied to account for surface waterand warm, light colored surfaces such as sand. These products are referred toas M0D10A1 for data from the Terra satellite, and MYD10A1 for data from theAqua satellite. Due to a band failure on the MODIS instrument on board the Aquasatellite, the Normalized Difference Snow Index (NDSI) for the MYD10A1 usesb7 (2.1 µm) in place of b6. Differentiation between snow and cloud is the primarypurpose of the inclusion of the near infrared bands in the NDSI. The NDSI takesadvantage of the fact that snow has a high reflectance in the visible range (b4) anda low reflectance in the near infrared when compared with clouds, which appearsimilar in the visible range (Dozier, 1989). In the algorithm, an NDSI > 0.4 resultsin a pixel being classified as snow cover. MODIS snow cover makes adjustmentsto the classification method to improve mapping of SCA under forest cover, withadjustments in the algorithm made when a pixel has a Normalized Difference Veg-25etation Index (NDVI) > 0.1 using the algorithm developed by Klein et al. (1998).In addition to the M0D10A1 and MYD10A1 products, individual band sur-face reflectances (known as MOD09GA for Terra and MYD09GA for Aqua) wereobtained and used to calculate the NDSI value for each pixel (calculated as in Eq.2.1). Binary snow data (MOD10A1 and MYD10A1) were obtained via file transferprotocol (FTP) from the National Snow and Ice Data Center 1, and the surface re-flectance data (MOD09GA and MYD09GA) were obtained from the United StatesGeological Survey 2.Landsat Along with satellite imagery from the MODIS instruments, imageryfrom the Landsat 8 Operational Land Imager (OLI) was also obtained. The Land-sat 8 OLI has a higher spatial resolution, but poorer temporal resolution than theMODIS instrument, and was used as a high resolution check on potential regionalspecific re-processing of the MODIS output. Landsat 8 produces images (referredto as tiles) containing 11 bands, ranging from 0.43 to 12.51 µm and with resolutionranging from 15 to 100 m 3. Landsat 8 images the entire earth once every 16 days,and each tile represents an area of approximately 185 km (cross track) by 180 km(along track) 4. The Landsat 8 satellite launched on February 11, 2013.In the south coast region, 66 tiles (representing 34 individual days) with ≤10% cloud cover were available from April 2013 through February 2015. Thesetiles were downloaded and processed using the Google Earth Engine applicationprogramming interface (API) for the Python programming language. Snow coverwas delineated at 30 m resolution for each tile using the algorithm developed byKlein et al. (1998). The Klein et al. (1998) algorithm became the basis of snow de-lineation in the MOD10A1 and MYD10A1 binary snow cover products producedby the NSIDC (Hall et al., 1998; Hall and Riggs, 2007) and uses a combinationof the NDSI and NDVI to assist with delineating snow under tree cover (describedpreviously in paragraph 2.2.2). Salomonson and Appel (2004) used the Klein et al.(1998) method to delineate snow in Landsat for comparison with MODIS snow1 input to the Klein et al. (1998) algorithm, the NDSI was calculated asfollows:NDSI =b3−b6b3+b6(2.2)with band 3 (b3, green, 0.53-0.59 µm) and band 6 (b6, shortwave infrared, 1.57-1.65 µm). The NDVI was calculated as follows:NDVI =b5−b4b5+b4(2.3)with band 4 (b4, red, 0.64-0.67 µm) and band 5 (b5, near infrared, 0.85-0.88 µm).Classification boundary limits in the ’Proposed additional snow field’ in Figure 6of Klein et al. (1998) were estimated as functions by Chokmani et al. (2010).Ground-based photographs of snow cover Ground-based photos of snow cover-age in the Cheakamus River at Daisy Lake catchment were collected over 10 days,shown in Table 2.5 using a Canon 60D digital single lens reflex (DSLR) camerawith a fixed focal length lens from the peak of Whistler Mountain and in multiplelocations in the Callaghan valley, on the Western side of the Cheakamus catchment.Photo locations were chosen after performing a viewshed analysis to determine lo-cations that had the most visibility of the catchment. An example viewshed, fromthe peak of Whistler Mountain, is shown in Figure 2.4. Five viewpoints were usedfor most survey days. The exceptions are 24 July, 2013, when only three view-points were used, and 9 July, 2014, when seven viewpoints were used. Viewpointsare listed in Table 2.6. Differences in total numbers of photographs between dayswith the same number of viewpoints are due to different levels of overlap betweenphotographs. A 50 mm fixed focal length lens was used for viewpoints 3 and 6 inTable 2.6. For all other viewpoints, a 28 mm lens was used.Software for projection of photos onto a DEM, from Corripio (2004), was usedto project the snow photographs onto the DEM for the Cheakamus River at DaisyLake. The software described in Corripio (2004) uses a perspective projection andcreates a linking function between the pixels of the digital photo and the pixels of27Table 2.5: Dates of ground based photo georeferencing and snow delineationin the Cheakamus catchment.Date Number of projected photos7 July, 2012 146 August, 2012 1720 August, 2012 1014 September, 2012 117 October, 2012 165 June, 2013 123 July, 2013 1124 July, 2013 810 September, 2013 129 July, 2014 14Table 2.6: Viewpoints locations used for ground based photo georeferencing.All survey days except 4 July, 2013 (which only used 1-3) used view-points 1-5. The survey from 9 July, 2014 used all 7 viewpoints.Viewpoint Number Latitude (◦N) Longitude (◦W) Elevation (m)1 50.0584 122.9579 21752 50.0592 122.9588 21683 50.1058 123.1200 7034 50.1357 123.1258 8305 50.1355 123.1260 8236 50.1931 123.1836 11857 50.1821 123.1678 1243the DEM. The DEM pixels overlaid upon the image are shown in an example pho-tograph taken at the peak of Whistler Mountain on 6 August, 2012, in Figure 2.5.This linking function is also used in the opposite direction in the software fromCorripio (2004) to map the photograph pixels onto the DEM. After mapping allof the individual photographs onto the DEM using the Corripio (2004) software,k-means cluster analysis was used to delineate snow and non-snow regions on thereprojected photographs. Manual inspection indicated that this method was un-likely to identify snow under full tree cover; hence, this method was only usedduring the summer months. An example result from the k-means cluster analysis,for all of the area of the Cheakamus River at Daisy Lake watershed visible fromthe viewpoints, is shown in Figure 2.6.28Figure 2.4: Visibility of the Cheakamus River at Daisy Lake catchment us-ing the 25-m SAMM-DEM from the viewpoint at the peak of WhistlerMountain (red circle). Visible areas are violet, non-visible areas aregreen.Two photo survey dates coincided (within four days) with Landsat 8 satelliteflyovers and cloudless weather, and allowed for accuracy assessment of the georef-erencing and snow delineation method. Because the photographs were projected atan equivalent resolution of 25 m, they were resampled to 30 m for direct compari-son to Landsat 8 on a pixel-by-pixel basis. The photo survey from 10 September,2013, coincided with a Landsat 8 flyover on 12 September, 2013. Results for areasvisible from the georeferencing viewpoints indicated 94% agreement. The photo29Figure 2.5: The 25 m SAMM-DEM (Rosin, 2010) projected onto a phototaken from Whistler Peak on August 6, 2012, using methods from Cor-ripio (2004). Red dots indicate centroids from the DEM pixels.survey from 9 July, 2014, coincided with a Landsat 8 flyover on 13 July, 2014;results indicated 90% agreement.2.2.3 Stream dischargeDaily discharge was obtained from the WSC HYDAT database for all watershedsshown in Figure 2.1, except for the Cheakamus River at Daisy Lake. CheakamusRiver at Daisy Lake is a BC Hydro operational watershed with power generatedat the Daisy Lake dam. BC Hydro calculates inflows to the Daisy Lake reservoiron the basis of water levels in Daisy Lake, electrical generation, and estimates ofthe amount of water leaving the reservoir (Fleming et al., 2010). Reservoir inflowvolumes for Cheakamus at Daisy Lake were obtained directly from BC Hydro. Allstreamflow data were quality controlled by either the WSC or BC Hydro. Datawere available through 2013 for Cheakamus at Daisy Lake. Data were availablethrough 2011 for the Cruickshank River, and all others were available through30Figure 2.6: Snow cover delineated via the method developed by Corripio(2004) and k-means cluster analysis in visible regions of the Cheaka-mus catchment on 6 August, 2012. White areas inside the CheakamusDaisy basin were not visible at any photo locations.2012.Gauge locations and WSC IDs are shown for all discharge locations in Ta-ble 2.7. The outlet of the Daisy Lake dam is listed as the gauge location forCheakamus River at Daisy Lake; however, as described previously, discharge forthis watershed is based on inflow to Daisy Lake.The annual hydrographs for the test basins are shown in Figure 2.7. The hy-drographs indicate that the Cheakamus River at Millar Creek and Daisy Lake, and31Table 2.7: Gauge locations for study watersheds.Name WSC ID Latitude (◦N) Longitude (◦W)Cheakamus Daisy NA 49.9797 123.1447Cheakamus Millar 08GA072 50.0797 123.0355Chemainus 08HA001 48.8791 123.7019Coquitlam 08MH141 49.4878 122.7925Cruickshank 08HB074 49.5786 125.2139Elaho River 08GA071 50.1139 123.4283Elaho River, are primarily dominated by the spring freshet, but experience somepeak flows in fall. The Coquitlam and Cruickshank rivers experience peak flowsdue to rain in late fall and early winter along with a spring freshet due to snowmelt.Chemainus River appears to be primarily dominated by fall/winter rainfall; how-ever, Fleming et al. (2007) found that Chemainus River exhibited a snowmelt com-ponent during ENSO cold-phase conditions. The hydrologic regimes of these wa-tersheds are consistent with what would be expected based on the BEC zone cov-erage for the watersheds (Table 2.4) and the conclusions from Trubilowicz et al.(2013). The association between Coastal Western Hemlock and and rainfall domi-nation appears to be applicable to the Chemainus River, which contains the highestpercentage of Coastal Western Hemlock (87%) of the study watersheds by nearlydouble (see Table 2.4).Halverson and Fleming (2015) applied a network theory based method for hy-drologic classification of watersheds in British Columbia and Yukon, using a com-munity detection algorithm to identify distinctive hydrologic communities in un-regulated Water Survey of Canada watersheds. Results from Halverson and Flem-ing (2015) indicated that Elaho and Cheakamus River at Millar Creek watershedswere members of the same network community, typified by significant influence ofsnowmelt and glacier melt, and relatively high elevations. The Chemainus, Coquit-lam and Cruickshank river watersheds were found to be members of the communitytypified by slightly lower elevations, and a significant influence of both rainfall andsnowmelt in the annual hydrograph. As it is a regulated watershed, CheakamusRiver at Daisy Lake was not included in the analysis of Halverson and Fleming(2015).32Cheakamus Daisy Cheakamus MillarChemainus CoquitlamCruickshank Elaho020040060005010015020001002003000255075010020003006009000 100 200 300 0 100 200 300Day of yearDischarge (m3 /s)Figure 2.7: Annual hydrographs for the period of record available for eachstudy watershed. Mean daily flow is shown in red. Each year’s data isshown as a grey line. Day of year 1 corresponds to January 12.2.4 Atmospheric river eventsA list of atmospheric river events that impacted the south coast of British Columbiawas supplied by Ruping Mo of Environment Canada. The atmospheric river eventswere selected from a list of events supplied to Ruping Mo by Paul Nieman, whoidentifed landfalling atmospheric rivers as areas over the Pacific Ocean with inte-grated water vapour depths greater than 2 cm, widths of less than 1000 km and33lengths of greater than 2000 km from SSMI satellite imagery in Neiman et al.(2008b) and Dettinger et al. (2011). Those events that impacted the BC coast werethen identified through examination of satellite imagery (Ruping Mo, pers. comm,2015).Suspected atmospheric river events from 2013 were included through exam-ination of BC Hydro weather records (noting days in which the meteorologicalforecaster noted atmospheric river conditions in the weather log), which were pro-vided by Tim Ashman of BC Hydro. Only atmospheric river events that overlappedwith the analysis period of the automated snow pillows (Oct. 2003 - Dec. 2013)were used. Atmospheric river events are listed in the following table (Table 2.8).Table 2.8: List of atmospheric river events identified for the BC south coastfrom 2003 to 2013. Events added via investigation of BC Hydro weatherrecords are indicated in bold.Start Date End Date1 2003-10-15 2003-10-182 2003-11-18 2003-11-183 2004-11-01 2004-11-014 2004-11-06 2004-11-065 2004-11-14 2004-11-156 2004-11-24 2004-11-247 2004-12-04 2004-12-048 2004-12-10 2004-12-109 2005-01-17 2005-01-2010 2005-09-29 2005-09-2911 2005-10-16 2005-10-1712 2005-12-23 2005-12-2313 2006-11-03 2006-11-0714 2007-01-02 2007-01-0315 2007-01-22 2007-01-2216 2007-03-10 2007-03-1117 2007-12-03 2007-12-0318 2009-01-07 2009-01-0719 2009-10-29 2009-10-3020 2009-11-15 2009-11-1621 2009-11-25 2009-11-2522 2010-01-11 2010-01-1523 2010-12-23 2010-12-2524 2012-01-03 2012-01-0425 2013-01-08 2013-01-0926 2013-02-28 2013-03-0227 2013-03-12 2013-03-1434Chapter 3Suitability of North AmericanRegional Reanalysis (NARR)output for hydrologic analysisand modelling3.1 IntroductionThe Predictions in Ungauged Basins (PUB) decade of 2003-2013 recognized a lackof observations as one of the primary obstacles to increasing our understanding ofhydrologic processes (Sivapalan et al., 2003). Great strides have been made inthe realms of remote sensing and alternative observation networks over this timeperiod (Hrachowitz et al., 2013), but high elevation meteorological observations re-main one of the major needs for hydrologic modelling and analysis in mountainousenvironments (Moore et al., 2013).In situations in which high-elevation data are lacking, observations from low-elevation stations are typically extrapolated to higher elevations, with orographiceffects addressed by application of vertical gradients and other topographic adjust-ments. Routines and software packages have been developed to estimate weathervariables at individual points (e.g Running et al., 1987) and in gridded spatial fields35(e.g. Liston and Elder, 2006; Thornton et al., 1997); however, fundamental assump-tions about the vertical gradients of variables must be made in these situations.Output from regional reanalysis models offers an alternative option with poten-tial for making up for the lack of high elevation meteorological stations in moun-tainous environments. With a 32-km spatial resolution, the North American Re-gional Reanalysis (NARR) model (Mesinger et al., 2006) is the highest resolutionreanalysis model for North America. NARR output has been used directly as in-put to macroscale hydrologic models (e.g. Choi et al., 2009; Woo and Thorne,2006), in regional water budget calculations (e.g. Luo et al., 2007; Ruane, 2010;Sheffield et al., 2012) and interpolated and downscaled for use in a glaciation model(Jarosch et al., 2012). Choi et al. (2009) and Keshta and Elshorbagy (2011) foundthat NARR output proved useful for hydrologic purposes, especially when otherobservations were unavailable. However, those two studies focused on the Cana-dian Prairies, where the topography is generally low-relief. The utility of NARRoutput for hydrologic purposes in mountainous regions has remained relatively un-explored.Although the horizontal resolution of NARR output is likely not detailed enoughto drive a hydrologic model for small- to medium-scale catchments (e.g., up to1000 km2 in area) without ground-based supporting observations in a mountain-ous region, it may be possible to use as a supplement to data from observationnetworks, particularly for extrapolating to higher elevations. For example, Holdenet al. (2015) created a 250 m resolution temperature dataset for the United StatesNorthern Rocky Mountains using a combination of modelled solar radiation andsoil moisture, geographic information systems (GIS) based topography and land-cover, and lapse rates derived from NARR temperature output, and found meanabsolute errors < 1.4 K when compared to surface observations for daily maxi-mum and minimum temperature. However, Holden et al. (2015) did not directlyevaluate the accuracy of NARR air temperatures for high elevations, nor did theyexamine variables other than air temperature. Indeed, the performance of NARRin replicating observational data at these locations has received little attention dueto the lack of high-elevation observations in most mountain regions (Jarosch et al.,2012).In this study, the primary goal was to assess the utility NARR output for hydro-36logic analysis and modelling in a mountainous region. This was done by (1) com-paring downscaled NARR output to surface observations at Automated WeatherStations (AWS) located at ridgetop locations (Shea and Moore, 2010), and (2) il-lustrating the potential for NARR output as high altitude intelligence for modellingsnow dynamics and streamflow response to an atmospheric river event, which areoften poorly forecast due to incorrectly specified storm freezing levels (Moore,1993).3.2 Methods3.2.1 Weather observationsAir temperature (Ta; ◦C), relative humidity (RH; %), wind speed (w; m s−1), in-coming shortwave radiation (K↓; Wm−2), and incoming longwave radiation (L↓;Wm−2) from the Bridge, Helm, Place and Weart glacier ridgetop automated weatherstations were used as a ground truth for testing of NARR interpolation methods.Station locations are shown in Figure 2.1 and variables recorded are summarizedin Table 2.1. See paragraph 2.2.1 and Shea and Moore (2010) for further details ondata collection.Post-processing was required to calculate variables that matched the NARRoutput. Vapour pressure (ea; hPa) was calculated as follows:ea = es · RH100 (3.1)where RH is relative humidity (%) and es is the saturation vapour pressure (kPa),computed as in Buck (1981):es =0.611 · exp(17.37·TaTa+238.9): Ta ≥ 00.611 · exp(22.45·TaTa+272.5): Ta < 0(3.2)3.2.2 NARR productsThe North American Regional Reanalysis model (Mesinger et al., 2006) uses sur-face, radar, and satellite observations as input for post processing weather model37output. Output is gridded at a resolution of 32 km, with a three-hour time stepat multiple pressure levels for more than 400 meteorological variables from 1979to near-present. Air temperature (Ta), specific humidity (qa), wind velocity fromeast-west (U) and north-south components (V), K↓, and L↓were extracted for com-parison to surface observations.For Ta, qa and wind velocity (U and V components), NARR output is given atmultiple pressure levels, with each pressure level having a corresponding dynamicgeopotential height. Thus, the interpolation was performed in three dimensions.For each time step, the two nearest heights to the elevation of a site were selectedat each of the four nearest NARR centroids (shown in yellow, Figure 2.1), andthe desired variable was linearly interpolated to the reference elevation of the site.Following this vertical interpolation, inverse distance weighting of the four nearestNARR centroids was used to interpolate horizontally to the station location. Thesemultilevel variables represent instantaneous measures at a three-hour time step.NARR shortwave and longwave radiation are only available as mono-level out-put at the earth’s surface as represented in the NARR digital elevation model. Thesevariables were only interpolated in two dimensions, using inverse distance weight-ing interpolation from the four nearest points. These mono-level variables repre-sent the mean irradiance over a three-hour time step, with the time stamp at thestart of the three-hour interval.The NARR specific humidity and U (westerly) and V (northerly) componentsof the wind vector were post-processed for direct comparison to surface obser-vations. Wind speed was computed from the vector winds with the Pythagoreantheorem and qa was converted to ea (kPa) as in Bolton (1980):ea = qaP0.0378 ·qa+0.0622 (3.3)where P is air pressure supplied directly from NARR (kPa) at the geopotentialheight of the station (m).NARR variables were extracted from the NOAA data server 1.1 Alternative methods for estimating incoming radiationSome commonly used methods for estimating incoming short- and longwave ra-diation at remote sites were applied in order to assess whether the NARR outputprovided any accuracy improvements relative to those methods.Potential direct incoming solar radiation (K↓P) was calculated at each site atboth the three hourly and daily timesteps using standard solar geometry equationsfrom Iqbal (1983) and the potential incoming shortwave solar radiation equationfrom Hock (1999), which was modified by Jost et al. (2012) as follows:K ↓P = I0 ·E0 · τp0cos(θs) · cos(θs) (3.4)p0 = e−z8434.5 (3.5)where I0 is the solar constant, E0 is the earth orbit eccentricity factor (Iqbal, 1983),τ is the mean atmospheric clear sky transmissitivity (set at 0.75 as in Hock (1999)),z is the site elevation, and θs is the solar zenith angle.In order to assess the potential benefits of using NARR output over more com-mon methods, incoming longwave radiation was estimated by the Stefan-Boltzmannlaw using both an assumed emissivity (ε) of 0.85 and an emissivity calculated usingthe Prata (1996) equation:ε = 1− (1+ξ )e−1.2+3.0ξ 0.5 (3.6)ξ =e0T0·46.5 (3.7)where T0 is the surface temperature (K) and e0 is the vapour pressure (hPA).3.2.4 Performance measuresThe NARR downscaling accuracy (and accuracy of alternative estimation meth-ods) was assessed using the mean bias error (MBE), the root mean squared er-ror (RMSE), the Nash-Sutcliffe model efficiency (Em) (Nash and Sutcliffe, 1970),and the coefficient of determination for a linear regression of NARR versus ob-39served (r2), as follows:MBE =1n∑(yi− xi) (3.8)RMSE =√1n∑(yi− xi)2 (3.9)Em = 1− ∑(yi− xi)2∑(xi− x¯)2 (3.10)r2 = 1− ∑(yˆi− y¯)2∑(yi− y¯)2 (3.11)where yi is the value predicted by NARR and xi is the observed value at time stepi, x¯ is the mean value of all observations, n is the number of observations, y¯ is themean of all NARR output, and yˆi is the value predicted by a linear model betweenNARR and observations. Performance measures for both daily mean and three-hourly downscaled data were calculated.3.3 Results3.3.1 Example time seriesSubsets of the downscaled and observed meteorological data at Bridge and Placeglacier sites illustrate the three hourly data (Figures 3.1 and 3.2). NARR-predictedair temperature follows the observed diel and multi-day variations at both sites.Vapour pressure follows well with observations at the Place Glacier site (Figure3.2) but NARR significantly overestimates ea at the Bridge Glacier site during thefirst three days (Figure 3.1). Wind speeds from NARR are of similar magnitudeto observations at the Bridge Glacier site, but do not appear to closely follow in-dividual observations (Figure 3.1). Downscaled incoming solar radiation tracksobservations during most clear days at the Bridge Glacier site, but not during peri-ods with cloud cover (Figure 3.1). Downscaled incoming longwave radiation is ofsimilar magnitude to observations at the Place Glacier site, but does not appear tocorrelate well with observations (Figure 3.2)40ll llllllll llllllll ll l lllll llllllll ll lllll l lllllll l ll ll lllllllllllllll lll l l l lll ll ll llll llll l ll l l lllll lll l llll lll ll l l lll l lll lllllll l lll l lll l l l lllllll l lllll llllll l l ll lllll llll lllllll l lllll l l llll l lll lllll ll lll l l l l ll lllllllllllllllllllll ll llllll l llll llll ll l llll lll lllll llllllllllll lll ll lll llllllll lllllllllll l lllll llllllllllllllllllll l lllllll l lllllll l lllllll l lllllll l lllllll l lll l llllllllllllllll lllllllllllllll llllll51015200. 02 Jul 04 Jul 06 Jul 08DateT a (°C)ea (kPa)w (m/s)K↓(W/m2 )l lNARR ObsFigure 3.1: Observed and downscaled variables for July 1 to 8, 2008, at theBridge Glacier site at a three-hour interval.41ll l ll l llll lllllll lll l l llll llllllll ll lllll llllllll l l lll l lllllllllllllllll llll ll lllll l llllll lllll llllllll l l ll ll lllllllllllllllll lll l ll lll llll llll l llllllll llll lll llllllll lllllll llllllllll lllll llllllll lllllll l llll llllll lllllllll ll llllllllllllllllllll l lllllllllllllllllllllllllllllll lll lllllllllllll llllllllll ll ll llllll051015200.60.81.0240270300330360Jul 02 Jul 04 Jul 06 Jul 08DateT a (°C)ea (kPa)L↓(W/m2 )l lNARR ObsFigure 3.2: Observed and downscaled variables for July 1 to 8, 2008, at thePlace Glacier site at a three-hour interval.3.3.2 Air temperatureThere were strong linear relations between observations and NARR downscaledmean daily air temperature at all sites (Figure 3.3). An over-prediction of 1-2K at temperatures below -10 ◦C can be seen at the Bridge and Weart sites. Theseasonality of error is shown in Figure 3.4. There is no strong seasonal componentto the error, but NARR-estimated temperatures at Bridge and Weart tend to begreater than observed in winter. As Place data were only available during summer,an investigation of seasonality was not possible. The downscaled NARR outputhad a slight negative bias at the Place site, whereas all other sites had a small,positive MBE. RMSE, Em, and r2 improved in all cases by integrating to dailyresolution from three-hour values (Table 3.1).42llllllllllllllllllllllllllllll lllllllll llllllllllllllllllllllllllllllllllllll llllllllll lll llllllllllllllllllllllll llllllllllllllllllllllllllll llllllllllllllllll llllllll llllllllllllllllllllllllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll llllllllllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllll llllllllllllllllllllllll lllllllllllllllllllllllllllllll lllllllllllllllllllllllllllllllllllllllll llllllll llllllllllllllllllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllllllllllllll l lllllllllllllllllllllllllllllllllllllllllBridgeHelmPlaceWeart−1001020−1001020−1001020−1001020−20 −10 0 10 20Observed Ta (°C)NARR Ta (°C)Figure 3.3: Comparison of NARR-derived and observed mean daily air tem-perature. The straight line indicates perfect agreement.43lllllllll lllllllllllllllllllllllllll llllllllllllllllllll lBridgeHelmPlaceWeart−50510−50510−50510−505101 2 3 4 5 6 7 8 9 10 11 12MonthError (°C)Figure 3.4: Monthly distribution of error between NARR-derived and ob-served mean daily air temperature.Table 3.1: Performance measures for air temperature at each site.3-Hour DailySite MBE (◦C) RMSE (K) Em r2 RMSE (K) Em r2Bridge 0.17 1.58 0.96 0.96 1.21 0.97 0.98Helm 0.43 1.67 0.95 0.95 1.07 0.98 0.98Place -0.13 1.31 0.94 0.94 0.80 0.98 0.98Weart 0.47 1.61 0.96 0.96 0.98 0.98 0.99443.3.3 Vapour pressureDownscaled values of mean daily ea are highly correlated with the observations(Figure 3.5), although more weakly than for Ta. Vapour pressure was over-predictedby NARR at the Bridge site, but under-predicted at the Helm site, with low bias forthe Place and Weart glacier sites. Figure 3.6 reveals this is mostly due to differ-ences in error during summer. Performance measures for vapour pressure (Table3.2) indicate good downscaling accuracy, although slightly less strong than for Ta.MBE was negative for Helm and positive for all other sites. As with Ta, dailyintegration improved RMSE, Em and r2 in all cases.Table 3.2: Performance measures for vapour pressure at each site.3-Hour DailySite MBE (hPa) RMSE (hPa) Em r2 RMSE (hPa) Em r2Bridge 0.09 0.13 0.58 0.84 0.11 0.62 0.89Helm -0.05 0.12 0.69 0.75 0.10 0.77 0.84Place 0.02 0.11 0.55 0.60 0.08 0.66 0.70Weart 0.02 0.08 0.85 0.86 0.06 0.91 0.9245lllllllllllllll llll llllll llllllllllllllllllllllllllllll lllllllll llllllllllllllllllllllllllllllllllll llll llllllllllllllllll llllllllll lllllllll lllllllllllllllllllllllllllllllllllllll llllllllll llllllllllllllllllllllllllllll llllllll lllll llllllllllllllll llllll llllll llllllllllllllll lllllllllll ll lllllll llllll lllllllll llllll lllllll lllll llllllllllll llllllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllll lllllll llllllllll lll llllllllllll llllllllll llllllllllll llllllllllllllllllllllllllllll lll ll lllll lllllllllllllllll llllllll lllllllllll llll ll lllllllllllllllllllllllllllllllllllllllllll lllllllll lll llllllll llllllllllllllll lllll llllllllll llllllllllllllllllllllllllllll lll llllll llllll lllllllllllllllllllllllllllll lllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllll lllll lllllllllllllllllllllllllllllllllllllllllllllll lllllllllllllll lll llllllllllllllllllllllllllllllllllllllllllllllllllll llllll llllllll lllllllllllll lll llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll llll lllllll lllll lllllllll lllllBridgeHelmPlaceWeart0. 0.25 0.50 0.75 1.00Observed ea (kPa)NARR ea (kPa)Figure 3.5: Comparison of NARR-derived and observed mean daily vapourpressure. The straight line indicates perfect agreement.46llll llll lllllllll lllllllllllll llllllllllll lllllllllllllllllllllllllllllllllllllllllllllllllllBridgeHelmPlaceWeart−−−− 2 3 4 5 6 7 8 9 10 11 12MonthError (kPa)Figure 3.6: Monthly variation of error between NARR-derived and observedmean daily vapour pressure.3.3.4 Wind speedDownscaled NARR output consistently and strongly over-predicted wind speed atHelm and Weart glaciers (Figure 3.7). The over-prediction was most prevalent inthe winter season, when the southern Coast Mountains experience frequent mid-latitude storms (Figure 3.8). The positive bias is most likely due to NARR being47designed to represent the free atmosphere, while surface friction acts to reducewind speeds measured on the ground.The bias in the NARR-derived wind speeds is reflected in the performanceindices (Table 3.3). In particular, the negative model efficiencies indicate that theobserved mean is a better predictor than the NARR values. RMSE was improvedin every case by resampling to daily resolution; however, Em was decreased forthe Weart site. As expected, MBE was positive in all cases. Although there wasa positive bias, there was a reasonable linear relation between NARR output andobserved wind speeds, particularly at the Weart Glacier site.Table 3.3: Performance measures for wind speed at each site.3-Hour DailySite MBE (m/s) RMSE (m/s) Em r2 RMSE (m/s) Em r2Bridge 0.75 2.60 -0.52 0.16 1.73 -0.19 0.38Helm 2.41 4.39 -1.87 0.34 4.33 -3.37 0.35Weart 4.54 5.84 -5.51 0.65 5.49 -7.18 0.8448lllllllllllllllllllll ll ll llllll llllllll lllllllllllllll lllllll lll l lllllllll llllll lll llllll ll llllll l llllllll llllll llllllllll lll ll lllll l lllllllllllll llllllllll llllllllll llllll llllll lll lllllllll lllllllllllllllll lllll lll lll lllll lllllll ll llllllll l lllll llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll lllllllllllllll lllllllllllllllllllllllllllll llllllllllllllllllllllllllllllll llllll lllllllllll lllllllllllll l lllllllllllllllllll llllllllllllll llllllllllllllllllllllllllllllllll llllllllll llllllllllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllll lllllllllllllllllllllllllllllll llllllllllllll lllllll lll lll lllllllllllllllllllllllllllllll lllll lllllllllllllllllllllllll llllllllll lllllll llllBridgeHelmWeart0102001020010200 5 10 15Observed u (m/s)NARR u (m/s)Figure 3.7: Comparison of NARR-derived and observed mean daily windspeed. The straight line indicates perfect agreement.49llllll l llllllllllllllll l llllllll lllBridgeHelmWeart0510150510150510151 2 3 4 5 6 7 8 9 10 11 12MonthError (m/s)Figure 3.8: Monthly distribution of error between NARR-derived and ob-served mean daily wind speed.3.3.5 Incoming shortwave radiationObserved and downscaled K↓ are correlated, although the spread is greater thanfor Ta, ea, or wind speed (Figure 3.9). Performance measures for K↓ indicatelarge improvements in RMSE when data are integrated to daily resolution (Table3.4). Em decreased when integrated to daily values in all cases, likely due to theremoval of the diurnal cycle of shortwave radiation, which is captured well byNARR (Figure 3.1). MBE was positive for all sites.The downscaled NARR output outperformed the modified Hock (1999) model,which yielded higher RMSE and lower Em and r2 (Table 3.5). The NARR outputand the Hock (1999) model had similar patterns of bias among the three sites.50l lllllllllllllllllll ll ll lll lllll llllllllllll lll lllllllll llll llll lllllllllllllllll lllllllll lll lllllll llllllllllllllllllllllllllllllllllllllllllllllll lllllllllllllllllllllllllllllll llllllllllllllll llllll llll lllllllllllllllllllllllllllll ll lllllllllllllll ll ll lllllllllll llll l llllllllllllll lllllllllllllllllllllllllllllllll lll lllllllllllllll ll llll lllllllllll ll llllll l llll lllll llllllllllllllllllll lll llllllll lllllllll llllllllll lllllllllll lllllllllllll lllllllllllllllllllllllllllllllll llllllllllllllllll lllllllllllll llll lllllllllllllllllllllllll lllllll llllllllllllllll lllllllll llllllllllllllllllll lllllllllllllllllll llllllllllll llllllllllll llll llllll lllllllllll lllllllllll lllllllllll lllllllllll lllllllllllllllllllllll llllllll ll llllllllllll llllllllllllllll lllllll l lllllll lllllll lll lllllllll lllllllll lll llllllllllll llllll llllllllllllll llllllllllllll llllllllllllll lllll lllllllllllllll lllllllllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllll llllllll llllll llllllll lllllllllllllllll llllllllll lllllll l ll lllll llllll llll lllllllllllllllll lll lllllllllll lllll llllllllll llllllllllllllllllllllllllllllll llllllllllllllllll llllllllll lll lllBridgeHelmWeart1002003004001002003004001002003004000 100 200 300 400Observed K ↓  (W/m2)NARR K ↓ (W/m2 )Figure 3.9: Comparison of NARR-derived and observed mean daily K↓. Thestraight line indicates perfect agreement.Table 3.4: Performance measures for K↓ at each site. Observations, MBEand RMSE are in units of Wm−2.3-Hour DailySite Obs. MBE RMSE Em r2 RMSE Em r2Bridge 171 6.3 84.2 0.89 0.89 42.44 0.87 0.87Helm 157 23.5 94.6 0.84 0.86 52.1 0.78 0.83Weart 168 18.4 92.8 0.86 0.87 46.0 0.82 0.8551Table 3.5: Performance measures for estimation of potential K↓ via the mod-ified Hock (1999) method. Observations, MBE and RMSE are in unitsof Wm−2.3-Hour DailySite Obs. MBE RMSE Em r2 RMSE Em r2Bridge 171 0.2 109.4 0.82 0.82 62.1 0.71 0.72Helm 157 28.0 122.3 0.74 0.80 71.6 0.59 0.67Weart 168 18.2 105.1 0.82 0.85 61.1 0.69 0.733.3.6 Incoming longwave radiationObserved and interpolated longwave radiation from NARR are poorly correlated(Figure 3.10), although the period of record was more limited than for the othervariables. With only one site, and only two summer’s data, it was more difficultto assess the performance of NARR for predicting L↓, or the seasonal variabilityin errors. Longwave radiation interpolated from NARR performed better than theempirically estimated alternative (Table 3.6). Although the Stefan-Boltzmann lawwith an emissivity of 0.85 yielded a slightly lower MBE than NARR, the NARRoutput outperformed the alternative approaches on all other performance measures(Table 3.6).llllllllll lllllllllllll llllllllllllllllllllll lllllllllllllllllllllllllllll lllllll lllllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllllllllPlace250300225 250 275 300 325 350Observed L ↓  (W/m2)NARR L↓ (W/m2 )Figure 3.10: Comparison of NARR-derived and observed mean daily L↓.The straight line indicates perfect agreement.52Table 3.6: Performance measures for L↓ at the Place site for the three testedmethods. Mean observed L↓ at the place site was 288 Wm−2. MBE andRMSE are in units of Wm−2.3-Hour DailyMethod MBE RMSE Em r2 RMSE Em r2NARR -5.7 30.8 0.15 0.35 19.2 0.47 0.54Stefan-Boltzmann 4.3 42.5 -0.61 0.01 37.0 -0.97 0.03Prata (1996) -33.2 51.0 -1.32 0.00 46.8 -2.14 0.003.4 Case study: use of NARR products to supplement airtemperature observations for modellingrain-on-snow events in mountain catchments3.4.1 BackgroundIn many applications of operational hydrologic models, air temperatures from val-ley bottom stations are distributed to higher elevations using vertical temperaturegradients calibrated by fitting modelled to observed runoff (e.g. Leavesley andStannard, 1995; Quick and Pipes, 1977). During rain-on-snow conditions asso-ciated with atmospheric rivers, which generate major flood events in the PacificNorthwest of the United States (Neiman et al., 2011) and coastal British Columbia(Spry et al., 2014), hydrologic forecasting can be difficult. One reason for this isthat environmental lapse rates can be substantially lower than the standard moistadiabatic lapse rate or modelled vertical temperature gradients calibrated for nor-mal conditions, leading to under-prediction of air temperature at higher elevationsand incorrect classification of precipitation as snow. As a consequence, the floodresponse can be under-predicted (Moore, 1993). In this section, the potential useof NARR air temperatures to improve hydrologic prediction during a major at-mospheric river event is explored. It is hypothesized that the under-prediction ofstreamflow response during atmospheric river events is caused in part by incor-rectly specified vertical temperature gradients and the associated problems withclassification of precipitation as snow over substantial portions of the catchment.533.4.2 Study catchment and modelling systemThis analysis focused on the 721 km2 Cheakamus River at Daisy Lake catchment(Cheakamus-Daisy), shown in Figure 2.1. Elevation ranges from 369 m to 2594m. The Cheakamus-Daisy catchment is 13% glacierized and approximately 50%of the total area is above treeline. Although the majority of annual runoff derivesfrom seasonal melting of snow and glacier ice, rain-on-snow events during autumnand winter can generate higher peak flows than those associated with snow andglacier melt.Cheakamus-Daisy has been regulated since 1957 for hydroelectric power gen-eration by BC Hydro, the primary power utility of British Columbia (BC Hydro,2002). BC Hydro calculates inflows to the Daisy Lake reservoir on a daily basis us-ing a water balance model based on the changes in reservoir levels and the amountof water released from the dam.The UBC Watershed Model (UBCWM) (Quick and Pipes, 1977) is the refer-ence hydrologic forecasting model for BC Hydro. UBCWM discretizes a water-shed into hydrologic response units within elevation bands. Input data include dailymaximum and minimum air temperature and precipitation, both from observationsand weather forecasts. Input data are distributed by elevation using a multi-segmentapproach in which the vertical temperature gradient varies with elevation (with aninflection point at 2000 m), the amount of daily precipitation, and the daily tem-perature range (Quick and Pipes, 1977). Temperature and precipitation gradientsare determined through calibration to match the observed hydrograph.The UBCWM configuration for predicting inflow into Daisy Lake uses a single,synthetic station as input, known as Cheakamus synthetic climate station (CMS).This station is a linear combination of three separate observation locations fallingwithin the watershed and has a representative virtual elevation of 827 m. The useof a single synthetic station and calibrated vertical temperature gradient, insteadof three stations separately, streamlines data processing, is more robust to missingdata from individual stations, reduces dimensionality of model optimization, limitspotentially unrealistic SWE distributions due to spatial heterogeneity of input data,provides for a single point for downscaling numerical weather predictions for shortterm hydrologic forecasts, and is conceptually consistent with other forecasting54methods (Fleming et al., 2010).The UBCWM configuration for Cheakamus-Daisy has a long-term Em (Nashand Sutcliffe, 1970) of 0.82, and daily forecasts are produced by hydrologic fore-casters to improve on this accuracy even further. However, the model often per-forms poorly during rain-on-snow events, typically under-predicting streamflowresponse.3.4.3 Using NARR to estimate vertical air temperature profilesHigh elevation NARR output was used to modify the standard UBCWM lapsemethod, in order to predict the elevation of the 0◦C isotherm, with as little changeto the standard method as possible. A second synthetic station at the Helm Glacierautomated weather station (AWS) site (NARR-Helm; 2192 m) was produced basedon NARR output, and the following rules were applied:• When there was no inversion and temperatures were below 0◦C at CMS, theUBCWM standard temperature gradient was applied with CMS as the basetemperature.• When temperatures were above 0◦C at NARR-Helm, or for inversions, whenNARR-Helm was warmer than CMS, the UBCWM standard temperaturegradient was applied with NARR-Helm as the base temperature.• For temperatures above 0◦C at CMS and below 0◦C at NARR-Helm, thestandard UBCWM temperature gradient was calculated with a linear inter-polation between the two synthetic stations.Following the estimation of the 0◦C isotherms using NARR guidance, the po-tential for NARR to assist in providing a more realistic vertical temperature gradi-ent during an atmospheric river event from January 2010 in the Cheakamus-Daisycatchment was investigated. The specific atmospheric river event occurred fromJanuary 11-15, 2010 (see Table 2.8). During this event, approximately 170 mmof rainfall was recorded at the CMS station. The atmospheric river event was pre-ceded by further significant rainfall and high temperatures, with an additional 60mm of rainfall occurring in the four days prior to the primary event. Prior to this,three days of relatively clear weather occurred. Temperatures remained high at the55CMS station for the duration of the event, reaching a maximum of 4.5◦C on Jan-uary 12, 2010. Daily minimum temperatures were above freezing from January 9through 13 at the CMS station.Following Moore (1993), the CMS input temperatures (maximum and mini-mum daily) were adjusted for the month of January 2010 to cause the UBCWM0◦C isotherms to match the 0◦C isotherms predicted by the combination of CMSand the NARR-Helm station. When the two 0◦C isotherms were predicted to begreater than 200 m apart, the temperature that would be required at CMS to repro-duce the elevation from the NARR/CMS combination was back calculated. TheUBCWM was then run, without any parameter adjustment or re-calibration, forthe Cheakamus-Daisy catchment to see how predicted reservoir inflows were al-tered when the UBCWM was operated using adjusted temperatures for January2010.Finally, the vertical temperature gradients were calculated, for the daily max-imum temperature, for select times (before, during, and after the peak flows wereobserved in the Cheakamus-Daisy catchment for the January 2010 event) using thestandard UBCWM method and the UBCWM method with temperatures at CMSadjusted to match the 0◦C isotherm predicted with guidance from the NARR-Helmstation, as described above. The full air temperature profile for the NARR gridpoint nearest to the center of the Cheakamus-Daisy catchment (shown on the south-west of the catchment boundary in Figure 2.1) was extracted for each of the selecteddays in order to compare temperature profiles produced by NARR output to thetwo alternative states of temperature profiles that were used within the UBCWMframework.3.4.4 ResultsThe atmospheric river event from January 7 to 15, 2010, would have been signifi-cantly under-forecasted if the watershed model had been run without interventionfrom BC Hydro’s forecasting team (UBCWM unadjusted, Panel 2, Figure 3.11).The standard UBCWM temperature gradient method produced a freezing level be-tween 800 and 1500 m during this event. As a result, the model predicted snowto be falling at the higher elevations of the watershed, when in reality only heavy56rains and rain-on-snow melt could cause the peak inflow rate of greater than 200m3 s−1, based on past experience.Prior to the atmospheric river event, nearly 40 mm of precipitation was recordedon January 1, with little runoff response (Figure 3.11). At that time, the standardUBCWM and NARR-Helm freezing levels were in general agreement. When theatmospheric river event began a few days later, the calculated freezing levels beganto diverge, and runoff was significantly under-predicted by the standard model.The elevation of the 0◦C isotherm predicted by interpolation between the NARR-Helm and BC Hydro’s standard CMS synthetic station data reached a maximumelevation of 2500 m, but the UBCWM standard temperature gradient produced amaximum freezing level of 1500 m. (Figure 3.11). The standard UBCWM methodwould never have more than 50% of the watershed above freezing, whereas incor-poration of NARR temperatures showed that there may have been more than 75%of the catchment receiving rainfall during the period of heaviest precipitation.Prior to the onset of the atmospheric river event (January 7), NARR output in-dicates the presence of an elevated inversion above 1600 m, which is not well rep-resented by the temperature profile assumed by UBCWM. Adjusting the UBCWMtemperature profile to match the NARR-predicted freezing level resulted in tem-peratures in the lower elevations of the catchment that were substantially higherthan both the UBCWM and NARR temperatures (Figure 3.12, upper panel), caus-ing UBCWM to predict apparently excessive snowmelt over a large fraction ofthe catchment and thus to over-predict reservoir inflows in the earlier part of theevent. For January 11, the wettest day of the AR, the NARR-adjusted temperatureprofile provides a closer match to the temperature profile from the raw NARR el-evation output (Figure 3.12, middle panel), producing an improved simulation ofpeak streamflow (Figure 3.11, panel 2). Later in the event, NARR output indicateda nonlinear temperature profile with an overall warmer air column compared to theUBCWM profile (Figure 3.12, lower panel).5702040608050100150200−5051000150020002500Jan 04 Jan 11 Jan 18 Jan 25 Feb 01DatePrecipitation (mm)Reservoir Inflow (m3 /s)Temp (°C)Elev. of 0°C (m)UBCWM NARR−adjusted UBCWM ObservedCMS NARR−HelmUBCWM lapse method Incorporation of NARRFigure 3.11: Summary of hydroclimatic conditions and inflow for BC Hy-dro’s Cheakamus River catchment during an atmospheric river eventduring January 2010. From top to bottom: precipitation, reservoir in-flow, range of daily maximum and minimum air temperature from boththe CMS and NARR-Helm sites, and elevation of the 0◦C isotherm.582010−01−072010−01−112010−01−18100020001000200010002000−10 0 10Air Temp (°C)Elevation (m)adjusted UBCWMNARRUBCWMFigure 3.12: Vertical temperature profiles as predicted by the UBCWM forthe standard CMS temperatures and temperatures at CMS adjusted tomatch the elevation of the 0◦C isotherm and temperature profile fromNARR output at the grid point nearest the centre of the Cheakamuscatchment.593.5 Discussion3.5.1 Potential for use of NARR output for hydrologic analysisThe NARR products are attractive for hydrologic analysis and modelling becausethey include all the variables required to drive physically based models of snowand ice melt, evapotranspiration, and geocryology. Due to the reasonable precisionand low bias, air temperature output from NARR appears to have strong poten-tial for hydrologic applications at high elevations in mountainous regions. Highelevation air temperature information can provide information about precipitationphase, snowmelt, and evapotranspiration in locations where no surface observa-tions are available. The use of NARR for air temperature can likely provide moreuseful information than extrapolation methods that assume static vertical tempera-ture profiles, especially in situations with isothermal or inverted temperature pro-files. The RMSE at the daily scale for NARR air temperatures (ranging from 0.98to 1.21 K) compare favorably to results from Fiddes and Gruber (2014), who founda daily RMSE of 1.95 K for downscaled climate model output in the Swiss Alps.Vapour pressure from NARR also appears reasonably accurate, and could havepotential for use in hydrologic applications in mountainous regions. Vapour pres-sure was relatively unbiased at Place and Weart glaciers, but exhibited a positivebias at the Bridge Glacier site and a negative bias at the Helm Glacier site. Thisindicates that NARR vapour pressure may exhibit a spatially variable bias within acatchment. Despite this, knowledge of vapour pressure where no observations areavailable could be potentially useful for quantifying the direction and magnitudeof latent heat exchange between the surface and the atmosphere, which can be animportant source of energy for snowmelt during major rain-on-snow events (Markset al., 1998).Downscaled wind speed from NARR has a positive bias, consistent with thefact that observed wind speeds are influenced by surface friction and results fromother studies (e.g. Abatzoglou, 2013; Stegall and Zhang, 2012). Fiddes and Gru-ber (2014) found similar RMSE (2.69 m s−1) but lower bias after applying thetopographic corrections from Liston and Sturm (1998) to ERA-Interim reanalysisoutput. The linear relations between downscaled and observed wind speeds, espe-60cially for Weart Glacier (Table 3.3), indicate that a bias correction could potentiallybe applied to more accurately downscale the NARR output to surface observations.Another possibility is that a wind record from a temporary gauge could be extendedby NARR through calibration.The Bridge Glacier site is located on a minor ridge below the main ridge crest.Consequently, wind speeds and directions are modified by anabatic and katabaticflow, surface friction, and steering of winds along the valley axis. This likely ex-plains the low r2 results at the 3-hr time step. Wind speed estimates from NARRare likely only useful as an estimator at the highest elevations in this region.A monolevel 10-m surface wind output is also available from the NARR model.This output was tested as an alternative to the pressure level winds. While therewas less positive bias, likely due to corrections for surface friction, the r2 valueswere similar to the results for the NARR pressure level winds for the Bridge site,and worse for the Place and Weart sites. Hence, it is argued that, because both prod-ucts would require bias correction, the downscaled pressure-level winds should bepreferred because they yielded higher r2 with observed wind speeds.On most sunny days, incoming shortwave radiation from NARR appears to bereasonably accurate (see Figure 3.1); however, cloud cover often does not appearto be captured well by NARR (e.g., July 3, 4 and 7 in Figure 3.1). The positiveMBE in summer (Table 3.4) may be due to small scale convective cloud formationthat is common on warm days in the southern Coast Mountains, which may not becaptured by the NARR model due to its spatial and temporal resolution. NARRdownscaled K↓ outperforms an estimate of potential incoming solar radiation, butnot by a wide margin. Considering the limited accuracy of shortwave radiationfrom NARR, and particularly its bias, it is recommended that caution should beexercised when using NARR shortwave radiation as input to hydrologic analysesand models. The positive bias of incoming solar radiation is consistent with otherstudies that found positive biases in NARR incoming shortwave radiation (e.g.Kumar and Merwade, 2011; Schroeder et al., 2009).Accuracy may be further degraded by the fact that incoming shortwave radia-tion is only available as mono-level output from NARR and hence may not be rep-resentative of the atmospheric transmissivity at a high elevation site. The NARRmono-level surface elevations are more than 500 m lower in elevation than the true61site elevation at the Helm and Weart glacier sites. At the Bridge Glacier site, theNARR mono-level elevation is within 10 m of the AWS elevation, and the BridgeGlacier site does display a lower MBE than either Helm or Weart (Table 3.4).Although fewer data were available to evaluate longwave radiation, it doesappear that NARR downscaled L↓ outperforms common techniques for estimatingincoming longwave radiation, including use of the Stefan-Boltzman law using botha constant emissivity of 0.85 and the Prata (1996) equation. Further evaluation isnecessary, especially during winter periods, to make a more definitive conclusionabout the utility of NARR downscaled L↓ for hydrology. NARR downscaled L↓also performed favourably compared to results from Fiddes and Gruber (2014),who found a RMSE of 29.4 W m−2 using the Stefan-Boltzman law modified witha cloud-cover component.Integrated daily values typically provided improved performance measures.This is thought to be due to the fine-scale time desynchronization errors associatedwith the three-hour timestep NARR output as compared to observations. NARRmodel output may not re-produce the timing of weather changes such as frontalpassages with accuracy on the order of hours; however, these fine-scale timingdiscrepancies are smoothed when values are integrated to a daily timestep.3.5.2 Use of NARR in hydrologic modellingFor the January 2010 atmospheric river event, using NARR output to define verticalair temperature gradients resulted in a higher snowline elevation and a substantiallylarger amount of the watershed experiencing rainfall than with the standard tem-perature gradient. This increase in rainfall caused an increase in predicted reser-voir inflows during this event. The incorporation of input temperatures that wereadjusted to match the NARR interpolated 0◦C isotherm predicted peak reservoirinflows that better matched observed inflows, but the inability of the static temper-ature gradient parameters of the UBCWM to capture inversion or isothermal con-ditions (even when the CMS temperatures were adjusted) is still a limitation andwill remain so unless dynamic temperature gradients can be implemented withinthe operational modelling framework. The adjustment of temperatures to matchelevated freezing levels during inversion or isothermal conditions resulted in tem-62peratures at lower elevations that were higher than in reality and likely resulted inunrealistic snowmelt (see panel 1, Figure 3.12)Given that operational hydrologic models are highly tuned to produce the mostaccurate forecasts possible with minimal input data, and the inputs to the UBCWMwere adjusted without any re-calibration, the results of this adjustment appearpromising. The over-prediction of reservoir inflows during the falling limb of thepeak event (approximately January 15 in 3.11) appears similar in shape to the orig-inal modelled falling limb, but transposed upwards following the increased peakflows predicted on January 11 and 12. Hence, the results may be improved withcalibration during atmospheric river events. The rapid onset of melt when temper-atures rise above freezing is another issue known to forecasters using the UBCWMand is primarily due to the lack of cold content tracking within UBCWM.Although NARR’s temporal latency limits its direct use for near-real-time fore-casting, NARR output could be valuable for assisting in model calibration, particu-larly for the parameters that control runoff generation during rain-on-snow events.Alternatively, switching to an ’atmospheric river’ mode, in which temperature pro-files are temporarily adjusted or perhaps set to isothermal, as suggested by Ro¨ssleret al. (2014) may assist operational forecasters in producing better runoff forecastswhen these meteorological conditions are expected.The ability of NARR temperature output to assist in providing more accuratestreamflow predictions is also an important piece of validation for the NARR airtemperature output itself. Streamflow response is an integrated snapshot of all ofthe processes in action within a watershed. The fact that incorporation of NARRtemperature output from high elevation into a hydrologic model improved accu-racy of discharge estimates while holding all other aspects equal indicates that theNARR output can provide an important benefit beyond traditional vertical temper-ature extrapolation within a mountain environment.Given the errors in prediction shown in this study, the use of NARR output di-rectly to drive energy balance snow models (through the use of wind speed, vapourpressure and radiation) may not generate accurate predictions, especially duringrain-on-snow events. Considering the errors associated with most of the NARRvariables, particularly shortwave and longwave radiation, it is unclear whether theiruse as input to an energy-balance snowmelt model would provide superior perfor-63mance compared to a temperature-index-based model. It would be a useful exer-cise to compare the relative performance of different complexities of snow modelsdriven by NARR output.3.6 ConclusionsThe comparison of NARR output with observations from four high-elevation sitesin the southern Coast Mountains of BC demonstrates that air temperature andvapour pressure are predicted with a high level of accuracy. The air temperaturefields, in particular, appear to be of sufficient accuracy to assist in predicting stormfreezing levels, and are likely to be of value in driving temperature-index snowmeltmodels. Wind speeds at the ridge crest sites were consistently lower than predictedby NARR, as expected, but the reasonable correlation at some sites suggests thatthe NARR values have potential value as indices, subject to bias correction. Short-wave radiation was predicted more poorly by NARR, and may not provide enoughof an improvement over standard models of potential incoming shortwave radiationto be of use in hydrologic applications. Longwave radiation requires further testingduring winter months, though it does provide increased accuracy over commonlyused estimation methods.64Chapter 4Point-scale analysis of the role ofthe snowpack in generating wateravailable for runoff duringrain-on-snow events4.1 IntroductionThe amount of water available for runoff (WAR) during a rain-on-snow event de-pends on (1) the capacity of the snowpack to store rainwater, which is, in turn, afunction of its water content and internal energy content, (2) the amount of rainduring the event, and (3) the availability of energy to melt snow. For example, if asnowpack is initially dry (no liquid water content) and has sub-zero temperatures,some of the energy available to melt snow will be consumed by warming of thesnowpack to 0◦C, and rainwater and meltwater may freeze into the snowpack, re-leasing latent heat and warming the snowpack. Once a snowpack layer is brought tothe melting point, additional liquid water will be retained to satisfy the snow’s wa-ter retention capacity. Note, however, that vertical preferred pathways can conveyliquid water to the ground surface even before the snowpack has become entirelyisothermal at its water-retention capacity (Conway and Benedict, 1994).65A major concern is the extent to which snowmelt during rain-on-snow aug-ments rainfall, thus increasing the magnitude of streamflow response. The UnitedStates Army Corps of Engineers (USACE) snow investigations of the 1940s and1950s concluded that rainfall typically supplies between 70 and 90% of the totalWAR during rain-on-snow (Harr, 1981). In extreme rain-on-snow events this rangecan be exceeded; for a record event in the Oregon Cascades in 1996, snowmeltcomprised up to 57% of WAR (Marks et al., 1998). A primary focus of rain-on-snow research has focused on these large, destructive rain-on-snow events thatresult in significant snowmelt (e.g. Garvelmann et al., 2015; Jennings and Jones,2015; Jones and Perkins, 2010; Marks et al., 1998, 2001). During large rain-on-snow events, the sensible and latent heat fluxes are generally found to be the pri-mary drivers for snowmelt (Berris and Harr, 1987; Marks et al., 1998; van Heeswijket al., 1996). In a modelling study, van Heeswijk et al. (1996) found that the com-bination of high winds and a saturated isothermal snowpack were most likely toproduce rain-on-snow events with greater snowmelt than rainfall in winter months.Although a significant number of smaller rain-on-snow events occur, relativelylittle attention has focused on their ability to generate (or not generate) WAR. Ina study that did investigate rain-on-snow over a magnitude of sizes, Mazurkiewiczet al. (2008) modelled the snowpack energy and mass balance during rain-on-snowevents at three sites in Oregon and found that net all-wave radiation, not the tur-bulent heat fluxes, was the dominant energy source for rain-on-snow throughout aseason. Supporting the USACE empirical results, Mazurkiewicz et al. (2008) foundthat percolation of rainwater through the snowpack (as opposed to snowmelt) wasthe dominant contributor to WAR generation during rain-on-snow on an annual ba-sis. Mazurkiewicz et al. (2008) primarily focused on the analysis of model output atthe seasonal timescale, and did not attempt to differentiate WAR sources betweenlarge, potentially flood generating rain-on-snow events, and smaller events.The advective energy from the rain that falls during rain on snow has generallybeen thought to have a minor contribution to the energy consumed for snowmelt.At a seasonal scale, Mazurkiewicz et al. (2008) found that the advective energy dueto rainfall was less than 3% of the total annual rain-on-snow energy budget; how-ever, the contribution may be quite variable depending on the event. van Heeswijket al. (1996) found that the advective energy due to rain contributed between 1 and666% of the energy budget for three events in the Oregon Cascades. Kormos et al.(2014) found that advective energy due to rainfall accounted for more than 17% oftotal energy fluxes during two mid-winter rain-on-snow events in a small watershedin Idaho, U.S.A., and Jennings and Jones (2015) found that advective energy dueto rainfall accounted for 29-44% of the total energy budget during rain-on-snowevents where persistent snowmelt occurred.The internal conditions of the snowpack prior to the onset of rain-on-snow canimpact the ability of the snowpack to either store (through refreezing or storage ofliquid water) or release rainwater (through percolation) and influence the tendencyfor a snowpack to begin melting (Kroczynski, 2004). The ability of a snowpack tostore water is dependent on the depth, density, and internal structure of the snow(e.g. grain size, internal ice lenses, and macropores), the internal energy contentof the snow, and the rainfall temperature and intensity (Singh et al., 1997). Anunderstanding of how the internal conditions of the snowpack vary at the onset of arange of different rain-on-snow events can potentially help to differentiate betweenevents with substantial WAR generation and those without.A primary limitation to investigating rain-on-snow events over a wide rangeof event sizes is data availability. Though rain-on-snow is relatively common ata range of elevations in coastal regions, measurement of enough rain-on-snowevents to capture and differentiate between large, potentially destructive eventsand small events requires a long term data collection program. Unfortunately, de-tailed studies including sufficient weather data to compute the energy budget, orinvolving lysimeters to assist in quantifying the water budget, are typically short-term, spanning at most a few years. Along with the previously mentioned study byMazurkiewicz et al. (2008), notable exceptions include Jennings and Jones (2015),who drew upon 10 years of lysimeter data supported by comprehensive meteoro-logical data, and Berg et al. (1991), who collected snowmelt lysimeter data for 20events over six years in the Sierra Nevada mountains.Automated snow pillow sites are an alternative source of potentially valuableinformation on rain-on-snow. These sites record variations in snowpack waterequivalent and are usually collocated with instrumentation to measure air tempera-ture and precipitation. In this chapter, data from five snow pillow sites are analysedto generate estimates of the snowpack water balance and water available for runoff67for a range of events over multiple years of record. Based on previous literature andprocess-based considerations, the analysis was guided by the general hypothesisthat the magnitude of WAR generated during rain-on-snow, and the relative contri-butions of rainfall and snowmelt, vary with rainfall amount, antecedent snowpackcondition, air mass characteristics, and season. More specific hypotheses are (1)that the effect of antecedent conditions should be greatest for events with low ratiosof total rainfall to snowpack water equivalent, and (2) that air mass characteristicsmay pose an upper limit to snowmelt rates during rain-on-snow, and thus snowmeltcontributions to WAR generation, for the largest rain-on-snow events at these sites.4.2 Methods4.2.1 Data sourcesThis study was based on data from five automated snow pillows in the south coastregion of British Columbia. Three are operated by the BC River Forecast Cen-tre (Chilliwack River, Jump Creek and Upper Squamish River), one is operatedby BC Hydro (Wolf River), and one is operated by Metro Vancouver (Disappoint-ment Lake). Air temperature (Ta, ◦C), total precipitation (P, mm), and snow waterequivalent (SWE) depth (mm) at an hourly timestep were used from these sites.Additionally, hourly Ta and manual weather observations from the Whistler Envi-ronment Canada station, and hourly Ta and P from the Daisy Lake DCP, operatedby BC Hydro, were used. Data sources are described in further detail in Chap-ter 2. Air temperature (Tn, K) and specific humidity (qa, kg kg−1) from the NorthAmerican Regional Reanalysis (NARR) model were downscaled to the snow pil-low locations using methods described in Chapter Discrimination between rain and snowAt the snow pillow sites, air temperature and precipitation were the only meteoro-logical observations available. Therefore, it was necessary to partition precipita-tion between rain and snow at the sites using only air temperature. The partitioningequation from the UBC Watershed Model (Quick and Pipes, 1977) was used:68Fs =0 Ta ≥ Tt +0.5∆T0.5+ Tt−Ta∆T Tt −0.5∆T < Ta < Tt +0.5∆T1 Ta ≤ Tt −0.5∆T(4.1)where Ta is the observed air temperature (◦C), Fs is the fraction of precipitationfalling as snow, Tt is the transition temperature parameter (◦C), and ∆T is the tem-perature range parameter (◦C). Since eq. 4.1 was applied at individual points in-stead of within a watershed model, as is more typical, an alternative calibrationapproach was developed. The values of Tt and ∆T were calibrated on the basisof manual hourly observations of precipitation phase at the Environment CanadaWhistler meteorological station (station ID Whistler; see Figure 2.1 for location),which is the only station in the region with observations of precipitation phase thatis considered here. The observed fraction of snow was computed as:Fs =SwxRwx+Swx(4.2)where Rwx and Swx are the manually observed rainfall and snowfall rates, respec-tively, with assigned values of 0 (none), 1 (light), 2 (moderate) and 3 (heavy) (Me-teorologic Service of Canada, 2015). The values of Tt and ∆T were optimized usingthe rgenoud package (Mebane and Sekhon, 2011) within the R software package (RDevelopment Core Team, 2015). The calibration minimized the root mean squarederror between observed and predicted Fs. Snowfall (S) and rainfall (R) were thencalculated as:S = Fs ·P (4.3)R = (1−Fs) ·P (4.4)where P is the total precipitation for each hourly time step.694.2.3 Identification of rain-on-snow eventsMazurkiewicz et al. (2008) identified a rain-on-snow event as a weather event thathad measurable rain occurring over every 3-hr period for at least 24 hours withsnow on the ground, continuing until three hours occurred with no rainfall. The cri-teria used in this analysis were a modified version of those used by Mazurkiewiczet al. (2008). Precipitation was integrated to three hour totals of rainfall and snow-fall. The following criteria for identifying a rain-on-snow event were then applied:• Measurable rainfall must have occurred over every three hour period for atleast 12 hours. This criterion was shortened to include shorter events, andbecause the precipitation gauges used at the snow pillow sites were less sen-sitive than those used by Mazurkiewicz et al. (2008).• At least 50% of the total event precipitation and 50% of the precipitationthat occurred over each 3-hr period must have been rain. This criterion wasincluded to account for error in the discrimination model for rain and snow.• At least 5 mm of rain must have fallen. This criterion was included to ac-count for potential errors in precipitation gauge measurements and the dis-crimination model for rain and snow.• At least 10 mm of SWE must have been present on the snow pillow at somepoint during the event, and the SWE must have changed during the event(i.e. not a flat line). This criterion was included to account for errors in thesnow pillow measurements.Following the automated application of these criteria, time series of SWE, pre-cipitation, and air temperature during each identified event were investigated man-ually to remove any events with suspected data errors (e.g. potential bridging onthe snow pillow causing no response during rainfall, or rain so light that it was notdiscernable from gauge noise).4.2.4 Calculation of the mass balance at a snow pillowIn order to characterize the behaviour of the snowpack during rain-on-snow, asnowpack water balance was calculated for each identified rain-on-snow event.70The water balance for a snowpack can be represented as follows:∆SWE = ∆W +∆I (4.5)where W is the liquid water storage (mm WE), and I is the water storage as ice(mm). The terms on the right hand side can be expressed as:∆W = R+M−OF−E (4.6)and:∆I = S+Wt −M−Su (4.7)where R is rainfall (mm), M is melt (mm), OF is outflow of liquid water (mm;which becomes WAR), E is evaporation/condensation (mm), S is snowfall (mm),Wt is net wind transport (mm; positive for deposition, negative for scouring), andSu is sublimation (mm). In this analysis, E and Su were assumed to be negligiblein terms of the water balance, although the associated energy exchanges could besignificant terms in the energy balance. As the focus of this study was rain-on-snow, where rainfall on the snowpack should minimize wind transport, the netwind transport was also assumed negligible. Combining equations 4.5 to 4.7 andre-arranging yields:∆SWE = R+M−OF +S−M = R+S−OF (4.8)For an entire event, one can compute WAR as follows:WAR = OF = R+SF−∆SWE (4.9)Separating contributions to WAR between through-flowing rain and melt is themain source of uncertainty in the water balance calculation. Two end-member con-ditions can be considered. First, if it is assumed that the snowpack is ”saturated”(used here to mean that the water retention capacity is filled) at the beginning andend of the event, then the change in SWE can be given as:71∆SWE = ∆I+ cwr ·∆I = (1+ cwr) ·∆I (4.10)so that:∆I =∆SWE1+ cwr(4.11)and:M = S− ∆SWE1+ cwr(4.12)where cwr is the water retention capacity of the snow, assumed to be 0.05, followingMarks and Dozier (1992) and Jost et al. (2012), for example. The alternative end-member condition is that the snowpack is dry at the start of the event. In this case,some of the melt will have been retained and not become WAR, and M will equal:M = S+SWE(t0)− SWE(te)1+ cwr (4.13)where t0 and te are the start and end times of the event. The difference betweenthe two conditions (which can be thought of as the maximum uncertainty in Massociated with lack of knowledge of initial snowpack water content) is:e = SWE(t0)cwr1+ cwr(4.14)4.2.5 Modelling of the cold content and liquid water content of thesnow packIn order to estimate the cold content (the required energy to raise the temperatureof the snowpack to 0◦C) and liquid water content of the snowpack prior to the onsetof rain-on-snow, the two-layer snowpack model from Jost et al. (2012) was run foreach snow pillow site at a 1-hr timestep. Internal snowpack mass and energy wascalculated as in the ’temperature index with cold content’ melt routine from theJost et al. (2012) model, which was ported from the C# language to C++ for use atthe point scale in this analysis. When air temperature (Ta) is above 0oC, the rate ofenergy consumption for snowmelt (Mr, MJ h−1 m−2) is computed as:72Mr = m f ·Ta (4.15)where m f is a snowmelt factor (MJ h−1 m−2 oC−1). When air temperature is below0oC but greater than the temperature of the upper snowpack layer (Ts), the energytransfer to the snowpack used for warming (Mr, MJ h−1 m−2) is computed as:Mr = m f · (Ta−Ts) (4.16)where Ts is the temperature of the upper snowpack layer (◦C). When air tempera-ture is below 0oC and Ta < Ts, the rate of energy loss from the snowpack (Mr, MJh−1 m−2) is computed as:Mr = cc f · (Ta−Ts) (4.17)where cc f is a cold content factor (MJ h−1 m−2 oC−1). Jost et al. (2012) definedm f as a melt factor specific to forested watersheds, with two parameters to adjustmelt by crown cover. In this analysis, since the snow model was run on a pointscale, m f was not adjusted by the crown closure.It is recognized that turbulent heat exchange, which is not tracked by thismodel, is a primary control on snowmelt during large rain-on-snow events (e.g.Harr, 1981; Marks et al., 1998; U.S. Army Corps of Engineers, 1956). The goalof running the Jost et al. (2012) model was only to provide a simplified estimateof cold content and snowpack liquid water content prior to the onset of rain-on-snow in order to evaluate the effect of antecedent conditions on observed WARgeneration, not to model WAR generation during rain-on-snow events themselves.The use of the Jost et al. (2012) model to estimate the internal liquid watercontent of the snow allowed for an estimate of the contribution of snowmelt (M) toWAR. Hence, the snowmelt (M) that occurred during the event could be estimatedas:M = S+SWE(t0)1+ ewr− SWE(te)1+ cwr(4.18)where ewr is the estimated liquid water retained in the snow (between 0 and 0.05),as calculated at the start of the event by the Jost et al. (2012) model. In eq. 4.18,73it is assumed that the snowpack is at its maximum water retention capacity at theend of the event.The snow model from Jost et al. (2012) cannot operate with missing data,which existed at all sites. Missing data (Ta and P) from each snow pillow wereinfilled by fitting regressions with the Ta and P data from the Daisy Lake DCP sta-tion, which had the least missing data of any climate station used in this analysis.Infilled data are described further in Appendix B. Infilled data were only used forrunning the Jost et al. (2012) model; all other analyses were performed using onlyobserved data.Due to the use of infilled data to operate the model, and the goal of trackingthe internal mass and energy (as opposed to accurate reproduction of SWE), theJost et al. (2012) model was not calibrated for the snow pillow sites. The valuesof m f (A f in Jost et al. (2012) paper) and cc f parameters fit by Jost et al. (2012)– 0.0262 and 0.0556, respectively – were used. The two layer snowpack modelwas specified with a maximum depth of the surface (active) layer of 100 mm, anda maximum liquid water content (cwr) of 0.05. Error statistics between simulatedand observed SWE were calculated for each site. The mean bias error (MBE), rootmean square error (RMSE), and the Nash-Sutcliffe model efficiency (Em) (Nashand Sutcliffe, 1970), as described in Section 3.2.4 were calculated for each site.4.2.6 Calculation of the influence of the energy content of rainfall onmeltFor events where snowmelt occurred, the energy consumed to produce snowmelt(Qm; MJ m−2) was calculated from the snow melt estimated via the mass balanceof the snow pillow as follows:Qm = M ·ρ ·h f ·B (4.19)where M is the snowmelt over the event (m water equivalent), ρ is the density ofwater (kg m−3), h f is the latent heat of fusion (MJ kg−1) and B is the fraction of icein a unit mass of snow (kg kg−1). Gray and Landine (1988) noted that B typicallyfalls between 0.95 and 0.97; here B was assumed to be 0.95. The heat content ofthe rainfall transferred to the snowpack summed over an entire event (Qr; MJ m−2)74was calculated as:Qr =Cv ·∑R · (Tw−Ts) (4.20)where R is precipitation in the form of rain for each timestep (mm h−1), Tw is thewet bulb temperature (◦C), used to estimate the temperature of falling rain at eachtimestep, and Cv is the volumetric heat capacity of water (MJ m−2 mm−1 K−1).The wet bulb temperature was calculated using the following empirical equationfrom Stull (2011):Tw = Ta ·atan[0.151977(RH +8.313659)1/2]+atan(Ta+RH)−atan(RH−1.676331)+0.00391838(RH)3/2 ·atan(0.023101 ·RH)−4.686035(4.21)where RH is the relative humidity (%, from downscaled NARR output). The ratioof Qr to Qm was calculated for each event to estimate the influence of the heatcontent of rainfall on snowmelt.4.2.7 Analysis of rain-on-snow eventsThe summary variables described in Table 4.1 were calculated for each identifiedrain-on-snow event at each site. The events were categorized on the basis of theratio of total rainfall to SWE at the start of the event. Events with R/SWE ≤ 0.05(the assumed cwr) were considered ’small’ events; events with R/SWE > 0.05 wereconsidered ’large’. This differentiation was made to separate between events thatmay not be large enough to cause full saturation of the snowpack (where rainfall isless than 5% of the antecedent SWE) in the event that the snowpack did not haveantecedent liquid water stored, and those events that were expected to cause fullsaturation of the snowpack, whether or not the snowpack had antecedent storedliquid water prior to the event onset.75Table 4.1: Summary variables calculated for each rain-on-snow event at eachsite.Variable Unit SourceEvent length hr MeasuredTotal precipitation mm MeasuredMean Ta oC MeasuredCum. degree hours > 0◦ ◦C MeasuredStarting SWE mm Measured∆SWE mm MeasuredTotal R mm eqs. 4.1 and 4.4Total S mm eqs. 4.1 and 4.4Total snowmelt (using ce) mm eq. 4.18Total snowmelt (assuming saturation) mm eq. 4.12Total snowmelt (assuming dry snowpack) mm eq. 4.13Maximum e in snowmelt estimation mm eq. 4.14Total WAR mm eq. 4.9Init. Ts oC (Jost et al., 2012)Init. snow saturation % (Jost et al., 2012)Init. snowpack cold content MJ h−1 m−2 (Jost et al., 2012)Total Qm MJ m−2 eq. 4.19Total Qr MJ m−2 eq. Recursive partitioning to predict WAR productionThe primary drivers for the production of WAR during rain-on-snow were investi-gated using recursive partitioning to create regression trees via the ’rpart’ packagein R (Therneau et al., 2015). This approach is more flexible for characterizingnonlinear dependencies (such as thresholds and interactions) than linear regressionmethods. Regression trees successively split the response variable into two cate-gories based on a threshold value of a predictor that maximizes the reduction inan error metric. When the error metric can no longer be reduced by any predictorvariables, the end is referred to as a leaf. Each split is referred to as a node. Themean of the predictand for all members is displayed at the end of each leaf. In thisanalysis, total R (mm), total S (mm), initial snow saturation (%), initial cold con-tent (MJ h−1 m−2), starting SWE level (mm), Qr (MJ m−2), and the accumulateddegree-hours above 0◦C were used as predictors for total WAR generation (mm).76A potential issue with regression trees is the presence of collinearity of inputvariables. If two predictor variables are highly correlated, either variable could bechosen as a split to the response variable. While this does not necessarily affectthe predictive ability of the regression tree, it does complicate the interpretationof the tree. Collinearity was a possibility in this analysis. In particular, Qr iscomputed from both Ta and rainfall; hence it is likely closely related to both ofthese variables. Additionally the snow that fell during an event is likely inverselyrelated to the accumulated degree-hours above 0◦C. Another potential issue withregression trees is over-fitting. In this analysis, each regression tree was pruned tothe minimum cross-validation error to help avoid overfitting. Regression trees werefit separately for the small (R/SWE ≤ 0.05) and large (R/SWE > 0.05) events toinvestigate how the dominant controls varied between event magnitudes.4.3 Results4.3.1 Discrimination between rain and snowResults of the fitting of equation 4.1 to the precipitation observations are shown inFigure 4.1, where points are observations made by Environment Canada weatherobservers, and the red line indicates the optimized Equation. The best fit parame-ters are Tt = 1.11 ◦C and ∆T = 4.39 ◦C.4.3.2 Point modelling at snow pillow sitesObserved and simulated SWE are shown in Figure 4.2, and performance measuresare summarised in Table 4.2. Overall, the model under-predicted SWE at all sitesbut Jump Creek. Performance was worst at the Chilliwack River site, which maycorrespond with the high levels of wind exposure at this site resulting in precip-itation gauge undercatch (Doris Leong, BC River Forecast Centre, pers. comm,2015). The Jost et al. (2012) model predicted that the majority of the identifiedrain-on-snow events had a saturated snowpack and no cold content present at theevent onset. Only 55 of the 286 total identified events were predicted to be un-saturated at the event onset, and of those, only seven were predicted to have coldcontent present.77ll llll lllll l lllllllll ll lllllllll l l llllllll ll l lll lll llll l ll l lllllllllll llllll llllllll lllllllllllll lll l lll l lllllllll l l l l lllll llllllll0.000.250.500.751.00−4 0 4Ta (°C)F snFigure 4.1: Rain-snow discrimination equation from the UBC WatershedModel, fit to data from the Environment Canada Whistler station. Man-ual observations are in black, the optimized model is in red.4.3.3 Event scale analysisThe 286 rain-on-snow events that were identified are summarized by snow pillowsite in Table 4.3. Figure 4.3 further breaks down the identified rain-on-snow eventson a site by site basis. Note that the identified events are not a comprehensive list,as each site had missing data at different times.Figure 4.4 shows the amount of WAR generated for each event in relation tothe amount of rainfall during the event, coloured by the R/SWE ratio. Symbolsindicate whether an event was associated with an atmospheric river event (as indi-cated in Table 2.8). For events with rainfall in excess of 0.05 of SWE, WAR variesapproximately linearly with event rainfall, and the majority of large events wereassociated with atmospheric river occurrence. For events with rainfall in excess of78CHILDISJUMPUSRWOLF0500100015002000250005001000150020000500100015002000050010001500200005001000150020002004 2006 2008 2010 2012 2014YearSWE (mm)SimulatedObservedFigure 4.2: Observed and simulated SWE (via the model from Jost et al.(2012)) at each snow pillow site.approximately 40 mm, the snowpack consistently enhanced the amount of WARgeneration through snowmelt. The following linear model was fit for events with> 40 mm:WAR = 1.26 ·R+17.93 (4.22)indicating that snowmelt typically enhanced WAR generation by over 25% aboverainfall alone. The linear model had an adjusted R2 of 0.87, standard error of 27.6mm, and p < 2.2 ·10−16The influence of site specific differences on the relationship between rainfall79Table 4.2: Error statistics for the application of the Jost et al. (2012) model ateach snow pillow site.Site MBE (mm) RMSE (mm) EmCHIL -426 568 0.31DIS -237 396 0.69JUMP 39 268 0.80USR -99 249 0.86WOLF -75 276 0.80and WAR for large events was investigated. Figure 4.5 shows all events with > 40mm of rainfall, coloured by ASP site. Linear models were fit for each individualsite. The slopes of linear models are similar across all sites except the ChilliwackRiver site, ranging from 1.19 (Wolf River) to 1.28 (Disappointment Lake). TheChilliwack River site linear model had a slope of 1.47, potentially indicating agreater enhancement via melt from WAR generation; however, this slope was in-fluenced by only one event with greater than 100 mm of rainfall at the site. The Up-per Squamish River and Disappointment Lake were the only sites with any heavyrainfall events in which the snowpack damped WAR generation.For events in which SWE decreased, the ratio of snowmelt to WAR variedwith the time of year (Figure 4.6). During mid-winter rain-on-snow, snowmelt wastypically a smaller contributor to WAR to than other sources (e.g. through flowingrain), whereas during the spring and early summer (May and June) snowmelt wasthe larger contributor to WAR. As an illustration, a very large (the largest recordedat any site in this analysis) mid-winter event is shown in the left panel of Figure 4.7.This event recorded 343 mm of total rainfall (and 5 mm of snowfall). Duringthis event, 449 mm of WAR was generated and snowmelt was estimated as 108mm, resulting in a snowmelt/WAR ratio of 0.24. A moderately sized late seasonevent, from the Jump Creek snow pillow, is shown in the right panel of Figure 4.7.Moderate rainfall (12 mm), and snowmelt of 31 mm resulted in WAR generationof 45 mm and a snowmelt/WAR ratio of 0.68.80Table 4.3: Summary of the number of rain-on-snow events identified at eachsnow pillow site.R/SWESite Name Abbreviation ≤ 0.05 > 0.05 TotalChilliwack River CHIL 37 22 59Disappointment Lake DIS 20 17 37Jump Creek JUMP 43 60 103Upper Squamish River USR 27 20 47Wolf River WOLF 18 22 404.3.4 Snowmelt rate and advective energy exchangeThe contribution of energy advected via rainfall to the total energy consumed bysnowmelt is shown in Figure 4.8 for all events in which snowmelt occurred. Typ-ically, Qr contributed less than 10% of the total energy required for melt, and theevents with the largest contribution were the smallest events (with less than 10 mmof snowmelt). The highest Qr/Qm ratio for an event with more than 10 mm ofsnowmelt was for a large summer event at the Disappointment Lake ASP, begin-ning on July 20, 2007. This event recorded 160 mm of rainfall and had an estimated77 mm of melt, with Qr contributing 14% of the total energy consumed for melt.4.3.5 Regression tree analysisIn the regression tree for large events (R/SWE > 0.05) (Figure 4.9), the primarysplit, and all successive splits except one, are based on storm rainfall, suggestingthat the variation of snowmelt among events was a second-order control on WARgeneration. The accumulated degree-hours above freezing, which integrated theeffects of temperature and the length of the event, was included for the smallest ofthe ’large’ events.For events with R/SWE≤ 0.05 (Figure 4.10), the regression tree is more com-plex, indicating that antecedent factors may may have a greater influence on WARgeneration than for larger rain events. The primary split is based on Qr, whichdepends on rainfall, air temperature and humidity, and thus integrates three of thekey meteorological controls on WAR generation via both rainfall and snowmelt.For events with little energy input due to rain, the saturation level of the snow ap-81lllllllllllll255075100CHIL DIS JUMP USR WOLFSiteEvent Length (hr)lllllllllllllllllllllllllllll0100200300CHIL DIS JUMP USR WOLFSiteEvent Precipitation (mm)snowrainlllllllllllllllll04080120CHIL DIS JUMP USR WOLFSiteEvent Melt (mm)lllllllllllllllll0100200300400CHIL DIS JUMP USR WOLFSiteWAR (mm)Figure 4.3: Summary of identified rain on snow events on a site by site ba-sis. See Table 4.3 for definitions of site name abbreviations. Note thatnegative melt indicates net water storage.peared to affect whether or not the snowpack enhanced WAR generation. If thesnow was unsaturated, the events showed a tendency toward low WAR production;if the snowpack was saturated, the SWE depth at the start of the event appeared tohave an effect. For events with more energy input from rainfall, the total rainfallwas the dominant predictor for WAR production.Though the regression tree for small events includes more variables (indicatingthat the factors dictating WAR generation for events with R/SWE≤ 0.05 are more8201002003004000 100 200 300Event Rainfall (mm)WAR produced (mm)Rain/SWE≤ 0.05> 0.05ARNon-ARFigure 4.4: The relationship between the amount of WAR generated and thetotal rainfall for an event. Points are coloured by the ratio of rainfall toSWE. A linear model, fit for events with≥ 40 mm rainfall is shown as adashed black line. A 1:1 line is shown in solid black. Symbols indicatedwhether or not the event occurred during a known atmospheric river.nuanced), minimum cross validation error was higher in the small event tree (0.54)than in the large event tree (0.25), indicating that there is substantial variability thatis not explained in the data available in this analysis for small events. Additionally,cross validation depends on random sorting of events, and since the predictive abil-ity of the regression tree for events with R/SWE ≤ 0.05 was not strong, pruningvia cross validation could result in anywhere between three and eight splits in theregression tree, depending on the random seed for cross validation. Thus, only thefirst few splits in the regression tree for events should be viewed as robust; splits83100200300400100 200 300Event Rainfall (mm)WAR produced (mm) ASPCHILDISJUMPUSRWOLFARNon-ARFigure 4.5: The relationship between the amount of WAR generated and thetotal rainfall for events with > 40 mm of rainfall, coloured by site. Alinear model is fit for each site. A 1:1 line is shown in solid black.Symbols indicated whether or not the event occured during a knownatmospheric the ends of the branches should not be over-interpreted. The primary conclu-sion that can be drawn from Figure 4.10 is that antecedent factors likely do havean effect on WAR generation during these events, whereas antecedent effects areoverwhelmed during large events.84Jan FebMar AprMay JunJul OctNov Dec024602460246024602460.25 0.5 0.75 0.25 0.5 0.75Snowmelt/WARCount Rain/SWE≤ 0.05> 0.05Figure 4.6: The distribution of snowmelt /WAR ratios for all events with pos-itive melt, separated by month. The black vertical line indicates a ratioof Discussion4.4.1 WAR generation during large eventsThis analysis demonstrated that rainfall was consistently enhanced by 25 to 30%during events with more than 40 mm of rainfall. This enhancement of heavy rain-fall by snowmelt can often make the difference in the generation of large floods,both due to its augmentation of precipitation and ability to synchronize peak dis-charge over multiple sub-basins (Jones and Perkins, 2010). Jones and Perkins85lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll ll lllll lll lllllll ll l l ll l ll ll lllll ll12345660901201500. 00:00 18 00:00 19 00:00 20 00:00 22 03:00 22 06:00 22 09:00 22 12:00 22 15:00 22 18:00Date Time (January 2005) Date Time (May 2006)T a (°C)SWE (mm)Precip. (mm) SnowRainFigure 4.7: Example rain-on-snow events from the Jump Creek ASP. Theleft panel shows results for a large mid-winter event from January 2005;the right panel illustrates a moderately sized late season event from May2006.(2010) also found that events with large amounts of rainfall resulted in enhance-ment of WAR through melt.Though WAR was enhanced by snowmelt during major rain-on-snow eventsat the ASP sites, Figures 4.4 and 4.9 demonstrate that the total rainfall occurringat a site was the dominant factor influencing WAR generation during large rain-on-snow events that have the greatest potential for flood generation in the southcoast region. Figure 4.4 also shows that these large high-elevation rain events arecommonly associated with atmospheric river occurrence. In order to know whenrain is falling, it is necessary to have accurate measurements of temperature at thesites of interest. This analysis reinforces the importance of understanding whenand how much rain is falling at high elevations within mountainous regions suchas southwestern BC. The ’low hanging fruit’ for potential improvement of pre-86lllllll llll l l llll llll lllll llllllllllllllllllllllllllllll lllllllllllllllllll lllllllllllllllll llllllllllllllll lllllllllllllllllllllllllllllllllllllllllllllllllllll llll llllllllllllllllllllllllllllllllllllllllllllllllllllll0. 30 60 90 120Total snowmelt (mm)Q rQ mFigure 4.8: The contribution of advective heat from rainfall to total energyconsumed by snowmelt.diction of response during rain-on-snow at the watershed scale may simply moreobservations of precipitation and temperature at high elevations, as has been calledfor on multiple occasions (e.g Moore et al., 2013; Sivapalan et al., 2003), and wasdemonstrated through incorporation of high elevation temperature information intoan operational hydrologic model for forecasting runoff during an atmospheric riverin Chapter 3. Knowledge of the amount of rain falling at high elevations, alongwith the finding of approximately 25 to 30% enhancement via snowmelt duringmajor rain-on-snow storms, such as winter atmospheric river events, should assisthydrologic forecasters in predicting runoff response during expected rain-on-snowevents.4.4.2 Seasonal variation of WAR compositionMost research into rain-on-snow has focused on mid-winter events, but when smaller,high frequency rain-on-snow events were included in the analysis, late seasonevents represented a significant number of the total identified events. Of the 286 to-87|rain< 77.09rain< 38.18degreehrs< 152.9rain< 148.332.8n=6772.1n=783.1n=38161n=21294n=8Figure 4.9: Regression tree fit for large rain-on-snow events (R/SWE >0.05). Total WAR generation for the event (mm) and the number ofmembers (n) are indicated at the end of each leaf. rain = total eventrainfall (mm), degreehrs = accumulated degree hours above freezing(◦C hr −1).tal events identified, 123 occurred from May to July, with June and May being thefirst and second most common months in which rain-on-snow events occurred inthis analysis. Most large events (R/SWE > 0.05) occurred between November andJanuary, months that are associated with shallower snowpacks and higher rainfalltotals compared to the late-season events.During mid-winter rain-on-snow, rainfall percolation was the largest contribu-88|sum.qr< 0.1047start.sat< 100start.SWE>=1446snow< 2.818rain< 30.46rain< 19.02degreehrs< 120snow>=10.9714n=1621n=19 30n=3243n=1742n=21 50n=1269n=868n=787n=13Figure 4.10: Regression tree fit for small rain-on-snow events (R/SWE ≤0.05). Total WAR generation for the event (mm) and the number ofmembers (n) are indicated at the end of each leaf. sum.qr = Qr for theevent (MJ m−2), rain = total event rainfall (R, mm), start.sat = initialsnowpack saturation (%), start.SWE = initial snow water equivalent(mm), snow= event snowfall (mm), degreehrs= accumulated degreehours above freezing (◦C hr −1).tor to WAR, whereas during the spring and early summer (May, June, July) meltwas the larger contributor to WAR (Figure 4.6). This seasonal shift in the dom-inant source of WAR likely reflects three factors. First, the snowpack in May toJuly is likely already saturated, devoid of cold content, and potentially already ab-89lating. Second, air temperatures are typically higher during the late season, leadingto greater availability of energy for melt. Third, very large rain events, such as the’Pineapple Express’ storms that affect the south coast region are more commonlylate fall to mid-winter storm events (Spry et al., 2014). In relation to the secondfactor, Kormos et al. (2014) found that net all-wave radiation played a greater rolein the snowpack energy balance during rain-on-snow events in the spring and earlysummer over the winter of 2010/2011 in southwest Idaho, USA.4.4.3 Effects of antecedent conditionsDuring events with R/SWE ≤ 0.05, the presence of an unsaturated snowpack (aspredicted by the Jost et al. (2012) model) did appear to reduce WAR generationwhen the contribution of advective energy due to rainfall was low (see Figure 4.10).The analyses did not detect any additional influence of cold content over and abovethe effects of saturation levels. However, only seven out of a total of 286 eventswere predicted to have any cold content prior to the onset of rain-on-snow. Thislack of cold content is not unexpected considering the strong maritime influencesin the study region. The dominance of isothermal, saturated conditions at the onsetof rain-on-snow is reflected in the lack of an observed tendency for the snowpackto retain a net amount of water (negative melt in Figure 4.3) during an event, aswas noted by van Heeswijk et al. (1996).4.4.4 Controls on melt generationIn this study, snowmelt was typically about 25 to 30% of total event rainfall duringmid-winter events with more than 40 mm of rainfall. Note, however, that rainfallwas likely under-measured due to wind effects, such that the estimated snowmeltcontributions will tend to be over-estimates. None of the very large events hada snowmelt (or excess runoff) contribution similar to rainfall, as has been docu-mented during some large events at open sites (e.g Marks et al., 1998). Investiga-tion of this relationship by ASP site (Figure 4.5) showed that snowmelt enhance-ment (and lack of very large events with large amounts of snowmelt) was fairlyconsistent across all sites. This likely reflects the nature of the snow pillow sites,which are typically located in small gaps in areas with open forest, and would thus90be somewhat sheltered from wind. Consequently, the turbulent exchanges of sen-sible and latent heat would be limited by wind speeds, relative to their potentialvalues at more exposed sites. The Chilliwack River ASP did show a tendency to-wards more enhancement via melt. Though this relationship should be interpretedwith caution due to the limited number of large events recorded at the site, it maybe an indication of greater enhancement of melt due to turbulent fluxes at this site.The Chilliwack River site is potentially the most wind exposed of the ASP sitesused in this analysis (Doris Leong, BC River Forecast Centre, pers. comm, 2015).Jones and Perkins (2010), who also found that events with more than 40 mm ofrainfall always resulted in enhancement of WAR through melt, attributed smallerrain-on-snow events where WAR was equal to or less than total rainfall to bothmixed rain and snowfall, and to the absorption of water by the snowpack.The advective heat due to rainfall had a relatively small role in the total energyrequired to produce snowmelt (Figure 4.8). The ratio of Qr to Qm was typicallyless than 0.1. Situations when Qr did contribute a larger percentage of Qm wereprimarily events where small amounts of snow melt (<10 mm) occured. The ad-vective heat due to rain was the primary split in the regression tree for small events,however (Figure 4.10). Though this inclusion indicates the importance of Qr, it iscorrelated with both rainfall and air temperature; hence, this inclusion does notnecessarily indicate the importance of the heat content of the rainfall itself.4.4.5 Sources of uncertaintyThe rain-on-snow events identified at these sites are not an exhaustive list of events.The identification criteria were utilized to provide the best chance that an event wasindeed rain-on-snow, but could have resulted in very light rain events and eventsthat were near the threshold between rain and snow being excluded from the anal-ysis. Additionally, missing data at some sites may result in particular events beingmissed at some locations. The analysis relied on the separation of rain and snowbased on air temperature. As described in Section 4.2.3, a threshold of >50% rain-fall for each time step was used to confirm that the ASP was actually experiencingrain during the event. However, this threshold also eliminated some of the small-est rain-on-snow from being included in this analysis. Although the identification91of small events and estimation of the water balance are potentially compromisedby these sources of uncertainty, the magnitudes of rainfall and snowmelt in thelarger events are sufficiently large that the results for the larger events should bereasonably robust.The mass balance equation assumed that all snowpack losses and gains wereassociated with rainfall, snowfall, outflow, and melt. This assumption is not alwaystrue; however, the specific conditions of rain-on-snow make this assumption moreappropriate. Evaporation of liquid water from the snowpack is unlikely during rain-on-snow, as the air is at or near 100% humidity during these situations, and the airtemperature is above 0◦C. In a mass and energy balance modelling study duringrain-on-snow, van Heeswijk et al. (1996) found that losses (or gains) in the snowmass balance due to evaporation and condensation were negligible. Additionally,sublimation, snow deposition and scour should be minimal during rain-on-snow,when the snow is warm and wet. Li and Pomeroy (1997) found that snow transportby wind was rare for snow that had experienced air temperatures greater than 5◦C,which was common for the events in this analysis.The performance of the snow accumulation and melt model introduces furtheruncertainty. Though the purpose of the snowpack model was to provide rudimen-tary tracking of the internal state of the snowpack, not to predict SWE, accuratesimulation is usually an indication of appropriate representation of internal con-ditions. The generally poor performance of the model in prediction of the annualcycle of SWE is thought to be due primarily to: (1) the use of a significant amountof infilled driving data; and (2) the use of precipitation gauges without windscreensat all sites, which can lead to an under-catch of 50% or more during snowfall (Ding-man, 2015).To investigate the sensitivity of the snow model to these issues, a test wasperformed to determine whether the internal state estimates (liquid water content,cold content, and snow surface temperature) of the Jost et al. (2012) model arerobust to issues with the driving data. Precipitation, during time steps where themajority of precipitation was thought to be snow, was adjusted at each site by theratio of annual peak simulated to modelled SWE for each year. The model wasre-run using the adjusted data. If the peak observed SWE was not available fora particular year, the average ratio for all available years was used. Re-analysis92using these modified model outputs (not shown) indicated an improvement in thesimulation of peak SWE but little change to the prediction of pre-event saturationand cold content.Though the turbulent exchanges of sensible and latent heat have been recog-nized as a significant contributor to snowmelt during rain-on-snow, direct compu-tation requires values of humidity and wind speed, in addition to air temperature,which are not readily available at automated snow pillow sites (or in most oper-ational monitoring networks). Wind speed has been shown to be the dominantfactor in turbulent heat exchange (Marks et al., 1998). As shown in Chapter 3,NARR output provides reasonable estimates of air temperature and vapour pres-sure, but (at best) systematically overestimates wind speed, and therefore cannotbe used directly to compute turbulent heat exchange. Furthermore, wind speedsfrom NARR are likely to be even more poorly representative of wind speeds atthe more sheltered ASP sites than the ridgetop locations examined in Chapter 3.The inaccuracy of NARR downscaled wind speed meant that it was not possibleto accurately estimate the turbulent energy exchange at the ASP sites. Whereas itmay be assumed that the very large rain events were usually associated with highwinds (since they were also often associated with atmospheric rivers), a lack ofwind speed observation was likely a factor in the relative inability to accuratelypredict WAR generation for smaller events with less than 40 mm of rainfall.4.5 ConclusionsThe primary conclusions that can be made from this analysis are (1) the totalamount of rainfall is the most important control on WAR generation for largeevents; (2) for large rainfall events (> 40 mm), WAR generation is consistentlyenhanced by approximately 25 to 30% over rainfall alone; (3) though it will beoverwhelmed during large events, the snowpack water content is likely the mostuseful internal snowpack characteristic for anticipating the potential snowmelt re-sponse during small rain-on-snow events, particularly in maritime environmentswhere a substantial snowpack cold content is uncommon; and (4) late season rain-on-snow results in a larger contribution of snowmelt to WAR production than formid-winter due to its occurrence on ripe, saturated snowpacks, and the generally93higher air temperatures that would provide greater energy input.The present study has focused on rain-on-snow response at snow pillow sites,which are typically located in small openings in open forest sites near or belowtreeline. However, these results may not represent conditions at open sites, whichare more exposed to wind, or sites with greater canopy cover, which would beless exposed to wind. In addition to effects of forest canopy, elevation can playa major role in controlling rain-on-snow occurrence and WAR generation. Fur-ther work, including long term, detailed tracking of meteorological conditions, andenergy and mass balance of the snowpacks at multiple sites with a range of differ-ent elevations and characteristics are necessary to continue towards understandingsnowmelt during both large and small rain-on-snow events.94Chapter 5Probabilistic estimation of snowcover from MODIS imagery5.1 IntroductionAlong with better knowledge of point scale processes involved in rain-on-snowphenomena, an ability to observe the spatial dynamics of snow is a major steptowards increasing knowledge of the processes governing snow dynamics in hy-brid rain-snow basins (McCabe et al., 2007). As described in Section 1.1.3, anultimate goal of remote sensing of snow for water resources is to measure SWE;however, this has not yet fully been achieved. Therefore, optical remote sensing ofSCA remains an important component for observing large scale spatial dynamicsof snow. Optical remote sensing of snow has been accomplished using a range ofsatellites and at a range of temporal and spatial resolutions. Geostationary observa-tions of snow from the Geostationary Operational Environmental Satellite (GOES)satellites provide near-continuous temporal resolution, but at a 4-km spatial reso-lution (for the GOES-West satellite, which images western North America) 1. Atthe other end of the spectrum, images from the Landsat 8 operational light im-ager are 30-m resolution, but have a 16-day return period, limiting their ability tocapture snow dynamics, particularly in commonly cloud obscured areas such as1 south coast region of BC. The Moderate Resolution Imaging Spectroradiome-ter (MODIS) instrument, onboard NASA’s Aqua and Terra satellites since 2002,is currently a primary operational source for satellite SCA observation due to itscombination of high temporal resolution (twice daily) and relatively high spatialresolution of 500 m (Hall et al., 2002). The normalized difference snow index(NDSI), detailed in section 2.2.2, is the primary basis for operational snow delin-eation with the MODIS instrumentation. The MODIS binary snow classificationhas been found to correctly classify a pixel as either snow or non-snow covered(when compared with ground observations) approximately 93% of the time (Halland Riggs, 2007); however, specific shortcomings in the challenging environmentof mountainous maritime regions have limited the utility for improved understand-ing of the spatial dynamics of snow in the rain-on-snow environment.Errors with binary, threshold based classification of snow from satellite obser-vations are greatest under heavy tree cover, and at times when the snowpack isthin and rapidly changing (Hall and Riggs, 2007). During accumulation (ablation)periods, snow delineation can switch rapidly from under- to over- (over- to under-) representation as pixels approach and cross the delineation threshold (G. Jost,Senior Hydrologic Modeller, BC Hydro, pers. comm, 2015).One option to address this shortcoming is by estimating the fractional snowcover to account for sub-pixel variability within 500 m MODIS pixels. The mostsuccessful of these methods are based on the use of spectral unmixing modelswhich use pre-defined spectral signatures of unmixed pixels to estimate the fractionof a pixel that is snow covered (Painter et al., 2009; Rittger et al., 2013).A relatively unexplored alternative to estimating the fraction of snow coverwithin a pixel is to make probabilistic estimates of snow cover. In one example,Musial et al. (2014), using the binary MODIS snow cover products as the ’groundtruth’ for a training algorithm, created a multi-tiered probabilistic classificationsystem based on artificial neural networks to classify output from the AdvancedVery High Resolution Radiometer (AVHRR) between cloudy, snow and clear skyconditions. However, as described previously, the binary MODIS snow cover prod-uct has uncertainties in delineating snow in mountainous maritime regions; hence,higher resolution training data may be more appropriate. Additionally, the primaryfocus of the Musial et al. (2014) study was differentiating between snow and clouds96instead of snow and bare land.The difficulties of delineating snow under heavy forest cover affect all opticalremote sensing methods due to direct blockage of the ground by tree cover. Thisissue is dealt with differently, depending on the method. NDSI-based methodsmake assumptions about snow cover underneath a canopy based on adjustments tothe classification threshold (Klein et al., 1998), whereas unmixing methods makeno assumptions and directly represent the visible amount of snow within a pixel(Rittger et al., 2013).Accuracy of snow-cover discrimination varies among times of observation andspatially within an image (Hall and Riggs, 2007; Parajka and Blo¨schl, 2006; Rittgeret al., 2013), but statements of classification accuracy are typically summary val-ues for a specific method and do not account for spatio-temporal variability. Whileobtaining more data via satellite snow observations is an important goal, the dataare only useful if the information obtained from these observations is accurate or,at least, are associated with realistic assessments of their uncertainties. It should berecognized that satellite observations are not direct observations, and hence have ahigher level of uncertainty associated with them than more typical hydrologic ob-servations. As Beven and Westerberg (2011) cautioned, taking observations withpotentially high levels of uncertainty as fact increases the risk of using disinforma-tion to inform hydrologic analysis.Cloud obscuration is a primary issue in snow cover delineation using opticalsensors (Dietz et al., 2012; Parajka and Blo¨schl, 2006). A range of algorithmsand methods have been developed to estimate snow cover under cloud obscura-tion, commonly following a tiered approach, whereby increasing amounts of cloudremoval are traded for decreasing levels of accuracy (e.g. Gafurov and Ba´rdossy,2009; Parajka and Blo¨schl, 2008). As with direct observation of snow cover viasatellite, it is rare for estimates of uncertainty to be included with cloud infillingmethods. In one exception to this, Lo´pez-Burgos et al. (2013) created a probability-of-snow algorithm for infilling cloud obscuration, based on locally weighted logis-tic regression in central Arizona. Though it was tested in a region with less cloudobscuration than the south coast of British Columbia, and verified with point data(as opposed to other pixel-scale data), Lo´pez-Burgos et al. (2013) found that thismethod was effective in reducing cloud obscuration, but appeared biased towards97over-predicting snow cover.The purpose of the current analysis was to develop and test a procedure forincreasing the amount of information conveyed about snow dynamics from thesatellite snow observations provided by the MODIS instruments, and to provide anassessment of confidence in the observations in order to maximize the transfer ofuseful information. This study had two specific objectives: (1) to develop a prob-abilistic snow cover classification model that improves accuracy and increases theinformation conveyed about snow cover in mountainous maritime regions such asthe BC south coast when compared to stock MODIS binary snow cover classifi-cation methods; and (2) to develop and test a probabilistic approach to infilling ofcloud-obscured pixels that maximizes information gain while avoiding the incor-poration of disinformation from erroneous predictions.5.2 MethodsThe methods for maximizing information gain from MODIS snow observationscan be split into two steps. First is the application to visible pixels. In this step,stock MODIS snow products and individual MODIS band values were convertedinto estimates of the probability that a pixel contained ≥ 50% snow cover (P50).Second is the estimation of P50 for cloud-obscured pixels. In this step, P50 esti-mations from the surrounding pixels, and from the same location on previous days(both created using methods described in step 1), are used to estimate P50 for thepixel in question. The summary flow chart in Figure 5.1 outlines the approach usedto assign a P50 value to each pixel. Individual steps will be described in detail inthe following sections.5.2.1 Data sourcesTwo types of data from the MODIS satellites (both Aqua and Terra) were used todevelop a probabilistic model for estimating P50. Only MODIS pixels that wereclassified as either snow covered or land covered by either Aqua or Terra wereused (i.e. no clouds, open water, ocean, etc.). Predictor variables for the prob-abilistic model (described in paragraph 2.2.2) were NDSI values calculated fromthe MOD09GA and MYD09GA individual MODIS bands.98Figure 5.1: Flow chart describing the calculation of P50 at each pixel.99Training data (predictands) were calculated from high resolution snow coverdata in the south coast region. As described in Chapter 2, high resolution snowcover data from two sources were available. The first was snow cover classifiedfrom 66 Landsat 8 tiles (representing 36 individual days between March 2013 andFebruary 2015) at 30 m resolution (paragraph 2.2.2). The second was data from 10days of ground based georeferenced photos, at an equivalent resolution of 25 m,within the Cheakamus Daisy catchment (Figure 2.1) between July 2012 and July2014 (paragraph 2.2.2).Data from both high resolution sources were resampled to the equivalent 500m resolution of the MODIS predictor data to produce a binary snow cover classi-fication. If a 500 m pixel contained (did not contain) > 50% snow cover, it wasclassified as snow covered (non-snow covered). A MODIS pixel was only usedwith the ground based photo-georeferencing data set if it had at least 50% visibil-ity of the entire 500 m pixel. The combined data provided 2,478,588 total visibleMODIS pixels (from at least one of Aqua or Terra) with a high resolution source.Of these pixels, the high resolution source was from Landsat 8 for 2,473,328 ofthe pixels and the georeferenced photos for 5260 pixels. Terra had 2.37 millionvisible MODIS pixels and Aqua had 2.33 million visible pixels. Aqua and Terraoverlapped in viewing 2.22 million pixels on the same day.5.2.2 Uncertainty model for unobscured pixelsNaive Bayes classificationNaive Bayes classification (Zhang, 2005) was used to probabilistically classifyMODIS data based on NDSI observations. Naive Bayes classification applies theprinciples of Bayes’ theorem to any number of predictors with the ’naive’ assump-tion of independence between predictors. This assumption is considered naivebecause in most cases the predictors are not independent; however, naive Bayesclassification has been found to perform well and the independence assumptionspeeds up and simplifies the classification algorithm (Mitchell, 1997). Since theclassification described here only used a single predictor, dependence among pre-dictors was not an issue.100For naive Bayes classification with one predictor variable, Bayes’ rule can bewritten as:P(Ck|x) = P(Ck) ·P(x|Ck)P(x) (5.1)where Ck is an event with k possible outcomes, x is a continuous variable, and P(.)is a probability. In this case Ck is the snow cover on a MODIS pixel as determinedfrom high resolution data with two possible outcomes, greater or less than 50%coverage. The MODIS NDSI is represented by x. In naive Bayes classification,P(x) is a constant, as it is independent of the outcome, and values of the predictorx will be known.A naive Bayes classifier model using continuous data is fit by computing thedistribution parameters of the predictor variable for all k classes. In this case,which is the most common assumption, a normal distribution of NDSI values wasassumed for all classes (Mitchell, 1997). The naive Bayes classifier model theninvolves estimating the mean and standard deviation for each class and the priorprobabilities of each outcome (P(Ck)).To apply the classifier model, the probability density is calculated for each out-come given an NDSI value and multiplied by the prior probability of each outcometo calculate probability scores. The probability for an individual category thenequals the score divided by the sum of all scores.The naive Bayes algorithm from the ’e1071’ package (Meyer et al., 2015) inthe R computing language was used to train a naive Bayes classifier separately forAqua and Terra. This algorithm performs calculations in log space for numericalstability (Meyer et al., 2015). The resulting probability (between 0 and 1) repre-sents the probability that a pixel has > 50% snow cover (P50).Application of the naive Bayes classifierThe full data set of ∼2.5 million pixels was randomly split into training (80%)and testing (20%) data sets, based on sampling from a uniform distribution. NaiveBayes models were fit independently for the Aqua and Terra training data on thepixels on the 1.9 (1.8) million pixels with observations from Terra (Aqua). Directaccuracy of binary classification via naive Bayes was compared to the accuracy101of the MOD10A1 and MYD10A1 binary snow classification products. For thiscomparison, if a pixel had a P50 > 0.5, it was classified as snow covered, andif P50 was ≤ 0.5, it was classified as non-snow covered. Further, probabilisticperformance was assessed by comparing the predicted P50, via the naive Bayesmodel, to the actual percent of pixels with > 50% snow cover (Pobs, based on thehigh resolution observations) in each of four binned P50 groups (0 to≤ 0.25, >0.25to ≤ 0.5, >0.5 to ≤ 0.75, >0.75 to ≤ 1) as follows:Pobs =NsNs+Nn(5.2)where Ns and Nn are the actual number of test 500 m pixels within the P50 bin tocontain and not contain > 50% snow cover, based on the high resolution obser-vations, respectively. Performance of the naive Bayes classification methods wastested individually for the Aqua and Terra models and for the mean of the predictedAqua and Terra probabilities when both observations were available.For use in the second part of the uncertainty model (predictions in unobservedpixels, Figure 5.1) and in Chapter 6, naive Bayes models were fit for both Aquaand Terra using all available data. When pixels had unobscured observations fromboth Aqua and Terra, the probabilities were averaged.5.2.3 Uncertainty model for cloud-obscured pixelsAs a first step towards reducing loss of information due to cloudcover, Aqua andTerra snow observations (produced in the previous steps) were combined into asingle image, as in Parajka and Blo¨schl (2008) and Gafurov and Ba´rdossy (2009).If a pixel was mapped as cloud-obscured by Terra, but not by Aqua (or vice versa),the unobscured pixel was used. If both pixels were obscured, the probability ofsnow coverage in cloud-obscured pixels was estimated first on the basis of visiblepixels in the current image and second on the basis of information from previousdays.Stratified neighbourhood selectionSpatial information was incorporated through the use a stratified moving windowfilter centred on each obscured pixel. A window filtering method assumes that the102MODIS snow pixels (in this case the combined Aqua/Terra probabilistic estimatesof a pixel having > 50% snow cover from the previous section) are most similar tothe nearby pixels, referred to as the pixel’s neighbourhood. In its most simple form,a ’window’ for an obscured pixel defines the neighbourhood as all of the visiblepixels within a certain distance of the obscured pixel. In this application, the dis-tance based neighbourhood was extended to include elevation and landcover (fromthe Globcover 2009 dataset, simplified into five categories, see paragraph 2.2.2) inorder to maximize the homogeneity of the neighbourhood.After simplifying to only pixels with the same landcover classification, theneighborhood of a given obscured pixel is defined as all cells within a distancetolerance (dtol , km) in both the x and y directions of the target cell and lying withinan elevation band defined as zi +/- ztol , where zi is the elevation of the target cell (m)and ztol is a defined elevation tolerance (m). Once the neighbourhood was defined,the P50 of the target cell was defined as the mean P50 of all visible cells within theneighbourhood.Smaller dtol and ztol values should provide more accurate predictions, whilelarger dtol and ztol values should provide for more total predictions. The dtol andztol were selected manually, with a goal of maximizing accuracy while includinga large enough neighbourhood that a decision could be produced even with signif-icant cloud cover. This classification was only applied if n ≥ 50. If an obscuredpixel was classified in the ’Permanent Snow/Ice’ category by the Globcover 2009dataset, it was assumed to have P50 = 0.99.Information from previous daysWhen n < 50, the most recent P50 at the obscured pixel location (up to seven daysprior) was extracted and exponentially decayed towards a probability of 0.5 asfollows:N(t) = (N0−0.5) ·2−tt0.5 +0.5 (5.3)where t is the number of days prior to the infilling image that the observation oc-cured, t0.5 is the half life, N0 was the observed P50, and N(t) was the decayedprobability. The calculated N(t) was then assigned as P50 for the obscured pixel.103Exponential decay was applied to decrease confidence in the observation (by push-ing it closer and closer to 0.5) as the observation became older. A half life of fourdays was chosen in order to remove the majority of the confidence in the previousobservation by the time eight days had passed since an observation. Eight daysis the return time of the MODIS satellite from which it observes a pixel from anidentical zenith angle. When a pixel did not have any observations in the previouseight days, or more than 50 members in its neighbourhood, or a classification aspermanent snow/ice by the Globcover 2009 dataset, P50 was set to 0.5.Benchmark modelThe performance of this cloud infilling method was compared with the well knownSNOWL filter from Parajka et al. (2010), which uses the median elevation of snowcovered and non-snow covered pixels within a region to estimate snow cover inobscured pixels. If an obscured pixel was above the median elevation of snowcovered pixels, it was classified as snow covered, and if a pixel was below the me-dian elevation of non-snow covered pixels, it was classified as non-snow covered.Pixels in between these two elevations were classified as ’partially snow covered’,though this may be equally well considered ’maybe snow covered’ (similar to apixel having a P50 value near 0.5). The SNOWL filter was not developed for use onprobabilistic observations; hence it was tested on the binary MODIS Aqua/Terracomposite binary data and had to be interpreted slightly differently than the initialintent from Parajka et al. (2010). The ’partial snow cover’ category was interpretedas a no-decision or low confidence pixel and the snow-covered and non-snow cov-ered were interpreted high confidence classifications for better comparison with theprobabilistic methods.Performance AssessmentPerformance of the probabilistic cloud infilling and SNOWL methods was assessedby applying cloud masks from different days to an image from a relatively clearday, and attempting to predict the synthetically obscured pixels. Random samplesof 1000 synthetically obscured points were taken on 12 different days. Four dif-ferent cloud masks, with four different levels of obscuration, were applied to the10412 test days. The probabilistic and SNOWL methods were used to predict valuesin obscured pixels. Since the two methods do not predict the same quantity, directcomparison between probabilistic and SNOWL methods was not possible. How-ever, it was possible to compare the percent of total high confidence predictionsthat were successful. For the probabilistic model, success was defined as occurringwhen the infilling model predicted P50 > 0.75 (or less than 0.25) for pixels thatactually had P50 > 0.75 (or less than 0.25) based on the pixel’s NDSI value andthe naive Bayes classifier. For the SNOWL filter, success was defined as occur-ring when the filter classified a pixel as snow (or non-snow) covered when it wasactually snow (or non-snow) covered in the Aqua/Terra composite image.Although making correct predictions is the key objective, increasing the num-ber of high confidence predictions (provided they are correct) would contributemore total information. Therefore, the total number of high confidence predictionsproduced by both models was also shown for each cloud mask/test day combina-tion.5.3 Results5.3.1 Probabilistic classification of observed pixelsFor the 474,000 pixels in the test data set for the Terra data, the standard bi-nary classification of the MOD10A1 data product provided the correct classifi-cation 90% of the time. Fitting of the naive Bayes classification model with theMOD10A1 and Terra NDSI improved this classification accuracy to 95%. Theprobabilistic prediction results for the Terra naive Bayes model are shown in Ta-ble 5.1.For the 466,000 pixels in the test data set for the Aqua data, the binary clas-sification of the MYD10A1 data product provided the correct classification 87%of the time. Fitting of the naive Bayes classification model with the MYD10A1and Aqua NDSI improved this classification accuracy to 94%. The probabilisticprediction results for the Aqua naive Bayes model are shown in Table 5.2.There were 445,000 pixels in the test data set that had observations from bothAqua and Terra. For these pixels, the Terra (MOD10A1) binary snow product105Table 5.1: Binned results for test data from the Terra naive Bayes classifica-tion model to predict the probability that a pixel has > 50% snow cover.Pobs is calculated as in eq. 5.2.Predicted probability Pobs Total # of pixels0 to ≤ 0.25 0.016 397889>0.25 to ≤0.50 0.342 7854>0.50 to ≤0.75 0.447 7925>0.75 to ≤1.0 0.841 61649Table 5.2: Binned results for test data from the Aqua naive Bayes classifica-tion model to predict the probability that a pixel has > 50% snow cover.”Actual probability” represents the fraction of the total group that had anobserved snow cover > 50%.Predicted probability Actual probability Total # of pixels0 to ≤ 0.25 0.019 391127>0.25 to ≤0.50 0.317 11132>0.50 to ≤0.75 0.459 11744>0.75 to ≤1.0 0.826 52137correctly classified whether a pixel had > 50% snow cover 90% of the time, andthe Aqua (MYD10A1) correctly classified whether a pixel had > 50% snow cover87% of the time. The mean of the two probabilities (re-classified as binary data)provided the correct classification 96% of the time. The probabilistic predictionresults for the combined Aqua/Terra model are shown in Table 5.3.Table 5.3: Binned results for test data for the combined Aqua and Terra naiveBayes classification models to predict the probability that a pixel has >50% snow cover. Actual probability represents the fraction of the totalgroup that had an observed snow cover > 50%.Predicted probability Actual probability Total # of pixels0 to ≤ 0.25 0.012 375757>0.25 to ≤0.50 0.259 16394>0.50 to ≤0.75 0.527 10228>0.75 to ≤1.0 0.897 434491065.3.2 Infilling of cloud-obscured pixelsData usedThe days used for prediction are summarized in Table 5.4. The SNOWL methodonly used the data from the specific day shown in Table 5.4, whereas the proba-bilistic method could use data from up to seven days prior to the specific day. Asub-sample of the days used for prediction with the probabilistic model is shownin Figure 5.2.Table 5.4: Summary of days used for testing cloud infilling methods. %Clouds represents the cloud coverage for the visible land pixels in theimage. % Non-Snow and % Snow represent the coverage only for non-obscured land pixels.Date % Clouds % Land % Snow2009-12-27 27 57 432010-01-16 8 45 552011-12-22 0 97 32012-01-10 3 25 752012-03-06 12 73 272012-09-07 7 44 562012-11-08 15 28 722012-11-10 3 46 542013-01-21 45 25 752013-04-01 8 35 652013-05-09 12 45 552013-06-04 4 76 24Cloud masks from alternate days (used to artificially obscure visible pixelsfor algorithm testing) were taken from the MOD10A1 (Terra) binary data and areshown in Figure 5.3. The cloud coverage in each is summarized in Table 5.5.Percent cloud cover represents the percent of visible land pixels in the image.Cloud infilling resultsA dtol of 100 km and a ztol of 100 m were determined to offer the best compromisebetween accuracy of predictions and total number of predictions made. Results of1072012−11−082013−01−212013−05−095e+056e+055e+056e+055e+056e+051000000 1100000 1200000 1300000EastingNorthing0.000.250.500.751.00P50Figure 5.2: Examples of the days used for testing of the method for proba-bilistic infilling of cloud obscuration. Only data available as information(visible data) is shown.108Oct 2, 2012 March 29, 2013Oct 17, 2012 March 13, 20135e+056e+055e+056e+051000000 1100000 1200000 1300000 1000000 1100000 1200000 1300000EastingNorthingNo Snow Snow CloudFigure 5.3: Spatial patterns for cloud masks that were overlaid upon test im-ages for cloud infilling.infilling of 1000 random points using the landcover and elevation based probabilis-tic model, with four different cloud masks and three different days, are shown in asubsample of days in Figure 5.4. The predicted probabilities are grouped into fivebins, with a bin width of 0.25. The high confidence predictions of P50 (either <0.25 or > 0.75) were largely correct, even with significant cloud obscuration. Thelow confidence predictions of P50 had a significantly larger spread, but generallycorresponded well with the observed range.The accuracy of high confidence predictions of P50 (either < 0.25 or > 0.75)using both models is summarized in Figure 5.5. The total numbers of predictionsmade with high confidence (out of 1000) are summarized in Figure 5.6. The prob-109Table 5.5: Percent coverage for masks that were overlaid upon test imagesfor cloud infilling.Date % Cloud Cover2012-10-02 252013-03-29 472012-10-17 822013-03-29 99abilistic method generally outperformed the SNOWL method in accuracy of highconfidence predictions for fall or winter days (with March 6, 2012, being the ex-ception). For spring and summer days, the SNOWL and probabilistic methodsproduced similar levels of accuracy.Figure 5.6 shows that the probabilistic and SNOWL methods produced similartotal numbers of high confidence predictions for fall and winter test days for lowerlevels of cloud cover; however, when cloud coverage increased, the number ofhigh confidence predictions generally decreased more rapidly for the probabilisticmethod than the SNOWL method (particularly January 16, 2010, and January 10,2012). For spring and summer days (particularly September 7, 2012, and May 9,2013), the probabilistic method produced more high confidence predictions thanthe SNOWL method.5.4 Discussion5.4.1 Regional accuracy and adjustmentDue to the global scope of MODIS, it is understandable that a globally optimizedclassification algorithm would be used for the stock Aqua and Terra binary snowcover data. However, results from sensitivity testing of the original algorithm inHall et al. (1995) indicated that there were regional differences that could be ac-counted for if a regionally optimized product was desired. Because of the signif-icant amount of training data that was obtained through Google Earth Engine, itwas possible to fit a classification algorithm separately for Aqua and Terra that wasbased directly on NDSI values. When forced to produce a binary estimate, the1102012−11−08 2013−01−21 2013−05−09lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll0.000.250.500.751. Clouds47% Clouds82% Clouds99% Clouds0:0.25 0.25:0.5 0.5:0.75 0.75:1 0:0.25 0.25:0.5 0.5:0.75 0.75:1 0:0.25 0.25:0.5 0.5:0.75 0.75:1Predicted from cloud obscuration modelCalculated in naive Bayes modelFigure 5.4: Examples of P50 prediction results from infilling 1000 syntheti-cally obscured pixels on three days with four levels of cloud masking.trained naive Bayes models improved prediction accuracy by approximately 5-7%for both the Terra and Aqua models as compared to the stock binary products, sup-porting the hypothesis that regional improvement in classification accuracy couldbe found. The classification accuracy of the binary data product was slightly higherfor Terra data than for Aqua, which may be expected due to Aqua’s modificationof the NDSI algorithm due to band failure (see paragraph 2.2.2) and lack of anadjustment to classification thresholds via the NDVI (Hall and Riggs, 2007). Thus,the snow cover delineation from Aqua has typically been considered as the secondchoice data source. The fitting of the naive Bayes probabilistic model to the Aquadata resulted in accuracy that was very similar to the Terra model.1112009−12−27 2010−01−162011−12−22 2012−01−102012−03−06 2012−09−072012−11−08 2012−11−102013−01−21 2013−04−012013−05−09 2013−06−0402550751000255075100025507510002550751000255075100025507510025% 47% 82% 99% 25% 47% 82% 99%Simulated cloud coveragePercent correctModelProbabilisticSNOWLFigure 5.5: Accuracy of high confidence predictions of P50 (either < 0.25 or> 0.75) for the probabilistic and SNOWL models.All accuracy numbers rely on the assumption of accuracy in the higher res-olution training data that were used. As with the MODIS pixels, these higherresolution training data are still remote sensing data, and hence have uncertaintyassociated with them as well. However, Hall and Riggs (2007) considered Landsatdata of high enough resolution to be used for ’absolute validation’ of MODIS data,and Landsat (from a range of different satellites in the lineage) has been used asthe test source for validation of MODIS snow classifications in many studies (Hall1122009−12−27 2010−01−162011−12−22 2012−01−102012−03−06 2012−09−072012−11−08 2012−11−102013−01−21 2013−04−012013−05−09 2013−06−0402505007501000025050075010000250500750100002505007501000025050075010000250500750100025% 47% 82% 99% 25% 47% 82% 99%Simulated cloud coverage# of high confidence decisionsModelProbabilisticSNOWLFigure 5.6: Total number of high confidence predictions made using the prob-abilistic and SNOWL al., 1995; Klein et al., 1998; Painter et al., 2009; Rittger et al., 2013; Salomonsonand Appel, 2004). The issue of mixed landcover pixels should be less importantat 30-m resolution than at 500-m, and the higher resolution can identify clearingareas within a 500-m pixel. However, issues of snow delineation under closed for-est canopies are still present at the finer scale resolution, and likely will remain forany snow cover product based on optical remote sensing imagery. At best, opti-cal remote sensing can only view the projected snow cover that is visible through113gaps in the forest canopy, and thus assumptions are required regarding snow coverunderneath the canopy (Rittger et al., 2013).5.4.2 Uncertainty in observed pixelsWhile an improvement in classification accuracy is welcome, the primary purposeof using the naive Bayes classifier was to produce a probabilistic prediction of >50% snow in a pixel (P50). As shown in Tables 5.1 and 5.2, the majority of pixelswere classified with high levels of confidence (with P50 either near 0 or near 1).The P50 output is primarily beneficial when the pixel in question comes close tothe binary switch between snow and non-snow. Binary data, even if correct thevast majority of the time, do not indicate when the classification is becoming lesscertain. The use of a binary threshold can lead to abrupt switches in delineatedsnow cover, which may not represent reality.While the shortcomings of binary classification of 500 m pixels have been rec-ognized previously (Dietz et al., 2012), most work to address this has focused onmaking sub-pixel estimations of fractional snow cover (e.g. Moosavi et al., 2014;Painter et al., 2009; Salomonson and Appel, 2004; Vikhamar and Solberg, 2003)instead of probabilistic estimates of snow cover. Fractional estimates can increaseaccuracy and smooth transitions during accumulation and ablation, but in thesemodels, data are still typically presented without associated uncertainty.The results in Tables 5.1, 5.2, and 5.3 indicate that, for the most part, thesesimple probabilistic estimates were similar to the actual probability that a pixel ineach bin contained 50% or more snow cover (based on the high resolution data).The probabilities became somewhat less successful the closer the estimate was to0.5, with the actual probabilities falling on the low end (or outside of) the predictedrange. While the exact probabilities in this middle, low confidence region are notpredicted as accurately as at the upper and lower ends, they do achieve the purposeof indicating when there is low confidence in the classification and when thereis not. Additionally, probabilistic classification of snow cover lends itself welltowards probabilistic estimates of snow under cloud cover, which benefits evenmore greatly from having a summary of the confidence in an estimation.1145.4.3 Infilling of cloud-obscured pixelsCloud obscuration is a significant factor in the utility of MODIS snow cover datafor assessing snow dynamics before, during, and after rain-on-snow in maritimeregions such as southwestern British Columbia, and hence the loss of informationmust be reduced as much as possible. Many common infilling methods produceestimates that are presented as actual observed data, and do not convey the addi-tional uncertainty that may be associated with the infilling method. The goal ofthe method developed here was to gain as much useful information about snowcover without adding disinformation. Combining Aqua and Terra (which observethe same location at different times of day) observations has been shown as anobvious first step towards reducing cloud obscuration, as accuracy should not besignificantly sacrificed (Gafurov and Ba´rdossy, 2009; Parajka and Blo¨schl, 2008).Parajka and Blo¨schl (2008) and Gafurov and Ba´rdossy (2009) both followed thiscompositing with tiered levels of infilling, with a trade off between the number ofpixels reduced and the accuracy of the classification. Since the accuracy of eachtiered level is estimated as a whole, equal accuracy for each infilled pixel is im-plied. The SNOWL method (Parajka et al., 2010) is one method that does includea level of uncertainty, in that some pixels are essentially no-decision pixels, eventhough they are referred to as ’partial snow’ pixels.The method described here instead employs probabilistic estimations of P50.When the evidence of the surrounding pixels is conclusive enough, the estimationwill be presented at similar levels of confidence as an observed pixel, whereas if theevidence is unclear, the probabilistic estimate approaches 0.5, as would be desired.This method also employs temporal information when the evidence from the pixelneighbourhood is minimal (a neighbourhood of less than 50 members). In thisanalysis, a larger regional scale window (100 km x 100 km) was used to determinethe obscured pixel’s neighbourhood. This approach has allowed for a significantamount of infilling, even in situations with 99% cloud cover over the full region.In comparison with the SNOWL filter, the performance of this moving windowmethod appears to be quite strong and more stable, with consistent accuracy levelsfor the entire year and all levels of obscuration. The SNOWL method outperformedthe window method by a small margin for the spring and summer test days. It is115suspected that this is because spring conditions are associated with a more definedsnowline, which is exactly what the SNOWL method relies upon. Nevertheless,accuracy of the probabilistic method during spring and summer is quite strong,and produces more high confidence estimates than the SNOWL method in thoseseasons. A significant advantage of the probabilistic method over the SNOWLmethod is the relatively stable accuracy levels as cloud coverage increases. Thisis due to the inherent use of evidence to produce the P50 estimate; if evidencefrom the neighbourhood is inconclusive (a wide range in P50 values within theneighbourhood or too small of a neighbourhood), then the estimated P50 will benear 0.5.For lower confidence predictions (between 0.25 and 0.75), predictions of theprobabilistic model tended to be near the probabilities from observed data (Fig-ure 5.4), especially at lower levels of cloud obscuration. As cloud obscurationincreased, the cloud infilling model tended towards conservative predictions, withthe infilling model making more low confidence predictions; however, a conserva-tive infilling model is preferable to a model that makes incorrect high confidencepredictions.This approach for classifying and infilling MODIS snow cover imagery couldbe useful for data assimilation to improve accuracy for hydrologic forecasting. TheEnsemble Kalman Filter, a commonly employed data assimilation technique, usesMonte-Carlo samples of model states to propagate errors throughout the system(Evensen, 2003) and hence allows for uncertainty in both observations and output.The Ensemble Kalman Filter has been used for assimilation of remotely sensedsnow cover into hydrologic forecasting models (Andreadis and Lettenmaier, 2006;Slater and Clark, 2006). Observations that implicitly include uncertainty, such asthe output produced in this method, could help to avoid over-correction in hydro-logic forecasts. Probabilistic maps of SCA could also be used to directly adjust theinternal snow state of a hydrologic model via a Monte Carlo sampling approach.Alternate model snow states could be sampled from the probabilistic estimates de-scribed in this chapter to represent potential ’true’ watershed snow states.In a less formalized approach, probabilistic maps of SCA can be used for qual-itative guidance to hydrologic forecasters, who improve forecasts using manualadjustments. In situations where confidence is high in snow observations, this116qualitative information can be taken as stronger evidence.Along with being used to improve discharge estimates (e.g. Finger et al., 2011;Franz and Karsten, 2013; Immerzeel et al., 2009; Parajka and Blo¨schl, 2008), satel-lite observations of snow cover have been used to specifically improve representa-tions of distributed snow parameters in land surface models (e.g Kuchment et al.,2010; Liston and Hiemstra, 2008). Finger et al. (2011) used spatial maps of snow-cover to assist in choosing parameter sets generated via Monte Carlo simulationin a hydrologic model in order to help minimize equifinality. Use of probabilisticestimates for these methods may help in ensuring that a model is not overfit todisinformation.5.5 ConclusionsThe primary conclusions from this analysis are: (1) accuracy of the MODIS binarysnow cover in southwest BC can be improved upon through a simple probabilisticclassification algorithm; (2) information from the surrounding (and previous days’)pixels can be used to infill obscured pixels with relatively high accuracy; and (3)the use of probabilistic estimates (both from observations and infilling of cloud ob-scuration) provides information on the confidence of predictions that are useful forlimiting the transfer of disinformation that can occur when estimates are presentedas true data.Future work should focus on the problem of detection of snow under dense for-est canopies. As described in Section 5.1, this issue has been dealt with differentlydepending on the snow delineation method. Some methods implicitly assume thatsnow exists under tree cover within mixed pixels (e.g Klein et al., 1998), which maycause over-representation in certain situations, and some directly represent only thevisible ’plan view’ snow cover within a pixel (e.g Painter et al., 2009), which mayunder-represent snow cover in certain situations. A probabilistic approach may bean appropriate way to deal with these situations; however, as always, appropriatedata to test against will be an issue, and optical sensors will not likely be adequate.117Chapter 6The influence of antecedent snowand high elevation rainfall onrunoff generation in atmosphericrivers6.1 IntroductionRain-on-snow events are considered one of the five main causative mechanismsof floods (Merz and Blo¨schl, 2003). In the western United States, rain-on-snowevents have occurred at nearly all sites that receive some winter precipitation assnowfall, and occur most notably in the Pacific Northwest (McCabe et al., 2007).The occurrence of rain-on-snow is increasing in central Europe (Freudiger et al.,2014); however, it may actually be decreasing in western North America due to lessconsistent snow accumulation (McCabe et al., 2007). Rain-on-snow events are alsoa significant source of flooding in eastern North America; Pradhanang et al. (2013)found that the majority of peak discharges occurring within the study period of2003-2013 in New York were due to rain-on-snow.As described in Section 1.1.5, a major source of mid-winter rain events, andby extension rain-on-snow floods, in southwest British Columbia and other mar-118itime regions, are atmospheric river events. Though atmospheric rivers have beenattributed as the cause of many flood events in the Pacific Northwest and BritishColumbia (e.g. Marks et al., 1998; Neiman et al., 2011; Spry et al., 2014), andare likely to be associated with rain-on-snow at some elevations due simply to thepersistent snow cover at elevations above 1000 m in the region, the importance ofsnow cover in flood generation during these events is unclear.Antecedent snow can serve to either damp or enhance runoff response, andits impact on runoff generation depends in part on the internal conditions of thesnowpack (Garvelmann et al., 2015; Kroczynski, 2004). Though McCabe et al.(2007) stated that the spatial snow coverage of a watershed is an important factorin understanding the potential response to rain-on-snow, the spatial contributionof snowmelt to stream discharge is not well understood at the watershed scale, andhas been primarily inferred from plot scale studies, which Perkins and Jones (2008)noted must be extrapolated with caution. Marks et al. (1998) found that the pres-ence of snow at low elevations from a prior storm was a major contributor to floodgeneration in a large rain-on-snow event in Oregon, USA. Garvelmann et al. (2015)found that the contribution to runoff was the greatest from upper elevations, whererainfall was highest and snowmelt was the largest in a mountainous catchment insouthwest Germany. Jones and Perkins (2010) found that the largest rain-on-snowfloods require synchronized snowmelt over a large percentage of the watershed toproduce their resultant peak flows, which Jennings and Jones (2015) noted is likelyenhanced further by the propagation of pressure waves of precipitation through thesnowpack to the river.However, as described in Section 1.1.5, atmospheric rivers have caused sig-nificant floods across the entire west coast of North America, in regions withoutsignificant antecedent snow on the ground, such as the Russian River of north-ern California (Neiman et al., 2008a), and lower elevation watersheds of westernWashington (Neiman et al., 2011). Neiman et al. (2008a) considered snowmelta secondary factor in flood generation during atmospheric rivers. In an analysisof 200+ ASP sites in Oregon, Washington, and California, Neiman et al. (2008b)found that atmospheric river events generally caused increases in SWE during fal-l/winter, and decreases in SWE during spring. Increases in SWE were more likelyin the Sierra Nevada mountains than in the Pacific Northwest, illustrating the com-119bined danger (due to flooding) and importance (as a water resource) of atmosphericrivers in the state of California. In Chapter 4 it was found that generation of WARduring rain-on-snow events associated with atmospheric rivers (which constitutedsome of the largest rain-on-snow events in the analysis) was typically enhanced bysnowmelt, indicating a loss of SWE during the events. However, to the author’sknowledge, a spatial analysis of the loss of SCA during atmospheric rivers, and thepotential impact on flood generation, has not yet been undertaken.During the most severe atmospheric river events along the west coast of NorthAmerica, often called Pineapple Express storms (Spry et al., 2014), the high levelsof rain falling over a wide range of elevations can be a major source of runoffgeneration, independent of whether or not snow was previously on the ground. InChapter 3 it was found that a better representation of freezing levels (which inturn provided a more realistic representation of the spatial extent of rain) duringan atmospheric river event helped to improve runoff predictions in an operationalhydrologic model.This chapter builds upon analysis from previous chapters to assess the signif-icance of atmospheric rivers in generating floods in southwest British Columbia,and to help determine the role of snow cover and high elevation rainfall in runoffgeneration during atmospheric rivers. The primary questions to be addressed are:(1) what is the significance of atmospheric rivers to flood generation in southwestBritish Columbia; (2) what role does high elevation rainfall play in runoff genera-tion during atmospheric rivers; and (3) does snow cover play a significant role inthe generation of floods during atmospheric rivers?6.2 Methods6.2.1 Data sourcesThis chapter is based on discharge measurements from six watersheds in the southcoast region: Cheakamus River at Daisy Lake, Cheakamus River at Millar Creek,Elaho River, Chemainus River, Cruickshank River, and Coquitlam River (Fig-ure 2.1). The watersheds are described in further detail, including elevation ranges,landcover, and hydrologic regime, in Section 2.2. Mean daily discharge data, ob-120tained from the Water Survey of Canada and BC Hydro, were used in all analyses.In order to assess the relationship between discharge and high elevation precip-itation, each study watershed was associated with the nearest ASP, as shown inTable 6.1. Precipitation, partitioned into rainfall and snowfall as in Chapter 4,from each of these ASP sites was used.Table 6.1: The nearest ASP site to each study watershed.Watershed name ASP name ASP abbreviationCheakamus Daisy Upper Squamish River USRCheakamus Millar Upper Squamish River USRChemainus Jump Creek JUMPCoquitlam Disappointment Lake DISCruickshank Wolf River WOLFElaho Upper Squamish River USRThe list of atmospheric river events influencing British Columbia (Table 2.8)during the analysis period (where ASP and SCA data was available) of 2003-2013were used as the analysis storms for assessing individual storm characteristics.For investigating the influence of atmospheric rivers and Pineapple Expressstorms on historic annual flood generation (and only this purpose), the historicalrecord of Dettinger et al. (2011) was used to extend the atmospheric river record.Dettinger et al. (2011) identified atmospheric river events from 1998 to 2008 usingSSMI imagery methods of Neiman et al. (2008b). Though it was not possible toformally identify all landfalling atmospheric rivers without SSMI imagery (pre-1998), Dettinger et al. (2011) identified Pineapple Express storms (which are asubset of major atmospheric river events) by identifying moisture fluxes from nearHawaii to the west coast of North America using output from the NCEP weatherreanalysis model from 1948 to 2008. In the overlapping period of 1998 to 2008,Dettinger et al. (2011) found that nearly all storms identified as pineapple expressstorms in NCEP model output were also identified as atmospheric rivers in SSMIoutput.Finally, cloud-infilled snow cover, represented by the P50 value as describedin Chapter 5, was used to investigate changes in SCA within the study watershedsbefore and after atmospheric river events from Table 2.8.1216.2.2 Flood frequency analysisThe mixed Gumbel distribution from Waylen and Woo (1983) was used to charac-terize the flood frequency of the study watersheds. Waylen and Woo (1983) notedthat a single distribution was not adequate to describe the flood frequency distri-bution in watersheds where multiple processes can give rise to floods. The mixedhydrologic regimes that are present in the watersheds used in this analysis, wherefloods can be generated from winter rain or from the spring freshet are examplesof watersheds with mixed runoff generating processes. As in Waylen and Woo(1983), the parameters for the Gumbel distribution were computed separately fromthe sample of peak annual daily discharges from winter (day of year 1 to 90 and244 to 365), which were thought to be driven by rainfall and/or rain-on-snow pro-cesses, and summer (day-of-year 91 to 243), which should be driven primarily bythe spring freshet and snowmelt processes. The mixed Gumbel distribution is givenas:F(Q′ > x) = 1−(exp−exp(−αs(x−βs)))·(exp−exp(−αw(x−βw)))(6.1)where F(Q′ > x) is the exceedence probability of a discharge x, and α and β areparameters of the Gumbel distribution, which can be estimated as:α = 1.281/σ (6.2)β = µ−0.45σ (6.3)where µ is the mean and σ is the standard deviation, computed separately for theannual series of peak winter (w) and spring/summer (s) discharges. In eqs. 6.2 and6.3, α and β were estimated using the sample means and standard deviations ofthe observed flood series (i.e. by the method of moments).6.2.3 River discharge and rain-on-snow eventsIn order to assess the role of precipitation at upper elevations on flood occurrenceduring atmospheric rivers, the total rainfall and snowfall during each atmosphericriver event identified from Table 2.8 was extracted for each ASP in Table 6.1. The122peak daily discharge that occurred for each watershed, either during or one dayafter each atmospheric river event, was extracted in order to assess the relationshipbetween rainfall (at the nearest ASP site to each watershed) and discharge.6.2.4 Changes in SCA during atmospheric riversIn order to investigate the potential relationship between loss (or gain) of SCA andrunoff during atmospheric rivers, the change in snow cover before and after theatmospheric river events identified in Table 2.8 was estimated. Though the meth-ods outlined in Section 5.2.3 were shown to offer accurate predictions, even withsubstantial cloud cover, the amount of high confidence predictions decreased withincreasing cloud cover (see Figure 5.6). Hence, days with less cloud coverage weredesirable to obtain more high confidence estimates of snow coverage within water-sheds. In order to estimate the change in SCA for an atmospheric river event withhigh confidence, days with cloudcover of less than approximately 50%, prior to andafter each atmospheric river event, were identified through investigation of MODIScloudcover in Google Earth Engine. For two pairs of atmospheric river events (De-cember 4 and 10, 2004 and November 15-16 and 25, 2009) there were no cleardays available in between the two events, and these events had to be grouped to-gether. For the chosen before/after days, the MODIS Aqua and MODIS Terra datawere used to calculate the P50 for visible pixels, and obscured pixels within eachstudy watershed were infilled using the stratified neighborhood method describedin Section 5.2.3. For each study watershed, P50 pixels were classified into four cat-egories, outlined in Table 6.2, and the areal coverage of each class was calculatedfor each study watershed (shown in Figure 2.1).Table 6.2: Categorization of P50 output from the model outlined in Chapter 5.P50 range Classification0 < P50 ≤ 0.25 Strong non-snow0.25 < P50 ≤ 0.5 Weak non-snow0.5 < P50 ≤ 0.75 Weak snow0.75 < P50 ≤ 1.0 Strong snowFollowing the calculation of the areal coverage of each category, the fractionof high confidence snow cover ( fsca) was calculated for each watershed and each123day as:fsca =AsAs+An(6.4)where As and An are the areas covered by pixels predicted with high confidence tobe snow covered or non-snow covered, respectively (km2). The absolute change infsca, referred to as ∆ fsca, was then calculated for each watershed and each atmo-spheric river event. The peak stream discharge for each study watershed during, orwithin one day after, each identified rain-on-snow event (or between the consecu-tive events) was obtained to assess the relationship between ∆ fsca and discharge.6.2.5 Regional response to atmospheric riversThe similarities in regional response to atmospheric rivers were investigated at boththe watershed and point scale in the time period from 2003-2013. A correlation ma-trix of Spearman’s ρ was calculated between peak daily discharge measurementsin the atmospheric river events from Table 2.8 for the six study watersheds. Ad-ditionally, Spearman’s ρ between high altitude WAR generation (calculated as inSection 4.2.4 for each of the atmospheric river events) was calculated at the ASPall sites in Table 6.1 except the Disappointment Lake ASP, which was eliminateddue to missing data.6.2.6 Influence of rainfall and initial snow cover on dischargeIn order to assess the influence of rainfall and initial fsca on the maximum dailydischarge that occurred during each atmospheric river event in Table 2.8, linearregression models, with increasing levels of complexity, were fit and compared foreach watershed. The models followed the form:Q = b0+b1 ·R (6.5)Q = b0+b1 ·R+b2 · fsca (6.6)Q = b0+b1 ·R+b2 · fsca+b3 ·R · fsca (6.7)where Q is the peak mean daily discharge within or one day after the event (m3124s−1), R is the total rainfall during the event (mm), and fsca is the antecedent snowcover calculated via eq. 6.4 prior to the onset of the event. The maximum 24-hour rainfall during each event was also tested as a predictor, but did not providebetter results than total rainfall. For each watershed, models were assessed usinganalysis of variance (ANOVA) to determine whether inclusion of the influences ofantecedent snow cover resulted in significant improvements in predictive ability.6.3 Results6.3.1 Flood frequency analysis and atmospheric riversFlood frequency analysis results are shown in Figure 6.1. Individual annual floodsare coloured by their occurrence in either winter or spring/summer, as defined inWaylen and Woo (1983). The flood of August 30, 1991, which is the flood of recordfor Cheakamus River at Millar Creek and one of the largest floods on the Cheaka-mus River at Daisy Lake, the Elaho, Cruickshank, and Coquitlam rivers, was notincluded in the fit of the mixed Gumbel distribution (it is still shown in Figure 6.1as an individual point). This flood, which was generated via heavy rainfall in latesummer, in some cases onto exposed glacier ice (especially for Cheakamus andElaho Rivers), did not fit the separation between rain dominant events and spring/-summer freshet associated events as outlined in Waylen and Woo (1983), and hencewas considered to be drawn from a different distribution.Annual floods in the Chemainus, Coquitlam, and Cruickshank rivers weredominated by winter rain events, which were often likely associated with rain-on-snow over at least a portion of the watershed. The annual floods for CheakamusRiver at Daisy Lake, Cheakamus River at Millar Creek, and Elaho Rivers weremixed between spring/summer peak discharges and winter peak discharges. Forall watersheds, a substantial number of annual floods were determined to be asso-ciated with atmospheric rivers (or the more intense Pineapple Express subset). Thepercent of annual floods associated with atmospheric rivers is shown in Table 6.3.Cruickshank River had the highest association, with more than 50% of the annualfloods associated with atmospheric rivers. Cheakamus River at Daisy Lake had thelowest association, with 22% of annual floods associated with atmospheric rivers.125lllllllllllllllllllllllllllllllll2003004005006007008009001.0021.01 1.1111.333 2 5 10 20 40100 500Return period (yrs)Discharge (m3 /s)Cheakamus Daisyllllllllllllllllllllll4050607080901002003001.0021.01 1.1111.333 2 5 10 20 40100 500Return period (yrs)Discharge (m3 /s)Cheakamus Millarllllllllllllllllllllllllll607080901002003004005006007001.0021.01 1.1111.333 2 5 10 20 40100 500Return period (yrs)Discharge (m3 /s)Chemainusllllllllllllllllll ll4050607080901001.0021.01 1.1111.333 2 5 10 20 40100 500Return period (yrs)Discharge (m3 /s)Coquitlamllllllllll ll l4050607080901002003004005001.0021.01 1.1111.333 2 5 10 20 40100 500Return period (yrs)Discharge (m3 /s)Season l lSpring/Summer WinterCruickshankllllllllllllllllll ll30040050060070080090010001.0021.01 1.1111.333 2 5 10 2040100 500Return period (yrs)Discharge (m3 /s)lAR Non−ARElahoFigure 6.1: Flood frequency analysis results for daily discharges at study wa-tersheds for the south coast region (identified in Section 2.2.3) using themixed Gumbel distribution from Waylen and Woo (1983)6.3.2 High elevation rainfall during atmospheric riversThe rain and snowfall measured at the Jump Creek (JUMP), Upper Squamish River(USR), and Wolf River (WOLF) ASP sites (Disappointment Lake was not showndue to missing data) during atmospheric river events identified in Table 2.8 are126Table 6.3: The percentage of annual floods, as indicated by + symbols inFigure 6.1, that were associated with atmospheric rivers, as identified byBC Hydro and Dettinger et al. (2011)Name Record length (Years) % associated with ARCheakamus Daisy 55 22Cheakamus Millar 31 29Chemainus 59 46Coquitlam 31 35Cruickshank 30 57Elaho 32 38llllllllllllll0100200300400JUMP USR WOLFSiteAtmospheric River Precipitation (mm)snowrainFigure 6.2: The quantity and phase of precipitation recorded at select ASPsites during atmospheric river events identified by BC Hydro.shown in Figure 6.2. The Upper Squamish River ASP was more likely to ex-perience snowfall during an atmospheric river than the other sites, and typicallyexperienced more snowfall than rainfall during an atmospheric river.For the atmospheric river events identified in Table 2.8, rainfall from the near-est snow pillow (when available) for each of the study watersheds is shown inFigure 6.3. ASP rainfall was linearly related to discharge at the nearby watersheds127llllllllllllllllllllllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllCheakamus DaisyUSRCheakamus MillarUSRChemainusJUMPCoquitlamDISCruickshankWOLFElahoUSR02004006000501001502000100200300400500−1001020304001002003000250500750100012500 50 100 150 0 50 100 150 0 100 200 300 40020 40 60 80 0 100 200 0 50 100 150Rainfall at ASP (mm)Discharge (m3 /s) FractionFigure 6.3: The relationship between rainfall at the snow pillow sites andpeak discharge during atmospheric river events. Points are colouredby the fraction of total event precipitation that occurred as snow. Panelheaders indicate the watershed and the nearest snow pillow. Green linesindicate the 2-year flood level. A linear model with 95% confidencelimits is shown for each watershed.for all sites except the Coquitlam River, which is limited by missing data at theDisappointment Lake ASP. Atmospheric river events that included a significantamount of snowfall at the upper level sites (the ’Snow Fraction’ indicated in Fig-ure 6.3), which was most common at the Upper Squamish River ASP, were relatedto lower discharge levels at the accompanying watersheds.1286.3.3 Snow covered area during atmospheric riversFigure 6.4 shows the fsca calculated at the onset of each atmospheric river eventand the resulting discharge in the study watersheds for identified atmospheric riverevents with available discharge from 2003-2013. The watersheds for the Cheaka-mus River and the Elaho River were most commonly nearly fully snow coveredat the onset of each event each event, and had the fewest exceedences of the 2-year flood during atmospheric river events. The Chemainus River watershed wasmost commonly snow-free (or nearly snow free) at the onset of each atmosphericriver event. No strong linear relationships between antecedent fsca and dischargeexisted for any of the watersheds, though the relationship appeared to be weaklypositive for the Chemainus River, which has a greater rain component and less con-sistent snow cover, and weakly negative for the watersheds with more consistentsnow cover (Cheakamus River watersheds and the Elaho River). The Coquitlamand Cruickshank rivers had the widest range of starting fsca during the identifiedatmospheric river events, but did not display a clear relationship.6.3.4 Individual event investigationFigure 6.5 illustrates the change in snow cover before and after the Oct. 16-18,2003 atmospheric river event, which resulted in the highest peak discharge for anyof the 2003-2013 atmospheric river events in all watersheds except for the Chemai-nus River. All watersheds (including the Chemainus River) experienced peak dis-charges greater than 2-year return period during this event, with the Elaho Riverexceeding 25-year return period and both Cheakamus River watersheds exceedingthe 50-year return period. The Cheakamus River watersheds and the Elaho Riverall transitioned from having an fsca > 0.5 before the event to less than 0.5 after theevent. During this event, the Upper Squamish River snow pillow recorded 173 mmof rainfall and generated 242 mm of total WAR (with nearly complete snowmelt),which was the largest event, in terms of WAR generation, recorded at that site bymore than 80 mm.Flood levels were significant in the other basins, though the estimated returnperiod was not quite as high; the Chemainus, Cruickshank and Coquitlam riverspeaked at flows exceeding 3-, 6- and 4-year return periods, respectively. Chemai-129ll lllllllllllllllllllllll llllllllllllllllllllll lllllllllllllllllllllllllllllllllllllllll lllllllllllllllllllll lllllllllllllll llCheakamus Daisy Cheakamus Millar ChemainusCoquitlam Cruickshank Elaho02004006000501001502000100200300400255075010020003006009000.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00Starting fscaDischarge (m3 /s)Figure 6.4: Association between fsca at the start of an event and peak dis-charge during identified atmospheric river events from 2003-2013.Green lines indicate the 2-year return period discharge for each studywatershed. A linear model with 95% confidence limits is shown foreach watershed.nus River was not thought to have any snow cover before or after the event. TheCoquitlam and Cruickshank Rivers had fsca levels less than 0.3 before the eventand were estimated to be completely snow free after the event.6.3.5 Regional response to atmospheric riversThe correlation (Spearman’s ρ) in WAR generation between analyzed ASP siteswhere data were available is shown in Table 6.4. The response at the Jump Creeksite was the most dissimilar to the other sites; however, all sites were correlated.Spearman’s ρ for the maximum daily discharge between study watersheds is shown1302003−10−14 2003−10−25025507510002550751000255075100025507510002550751000255075100Cheakamus DaisyCheakamus MillarChemainusCoquitlamCruickshankElaho% coverageStrong Non−snow Weak Non−snow Weak Snow Strong SnowFigure 6.5: Changes in SCA, as represented by the calculated P50 before andafter the October 15-18, 2003, atmospheric river event.131in Table 6.5. Discharge was most correlated between the watersheds of the Cheaka-mus River and the Elaho River. As Cheakamus River at Millar Creek is locatedwithin Cheakamus River at Daisy Lake, discharge was highly correlated betweenthese two gauges.Table 6.4: Spearman’s ρ for WAR production between sites during atmo-spheric river events from 2003-2013 in the south coast regionWOLF USR JUMPWOLF - 0.77 0.67USR 0.77 - 0.54JUMP 0.67 0.54 -Table 6.5: Spearman’s ρ for discharge between study watersheds during at-mospheric river events in the south coast region. Daisy = CheakamusRiver at Daisy Lake, Millar = Cheakamus River at Millar Creek.Daisy Millar Chemainus Coquitlam Cruickshank ElahoDaisy - 0.96 0.51 0.72 0.80 0.92Millar 0.96 - 0.37 0.61 0.73 0.95Chemainus 0.51 0.37 - 0.83 0.68 0.25Coquitlam 0.72 0.61 0.83 - 0.73 0.54Cruickshank 0.80 0.73 0.68 0.73 - 0.58Elaho 0.92 0.95 0.25 0.54 0.58 -6.3.6 Influence of rainfall and initial snow cover on dischargeResults of the site and watershed specific multiple linear regression are shown inTable 6.6. The model that included only rainfall (eq. 6.5) was found to be signifi-cant at levels of 0.95 for all watersheds. Inclusion of the antecedent fsca resulted innominal increases in adjusted R2 for all sites; however, the addition was not foundto be significant at a level of 0.95 for any watershed. Inclusion of antecedent fsca(as in eq. 6.6) appeared to provide the most improvement for the Chemainus Riverwatershed (with significance at level of 0.9), which was also the least correlatedwith high elevation rainfall.132Table 6.6: Results of the multiple regression models fit to predict discharge based on rainfall at the nearest ASP andfsca prior to the onset of the event. Daisy = Cheakamus River at Daisy Lake, Millar = Cheakamus River at MillarCreek.Catchment Model b0 b1 b2 b3 adj. R2 se p d f F p(> F)Daisy eq. 6.5 53.15 3.10 0.83 62.48 0.00 20- eq. 6.6 -13.70 3.30 77.69 0.84 60.57 0.00 19 2.19 0.16- eq. 6.7 4.72 2.82 53.82 0.78 0.85 61.75 0.00 18 0.28 0.60Millar eq. 6.5 7.67 0.91 0.88 16.17 0.00 17- eq. 6.6 -0.76 0.93 9.65 0.88 16.50 0.00 16 0.39 0.54- eq. 6.7 -23.56 1.56 37.19 -0.88 0.90 15.50 0.00 15 3.13 0.10Chemainus eq. 6.5 85.44 0.66 0.39 92.45 0.00 19- eq. 6.6 59.75 0.66 128.65 0.49 86.54 0.00 18 3.50 0.08- eq. 6.7 63.00 0.62 109.85 0.24 0.50 88.80 0.01 17 0.09 0.76Cruickshank eq. 6.5 33.72 0.87 0.74 43.50 0.00 11- eq. 6.6 28.56 0.87 11.13 0.74 45.41 0.00 10 0.09 0.77- eq. 6.7 21.29 1.00 26.68 -0.28 0.75 46.99 0.00 9 0.34 0.57Elaho eq. 6.5 95.58 5.42 0.76 143.40 0.00 17- eq. 6.6 107.44 5.40 -13.02 0.76 147.79 0.00 16 0.01 0.94- eq. 6.7 52.83 6.63 48.94 -1.56 0.76 152.31 0.00 15 0.06 0.811336.4 Discussion6.4.1 Flood frequency analysisThe primary purpose of the flood frequency analysis was to determine the dis-charge levels that represented flood conditions (or near flood conditions) in thesubsequent analyses of individual atmospheric river events. However, the floodfrequency analysis did provide some insights in its own right. Cheakamus Riverat Daisy Lake and Millar Creek, and Elaho River exhibit a mixed distribution ofannual peaks, with annual floods during both spring freshet and winter. For theCoquitlam and Cruickshank rivers, inspection of the annual hydrograph ensemblein Figure 2.7 shows that they do typically experience a spring freshet; however,the annual peak discharges are dominated by winter events. Chemainus River dis-plays little contribution to the annual flood series from spring events. The mixedGumbel distribution still provided a satisfactory fit for watersheds that do not dis-play mixed flood generation processes in the annual flood series. In this scenario,where the winter peak discharge dominates over the summer peak discharges formost years (such as the Chemainus, Coquitlam and Cruickshank rivers), the re-sults of the mixed distribution equation (6.1) will trend towards a standard Gumbeldistribution (Waylen and Woo, 1983).As described in Section 6.3.1, the flood of August 30, 1991, which was thelargest flood on record at the Cheakamus River at Millar Creek and one of themajor floods on the Cheakamus River at Daisy Lake, the Elaho, Cruickshank, andCoquitlam rivers, did not fit the paradigm of Waylen and Woo (1983). Thoughit occurred in the spring/summer, it was driven by a rare major late summer rainevent, and was likely associated with (also rare) heavy rain onto exposed glacierice. This event, combined with the (albeit weak) negative relationship betweeninitial snow cover and discharge levels for the more snow dominant Elaho andCheakamus River basins (Figure 6.4), is an indication that the initial snow covermay not be a strong driver of runoff production during atmospheric rivers in thesewatersheds.It is appropriate to consider that the spring/summer peaks were likely snowmeltdriven, perhaps with an additional rainfall contribution, and the winter peaks were134likely rainfall driven (often with a rain-on-snow component, as Harr (1981) foundin the Willamette River in Oregon, U.S.A.). However, there may still be spring/-summer peaks in which rain-on-snow at least played some part in the flood gen-eration process. Kormos et al. (2014) found that snow melt during late seasonrain-on-snow was largely driven by net all-wave radiation. In Chapter 4, it wasshown that a substantial number of late season rain-on-snow events occurred from2003-2013, generating the possibility that these events could supplement runoffthat was already high due to the spring freshet. In most cases the spring eventsidentified from 2003-2013 in Chapter 4 were associated with light rainfall and didnot correspond with discharge that appeared on the annual flood series.However, there were exceptions, and the largest example is the event from July20 to 22, 2007, in which 160 mm of rainfall was recorded at the DisappointmentLake snow pillow (the Upper Squamish River snow pillow was not operating atthe time). This event resulted in the peak flood for 2007 for the Cheakamus Riverat both Millar Creek (80 m3 s−1) and Daisy Lake (293 m3 s−1), and at the ElahoRiver (528 m3 s−1). The event of July 22 occurred at the time of year when thespring freshet is also typically occurring within these watersheds (see Figure 2.7);hence the rainfall event was likely supplemented an already high discharge.Along with their contribution of an average of 90% of the global polewardtransfer of moisture (Zhu and Newell, 1998), atmospheric rivers appear to con-tribute a significant amount to the annual flood series in southwest British Columbia,regardless of whether the watershed has a large snowmelt component to the annualhydrograph (e.g. Elaho and Cheakamus river watersheds) or not. Chemainus River,which had the largest contribution of atmospheric rivers to the annual flood series(see Table 6.3), was the least likely to have snow at the onset of an atmosphericriver (Figure 6.4). Similarly, Neiman et al. (2011) found that a large majority ofannual floods on rain dominated rivers of western Washington were due to atmo-spheric river events, and unassociated with either antecedent snowcover or snowfalling during the event.1356.4.2 Floods and high elevation runoff generationAs a whole, the rainfall at high elevations, which is the primary contributor toWAR generation during large rain-on-snow events in the region (see Section 4.4.1),does appear to play a dominant role in flood generation in these watersheds. Thissupports the conclusions from Ro¨ssler et al. (2014), who hypothesized that under-representation of rainfall at high elevation was one of the primary causes of under-forecasting of a major rain-on-snow flood in Switzerland. Accurate monitoring ofhigh elevation precipitation and temperature and accurate representation of verti-cal temperature gradients may provide significant improvement of both the under-standing of processes that occur during large rain-on-snow events associated withatmospheric rivers, and improve the ability to model and forecast response at thecatchment scale during these events.The Upper Squamish River ASP experienced more snowfall during atmosphericriver events (see Figure 6.2) than either the Wolf River or Jump Creek ASP. How-ever, the atmospheric river events in which the majority of precipitation fell as rainwere typically associated with high discharge levels (see Figure 6.3) within thenearby watersheds. These three watersheds all display the outlier event of Octo-ber 2003, which was the largest rainfall and WAR generating event at the UpperSquamish River site by a wide margin during the analysis period of 2003-2013,and resulted in the largest floods during the analysis period.The influence of atmospheric rivers on runoff production displayed regionalcoherence across all of the study basins (Table 6.5). Particularly, the responses ofthe Cheakamus and Elaho rivers were very similar across all events. It appears thatsimilarities in hydrologic regime (i.e. if a watershed typically had seasonal snowcover or not), influenced the similarities in response more than spatial proximity.Response at Cruickshank River, which has some snowmelt influence, was morecorrelated with the responses at the Coquitlam, Elaho, and Cheakamus River thanthe Chemainus River, which is in closer proximity. The Chemainus River had themost dissimilar response during atmospheric river events.1366.4.3 Snow covered areaSynchronized snowmelt over large areas is considered a primary requirement forextreme rain-on-snow flooding in the transient snow zone (Jones and Perkins,2010). In general, the Elaho River and two Cheakamus River watersheds water-sheds lie entirely above the transient snow zone in most years. They contain alpineand glacier cover (see Table 2.3) and have a significant contribution of snowmelt tothe annual hydrograph (see Figure 2.7). Additionally, these watersheds commonlyhad nearly full snow coverage at the onset of atmospheric river events (Figure 6.4).The October 2003 appears to be an exceptional event for these watersheds par-tially due to its occurrence during the early fall after one of the first snowfalls ofthe season; snow coverage was thin across most of the region during the October2003 event, even at higher elevations. Melting of this thin, yet extensive snowcover likely contributed to the very high flood levels (> 50-year return period atthe Cheakamus River watersheds, and > 25-year year return period at the ElahoRiver). This result, along with extreme levels of high elevation rainfall at the Up-per Squamish River ASP, support the idea that upper level snowmelt and loss oflower elevation snow cover contributed to the generation of the high dischargesexperienced during this event, utilizing both of the primary flood generation mech-anisms suggested by Jones and Perkins (2010). However, Neiman et al. (2011)observed that this event also produced significant flooding in western Washingtonin watersheds that did not have any snowpack present. During this event, the alltime 24-hour rainfall total was broken by almost 50% in Seattle (128 mm from aprevious record of 87 mm) (Neiman et al., 2011).SCA did not change enough for other events to provide conclusive data on therole of snowmelt in runoff production during atmospheric rivers in these water-sheds. Since the October 2003 event occurred so early in the winter season, it cor-responded with a thin snow pack at a range of elevations in these catchments. Thesnow pack in these regions deepens quickly beginning in November, and the totalSCA typically remains relatively constant through winter until the spring freshet.The lack of a significant impact of antecedent fsca on the ability to predict dis-charge during atmospheric river events could reflect the lack of variability in fscaat the start of individual events.137For the Chemainus River, which has the most transient snow cover of any ofthe watersheds, snow cover prior to the onset of an atmospheric river event wasweakly associated with discharge response (Figure 6.4). If the Chemainus water-shed contained snow prior to the atmospheric river, it appears that higher flowswere more likely. The fitted multiple linear regression model (Table 6.6) foundthat the addition of antecedent fsca to the model for this watershed was weaklysignificant (p(> F) = 0.08).Snow cover differences may explain why flood levels at the Chemainus Riverwere not as large for the October 2003 event within this watershed. During theOctober 2003 event, snow was not thought to be present on the Chemainus Riverwatershed. Hence, this event likely constituted a pure rain event at the ChemainusRiver, and flows exceeded 3-year return period (as opposed to 25-50 year returnperiods at the Cheakamus and Elaho River gauges). However, this cannot be con-cluded with certainty; differences in local weather may have played a more signif-icant role. Atmospheric rivers are focused bands of moisture, and hence there canbe great differences in precipitation totals over short distances. Unfortunately, dataerrors at the Jump Creek ASP during this event meant that upper level precipitationduring the event cannot be compared to the other ASP sites.It is recognized that the use of days up to 8 days before or after an atmosphericriver event could result in inaccurate estimates of immediate pre- and post- stormsnow coverage. However, it was assumed that the benefits of using days with loweramounts of cloud obscuration outweighed this potential inaccuracy. Additionallyit was assumed that in most situations, the atmospheric river event would be thedominant weather event that occurred within that time period between clear (orrelatively clear) days.The sensitivity to using fsca that was calculated from only the strong snow ornon-snow results from the P50 output was investigated. It was found that when theP50 results were forced into a binary classification of snow (0.5 being the cutofffor snow), essentially allowing both strong and weak predictions of snow and non-snow, the results were similar. The likely reason for this is the use of days with lowlevels of cloud obscuration that resulted in primarily high confidence predictionsat all study watersheds.Based on the results in Table 6.6, it appears that accurate forecasts of high138elevation rainfall will have a more immediate impact on the ability to predict floodresponse during major atmospheric river events than knowledge of the initial snowcover within the watershed. This result is supported by Fleming et al. (2015), whodeveloped ensembles of artificial neural network models for short term hydrologicforecasts on the Englishman River, on Vancouver Island, BC. Fleming et al. (2015)found that addition of antecedent snow as a predictor variable provided little benefitin improved forecasting ability.The influence of snow cover on runoff production could not be detected withthe analysis undertaken here; however, lack of a statistically significant increase inpredictive ability does not suggest that snow cover is unimportant in runoff gener-ation during atmospheric rivers. Chapter 4 showed that large rain-on-snow eventswere typically enhanced by 25-30% through snowmelt at least at upper elevations.This influence may be obscured by within storm spatial variability in precipitation(the ASP sites are not necessarily representative of the precipitation of the full area)and differences in freezing levels between atmospheric river events.For the Chemainus River, which contained the most variable snow cover at theonset of atmospheric river events of the study watersheds, the starting fsca appearedto have some (albeit only weakly statistically significant) detectable influence onrunoff generation during atmospheric rivers.6.5 ConclusionsThe primary conclusions that can be drawn from this study are as follows: (1) theannual flood series in the south coast of BC is significantly influenced by the occur-rence of atmospheric rivers, and this result appears to be independent of whether ornot the watershed experiences a seasonal snow cover; (2) heavy rainfall at high ele-vations appears to be a significant factor in the runoff that is generated at the catch-ment scale, regardless of whether or not a watershed contains a seasonal snowmeltcomponent; (3) antecedent SCA was only found to slightly increase the ability topredict runoff response during atmospheric rivers in the study watershed with theleast consistent snowcover.139Chapter 7ConclusionsThis research has used a combination of historical data from a variety of sources,along with remotely sensed information on snow cover and weather reanalysis out-put, to address issues limiting the understanding of rain-on-snow at multiple scales.These issues focused on a lack of high elevation data for analysis, a lack of under-standing of the nature of rain-on-snow over a wide range of event sizes, uncertaintyand issues of cloud obscuration in remotely sensed snow cover information, and therole of snow and high elevation rain on flood generation during atmospheric rivers.7.1 Summary of key findingsIn Chapter 3, it was shown that among several potentially relevant output variablesfrom the North American Regional Reanalysis (NARR) model, air temperature andvapour pressure likely had the greatest potential use for hydrological applicationsin this region. Wind speed was systematically over-predicted and had variablelevels of correlation with observed data depending on the site. Shortwave radiationwas accurate during sunny days, but did not capture the effects of cloud coverwell. While longwave radiation out performed some more common methods ofmodelling longwave radiation, there were not enough observations to fully test theoutput.Additionally in Chapter 3, the utility of NARR downscaled air temperature wasemployed as high altitude guidance in an operational hydrologic forecasting model140for the Cheakamus River at Daisy Lake on a single atmospheric river event that wasunder-predicted by the standard operational hydrologic model. Results showed thatnudging air temperature model inputs during an atmospheric river event (which areoften associated with atypical vertical temperature gradients) to match the freezinglevel predicted with information from NARR air temperatures resulted in a greaterpercentage of the modelled catchment receiving rainfall instead of snowfall andpeak flows that better matched observations. Further work on approaches to incor-porate dynamic lapse rates in forecasting models is needed to avoid over-predictingsnow melt at low elevations using this method.Chapter 4 addressed the lack of understanding of point processes antecedentto and during rain-on-snow events over ten years of observations at five differentASP sites. Results showed that WAR is typically enhanced by 25-30% at the ASPsites during large rain-on-snow events. For smaller events, a number of antecedentand meteorologic factors likely influence WAR generation along with rainfall, in-cluding the antecedent liquid water content of the snowpack and air temperaturesduring the event; however, more detailed observations are required to determine theexact nature of these influences at the ASP sites. Additionally, spring and summerwere identified as the most common time of year for rain-on-snow in this region.These late season events were distinguished by a substantially larger portion ofWAR generation due to snow melt than mid-winter events due to their concurrentoccurrence with spring and summer snow ablation.Chapter 5 addressed difficulties in optical remote sensing of snow cover in re-gions that experience significant cloud cover via the MODIS instruments. A novelapproach was developed to probabilisitically classify snow cover in both visibleand cloud obscured pixels. Results showed that the classifier slightly improved ac-curacy over the stock binary classification methods when forced to make a decision(between a pixel being snow or non-snow covered). However, more importantly,the classification model implicitly provided an estimate of the certainty of the clas-sification that could be synthesized with estimates of the certainty of predictionswithin cloud obscured pixels. Accuracy of the cloud infilling method was found tocompare well with the SNOWL method from Parajka et al. (2010) for infilling bi-nary MODIS snow cover data, and accuracy did not diminish as cloud obscurationincreased. Instead, only the confidence in the predictions decreased, as would be141expected.Chapter 6 used information and methods developed in the previous chapters,along with flood frequency analysis, to investigate the drivers of flood generationduring atmospheric river events at six catchments in the south coast region of BC.A large portion of annual peak flow events were associated with atmospheric riversin this region, and this result was independent of whether or not a watershed hada large contribution of snowmelt to the annual hydrograph. High elevation rain-fall, measured at the ASP sites nearest to the study watersheds, was a significantpredictor of discharge during atmospheric river events, but the role of snow coverwas unclear. Antecedent snow covered area had the greatest potential predictivepower in the watershed with the least consistent snow cover, but the results werenot statistically significant.A summary of key findings follows:• Air temperature and vapour pressure from NARR reanalysis output likelyhold the greatest hydrologic potential, but further work should focus on pre-cipitation• High altitude temperature nudging from NARR assisted in improving accu-racy of modelled runoff during an atmospheric river event in coastal BC.• At the point scale, rainfall percolation was the primary source of WAR gen-eration during a range of magnitudes of rain-on-snow events, but events with> 40 mm rainfall were typically enhanced 25 to 30% by snowmelt.• When considering events of all sizes, late spring and summer is the mostcommon time of year for rain-on-snow events in southwest BC.• A method for probabilistic classification and cloud infilling of satellite ob-served snowcover compared favourably with more common methods, andprovided the added benefit of directly supplying confidence estimation inthe snowcover classification.• Atmospheric rivers contribute to a significant portion of the peak flow eventsin watersheds of southwest BC, whether or not they are traditionally consid-ered rain dominant watersheds.142• Inclusion of antecedent snowcover, as delineated by satellite, did not signifi-cantly increase predictive ability of runoff response during atmospheric riverevents in southwest BC over rainfall alone. However, in the catchment withthe greatest variability in antecedent snow cover, that variable was a weaklysignificant predictor of peak flow magnitude.7.2 Significance and implications of the findingsThe results of this thesis has advanced our process based understanding of rain-on-snow at a range of scales, provided tools and methods that other researchers canuse to continue along this path, and generated new knowledge to help improve thehydrologic forecasting of rain-on-snow events.While many of the research findings in this thesis are specific to the region andthe specific sites and watersheds that were analyzed, these findings are relevant toother parts of the world that experience rain-on-snow, such as the mountains ofNew Zealand, Japan, and mainland Europe. Rain-on-snow is a global phenomenon(Merz and Blo¨schl, 2003), and in some regions it may be increasing due to highertemperatures (Freudiger et al., 2014), while at the same time decreasing in other re-gions due to a loss of snow cover (McCabe et al., 2007). These results suggest that,in the context of a warming climate, there may be a period of peak rain-on-snow in-fluence on flood generation. Indeed, Brunengo (2012) determined, through MonteCarlo analysis of synthetic storm characteristics and snow occurrence, where theoptimal elevation band of rain-on-snow is located in the Oregon Cascades, and thatthe location of these bands should be expected to be highly sensitive to climatechange.This thesis reinforced the notion that the range of different watershed typesin southwestern BC all experience peak flows due to rain-on-snow at one timeor another. Whitfield et al. (2002) found that climate change will likely impactthese different watershed types differently, with earlier snowmelt onset for snowdriven watersheds, increased winter peak flows in rain dominant watersheds, andtransition to rain dominance for hybrid watersheds.Schnorbus et al. (2014) found that the hydrology of the Campbell River basin,on Vancouver Island, BC, was the most sensitive of three study basins (within three143different climatic regions) in British Columbia to climate change. The other basinsthey assessed (Peace River and Upper Columbia River) were predicted to remainlargely snowmelt driven watersheds, whereas the Campbell River was expected tohave a significant decrease in snow accumulation and hence spring snowmelt. In-fluence of winter rainfall was expected to become a more significant runoff source.Additionally, Loukas et al. (2002) found that increasing winter rainfall and de-creasing snowmelt contribution due to climate change is expected to increase themagnitude of peak flows in the Upper Campbell River by 6%. The results of thisthesis reinforce the notion that the influence of increasing winter rainfall will likelydominate over the effects of decreasing winter snowcover within these regions, im-plying a potential net increase in winter peak flows within the study watersheds.Within watersheds that already experience most peak flows due to winter rainin BC, peak flows are also expected to increase due to climate change. Westonet al. (2003) predicted that peak annual floods would increase by approximately17% by 2080 on the Englishman River on Vancouver Island, BC. This increaseimplied that the current 20-year flood will represent the 10-year flood by 2080, andthat the 2080 bankfull flood is expected to be similar in magnitude to the current2-year flood.Appreciation of the significance of atmospheric rivers to flood generation inwestern North America in general has been reinforced in this thesis. Atmosphericriver occurrence is projected to increase with climate change (Radic´ et al., 2015),and this thesis has demonstrated the significant amount of annual floods that aredriven by atmospheric river events within this region. This thesis has shown thatthe antecedent snowcover during atmospheric rivers likely plays a secondary rolein flood generation during atmospheric rivers, hence the potential decrease in snowaccumulation due to climate change is unlikely to mitigate the potential increasesin atmospheric river occurrence.7.3 Recommendations for future researchDue to a lack of observations, it was not possible to test the accuracy of NARR forestimating surface precipitation. Future work should assess the accuracy of NARRprecipitation compared to surface observations in mountainous areas. However,144testing NARR precipitation output would be more challenging than for variablessuch as air temperature because (a) precipitation can exhibit greater sub-grid hor-izontal variability, posing greater uncertainty in downscaling from gridded valuesto points, and (b) measured precipitation at high-elevation sites can be subject toconsiderable undercatch due to high wind speed, especially for snow. If NARRprecipitation is found to be reasonably accurate, however, high elevation precipita-tion intelligence from NARR could prove useful for increasing both hydrologicalprocess understanding and model accuracy.In relation to operational modelling of rain-on-snow, an analysis of a singleatmospheric river event in Chapter 3 showed the potential utility of more accu-rate representation of the vertical temperature gradients in improving predictions.This result indicates that there may be attainable gains in hydrologic forecast ac-curacy during major rain-on-snow events that do not require either substantiallymore complex observations (i.e. full energy balance measurement) or substantialchanges to the model setup. Future work should focus on optimal methods forassimilating NARR air temperature into a forecast system. As operational hydro-logic models are highly tuned to specific watersheds, the most effective method forassimilation of NARR output will likely require extensive calibration in order toimprove runoff forecasts during non-standard situations while maintaining stabil-ity during standard forecasting situations. Though NARR’s temporal availabilitywould not allow it to be used in place of high altitude observations during real-time forecasting, assimilation and use of NARR air temperature for hind-castingof historic rain-on-snow events could provide guidance for model formulation andcalibration and be used to help determine the most useful locations for new AWSinstallation. Additionally, further work should test the utility of directly using nu-merical weather prediction model output, as opposed to reanalysis model output,in order to implement changes to lapse rate in real-time for hydrologic forecasting.The analysis in Chapter 4 showed that there were insights to be gained aboutthe rainfall component of rain-on-snow and the typical enhancement of snowmeltfor large events; however, long term tracking of the energy balance in order to ac-curately predict snow melt for smaller events requires more detailed observationsthan are available from operational monitoring networks. Additionally, uncertaintyin the partitioning between rain and snow likely led to the exclusion of the smallest145rain-on-snow events in this analysis. At the point scale, recommendations for fu-ture work center on the long term collection of more detailed data at sites similar tothe automated snow pillow sites and over a range of elevations and wind exposurelevels.In an ideal long-term monitoring program, the most important additional ob-servations that should be collected are: (1) on-site wind speed and relative humid-ity measurements to help estimate turbulent fluxes and allow for the use of wetbulb temperature to better predict rain/snow discrimination; (2) cameras or opticalprecipitation gauges to assist in determining if rain or snow is falling with greateraccuracy; (3) measured WAR via snowmelt lysimeters; (4) snow depth informationvia optical or sonic measurements (the Upper Squamish River ASP does measuresnow depth, however data were deemed too noisy during rain-on-snow at an eventscale to provide reliable insights); (5) measurement of long- and shortwave radi-ation fluxes; and (6) measurements of the temperature profile and stratigraphy ofthe snowpack. It is recognized that the irregular nature of rain-on-snow makes thislevel of data collection, which would have to be facilitated over the long term, alarge undertaking.Remote sensing of SWE would provide the greatest utility for helping to deter-mine in finer detail the spatial contribution of snow to runoff during atmosphericrivers as well as other rain-on-snow events. Alternative methods for obtaining re-motely sensed SWE, either from passive or active microwave methods, hold themost promise for achieving this goal, and are less influenced by cloud obscuration.However, passive microwave methods have been limited by saturation associatedwith deep snowpacks and overly coarse resolution, and active microwave methodsmay be affected during times with substantial atmospheric moisture such as dur-ing rain-on-snow (Nolin, 2010). LiDAR measurement of snow depth (the primarysource of spatial variability of SWE within a watershed) is only beginning to realizeits significant promise for measuring the spatial variability of snow (Deems et al.,2013). Currently, airborne LiDAR methods are hindered by the high cost of flights,but price may become less of an issue as the proliferation of unmanned aerial ve-hicles (drones) continues. Remote sensing technologies are evolving rapidly andthese methods will become more feasible over time as both new airborne and satel-lite based instruments are deployed and tested.146Further work at the point scale and also on charactering spatial snow distribu-tion will help to inform knowledge on rain-on-snow at the watershed scale. Gaug-ing of nested catchments with different elevation ranges would be valuable in gain-ing more insight into the elevational distribution of runoff generation during rain-on-snow. Given the proper driving data, spatially distributed hydrologic modellingcould be used as a heuristic tool to further explore the sensitivity of streamflow dur-ing atmospheric rivers to changing antecedent snowcover within these watersheds.147ReferencesAbatzoglou, J. T. 2013. Development of gridded surface meteorological data forecological applications and modelling. International Journal of Climatology,33(1):121–131. → pages 60Andreadis, K. M. and Lettenmaier, D. P. 2006. Assimilating remotely sensedsnow observations into a macroscale hydrology model. Advances in WaterResources, 29(6):872–886. → pages 3, 116Arino, O., Perez, J., Kalogirou, V., Defourny, P., and Achard, F. 2009.GLOBCOVER 2009. Technical report, European Space Agency, Frascati, Italy.→ pages 23, 164, 165Barry, R. G. 2008. Mountain Weather and Climate. Cambridge University Press.→ pages 5Bastola, S. and Misra, V. 2014. Evaluation of dynamically downscaled reanalysisprecipitation data for hydrological application. Hydrological Processes,28(4):1989–2002. → pages 5BC Hydro 2002. Bridge-Coastal Fish and Wildlife Restoration Program, StrategicPlan, Volume 2. Technical report, BC Hydro. Burnaby, BC. → pages 54Berg, N., Osterhuber, R., and Bergman, J. 1991. Rain-induced outflow from deepsnowpacks in the central Sierra Nevada, California. Hydrological SciencesJournal, 36(6):611–629. → pages 8, 67Berris, S. N. and Harr, R. D. 1987. Comparative snow accumulation and meltduring rainfall in forested and clear-cut plots in the Western Cascades ofOregon. Water Resources Research, 23(1):135. → pages 7, 66Beven, K. and Westerberg, I. 2011. On red herrings and real herrings:disinformation and information in hydrological inference. HydrologicalProcesses, 25(10):1676–1680. → pages 10, 97148Bolstad, P. V., Swift, L., Collins, F., and Re´gnie`re, J. 1998. Measured andpredicted air temperatures at basin to regional scales in the southernAppalachian mountains. Agricultural and Forest Meteorology, 91(3):161–176.→ pages 4Bolton, D. 1980. The Computation of Equivalent Potential Temperature. MonthlyWeather Review, 108(7):1046–1053. → pages 38Bristow, K. L. and Campbell, G. S. 1984. On the relationship between incomingsolar radiation and daily maximum and minimum temperature. Agriculturaland Forest Meteorology, 31(2):159–166. → pages 5Brunengo, M. J. 2012. Where is the rain-on-snow zone in the west-centralWashington Cascades? Monte Carlo simulation of large storms in theNorthwest. Phd thesis, Portland State University, Portland. → pages 143Buck, A. 1981. New equations for computing vapor pressure and enhancementfactor. Journal of Applied Meteorology, 20:1527–1532. → pages 37Carlson, T. N. 1980. Airflow through midlatitude cyclones and the comma cloudpattern. Monthly Weather Review, 108(10):1498–1509. → pages 13Choi, W., Kim, S. J., Rasmussen, P. F., and Moore, A. R. 2009. Use of the NorthAmerican Regional Reanalysis for hydrological modelling in Manitoba.Canadian Water Resources Journal, 34(1):17–36. → pages 36Chokmani, K., Dever, K., Bernier, M., Gauthier, Y., and Paquet, L. M. 2010.Adaptation of the SNOWMAP algorithm for snow mapping over easternCanada using Landsat-TM imagery. Hydrological Sciences Journal,55(4):649–660. → pages 27Conway, H. and Benedict, R. 1994. Infiltration of water into snow. WaterResources Research, 30(3):641–649. → pages 9, 65Corripio, J. G. 2004. Snow surface albedo estimation using terrestrialphotography. International Journal of Remote Sensing, 25(24):5705–5729. →pages 27, 28, 30, 31Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S.,Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars,A., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes,M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Ho´lm, E. V.,Isaksen, L., Ka˚llberg, P., Ko¨hler, M., Matricardi, M., McNally, A. P.,149Monge-Sanz, B. M., Morcrette, J. J., Park, B. K., Peubey, C., de Rosnay, P.,Tavolato, C., The´paut, J. N., and Vitart, F. 2011. The ERA-Interim reanalysis:Configuration and performance of the data assimilation system. QuarterlyJournal of the Royal Meteorological Society, 137(656):553–597. → pages 6Deems, J. S., Painter, T. H., and Finnegan, D. C. 2013. Lidar measurement ofsnow depth: A review. Journal of Glaciology, 59(215):467–479. → pages 146Deng, X. and Stull, R. 2007. Assimilating surface weather observations fromcomplex terrain into a high-resolution numerical weather prediction model.Monthly Weather Review, 135(3):1037–1054. → pages 6Dettinger, M. 2004. Fifty-two years of pineapple-express storms across the WestCoast of North America. Technical report, U.S. Geological Survey, ScrippsInstitution of Oceanography for the California Energy Commission, PIEREnergy-Related Environmental Research, La Jolla, CA. → pages 14Dettinger, M. and Ingram, B. 2013. The coming megafloods. Scientific American,308:64–71. → pages 14Dettinger, M. D., Ralph, F. M., Das, T., Neiman, P. J., and Cayan, D. R. 2011.Atmospheric rivers, floods and the water resources of California. Water,3(2):445–478. → pages 14, 34, 121, 127Dietz, A. J., Kuenzer, C., Gessner, U., and Dech, S. 2012. Remote sensing ofsnow: a review of available methods. International Journal of Remote Sensing,33(13):4094–4134. → pages 9, 97, 114Dingman, S. L. 2015. Physical Hydrology. Waveland Press, Inc., Long Grove, IL,third edition. → pages 92Dozier, J. 1989. Spectral signature of alpine snow cover from the LandsatThematic Mapper. Remote Sensing of Environment, 28:9–22. → pages 25Dozier, J., Painter, T. H., Rittger, K., and Frew, J. E. 2008. Time–space continuityof daily maps of fractional snow cover and albedo from MODIS. Advances inWater Resources, 31(11):1515–1526. → pages 10Eaton, B., Moore, R. D., Spittlehouse, D. L., Whitfield, P. H., and Stahl, K. 2010.Regional Hydrology. In Pike, R. G., Redding, T. E., Moore, R. D., Winkler,R. D., and Bladon, K. D., editors, Compendium of Forest Hydrology andGeomorphology in British Columbia, chapter 4, pages 85–110. BC Ministry ofForests and Range, Forest Science Program. → pages 2, 18150Evensen, G. 2003. The ensemble Kalman filter: Theoretical formulation andpractical implementation. Ocean Dynamics, 53(4):343–367. → pages 116Fiddes, J. and Gruber, S. 2014. TopoSCALE v.1.0: downscaling gridded climatedata in complex terrain. Geoscientific Model Development, 7(1):387–405. →pages 6, 60, 62Finger, D., Pellicciotti, F., Konz, M., Rimkus, S., and Burlando, P. 2011. Thevalue of glacier mass balance, satellite snow cover images, and hourlydischarge for improving the performance of a physically based distributedhydrological model. Water Resources Research, 47(7):14. → pages 117Fitzharris, B., Lawson, W., and Owens, I. 1999. Research on glaciers and snow inNew Zealand. Progress in Physical Geography, 23(4):469–500. → pages 2Fleming, S. W., Bourdin, D. R., Campbell, D., Stull, R. B., and Gardner, T. 2015.Development and operational testing of a super-ensemble artificial intelligenceflood-forecast model for a Pacific Northwest river. JAWRA Journal of theAmerican Water Resources Association, 51(2):502–512. → pages 139Fleming, S. W., Weber, F., and Weston, S. 2010. Multiobjective, manifoldlyconstrained Monte Carlo optimization and uncertainty estimation for anoperational hydrologic forecast model. In American Meteorological Society90th Annual Meeting, page 30, Atlanta, GA. → pages 30, 55Fleming, S. W., Whitfield, P. H., Moore, R. D., and Quilty, E. J. 2007.Regime-dependent streamflow sensitivities to Pacific climate modes cross theGeorgia-Puget transboundary ecoregion. Hydrological Processes,21:3264–3287. → pages 1, 18, 32Franz, K. J. and Karsten, L. R. 2013. Calibration of a distributed snow modelusing MODIS snow covered area data. Journal of Hydrology, 494:160–175. →pages 117Freudiger, D., Kohn, I., Stahl, K., and Weiler, M. 2014. Large-scale analysis ofchanging frequencies of rain-on-snow events with flood-generation potential.Hydrology and Earth System Sciences, 18(7):2695–2709. → pages 118, 143Gafurov, A. and Ba´rdossy, A. 2009. Cloud removal methodology from MODISsnow cover product. Hydrology and Earth System Sciences, 13(7):1361–1373.→ pages 10, 97, 102, 115151Garen, D. C. and Marks, D. 2005. Spatially distributed energy balance snowmeltmodelling in a mountainous river basin: estimation of meteorological inputsand verification of model results. Journal of Hydrology, 315(1-4):126–153. →pages 4, 5Garvelmann, J., Pohl, S., and Weiler, M. 2015. Spatio-temporal controls ofsnowmelt and runoff generation during rain-on-snow events in a mid-latitudemountain catchment. Hydrological Processes, 29:3649–3664. → pages 2, 3,66, 119Gerdel, R. W. 1954. The transmission of water through snow. AmericanGeophysical Union, Transactions, 35:475–485. → pages 8Gray, D. M. and Landine, P. G. 1988. An energy-budget snowmelt model for theCanadian Prairies. Canadian Journal of Earth Sciences, 25(8):1292–1303. →pages 74Hall, D. K., Foster, J. L., Verbyla, D. L., Klein, A. G., and Benson, C. S. 1998.Assessment of snow-cover mapping accuracy in a variety of vegetation-coverdensities in central Alaska. Remote Sensing of Environment, 66(2):129–137. →pages 26Hall, D. K. and Riggs, G. A. 2007. Accuracy assessment of the MODIS snowproducts. Hydrological Processes, 21(12):1534–1547. → pages 26, 96, 97,111, 112Hall, D. K., Riggs, G. A., and Salomonson, V. V. 1995. Development of methodsfor mapping global snow cover using moderate resolution imagingspectroradiometer data. Remote Sensing of Environment, 54(2):127–140. →pages 110, 112Hall, D. K., Riggs, G. A., Salomonson, V. V., DiGirolamo, N. E., and Bayr, K. J.2002. MODIS snow-cover products. Remote Sensing of Environment,83(1-2):181–194. → pages 96Halverson, M. J. and Fleming, S. W. 2015. Complex network theory, streamflow,and hydrometric monitoring system design. Hydrology and Earth SystemSciences, 19(7):3301–3318. → pages 32Harr, R. D. 1981. Some characteristics and consequences of snowmelt duringrainfall in western Oregon. Journal of Hydrology, 53(3):277–304. → pages 7,11, 12, 66, 73, 135152Harr, R. D. 1986. Effects of clearcutting on rain-on-snow runoff in westernOregon: A new look at old studies. Water Resources Research,22(7):1095–1100. → pages 7, 12Hasenauer, H., Merganicova, K., Petritsch, R., Pietsch, S. A., and Thornton, P. E.2003. Validating daily climate interpolations over complex terrain in Austria.Agricultural and Forest Meteorology, 119(1-2):87–107. → pages 4Hock, R. 1999. A distributed temperature-index ice- and snowmelt modelincluding potential direct solar radiation. Journal of Glaciology,45(149):101–111. → pages 39, 50, 52Holden, Z. A., Swanson, A., Klene, A. E., Abatzoglou, J. T., Dobrowski, S. Z.,Cushman, S. A., Squires, J., Moisen, G. G., and Oyler, J. W. 2015.Development of high-resolution (250 m) historical daily gridded airtemperature data using reanalysis and distributed sensor networks for the USNorthern Rocky Mountains. International Journal of Climatology. → pages 36Hollinger, J., Peirce, J., and Poe, G. 1990. SSM/I instrument evaluation. IEEETransactions on Geoscience and Remote Sensing, 28(5):781–790. → pages 13Hopkinson, C., Collins, T., Anderson, A., Pomeroy, J., and Spooner, I. 2012.Spatial snow depth assessment using LiDAR transect samples and public GISdata layers in the Elbow River Watershed, Alberta. Canadian Water ResourcesJournal, 37(2):69–87. → pages 11Hrachowitz, M., Savenije, H. H. G., Bloschl, G., McDonnell, J. J., Sivapalan, M.,Pomeroy, J. W., Arheimer, B., Blume, T., Clark, M. P., Ehret, U., Fenicia, F.,Freer, J. E., Gelfan, A., Gupta, H., Hughes, D., Hut, R., Montanari, A., Pande,S., Tetzlaff, D., Troch, P. a., Uhlenbrook, S., Wagener, T., Winsemius, H.,Woods, R., Zehe, E., and Cudennec, C. 2013. A decade of Predictions inUngauged Basins (PUB): a review. Hydrological Sciences Journal,58(6):1198–1255. → pages 35Immerzeel, W., Droogers, P., de Jong, S., and Bierkens, M. 2009. Large-scalemonitoring of snow cover and runoff simulation in Himalayan river basinsusing remote sensing. Remote Sensing of Environment, 113(1):40–49. → pages117Iqbal, M. 1983. An Introduction to Solar Radiation. Elsevier, Don Mills, Ontario.→ pages 39153Jarosch, A. H., Anslow, F. S., and Clarke, G. K. C. 2012. High-resolutionprecipitation and temperature downscaling for glacier models. ClimateDynamics, 38(1-2):391–409. → pages 36Jarvis, C. H. and Stuart, N. 2001. A comparison among strategies forinterpolating maximum and minimum daily air temperatures. Part I: Theselection of ”guiding” topographic and land cover variables. Journal of AppliedMeteorology, 40:1060–1074. → pages 4Jennings, K. and Jones, J. 2015. Precipitation-snowmelt timing and snowmeltaugmentation of large peak flow events, western Cascades, Oregon. WaterResources Research, 51:7649–7661. → pages 3, 8, 66, 67, 119Jones, J. and Perkins, R. 2010. Extreme flood sensitivity to snow and forestharvest, western Cascades, Oregon, United States. Water Resources Research,46:21. → pages 2, 3, 11, 12, 66, 85, 91, 119, 137Jost, G., Moore, R. D., Smith, R., and Gluns, D. R. 2012. Distributedtemperature-index snowmelt modelling for forested catchments. Journal ofHydrology, 420-421:87–101. → pages 39, 72, 73, 74, 76, 77, 79, 80, 90, 92Kattelmann, R. 1986. Measurements of snow layer water retention. InProceedings of the Symposium: Cold Regions Hydrology, pages 377–386,Bethesda, Maryland. American Water Resources Association. → pages 8, 9Kattelmann, R. 1987. Water release from a forested snowpack during rainfall.Forest Hydrology and Watershed Management - IAHS Publication,167:265–272. → pages 7Keshta, N. and Elshorbagy, A. 2011. Utilizing North American RegionalReanalysis for modeling soil moisture and evapotranspiration in reconstructedwatersheds. Physics and Chemistry of the Earth, Parts A/B/C, 36(1):31–41. →pages 36Kimball, J., Running, S., and Nemani, R. 1997. An improved method forestimating surface humidity from daily minimum temperature. Agriculturaland Forest Meteorology, 85(1-2):87–98. → pages 5Klein, A. G., Hall, D. K., and Riggs, G. A. 1998. Improving snow cover mappingin forests through the use of a canopy reflectance model. HydrologicalProcesses, 12(10-11):1723–1744. → pages 26, 27, 97, 113, 117154Klinka, K., Pojar, J., and Meidinger, D. V. 1991. Revision of BiogeoclimaticUnits of coastal British Columbia. Northwest Science, 65(1):32–47. → pages18, 23, 24Kormos, P. R., Marks, D., McNamara, J. P., Marshall, H., Winstral, A., andFlores, A. N. 2014. Snow distribution, melt and surface water inputs to the soilin the mountain rainsnow transition zone. Journal of Hydrology, 519:190–204.→ pages 7, 67, 90, 135Kroczynski, S. 2004. A comparison of two rain-on-snow events and thesubsequent hydrologic responses in three small river basins in centralPennsylvania. Technical report, NOAA/National Weather Service, StateCollege, PA. → pages 2, 3, 11, 12, 13, 67, 119Kuchment, L. S., Romanov, P., Gelfan, a. N., and Demidov, V. N. 2010. Use ofsatellite-derived data for characterization of snow cover and simulation ofsnowmelt runoff through a distributed physically based model of runoffgeneration. Hydrology and Earth System Sciences, 14(2):339–350. → pages11, 117Kumar, S. and Merwade, V. 2011. Evaluation of NARR and CLM3.5 outputs forsurface water and energy budgets in the Mississippi River Basin. Journal ofGeophysical Research, 116:21. → pages 61Lavalle´e, S., Brissette, F. P., Leconte, R., Larouche, B., and Others 2006.Monitoring snow-cover depletion by coupling satellite imagery with adistributed snowmelt model. Journal of water resources planning andmanagement, 132:71. → pages 3Leathers, D. J., Kluck, D. R., and Kroczynski, S. 1998. The severe flooding eventof January 1996 across North-Central Pennsylvania. Bulletin of the AmericanMeteorological Society, 79(5):785–797. → pages 2, 3Leavesley, G. H. and Stannard, L. G. 1995. The precipitation-runoff modelingsystem-PRMS. In Singh, V. P., editor, Computer models of watershedhydrology, pages 281–310. Water Resources Publications, Fort Collins, CO. →pages 53Li, L. and Pomeroy, J. W. 1997. Estimates of threshold wind speeds for snowtransport using meteorological data. Journal of Applied Meteorology,36(3):205–213. → pages 92155Liston, G. E. and Elder, K. 2006. A meteorological distribution system forhigh-resolution terrestrial modeling (MicroMet). Journal of Hydrometeorology,7(2):217–234. → pages 4, 36Liston, G. E. and Hiemstra, C. a. 2008. A simple data assimilation system forcomplex snow distributions (SnowAssim). Journal of Hydrometeorology,9(5):989–1004. → pages 117Liston, G. E. and Sturm, M. 1998. A snow-transport model for complex terrain.Journal of Glaciology, 44(148):498–516. → pages 60Lo´pez-Burgos, V., Gupta, H. V., Clark, M., Lopez-Burgos, V., Gupta, H. V., andClark, M. 2013. Reducing cloud obscuration of MODIS snow cover areaproducts by combining spatio-temporal techniques with a probability of snowapproach. Hydrology and Earth System Sciences, 17(5):1809–1823. → pages97Loukas, A., Vasiliades, L., and Dalezios, N. 2002. Potential climate changeimpacts on flood producing mechanisms in southern British Columbia, Canadausing the CGCMA1 simulation results. Journal of Hydrology, 259:163–188. →pages 144Loukas, A., Vasiliades, L., and Dalezios, N. R. 2000. Flood producingmechanisms identification in southern British Columbia, Canada. Journal ofHydrology, 227(1-4):218–235. → pages 2, 3, 18Luo, Y., Berbery, E. H., Mitchell, K. E., and Betts, A. K. 2007. Relationshipsbetween land surface and near-surface atmospheric variables in the NCEPNorth American Regional Reanalysis. Journal of Hydrometeorology,8(6):1184–1203. → pages 36Marks, D. and Dozier, J. 1992. Climate and energy exchange at the snow surfacein the alpine region of the Sierra Nevada. Water Resources Research,28(11):3043–3054. → pages 72Marks, D., Kimball, J., Tingey, D., and Link, T. 1998. The sensitivity of snowmeltprocesses to climate conditions and forest cover during rain on snow: a casestudy of the 1996 Pacific Northwest flood. Hydrological Processes,12:1569–1587. → pages 2, 3, 7, 8, 12, 15, 60, 66, 73, 90, 93, 119Marks, D., Link, T., Winstral, A., and Garen, D. 2001. Simulating snowmeltprocesses during rain-on-snow over a semi-arid mountain basin. Annals ofGlaciology, 32(1):195–202. → pages 3, 7, 66156Mastin, M., Gendaszek, A., and Barnas, C. 2010. Magnitude and Extent ofFlooding at Selected River Reaches in Western Washington , January 2009.Technical report, United States Geological Survey, Reston, VA. → pages 14Mazurkiewicz, A. B., Callery, D. G., and McDonnell, J. J. 2008. Assessing thecontrols of the snow energy balanc and water available for for runoff in arain-on-snow environment. Journal of Hydrology, 354(1-4):1–14. → pages 3,8, 66, 67, 70McCabe, G. J., Clark, M. P., and Hay, L. E. 2007. Rain-on-snow events in thewestern United States. Bulletin of the American Meteorological Society,88(3):319–328. → pages 2, 95, 118, 119, 143Mebane, W. R. and Sekhon, J. S. 2011. Genetic optimization using derivatives:the rgenoud package for R. Journal of Statistical Software, 42(11):1–26. →pages 69Merz, R. and Blo¨schl, G. 2003. A process typology of regional floods. WaterResources Research, 39(12):20. → pages 2, 118, 143Mesinger, F., DiMego, G., Kalnay, E., Mitchell, K., Shafran, P. C., Ebisuzaki, W.,Jovic, D., Woollen, J., Rogers, E., and Berbery, E. H. 2006. North AmericanRegional Reanalysis. Bulletin of the American Meteorological Society,87(3):343–360. → pages 6, 36, 37Meteorologic Service of Canada 2015. MANOBS - Manual of Surface WeatherObservations. Technical report, Environment Canada, Gatineau, QC. → pages69Meyer, D., Dimitriadou, E., Hornik, K., Weingessel, A., Leisch, F., Chang, C., andLin, C. 2015. e1071: Misc Functions of the Department of Statistics,Probability Theory Group (Formerly: E1071), TU Wien. → pages 101Mitchell, T. 1997. Machine Learning. New York, NY. → pages 100, 101Moore, R. and Prowse, T. 1988. Snow hydrology of the Waimakariri catchment,South Island, New Zealand. Journal of Hydrology, 27(1):44–68. → pages 15Moore, R. D. 1993. Application of a conceptual streamflow model in a glacierizeddrainage basin. Journal of Hydrology, 150(1):151–168. → pages 37, 53, 56Moore, R. D., Woods, R. W., and Boyle, D. P. 2013. Putting PUB into practice inmountainous regions. Streamline Watershed Management Bulletin,15(2):12–21. → pages 2, 4, 35, 87157Moosavi, V., Malekinezhad, H., and Shirmohammadi, B. 2014. Fractional snowcover mapping from MODIS data using wavelet-artificial intelligence hybridmodels. Journal of Hydrology, 511:160–170. → pages 114Mott, R. and Lehning, M. 2010. Meteorological modeling of very high-resolutionwind fields and snow deposition for mountains. Journal of Hydrometeorology,11(4):934–949. → pages 5Musial, J. P., Hu¨sler, F., Su¨tterlin, M., Neuhaus, C., and Wunderle, S. 2014.Probabilistic approach to cloud and snow detection on advanced very highresolution radiometer (AVHRR) imagery. Atmospheric MeasurementTechniques, (2005):799–822. → pages 96Nash, J. E. and Sutcliffe, J. V. 1970. River flow forecasting through conceptualmodels part 1 - A discussion of principals. Journal of Hydrology,10(3):282–290. → pages 39, 55, 74Neiman, P. J., Ralph, F. M., Wick, G. A., Kuo, Y.-H., Wee, T.-K., Ma, Z., Taylor,G. H., and Dettinger, M. D. 2008a. Diagnosis of an intense atmospheric riverimpacting the Pacific Northwest: storm summary and offshore vertical structureobserved with COSMIC satellite retrievals. Monthly Weather Review,136(11):4398–4420. → pages 14, 119Neiman, P. J., Ralph, F. M., Wick, G. A., Lundquist, J. D., and Dettinger, M. D.2008b. Meteorological characteristics and overland precipitation impacts ofatmospheric rivers affecting the West Coast of North America based on eightyears of SSM/I satellite observations. Journal of Hydrometeorology,9(1):22–47. → pages 13, 18, 34, 119, 121Neiman, P. J., Schick, L. J., Ralph, F. M., Hughes, M., and Wick, G. A. 2011.Flooding in western Washington: The connection to atmospheric rivers.Journal of Hydrometeorology, 12(6):1337–1358. → pages 14, 15, 53, 119, 135,137Newell, R., Newell, N., Zhu, Y., and Scott, C. 1992. Tropospheric rivers?-A pilotstudy. Geophysical Research Letters, 19(24):2401–2404. → pages 13Nolin, A. W. 2010. Recent advances in remote sensing of seasonal snow. Journalof Glaciology, 56(200):1141–1150. → pages 11, 146Painter, T. H., Rittger, K., McKenzie, C., Slaughter, P., Davis, R. E., and Dozier, J.2009. Retrieval of subpixel snow covered area, grain size, and albedo fromMODIS. Remote Sensing of Environment, 113(4):868–879. → pages 10, 96,113, 114, 117158Parajka, J. and Blo¨schl, G. 2006. Validation of MODIS snow cover images overAustria. Hydrology and Earth System Sciences, 10:679–689. → pages 9, 97Parajka, J. and Blo¨schl, G. 2008. Spatio-temporal combination of MODISimages–potential for snow cover mapping. Water Resources Research,44(3):W03406. → pages 10, 97, 102, 115, 117Parajka, J., Pepe, M., Rampini, a., Rossi, S., and Blo¨schl, G. 2010. A regionalsnow-line method for estimating snow cover from MODIS during cloud cover.Journal of Hydrology, 381(3):203–212. → pages 10, 104, 115, 141Perkins, R. M. and Jones, J. A. 2008. Climate variability, snow, andphysiographic controls on storm hydrographs in small forested basins, westernCascades, Oregon. Hydrological Processes, 22(25):4949–4964. → pages 2, 11,12, 15, 119Pradhanang, S. M., Frei, A., Zion, M., Schneiderman, E. M., Steenhuis, T. S., andPierson, D. 2013. Rain-on-snow runoff events in New York. HydrologicalProcesses, 27(21):3035–3049. → pages 2, 12, 118Prata, A. J. 1996. A new long-wave formula for estimating downward clear-skyradiation at the surface. Quarterly Journal of the Royal Meteorological Society,122(533):1127–1151. → pages 5, 39, 53, 62Quick, M. C. and Pipes, A. 1977. U.B.C. Watershed Model. HydrologicalSciences Journal, 22(1):153–162. → pages 53, 54, 68R Development Core Team 2015. R: a language and environment for statisticalcomputing. → pages 69Radic´, V., Cannon, A. J., Menounos, B., and Gi, N. 2015. Future changes inautumn atmospheric river events in British Columbia, Canada, as projected byCMIP5 global climate models. Journal of Geophysical Research: Atmospheres,120(18):9279–9302. → pages 144Ralph, F. M., Neiman, P. J., and Rotunno, R. 2005. Dropsonde Observations inLow-Level Jets over the Northeastern Pacific Ocean from CALJET-1998 andPACJET-2001: Mean Vertical-Profile and Atmospheric-River Characteristics.Monthly Weather Review, 133(4):889–910. → pages 13Rittger, K., Painter, T. H., and Dozier, J. 2013. Assessment of methods formapping snow cover from MODIS. Advances in Water Resources, 51:367–380.→ pages 96, 97, 113, 114159Rosin, K. 2010. Development, evaluation and application of dominant runoffgeneration processes in hydrological modeling. PhD thesis, University ofBritish Columbia, Vancouver. → pages 23, 30Ro¨ssler, O., Froidevaux, P., Bo¨rst, U., Rickli, R., Martius, O., and Weingartner, R.2014. Retrospective analysis of a nonforecasted rain-on-snow flood in the Alps:a matter of model limitations or unpredictable nature? Hydrology and EarthSystem Sciences, 18(6):2265–2285. → pages 2, 3, 4, 11, 12, 15, 63, 136Ruane, A. C. 2010. NARR’s atmospheric water cycle components. Part I: 20-yearmean and annual interactions. Journal of Hydrometeorology, 11(6):1205–1219.→ pages 36Running, S. W., Nemani, R. R., and Hungerford, R. D. 1987. Extrapolation ofsynoptic meteorological data in mountainous terrain and its use for simulatingforest evapotranspiration and photosynthesis. Canadian Journal of ForestResearch, 17(6):472–483. → pages 4, 5, 35Salomonson, V. V. and Appel, I. 2004. Estimating fractional snow cover fromMODIS using the normalized difference snow index. Remote Sensing ofEnvironment, 89(3):351–360. → pages 26, 113, 114Schnorbus, M., Werner, A., and Bennett, K. 2014. Impacts of climate change inthree hydrologic regimes in British Columbia, Canada. HydrologicalProcesses, 28(3):1170–1189. → pages 1, 143Schroeder, T. A., Hember, R., Coops, N. C., and Liang, S. 2009. Validation ofsolar radiation surfaces from MODIS and reanalysis data over topographicallycomplex terrain. Journal of Applied Meteorology and Climatology,48(12):2441–2458. → pages 61Shea, J. M. 2010. Regional-scale distributed modelling of glacier meteorologyand melt, southern Coast Mountains, Canada. PhD thesis, University of BritishColumbia, Vancouver. → pages 21, 23Shea, J. M. and Moore, R. D. 2010. Prediction of spatially distributedregional-scale fields of air temperature and vapor pressure over mountainglaciers. Journal of Geophysical Research, 115(D23):15. → pages 23, 37Sheffield, J., Livneh, B., and Wood, E. F. 2012. Representation of terrestrialhydrology and large-scale drought of the continental United States from theNorth American Regional Reanalysis. Journal of Hydrometeorology,13(3):856–876. → pages 36160Singh, P., Spitzbart, G., Hubl, H., Weinmeister, H. W., Hu¨bl, H., Weinmeister,H. W., Hubl, H., and Weinmeister, H. W. 1997. Hydrological response ofsnowpack under rain-on-snow events: A field study. Journal of Hydrology,202(1-4):1–20. → pages 2, 8, 9, 12, 67Sivapalan, M., Takeuchi, K., Franks, S. W., Gupta, V. K., Karambiri, H., Lakshmi,V., Liang, X., McDonnell, J. J., Mendiondo, E. M., O’Connell, P. E., Oki, T.,Pomeroy, J. W., Schertzer, D., Uhlenbrook, S., and Zehe, E. 2003. IAHSDecade on Predictions in Ungauged Basins (PUB), 2003-2012: Shaping anexciting future for the hydrological sciences. Hydrologic Sciences Journal,48(6):857–880. → pages 35, 87Slater, A. G. and Clark, M. P. 2006. Snow data assimilation via an ensembleKalman filter. Journal of Hydrometeorology, 7(3):478–493. → pages 116Spry, C. M., Kohfeld, K. E., Allen, D. M., Dunkley, D., and Lertzman, K. 2014.Characterizing pineapple express storms in the lower mainland of BritishColumbia, Canada. Canadian Water Resources Journal, 39(3):302–323. →pages 14, 15, 53, 90, 119, 120Stahl, K., Moore, R. D., Floyer, J. A., Asplin, M. G., and McKendry, I. G. 2006.Comparison of approaches for spatial interpolation of daily air temperature in alarge region with complex topography and highly variable station density.Agricultural and Forest Meteorology, 139(3-4):224–236. → pages 4Stegall, S. T. and Zhang, J. 2012. Wind field climatology, changes, and extremesin the Chukchi - Beaufort Seas and Alaska North Slope during 1979 - 2009.Journal of Climate, 25(23):8075–8089. → pages 60Stull, R. 2011. Wet-bulb temperature from relative humidity and air temperature.Journal of Applied Meteorology and Climatology, 50(11):2267–2269. → pages75Sui, J. and Koehler, G. 2001. Rain-on-snow induced flood events in southernGermany. Journal of Hydrology, 252(1-4):205–220. → pages 2Tedesco, M. and Narvekar, P. S. 2010. Assessment of the NASA AMSR-E SWEProduct. Selected Topics in Applied Earth Observations and Remote Sensing,IEEE Journal of, 3(1):141–159. → pages 11Therneau, T., Atkinson, B., and Ripley, B. 2015. rpart: Recursive partitioning andregression trees. → pages 76161Thornton, P. E., Running, S. W., and White, M. A. 1997. Generating surfaces ofdaily meteorological variables over large regions of complex terrain. Journal ofHydrology, 190(3-4):214–251. → pages 4, 36Trubilowicz, J. W., Moore, R. D., and Buttle, J. M. 2013. Prediction ofstream-flow regime using ecological classification zones. HydrologicalProcesses, 27(13):1935–1944. → pages 24, 32U.S. Army Corps of Engineers 1956. Snow hydrology: summary report of thesnow investigations. Technical report, U.S. Army Corps of Engineers, Portland,OR. → pages 7, 73van Heeswijk, M., Kimball, J. S., and Marks, D. G. 1996. Simulation of wateravailable for runoff in clearcut forest openings during rain-on-snow events inthe western Cascade Range of Oregon and Washington. Technical report,United States Geological Survey, Tacoma, Washington. → pages 3, 7, 8, 66,90, 92Vikhamar, D. and Solberg, R. 2003. Snow-cover mapping in forests byconstrained linear spectral unmixing of MODIS data. Remote Sensing ofEnvironment, 88(3):309–323. → pages 114Wade, N. L., Martin, J., and Whitfield, P. H. 2001. Hydrologic and climaticzonation of Georgia Basin, British Columbia. Canadian Water ResourcesJournal, 26(1):43–70. → pages 1, 18Way, R. G. and Bonnaventure, P. P. 2015. Testing a reanalysis-based infillingmethod for areas with sparse discontinuous air temperature data in northeasternCanada. Atmospheric Science Letters, 16(3):398–407. → pages 5Waylen, P. and Woo, M.-K. 1983. Annual floods in southwestern BritishColumbia, Canada. Journal of Hydrology, 62(1-4):95–105. → pages 122, 125,126, 134Weber, F., Garen, D., and Gobena, A. 2012. Invited commentary: themes andissues from the workshop Operational river flow and water supply forecasting.Canadian Water Resources Journal, 37(3):151–161. → pages 4Weston, S., Guthrie, R., and McTaggart-Cowan, J. 2003. The vulnerability ofLower Englishman River to modelled climate change. Canadian WaterResources Journal, 28(4):657–669. → pages 144162Whitaker, A. C. and Sugiyama, H. 2005. Seasonal snowpack dynamics and runoffin a cool temperate forest: lysimeter experiment in Niigata, Japan.Hydrological Processes, 19(20):4179–4200. → pages 2Whitfield, P. H., Cannon, A. J., and Reynolds, C. J. 2002. Modelling streamflowin present and future climates: examples from the Georgia Basin, BritishColumbia. Canadian Water Resources Journal, 27(4):427–456. → pages 1, 143Wilby, R. L., Hay, L. E., Gutowski, W. J., Arritt, R. W., Takle, E. S., Pan, Z.,Leavesley, G. H., and Clark, M. P. 2000. Hydrological responses todynamically and statistically downscaled climate model output. GeophysicalResearch Letters, 27(8):1199–1202. → pages 6Winkler, R. D., Moore, R. D., Redding, T. E., Spittlehouse, D. L., Carlyle-Moses,D. E., and Smerdon, B. D. 2010. Hydrologic Processes and WatershedResponse. In Compendium of Forest Hydrology and Geomorphology in BritishColumbia, chapter 6, pages 133–177. Ministry of Forests and Range ForestScience Program. → pages 18Winstral, A. and Marks, D. 2002. Simulating wind fields and snow redistributionusing terrain-based parameters to model snow accumulation and melt over asemi-arid mountain catchment. Hydrological Processes, 16(18):3585–3603. →pages 5Woo, M. K. and Thorne, R. 2006. Snowmelt contribution to discharge from alarge mountainous catchment in subarctic Canada. Hydrological Processes,20(10):2129–2139. → pages 36Zhang, H. 2005. Exploring conditions for the optimality of naı¨ve bayes.International Journal of Pattern Recognition and Artificial Intelligence,19(02):183–198. → pages 100Zhu, Y. and Newell, R. E. 1998. A proposed algorithm for moisture fluxes fromatmospheric rivers. Monthly Weather Review, 126(3):725–735. → pages 13,135163Appendix AGlobcover 2009The legend for the Globcover 2009 dataset from Arino et al. (2009) is shown inTable A.1. The accompanied simplified class, as described paragraph 2.2.2 areshown with each category.164Table A.1: Simplified Categories from the Globcover 2009 dataset (Arino et al., 2009)Value Label Simplified Class11 Post-flooding or irrigated croplands (or aquatic) 114 Rainfed croplands 120 Mosaic cropland (50-70%) / vegetation (grassland/shrubland/forest) (20-50%) 230 Mosaic vegetation (grassland/shrubland/forest) (50-70%) / cropland (20-50%) 240 Closed to open (>15%) broadleaved evergreen or semi-deciduous forest (>5m) 250 Closed (>40%) broadleaved deciduous forest (>5m) 360 Open (15-40%) broadleaved deciduous forest/woodland (>5m) 270 Closed (>40%) needleleaved evergreen forest (>5m) 390 Open (15-40%) needleleaved deciduous or evergreen forest (>5m) 2100 Closed to open (>15%) mixed broadleaved and needleleaved forest (>5m) 2110 Mosaic forest or shrubland (50-70%) / grassland (20-50%) 2120 Mosaic grassland (50-70%) / forest or shrubland (20-50%) 2130 Closed to open (>15%) (broadleaved or needleleaved evergreen or deciduous) shrubland (<5m) 2140 Closed to open (>15%) herbaceous vegetation (grassland savannas or lichens/mosses)” 2150 Sparse (<15%) vegetation 1160 Closed to open (>15%) broadleaved forest regularly flooded (semi-permanently or temporarily) - Fresh or brackish water 2170 Closed (>40%) broadleaved forest or shrubland permanently flooded - Saline or brackish water 3180 Closed to open (>15%) grassland or woody vegetation on regularly flooded or waterlogged soil - Fresh brackish or saline water 2190 Artificial surfaces and associated areas (Urban areas >50%) 6200 Bare areas 1210 Water bodies 5220 Permanent snow and ice 4230 No data (burnt areas clouds) NA165Appendix BInfilled time series at ASP SitesMissing data (P and Ta) at each ASP site was infilled using the Daisy Lake cli-mate station. For air temperature, missing data was infilled based on linear re-gression models, stratified by month, between the ASP sites and the Daisy Lakestation using functions provide by R.D. Moore (available at:∼rdmoore/Rcode.htm). For precipitation, a single linear regression, forced tohave a 0 intercept, was fit between monthly total precipitation (mm) at each ASPsite and the Daisy Lake site. Up to four days (96 hours) of missing precipitationdata was allowed for a monthly total to be included in the linear model fit for theChilliwack River, Jump Creek, Upper Squamish River and Wolf River sites. Dueto the substantial amount of missing data at the Disappointment lake ASP, thisthreshold was expanded to allow up to 8 days (192) hours of missing data for themonth for a month to be included in the regression. It is recognized that this couldlead to an under-representation of the precipitation gradient between the two sites;however, it was thought that the benefits of inclusion of more months in the regres-sion model out-weighed this. As can be seen in Figure B.9, data was not availableto meet this threshold for any winter months (December through March).166−10 0 10 20 30 40−200102030Daisy Ta (°C)ASP Ta (°C)lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll1 3 5 7 9 11−20−1001020Monthy − y^ (°C)Figure B.1: Correlation between hourly Ta at the Chilliwack River ASP andDaisy Lake climate station on the left pane. Residuals by month areshown in the right pane.167−15 −10 −5 0 5 10−20010Daisy Ta (°C)ASP Ta (°C) Jan−10 −5 0 5 10 15−20−10010Daisy Ta (°C)ASP Ta (°C) Feb−15 −5 0 5 10 15 20−20010Daisy Ta (°C)ASP Ta (°C) Mar−5 0 5 10 15 20 25−15−5515Daisy Ta (°C)ASP Ta (°C) Apr0 5 10 15 20 25 30−5515Daisy Ta (°C)ASP Ta (°C) May10 20 30 40−551525Daisy Ta (°C)ASP Ta (°C) Jun5 10 15 20 25 30 35 40051525Daisy Ta (°C)ASP Ta (°C) Jul5 10 15 20 25 30 35051525Daisy Ta (°C)ASP Ta (°C) Aug0 5 10 15 20 25 30051525Daisy Ta (°C)ASP Ta (°C) Sep−5 0 5 10 15 20 25−1001020Daisy Ta (°C)ASP Ta (°C) Oct−15 −10 −5 0 5 10−20010Daisy Ta (°C)ASP Ta (°C) Nov−15 −10 −5 0 5 10−20−10010Daisy Ta (°C)ASP Ta (°C) DecFigure B.2: Correlation between hourly Ta at the Chilliwack River ASP andDaisy Lake climate station, stratified by month.168−30−100102030ASP Ta (°C)2004 2005 2006 2007 2008 2009 2010 2011 2012 2013actual infilledFigure B.3: Observed and infilled Ta at the Chilliwack River ASP.llllllllllllllllllllllllllllllllllllll llll lllll0 100 200 300 400 500100300500Daisy monthly P (mm)ASP monthly P (mm)l1 3 5 7 9 11−100100300MonthPrecip (mm)Figure B.4: Correlation between monthly P at the Chilliwack River ASP andDaisy Lake climate station on the left pane. Residuals by month areshown in the right pane.169−200204060P (mm)2004 2005 2006 2007 2008 2009 2010 2011 2012 2013actual infilledFigure B.5: Observed and infilled hourly P (mm) at the Chilliwack RiverASP.170−10 0 10 20 30 40−100102030Daisy Ta (°C)ASP Ta (°C)llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll lll1 3 5 7 9 11−10051015Monthy − y^ (°C)Figure B.6: Correlation between hourly Ta at the Disappointment Lake ASPand Daisy Lake climate station on the left pane. Residuals by monthare shown in the right pane.171−15 −10 −5 0 5 10−15−505Daisy Ta (°C)ASP Ta (°C) Jan−10 −5 0 5 10 15−100515Daisy Ta (°C)ASP Ta (°C) Feb−15 −5 0 5 10 15 20−100515Daisy Ta (°C)ASP Ta (°C) Mar−5 0 5 10 15 20 25−100515Daisy Ta (°C)ASP Ta (°C) Apr0 5 10 15 20 25 30051525Daisy Ta (°C)ASP Ta (°C) May10 20 30 40051525Daisy Ta (°C)ASP Ta (°C) Jun5 10 15 20 25 30 35 4051525Daisy Ta (°C)ASP Ta (°C) Jul5 10 15 20 25 30 35051525Daisy Ta (°C)ASP Ta (°C) Aug0 5 10 15 20 25 3051020Daisy Ta (°C)ASP Ta (°C) Sep−5 0 5 10 15 20 25−551525Daisy Ta (°C)ASP Ta (°C) Oct−15 −10 −5 0 5 10−100510Daisy Ta (°C)ASP Ta (°C) Nov−15 −10 −5 0 5 10−505Daisy Ta (°C)ASP Ta (°C) DecFigure B.7: Correlation between hourly Ta at the Disappointment Lake ASPand Daisy Lake climate station, stratified by month.172−30−100102030ASP Ta (°C)2004 2005 2006 2007 2008 2009 2010 2011 2012 2013actual infilledFigure B.8: Observed and infilled Ta at the Disappointment Lake ASP.lllll llllllllll50 100 200 300100200300400500Daisy monthly P (mm)ASP monthly P (mm)4 5 6 7 8 9 11−50050100MonthPrecip (mm)Figure B.9: Correlation between monthly P at the Disappointment Lake ASPand Daisy Lake climate station on the left pane. Residuals by monthare shown in the right pane.173−200204060P (mm)2004 2005 2006 2007 2008 2009 2010 2011 2012 2013actual infilledFigure B.10: Observed and infilled hourly P (mm) at the DisappointmentLake ASP.174−10 0 10 20 30 40−200102030Daisy Ta (°C)ASP Ta (°C)lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll1 3 5 7 9 11−20−10010Monthy − y^ (°C)Figure B.11: Correlation between hourly Ta at the Jump Creek ASP andDaisy Lake climate station on the left pane. Residuals by month areshown in the right pane.175−15 −10 −5 0 5 10−15−5515Daisy Ta (°C)ASP Ta (°C) Jan−10 −5 0 5 10 15−100510Daisy Ta (°C)ASP Ta (°C) Feb−15 −5 0 5 10 15 20−100515Daisy Ta (°C)ASP Ta (°C) Mar−5 0 5 10 15 20 25−50510Daisy Ta (°C)ASP Ta (°C) Apr0 5 10 15 20 25 30−5515Daisy Ta (°C)ASP Ta (°C) May10 20 30 40051525Daisy Ta (°C)ASP Ta (°C) Jun5 10 15 20 25 30 35 4051525Daisy Ta (°C)ASP Ta (°C) Jul5 10 15 20 25 30 35051525Daisy Ta (°C)ASP Ta (°C) Aug0 5 10 15 20 25 30−2001020Daisy Ta (°C)ASP Ta (°C) Sep−5 0 5 10 15 20 25−1001020Daisy Ta (°C)ASP Ta (°C) Oct−15 −10 −5 0 5 10−15−5515Daisy Ta (°C)ASP Ta (°C) Nov−15 −10 −5 0 5 10−15−55Daisy Ta (°C)ASP Ta (°C) DecFigure B.12: Correlation between hourly Ta at the Jump Creek ASP andDaisy Lake climate station, stratified by month.176−30−100102030ASP Ta (°C)2004 2005 2006 2007 2008 2009 2010 2011 2012 2013actual infilledFigure B.13: Observed and infilled Ta at the Jump Creek ASP.lllllll lllllllllllllllllllllllllllllllllllllll llllllllll llllllllllllll0 100 200 300 4000200400600800Daisy monthly P (mm)ASP monthly P (mm)ll lllll1 3 5 7 9 11−2000200400MonthPrecip (mm)Figure B.14: Correlation between monthly P at the Jump Creek ASP andDaisy Lake climate station on the left pane. Residuals by month areshown in the right pane.177−200204060P (mm)2004 2005 2006 2007 2008 2009 2010 2011 2012 2013actual infilledFigure B.15: Observed and infilled hourly P (mm) at the Jump Creek ASP.178−10 0 10 20 30 40−200102030Daisy Ta (°C)ASP Ta (°C) lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll1 3 5 7 9 11−1001020Monthy − y^ (°C)Figure B.16: Correlation between hourly Ta at the Upper Squamish RiverASP and Daisy Lake climate station on the left pane. Residuals bymonth are shown in the right pane.179−15 −10 −5 0 5 10−20010Daisy Ta (°C)ASP Ta (°C) Jan−10 −5 0 5 10 15−20−10010Daisy Ta (°C)ASP Ta (°C) Feb−15 −5 0 5 10 15 20−20010Daisy Ta (°C)ASP Ta (°C) Mar−5 0 5 10 15 20 25−15−5515Daisy Ta (°C)ASP Ta (°C) Apr0 5 10 15 20 25 30−551525Daisy Ta (°C)ASP Ta (°C) May10 20 30 40051525Daisy Ta (°C)ASP Ta (°C) Jun5 10 15 20 25 30 35 40051525Daisy Ta (°C)ASP Ta (°C) Jul5 10 15 20 25 30 35051525Daisy Ta (°C)ASP Ta (°C) Aug0 5 10 15 20 25 30051525Daisy Ta (°C)ASP Ta (°C) Sep−5 0 5 10 15 20 25−1001020Daisy Ta (°C)ASP Ta (°C) Oct−15 −10 −5 0 5 10−20010Daisy Ta (°C)ASP Ta (°C) Nov−15 −10 −5 0 5 10−25−15−55Daisy Ta (°C)ASP Ta (°C) DecFigure B.17: Correlation between hourly Ta at the Upper Squamish RiverASP and Daisy Lake climate station, stratified by month.180−30−100102030ASP Ta (°C)2004 2005 2006 2007 2008 2009 2010 2011 2012 2013actual infilledFigure B.18: Observed and infilled Ta at the Upper Squamish River ASP.lllllllllllllllll llllllllllllllllllllllllllllllllllllll0 100 200 300 400 5000100300500Daisy monthly P (mm)ASP monthly P (mm)llll1 3 5 7 9 11−100050100MonthPrecip (mm)Figure B.19: Correlation between monthly P at the Upper Squamish RiverASP and Daisy Lake climate station on the left pane. Residuals bymonth are shown in the right pane.181−200204060P (mm)2004 2005 2006 2007 2008 2009 2010 2011 2012 2013actual infilledFigure B.20: Observed and infilled hourly P (mm) at the Upper SquamishRiver ASP.182−10 0 10 20 30 40−200102030Daisy Ta (°C)ASP Ta (°C)lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll1 3 5 7 9 11−20−1001020Monthy − y^ (°C)Figure B.21: Correlation between hourly Ta at the Wolf River ASP and DaisyLake climate station on the left pane. Residuals by month are shownin the right pane.183−15 −10 −5 0 5 10−20010Daisy Ta (°C)ASP Ta (°C) Jan−10 −5 0 5 10 15−2001020Daisy Ta (°C)ASP Ta (°C) Feb−15 −5 0 5 10 15 20−15−5515Daisy Ta (°C)ASP Ta (°C) Mar−5 0 5 10 15 20 25−1001020Daisy Ta (°C)ASP Ta (°C) Apr0 5 10 15 20 25 30−1001020Daisy Ta (°C)ASP Ta (°C) May0 10 20 30 40−551525Daisy Ta (°C)ASP Ta (°C) Jun5 10 15 20 25 30 35 40051525Daisy Ta (°C)ASP Ta (°C) Jul5 10 15 20 25 30 35051525Daisy Ta (°C)ASP Ta (°C) Aug0 5 10 15 20 25 30 35051525Daisy Ta (°C)ASP Ta (°C) Sep−5 0 5 10 15 20 25−551525Daisy Ta (°C)ASP Ta (°C) Oct−15 −5 0 5 10−2001020Daisy Ta (°C)ASP Ta (°C) Nov−15 −10 −5 0 5 10−20−1005Daisy Ta (°C)ASP Ta (°C) DecFigure B.22: Correlation between hourly Ta at the Wolf River ASP and DaisyLake climate station, stratified by month.184−30−100102030ASP Ta (°C)2004 2005 2006 2007 2008 2009 2010 2011 2012 2013actual infilledFigure B.23: Observed and infilled Ta at the Wolf River ASP.lllllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllllll0 100 200 300 400 5000100300500Daisy monthly P (mm)ASP monthly P (mm)llll1 3 5 7 9 11−100050100MonthPrecip (mm)Figure B.24: Correlation between monthly P at the Wolf River ASP andDaisy Lake climate station on the left pane. Residuals by month areshown in the right pane.185−200204060P (mm)2004 2005 2006 2007 2008 2009 2010 2011 2012 2013actual infilledFigure B.25: Observed and infilled hourly P (mm) at the Wolf River ASP.186


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items