MODELLING SNOWPACKS IN TOPOGRAPHICALLY COMPLEX TERRAIN DURING THE WINTER ACCUMULATION SEASON, COLUMBIA RIVER BASIN, BRITISH COLUMBIA, CANADA by Marzieh Mortezapour M.Sc., Warsaw University of Technology, 2014 M.Sc., Tarbiat Modares University, 2010 B.Sc., Hormozgan University, 2002 DISSERTATION SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY IN NATURAL RESOURCES AND ENVIRONMENTAL STUDIES UNIVERSITY OF NORTHERN BRITISH COLUMBIA December 2020 ©Marzieh Mortezapour, 2020 Abstract This dissertation investigates winter accumulation and snow cover change in the Columbia Mountains of British Columbia. In chapter 1, I start with an introduction that describes the study area, and then outlines the objectives and structure of this dissertation. In chapter 2, I examine the performance of two snow evolution models with different complexities (SnowModel and Alpine3D) at simulating winter glacier mass balance on four individual glaciers using two different forcing datasets, the Weather Research and Forecasting model (WRF) outputs and the North American Land Data Assimilation System (NLDAS). My results show that both models can simulate winter accumulation with less than 20% bias for each glacier, with SnowModel forced by WRF yielding the least overall bias. In chapter 3, I study the effect of wind on snow patterns to determine the impact of snow redistribution by wind in terms of erosion, deposition, and sublimation on winter mass balance estimation. The results demonstrate that modelled redistribution of snow by wind produces a visually realistic pattern of snow accumulation when compared to observed snow depth, but its impact on the glacier-averaged winter mass balance estimation is negligible (< 4%). The results also suggest that drifting snow sublimation is highly time and space dependent. Considering the model performance from previous chapters, in chapter 4, I analyzed the future snow cover change over the upper Columbia Basin under the Representative Concentration Pathway (RCP8.5) climate scenario by the end of the 21st century. I used downscaled climate projections of the Community Earth System Model (CESM1) by WRF, along with statistically downscaled data provided from the Pacific Climate Impacts Consortium (PCIC) to force SnowModel. The simulated snow maps represent a higher dynamically downscaled mean snow water equivalent (SWE) reduction – reaching up to 30% by the end of the century - than the ii statistically downscaled SWE reduction. While SWE reduction of more than 60% happens at lower and mid-elevations, altitudes higher than 2000 m are less vulnerable to climate change. I conclude this dissertation (Chapter 5) with a summary of the progress gained, study limitations, suggestions for future research, and research implications. iii Acknowledgements During my doctoral research at UNBC, I gained a lot of support from many people. First and foremost, I would like to express my sincere gratitude to my supervisors, Dr. Brian Menounos and Dr. Peter L. Jackson, for allowing me and supporting me to research this interesting subject. Their guidance, advice, and encouragement have contributed greatly to the completion of my dissertation. I also thank the other committee members, Dr. Stephen J. Déry, Dr. Valentina Radic, and Dr. Youmin Tang, for their contribution to this work. In addition, I am very grateful to Andre R. Erler for providing the main data of this research and Siraj ul Islam for his helpful comments and explanations. I am thankful to all the researchers involved in the ground-based measurements and LiDAR data collection and preparation, especially Ben M. Pelto and all his support. I thank Glen E. Liston for providing the SnowModel code to us and Mathias Bavay for his comments and explanations related to software errors and limitations with Alpine3D. I also thank Dr. You Qin Wang for continued support at the HPC lab. This journey would not have been possible without the support of my family, professors, staff and friends. My deepest thanks go to my family for encouraging me in all my pursuits and inspiring me to follow my dreams. I am especially grateful to my parents for their unconditional love and support. Last but not least, I would like to thank my husband, Hamid, who has been extremely supportive of me throughout this entire process and has made countless sacrifices to help me get to this point. This research was supported by Discovery Grants to Drs. Jackson and Menounos from the Natural Sciences and Engineering Research Council of Canada, the Canada Research Chairs Program and the Columbia Basin Trust. iv Acronyms ASWS Automated snow weather stations BCCA The Bias Correction/Constructed Analogues BCCAQ The Bias Correction/Constructed Analogues with Quantile mapping reordering approach CESM The Community Earth System Model CHM The Canadian Hydrological Model ClimateWNA The climate data for Western North America CRU Climatic Research Unit DEM Digital elevation model DD Dynamically downscaled GCM Global Climate Model GLC Global Land Cover database GPCC The Global Precipitation Climatology Centre HRM High Resolution Meteorological IPCC Intergovernmental Panel on Climate Change LiDAR Light Detection and Ranging data MAE Mean absolute error v NLDAS North American Land Data Assimilation System NSE Nash-Sutcliffe efficiency PBias Percent Bias (PB) PCIC Pacific Climate Impacts Consortium PRISM Parameter-Elevation Relationships on Independent Slopes Model QMAP Quantile mapping RCM Regional Climate Model RCP Representative Concentration Pathway RMSE Root mean squared error SD Statistically downscaled SRTM Shuttle Radar Topography Mission SWE Snow water equivalent WI Willmott Index of Agreement WRF Weather Research and Forecasting model vi Glossary1 Ablation All processes that reduce the mass of the glacier. The main processes of ablation are melting and calving (or, when the glacier nourishes an ice shelf, ice discharge across the grounding line). On some glaciers sublimation, loss of windborne snow and avalanching are significant processes of ablation. Ablation zone The part of the glacier where ablation exceeds accumulation in magnitude, that is, where the cumulative mass balance relative to the start of the mass-balance year is negative. Unless qualified, for example by giving a date within the year, references to the ablation zone refer to its extent at the end of the mass-balance year. The extent of the ablation zone can vary strongly from year to year. Accumulation The main process of accumulation is snowfall. Accumulation also includes deposition of hoar, freezing rain, solid precipitation in forms other than snow, gain of windborne snow, avalanching and basal accumulation (often beneath floating ice). Accumulation season A time span extending from a seasonal minimum of glacier mass to a seasonal maximum. 1 (Cogley et al., 2011) vii Accumulation zone The part of the glacier where accumulation exceeds ablation in magnitude, that is, where the cumulative mass balance relative to the start of the mass-balance year is positive. Albedo The ratio of the reflected flux density to the incident flux density. The albedos of glacier surfaces exceed 0.8 for freshly fallen snow, are less for aged snow and firn, and are significantly less for exposed glacier ice. Snow and ice that are sediment-laden or covered by debris can have albedos still lower. The difference between the albedos of snow and glacier ice is significant in the seasonal evolution of the energy balance and therefore of the rate of surface ablation. Alpine glacier See mountain glacier Altitude The vertical distance of a point above a datum. The vertical datum is usually an estimate of mean sea level. Annual Descriptive of a period equal or approximately equal in duration to a calendar year, such as a hydrological year or mass-balance year. Area-averaged Descriptive of a quantity that has been averaged over part or all of the area of the glacier. viii Avalanche A slide or flow of a mass of snow, firn or ice that becomes detached abruptly, often entraining additional material such as snow, debris and vegetation as it descends. Cirque glacier A glacier occupying a cirque. A cirque is a rounded recess with steep sides and back wall, formed on a mountainside by glacial erosion. Conservation of mass The principle that mass in a system is neither created nor destroyed, expressed by the relation: the rate of change of mass in an element of the system equals the rate at which mass enters the element minus the rate at which mass leaves the element. Cumulative Descriptive of a quantity that has been summed over a span of time. Density The ratio of the mass of any substance to the volume that it occupies. Density is expressed in kg m–3. Deposition The process by which a vapour changes phase directly into a solid. Drifting snow Snow entrained and transported above the surface by the wind. ix Elevation change Vertical change in glacier surface elevation (altitude), typically derived from two elevation measurements. Energy balance A relation describing the change in the amount of energy stored within a defined volume owing to flows of energy across the boundary of the volume. Equilibrium line The set of points on the surface of the glacier where the climatic mass balance is zero at a given moment. The equilibrium line separates the accumulation zone from the ablation zone. - Transient equilibrium line The equilibrium line at any instant, where cumulative ablation balances cumulative accumulation relative to the start of the mass-balance year. Freezing point The temperature, equal to 273.15 K (0 ºC) when the pressure is equal to a standard pressure of 101325 Pa, at which pure water freezes, releasing an amount of energy known as the latent heat of fusion. Geodetic method Any method for determining mass balance by repeated mapping of glacier surface elevations to estimate the volume balance. Glacierized Of a region or terrain, containing glaciers or covered by glacier ice today. x Glacier-wide Descriptive of a quantity that, whether or not it is expressed in specific units, has been measured or estimated over the entire glacier. Glaciological method A method of determining mass balance in-situ on the glacier surface by measurements of accumulation and ablation, generally including measurements at stakes and in snow pits. Surface hoar A surface deposit of interlocking ice crystals formed by the deposition of water vapour. Internal ablation Internal ablation can occur due to strain heating of temperate ice as the ice deforms. However, the largest heat sources for internal ablation are likely to be the potential energy released by downward motion of the ice and of meltwater. Internal accumulation Internal accumulation proceeds by the freezing of water that percolates early in the ablation season into firn that is still cold, heating the firn in the process, or by the freezing of retained pore water during the accumulation season, also releasing latent heat and thus slowing the downward advance of the winter cold wave. Inversion A layer of the atmosphere in which temperature increases with height. xi Latent heat The energy taken up or released per unit mass by a system in a reversible change of phase at constant temperature and pressure. Mass balance The change in the mass of a glacier, or part of the glacier, over a stated span of time; the term mass budget is a synonym. See mass-balance units for recommended units. The span of time is often a year or a season. A seasonal mass balance is nearly always either a winter balance or a summer balance. Point mass balance Mass balance at a particular location on the glacier, for example at an ablation stake or a snow pit. Meltwater The liquid resulting from melting of ice, firn or snow. Mountain glacier A glacier that is confined by surrounding mountain terrain, also called an alpine glacier. Precipitation Liquid or solid products of the condensation of water vapour that fall from clouds or are deposited from the air onto the surface. Retreat Decrease of the length of a flowline, measured from a fixed point. xii Snow depth In the firn area, the vertical distance between the glacier surface and the summer surface; outside the firn area, the vertical distance between the glacier surface and the ice surface (which may be superimposed ice or glacier ice) at the time of observation. Sublimation The process by which a solid changes phase directly into a vapour without melting. Water equivalent A unit, in full the “metre [of] water equivalent”, that is an extension of the SI for describing glacier mass in specific units as the thickness of an equal mass having the density of water. Windborne snow Blowing snow or drifting snow. Windborne snow may be redistributed from one part of the glacier to another. xiii Contents Abstract .................................................................................................................................... ii Acknowledgement ................................................................................................................... iv Acronyms .................................................................................................................................. v Glossary .................................................................................................................................. vii Contents ................................................................................................................................. xiv List of Tables .......................................................................................................................xviii List of Figures ....................................................................................................................... xxi Chapter 1: Introduction 1.1. Importance of glaciers and snow modelling ........................................................................ 1 1.1.1. Snow models ......................................................................................................... 4 1.1.2. Downscaling methods ............................................................................................ 7 1.2. Dissertation objectives .................................................................................................... 8 1.3. Study Region ........................................................................................................................... 10 1.4. Outline of chapters ........................................................................................................ 11 Chapter 2: The Role of Meteorological Forcing and Snow Model complexity in Winter Glacier Mass Balance Estimation, Columbia River Basin, Canada Abstract ............................................................................................................................... 14 2.1. Introduction ................................................................................................................... 15 2.2. Study areas and data used in this study ......................................................................... 17 xiv 2.2.1. Study area .......................................................................................................... 17 2.2.2. Input and evaluations data ................................................................................. 19 2.3. Snow models description and simulation ..................................................................... 24 2.3.1. SnowModel ........................................................................................................ 24 2.3.2. Alpine3D............................................................................................................. 26 2.3.3. Model simulations .............................................................................................. 27 2.4. Results .......................................................................................................................... 28 2.4.1. Evaluation of temperature and precipitation...................................................... 29 2.4.2. Evaluation of SWE and Snow depth ................................................................ 31 2.5. Sensitivity analysis and potential uncertainties ............................................................ 39 2.6. Discussion ..................................................................................................................... 41 2.7. Conclusion ................................................................................................................... 48 Chapter 3: Wind influenced snow redistribution and sublimation over alpine terrain using dynamically downscaled data Abstract ............................................................................................................................... 51 3.1. Introduction .................................................................................................................. 52 3.2. Study sites .................................................................................................................... 54 3.3. Snow evolution model description ............................................................................... 55 3.4. Meteorological data ..................................................................................................... 60 3.5.Model simulations ......................................................................................................... 64 xv 3.6. Results .......................................................................................................................... 66 3.6.1. Patterns of snow erosion and deposition ............................................................. 66 3.6.2. Surface and drifting snow sublimation ............................................................... 73 3.7. Discussion ................................................................................................................... 76 3.8.Conclusion ..................................................................................................................... 81 Chapter 4: Changes in Snow Cover over the Columbia Mountains, Canada, under a Future Climate Change Scenario using a Distributed Snow Model Abstract ................................................................................................................................ 84 4.1. Introduction ................................................................................................................... 85 4.2. Data and Methods ......................................................................................................... 88 4.2.1. Study region.......................................................................................................... 88 4.2.2. SnowModel........................................................................................................... 90 4.2.3. Input data ............................................................................................................. 90 4.2.4. Station data and forcing data correction .............................................................. 92 4.3. Future projections ......................................................................................................... 97 4.3.1. Temperature and precipitation .............................................................................. 97 4.3.2. Mean snow depth and SWE ................................................................................. 98 4.3.3. Interannual variability and elevation dependence ............................................. 103 4.4. Discussion ................................................................................................................... 105 4.4.1. Study limitations ............................................................................................... 109 4.5. Conclusion ............................................................................................................... 110 Chapter 5: Conclusion xvi 5.1. Dissertation summary ................................................................................................. 113 5.2. Suggestions for future research ................................................................................... 118 5.3. Research implications ................................................................................................. 121 Appendix 1. Supplementary information of Chapter 2 ................................................................. 123 2. Supplementary information of Chapter 3 ................................................................. 132 3. Supplementary information of Chapter 4 ................................................................. 134 References ............................................................................................................................. 136 xvii List of Tables Table 2.1. Automated Snow pillow Station characteristics. For Glacier-NMF station, temperature and precipitation derived from BC hydro network (station ID: FDL), while snow depth and SWE derived from Automated snow pillow network from Ministry of Environment (station ID: 2A17). ....................................................................................... 21 Table 2.2. Glaciological visits (number of locations, dates) and geodetic acquisition information (dates and point density) for each glacier. Field dates are the median date of glacier visit. (Derived from Pelto et al., 2019) ................................................................. 23 Table 2.3. The monthly temperature lapse rate and precipitation adjustment factors derived from PRISM, WRF and SnowModel default values (Kunkel, 1989). ............................... 28 Table 2.4. Statistical analysis of SWE time series (in meter water equivalent) directly out of WRF and simulated by SnowModel and ALPINE3D forced with WRF fields and NLDAS dataset for Winter 2013-2014. ........................................................................................... 31 Table 2.5. Comparison of integrated snow depth simulated by models over the whole domain for four glaciers with lidar measurement footprints using WRF and NLDAS forcing data. Winter 2016. ...................................................................................................................... 35 Table 2.6. Comparison of integrated WRF-forced and NLDAS-forced snow depth (m) with lidar measurements on the glaciated and in free glacier areas (same as Table 5 but separated for glaciated and glacier-free areas). ................................................................................. 37 Table 2.7. Snow density (kg m-3) comparison between simulated density by SnowModel and Alpine3D and observed data from snow pit measurements and density calculated from snow pillow data ................................................................................................................ 46 xviii Table 3.1. Snow pillow station characteristics. The observed mean temperature and mean monthly precipitation compared with simulated data from WRF and derived data from NLDAS for winter 2015-2016. .......................................................................................... 61 Table 3.2. Automated weather station characteristics that have been used for wind evaluation for accumulation period 2015-2016................................................................................... 62 Table 3.3. Observed and simulated 95 th percentile and means of wind speed (ms -1) by WRF and MicroMet and standard deviations at station locations............................................... 62 Table 3.4. Statistical analysis of simulated daily snow depth and SWE time series using SnowModel with and without wind for four snow pillows for accumulation period 20132014 ................................................................................................................................... 67 Table 3.5. Comparison of integrated snow depth forced by WRF and NLDAS over the whole domain of each glacier with LiDAR measurements and the related percent bias for accumulation period – Winter 2015-2016 ......................................................................... 68 Table 3.6. Winter mass balance comparison between simulations and observations (stake and LiDAR measurements) on the glaciers for winter 2015-2016........................................... 72 Table 4.1. Stations information used for meteorological evaluation of forcing projections and evaluation of simulated Snow depth by SnowModel for historical period (1979-1994) .. 93 Table 4.2. The NSE, Mean absolute error and correlation coefficients r (all significant at the p level of 0.01) and the relative biases for the (1979–1994) period for 10 specified stations ........................................................................................................................................... 96 Table S1.1. The difference between snow evolution modules of two snow models Alpine3D and SnowModel with more details .................................................................................. 129 xix Table S1.2. The percent bias of daily SWE for SnowModel performance using different lapse rates. User lapse rate derived from PRISM data and default lapse rates within the model ......................................................................................................................................... 130 Table S1.3. Statistics of simulated hourly SWE in respect to the observations at four station locations during winter 2015-2016 .................................................................................. 131 Table S1.4. Comparison of simulated density and SWE by SnowModel, averaged over Kokanee Glacier during winter accumulation 2015-2016 using one-layer and multilayer (6 layers) methods ................................................................................................................ 131 Table S2.1. Statistics of SWE (m.w.e) forced by WRF, with and without wind, against snow pillows observations for accumulation period 2015-2016............................................... 133 xx List of Figures Figure 2.1. Top) Study sites: Zillmer, Nordic, Conrad, and Kokanee glaciers and the nearest weather station locations to each glacier (Table 2.1). The blue outline represents the Canadian portion of Columbia River basin. Bottom) Topography of each individual glacier with LiDAR footprint boundaries (Yellow outlines) and main (ON) glacier boundaries (Red outlines)..................................................................................................................... 18 Figure 2.2. Measurement networks for the (a) Zillmer, (b) Nordic, (c) Kokanee, and (d) Conrad glaciers. Snow depth measurement locations, and snow pit/core locations are pictured. This figure has been adopted from Pelto et al. (2019) .............................................................. 22 Figure 2.3. a) Monthly mean daily temperature (˚C) and b) Monthly mean daily precipitation (mm) time series for the accumulation period 2015-2016 at the station locations. Temperature values derived from stations, WRF outputs, NLDAS dataset, and also from pre-processed modules (MeteoIo and MicroMet) of two snow models forced by WRF. Statistical values (R_squered, RMSE, and Bias) of WRF and NLDAS daily temperature and monthly mean daily precipitation against stations (refers to Table1) data are shown on the graphs ........................................................................................................................... 30 Figure 2.4. Simulated SWE time series by SnowModel and Alpine3D forced by NLADS and WRF data compared with observed snow pillow data at station locations for Winter season 2015-2016. SWE derived directly from WRF also is represented. ................................... 33 Figure 2.5. Boxplots comparing snow water equivalent [m] over terrain at the end of accumulation season 2016 for four glaciers. Black boxes - measurements, Blue boxes SnowModel forced by WRF, Green boxes - SnowModel forced by NLDAS, Olive boxes Alpine3D forced by WRF, and Gray boxes - Alpine3D forced by NLDAS. The boxes span xxi the interquartile range and the median is marked by a horizontal line inside the boxes. The whiskers are the dashed lines that extend to the highest and lowest values of each dataset. For the number of stake observations at each elevation band, please refer to Fig. 2 ........ 34 Figure 2.6. Distributed snow depth changes at the end of the accumulation period (end of April 2016) over four glaciers using SnowModel and Alpine3D with NLDAS and WRF forcing over LiDAR footprints. Black outlines represent the main glacier borders where probing of snow depth was completed. Note: The wind effects on snow pattern are neglected for the simulations ......................................................................................................................... 36 Figure 2.7. Comparison of simulated winter mass balance [m w.e.] by SnowModel and Alpine3D with stake and LiDAR measurements on the main glaciers for the accumulation period of 2015-2016. The notation ‘_B’ denotes the bias (absolute magnitude) of simulated winter mass balance in respect to the observations. ......................................................... 38 Figure 3.1. Four glaciers locations (Stars) over Columbia Mountains, British Columbia, Canada. The topography of each glacier also is displayed using contours over each simulation domain. Yellow outline shows the LiDAR footprint and glaciers with surface mass balance observations are represented by red outlines ............................................... 55 Figure 3.2. Wind roses of observations and WRF simulations, as well as down scaled WRF data by MicroMet at four station locations for the period of 1 October to 30 April 20152016. Colors indicate the 6-hourly mean wind velocity (m s-1), the frequency of occurrence is represented in percentage circles ................................................................................... 64 Figure 3.3. SnowModel flowchart with and without SnowTran3D. ..................................... 66 xxii Figure 3.4. Comparison of simulated snow depth forced by WRF considering wind and without wind effects and simulated snow pattern forced by NLDAS with wind with LiDAR footprint measurements at four glaciers for the accumulation period 1 Oct 2015- 30 Apr 2016. Last column shows the subtracted LiDAR measurements from WRF derived simulated snow depth with wind effects ............................................................................ 69 Figure 3.5. Measured and simulated WRF forced snow depth distribution with elevation over glacier free areas, surrounding each glacier for winter accumulation 2015-2016 ............. 70 Figure 3.6. Mean monthly drifting snow sublimation, Relative humidity and solid precipitation per time step (top); Accumulated transported snow, mean monthly wind speed, and monthly wind direction(bottom) averaged over the main glaciers (Fig.7) during accumulation period........................................................................................................... 72 Figure 3.7. Cumulative mean reduction of SWE (mm w.e.) averaged over the whole domain due to static sublimation at the surface (black, dashed), and drifting snow sublimation (black, solid) and the total sublimation (red, solid) during the complete simulation for four glacier locations ............................................................................................................. 73 Figure 3.8. Scatter plots of drifting snow sublimation rate against air temperature, relative humidity, wind speed and wind direction over main glaciers ........................................... 74 Figure 3.9. Summed static surface (up(a)) and drifting snow sublimation(bottom(b)) at four glaciers in (mm w.e.). Fluxes are cumulated over the simulation period (winter 2015-2016). The black and red outlines show the LiDAR footprint and main glaciers boundaries, respectively ....................................................................................................................... 75 xxiii Figure 3.10. Downscaled mean wind speed (m/s) and direction over each domain by MicroMet for winter 2015-2016 ......................................................................................................... 77 Figure 4.1. Digital elevation map of Columbia Mountains. Weather stations used for simulations are marked with points ................................................................................... 89 Figure 4.2. Mean monthly dynamically and statistically downscaled temperature and precipitation against observations for period 1979-1994 at station locations (See table 1 for stations information). The mean monthly offset used for precipitation and temperature adjustment for DD forced data during the whole period 1979 through 2099 .................... 94 Figure 4.3. Averaged mean monthly simulated temperature, precipitation and snow depth against observed data from 10 stations (Table.4.1) for reference period 1979-1994 ........ 95 Figure 4.4. Downscaled mean daily September to May precipitation rate (left) and temperature (right) by MicroMet for the reference and future periods over the whole area ................. 98 Figure 4.5. Monthly averaged SWE [m] over entire domain from CESM1-WRF forced into SnowModel for four ensemble members during historical period (1979-1992). Red line with circle marker shows the average ................................................................................ 99 Figure 4.6. The averaged solid precipitation rate change between future simulations and reference period for (a) mid-century and (b) end of the century using statistical downscaled (SD) and dynamical downscaled (DD) forcing ............................................................... 100 Figure 4.7. Averaged snow water equivalent (SWE) changes between future simulations and reference period for (a) mid-century and (b) end of century using statistical downscaled (SD) and dynamical downscaled (DD) forcing data........................................................ 102 xxiv Figure 4.8. Mean simulated snow depth and relative difference of snow depth (a) and mean simulated snow water equivalent and relative difference of SWE (b) for the reference period and two projections for mid-century and end of the century forced by statistical downscaled (SD) and dynamical downscaled (DD) forcing data .................................... 104 Figure 4.9. Relative changes (%) of the snow depth in the area for two forcing data in two different periods depend on the elevation. The blue represents lower and the red color shows the higher frequency of snow depth changes. The black line is the average line . 105 Figure S1.1. Scatter plots of daily mean temperature derived from WRF and NLDAS against observations for 2013-2014 (top) and 2015-2016 (bottom) at four station locations ...... 123 Figure S1.2. Topography (km) and outline of the outer at 7 km and inner WRF domains at 1 km. .................................................................................................................................. 124 Figure S1.3. Simulated time series of daily SWE for calibration period winter 2013-2014 (top) and snow depth for evaluation period for winter 2015-2016 (bottom) against observed data for four stations ................................................................................................................ 125 Figure S1.4. Comparison of simulated SWE [mm w.e.] by SnowModel averaged over domain at East-Creek station for winter 2013-2014. Light-blue line represents SWE forced by hourly data, and dashed-black line shows simulated SWE forced by 6-hourly data ....... 126 Figure S1.5. Comparison of simulated SWE (Left) and snow depth (Right) using rain/snow threshold of 1.3˚C (Light blue) and 2˚C (Dashed black) by Alpine3D at point scale for accumulation period winter 2015-2016 ........................................................................... 126 xxv Figure S1.6. Spatial scale comparison of simulated SWE [m w.e.] by Alpine3D, accumulated over accumulation period 2015-2016 at Kokanee Glacier using different thresholds of 1.3˚C (A) and 2˚C (B ....................................................................................................... 127 Figure S1.7. Spatial comparison of simulated snow density by SnowModel at Kokanee glacier using multilayer (A) and one-layer (B) methods and their difference (multilayer-onelayer) (C) for winter accumulation 2015-2016 .......................................................................... 127 Figure S1.8. Spatial comparison of simulated SWE [mm w. e.] by SnowModel with multilayer (A) and one-layer (B) averaged over accumulation period 2015-2016. Difference between both methods illustrates in (C)......................................................................................... 128 Figure S1.9. Time series of daily snow density [kg m-3] simulated by SnowModel (SM) and Alpine3D (ALP) for winter 2015-2016. ......................................................................... 128 Figure S2.1. Flowchart of model simulations with considering wind effects and without wind. ......................................................................................................................................... 132 Figure S3.1. Mean temperature of statistically downscaled (SD) and dynamical downscaled (DD) data from GESM1 for three 15-year periods (historical, mid-century and end of century) averaged over Columbia Basin based on RCP 8.5 ............................................ 134 Figure S3.2. Averaged snow depth time series over 10 stations compared with the DD-forced snow depth (mean of 4 ensembles) for the historical period 1979-1995. ....................... 134 Figure S3.3. Projected snow depth change (%) using statistically downscaled (SD) and dynamical downscaled (DD) data from CESM1 for mid-century and end of century over Columbia Basin ............................................................................................................... 135 xxvi Figure S3.4. Simulated incoming short wave (SW) (Right) and long wave (LW) (Left) radiations by SnowModel forced by SD and DD datasets for period 1979-1994 ........... 136 xxvii Chapter 1: Introduction This dissertation investigates the important factors affecting the simulation of snow accumulation in alpine terrain in terms of input data, model complexity and the effect of wind induced snow transport processes on winter Glacier Mass Balance. The dissertation presents results of high resolution physically-based snow model runs forced by High Resolution Meteorological (HRM) data applied on four glaciers over the Columbia Basin. Furthermore, the future changes in snow depth and SWE by the end of this century over this basin is presented using available statistically and dynamical downscaled data. 1.1. Importance of glaciers and snow modelling Surface mass balance is a direct response of a glacier to meteorological conditions that cause a glacier to accumulate or lose ice or snow over a specified time period. Surface accumulation occurs by snow and ice precipitation, refreezing liquid precipitation within the snowpacks, avalanches, and windblown snow. Processes that cause the glacier to lose mass include surface melt, surface meltwater runoff, sublimation, avalanching and snow transport by wind. The close relation of snow and ice melt with air temperature (e.g. Ohmura, 2001) results in sensitivity of the cryosphere to temperature fluctuations. According to the IPCC (2019), the global annual average of human-induced warming in 2017 reached about 1˚C above preindustrial levels and increases at a rate of 0.2˚C per decade, within an expected warming of 34˚C by 2100 (Rogelj et al., 2016). Alpine regions are particularly sensitive to climate fluctuations. Precipitation patterns and variability of winter temperatures can significantly affect spring and summer river flows in cold regions (Barnett et al., 2005; Adam et al., 2009). 1 The impacts of climate change are already evidenced by snowpack decline and the rapid volume losses of many alpine glaciers (e.g. Bauder et al., 2007; Larsen et al., 2007; Schiefer et al., 2007; Zemp et al., 2008; Bolch et al., 2010; Menounos et al., 2019). Glaciers and alpine snow cover provide an important water resource for Western Canada, and they too are affected by climate change. Expected changes in glaciers and snow cover that are affecting the runoff regime of alpine rivers result in a significant reduction in summer runoff in lowlands due to glacier retreat (Moyer et al., 2016; Shrestha et al., 2017; Brahney et al., 2017). Such changes can significantly influence the hydrology of glacierized basins, which affects water availability for human consumption, agriculture, tourism, and hydroelectricity generation, posing challenges for long-term water management (Khadka, 2013). Direct observations provide one approach to evaluate the links between meteorological conditions that favour glacier mass gain or loss (e.g. Brock et al., 2000; Hock and Holmgren, 2005; Andreassen et al., 2008). Prior work examines the relations between climate and winter accumulation using statistical (e.g. Trachsel and Nesje, 2015) and geodetic methods (e.g. Sold et al., 2013; Beedle et al., 2015; Pelto et al., 2019). In recent decades, scientists have developed energy-balance models to simulate glacier mass balance and snow cover forced by meteorological and climatological factors. Many studies model surface melt for alpine glaciers (e.g. Hock and Holmgren, 2005; Lott and Lundquist, 2008; DeBeer and Pomeroy 2009, 2017), whereas far fewer examine seasonal snowpack accumulation (e.g. Winstral et al., 2013; Reijmer et al., 2008; Franz et al., 2008; Liston et al., 2008; Mernild et al., 2005). Energybalance models are physically-based and improve simulation by considering the spatial and temporal contributions of energy into the snowpack (e.g. radiation, sensible and latent heat exchange, ground and rain heat fluxes) that result in melting and its response to climate changes 2 and surface conditions. Such studies are important for accurate prediction of summer glacier mass balance changes, especially in mountainous areas due to the complexity of topography and its influence on radiation (Arnold et al., 1996). However, the calculation of winter accumulation, which is influenced by wind and topography, is simplified within these models. The distribution of snow in mountainous regions is variable due to spatial differences in factors such as topography, shading and redistribution by wind that cause differences in accumulation of snow in alpine terrain (Winstral and Marks, 2002; Winstral et al., 2013), leading to variations in glacier mass balance (McGrath et al., 2018; Sold et al., 2013). Wind affects snow accumulation by transporting it to leeward and sheltered locations; wind also affects sublimation and evaporation processes (Gordon et al., 2006; Strasser et al., 2007). Large-scale orographic precipitation and local scale processes also play key roles in affecting broad scale patterns in snow accumulation (Houze, 2012; Mott et al., 2014). Input data uncertainties in terms of both spatial variability and uncertainties in parameterizations can lead to a significant loss of confidence in key simulation results related to the surface energy and mass budgets (Sauter and Obleitner, 2015). Coarse resolution meteorological model outputs or data sets, for example from Global Climate Models or model reanalysis datasets, are often unable to properly resolve the topographical details and processes important for glacier mass balance. Dynamical downscaling uses a physically-based atmospheric model that solves equations for conservation of mass, momentum and energy in the atmosphere and typically the land surface to provide a finer resolution of forcing data. Dynamical downscaling is capable of providing detailed information on the interaction of the land and the atmosphere, especially orographic precipitation, but is computationally expensive. Recently, meteorological data from regional atmospheric models have been used to directly 3 drive surface energy balance models, while surface fluxes can directly influence the atmospheric modelling when integrated in time (coupled models) (e.g. Mölg and Kaser, 2011; Gallée et al., 2012; Collier et al., 2013; Vionnet et al., 2014) or used indirectly by adding a simulated wind field to observations as input to the snow models (e.g. Mott et al., 2008, 2010). Coupled models are proper tools to study snowpack-atmosphere interactions; however, it is important that the modelling approach be consistent with the purpose of its application and required temporal and spatial resolution in the simulation (Kirnbauer et al., 1994; Hock and Jansson, 2005). 1.1.1. Snow models Many different snow models exist to calculate the vertical and spatial development of snow cover (e.g. Pomeroy et al., 1993; Liston and Sturm, 1998; Déry and Yau, 1999; Winstral and Marks, 2002; Lehning et al., 2006; Marsh et al., 2020; Vionnet et al., 2020), but only a few studies have applied these seasonal snow models to simulate glacier mass balance. These snow models vary in complexity and incorporate different approaches to model ablation, and accumulation with different requirements in terms of input and calibration data. According to Marsh (1999), Liston (2004) and Liston and Elder (2006), most of these models have physically-based representations of processes, which can be categorized as one-dimensional, two-dimensional, or three-dimensional (Bernhardt, 2008). A) One-dimensional models. These models are multi-layer snowpack models. In many cases these are used to assess avalanche or flood risk. These models are developed to accurately calculate the vertical stratification of the snowpack and the snow metamorphism and its association with mechanical properties such as thermal conductivity, viscosity and strength 4 (Spreitzhofer et al., 1993). Complex and physically-based formulations are necessary for a proper description of the above variables. A numerical solution is possible only if extensive input information is available. As such, their application is limited to well instrumented locations. Some of the 1-dimensional models are: CROCUS (Brun et al., 1989; Brun et al., 1992), SNOWPACK (Lehning et al., 2002; Bartelt and Lehning, 2002) and SNTHERM (Jordan, 1991), among others. Météo-France created the numerical model CROCUS to predict the evolution and stability of snowpack in relation to avalanches. To support avalanche warnings in Switzerland, SNOWPACK was created by the Swiss Federal Institute for Snow Research and Lawine Research (SLF). SNTHERM was created to predict runoff. B) Two-dimensional models are typically used as subroutines in atmospheric models and soil vegetation atmosphere transfer schemes, some of which are stand-alone models (e.g. the AMUNDSEN model; Strasser, 2007). Two-dimensional models conventionally have one layer and do not explicitly simulate the metamorphism of snow crystals. Lateral wind or gravity transport is often ignored, and a sub-scale snow distribution is based on the mean changes to albedo due to snow evolution (Liston, 2004). C) Three-dimensional models consist of at least two parts, one calculating the evolution of snow columns while the other simulates the process of snow transport. The complexity of each component determines whether a model can be used to simulate the seasonal snow cover or whether it only identifies a short-term snow transport event (e.g. during storm events). The ability of the models to simulate seasonal or event-based snow cover depends on the complexity of each component. The snow drift models are generally quite complicated and computationally intensive. They are often used for basic research, where they are split between event-based (more complicated) and seasonal models (simpler). For instance, a snow drift 5 routine was added in Lehning et al. (2002) to the Advanced Regional Projection System (Xue et al., 2000) atmospheric model. Compared to event-based models, seasonal snow-transport models usually have a medium level of complexity. They usually use a first order transport physics approximation and a single snow evolution model layer. The Prairie Blowing Snow Model (PBSM) (Pomeroy et al., 1997) is the first recognized seasonal and physical snow transport model. SnowTran-3D (Liston and Sturm, 1998; Liston et al., 2007), which is a part of SnowModel, is one of PBSM's successors. ALPINE-3D (Lehning et al., 2006) also is a three-dimensional model, which is a combination of SNOWPACK and snowdrift models. Recently, a 3‐D advection‐diffusion blowing snow model was presented by Marsh et al. (2020). A physically based mass balance model, iSNOBAL (Marks et al., 1999), was used for snow cover simulation and implemented in the Canadian Hydrological Model (CHM) (e.g. Marsh et al., 2020b; Vionnet et al., 2020). This model utilizes a variable resolution unstructured mesh that significantly reduces the total number of computational elements and produces heterogeneous snow covers without the need for calibration. The pattern of snow accumulation is regulated by (a) preferential deposition of snowfall, (b) snow redistribution by wind, and (c) other snow transports such as avalanches, and snow slides (Kuhn, 1995; Lehning et al., 2008). Models that simulate such processes are subject to large uncertainties due to the lack of accurate solid precipitation data, especially in mountainous areas (Sevruk, 1985), and makes it difficult to simulate snow accumulation (Dadic et al., 2010). The one-dimensional snow models are able to simulate snowpack evolution near weather stations (Etchevers et al., 2002). However, the results are less reliable for 2- or 3-dimensional spatially distributed models, especially in areas with complex or alpine terrain (e.g. Liston, 2006). This is thought to be due to errors and uncertainties in forcing data, as well as the 6 simplification of model formulas (e.g. Lehning et al., 2006). Distributed snow models compute snow accumulation using one of three methods: using precipitation data from nearby weather stations using a precipitation lapse rate to adjust data for all elevations; from gridded precipitation data sets; or from climate model outputs utilizing a threshold temperature to differentiate liquid and solid precipitation (Hock, 1999; Finger et al., 2012; Salzmann et al., 2012). In addition, some snow models also have a simple parameterization for snow redistribution processes based on terrain parameters (e.g. Huss et al., 2009; Mayr et al., 2013). In alpine environments, a high frequency of blowing snow conditions with wind speeds higher than 3 m s-1 at 10 m height above the surface (Liston and Sturm, 1998) and a lack of vegetation to stabilize snow cover (Pomeroy et al., 1997; Pohl et al., 2007) are the main reasons for large wind impacts on the snowpack. Due to the complexity of wind fields in steep terrain, drifting snow is implemented in the models mainly as a function of terrain parameters (Liston and Sturm, 1998; Winstral and Marks, 2002; Winstral et al., 2013). 1.1.2. Downscaling methods Many glaciers are located in remote areas, which makes it difficult to measure snow properties by direct observation. A modelling approach thus provides one method to estimate glacier mass balance. However, in data sparse regions it is challenging to find alternatives to stations for forcing data. Global Climate Models (GCMs) are too coarse to address highresolution weather conditions at a local scale, thus limiting their applicability for snowpack climate impact studies (Radić and Hock, 2006; Radić and Clarke, 2011). Statistical or dynamical downscaling approaches can be used to fill the gap between global climate model (GCM) resolution and the study area's regional scale. 7 Statistical downscaling uses estimates derived from empirical relationships between GCM simulated large-scale atmospheric variables and observed local-scale meteorological variables with an assumption that the relationship between variables remains unchanged over time. Statistical downscaling is computationally efficient. However, its main disadvantage is its assumption that climate dynamics of the past and present will be the same in the future. Furthermore, this approach can typically be generalized only to precipitation and air temperature due to the limited availability of observed data (Maraun et al., 2010; Teutschbein et al., 2011). Dynamical downscaling, on the other hand, uses a Regional Climate Model to produce a finer-scale climate projection. This method can represent detailed information about interactions between land and atmosphere, including orographic precipitation. The models used are generally numerical weather prediction models, often modified to allow for long-term climate simulations. This approach, however, does not require a stationarity assumption and accounts for regional effects on the climate system such as interaction of air masses with mountains. However, computationally, dynamical downscaling costs considerably more than statistical downscaling (Erler et al., 2017). Having downscaled data from climate models makes it possible to project future snow cover changes over the study area by the end of century. 1.2. Dissertation objectives In recent decades, snow models have been developed and improved in light of enhanced understanding of snow process physics and growth in computational power, leading to higher complexity of distributed physics-based snow models. Although snow models have been developed and applied in many regions of the globe, some aspects of modelling are not well 8 investigated. Melting processes are centred in most of the research, while winter accumulation affects summer ablation directly and has a similar significance. In addition, mostly the observed data from weather stations are used to force snow models. Here, not only I am using all high-resolution meteorological fields to force the snow models, I examine different models to assess the effect of model complexity on winter accumulation. I also investigate how wind affects glacier winter mass balance in this region, to determine whether or not the wind effect is important to consider given the fact that such processes require high computational and time demands. Furthermore, future snow cover changes over the Columbia Basin by the end of the 21st century are projected based on the RCP 8.5 scenario. Four individual glaciers over the Columbia Mountains of British Columbia are selected as test sites to simulate the winter accumulation and investigate the effect of wind on the winter mass balance. Furthermore, an assessment of global warming on snow cover and water resources in cold regions is a major challenge for understanding climate change on the regional scale, which is of great interest in modelling impacts and informing strategies for adaptation. Consequently, I project changes in future snow cover over the Columbia Basin by the end of the century. Specific objectives of this research are to: 1) Examine the importance of meteorological input data used to force the snow models, on winter accumulation in the absence of weather station data, 2) Investigate the importance of snow evolution model complexity on winter mass balance, 3) Determine the effect of wind on glacier winter mass balance in terms of snow erosion, deposition and sublimation, and 9 4) Assess changes in snow depth and SWE over the Columbia Basin by the end of the 21st century using a snow model driven by available dynamically and statistically downscaled data. 1.3. Study region In this study, I estimate winter mass balance of glaciers and snow change over the Canadian portion of the Columbia Basin (Fig. 2.1). The Canadian portion of the Columbia Mountains is located in the southern interior of British Columbia and is bounded by the Rocky Mountain Trench (east), the Columbia River (south), the Interior Plateau (west) and the Fraser River (north) with peaks above sea level (a.s.l.) of 3000 m with complex topographic features that influence its climate significantly. The Upper Columbia Basin comprises about 55,000 km2 of rivers or lakes, along with population centres of varying sizes. In late summer, snow melt induced river discharge decreases rapidly, which plays a key role in the region's water supply. The heavy winter snowfall is controlled by mid-latitude cyclones with wet and mild air masses from the Pacific Ocean, which are intercepted by the Columbia Mountains. Characteristics of this region include high annual precipitation, deep snow accumulations and relatively moderate winter temperatures. According to Bolch et al. (2010), glacier area (in 2005) of the British Columbia (BC) Southern Interior and Southern Rocky regions were reported to be 1910.4 and 1351.7 km 2, respectively. The Columbia Basin includes a part of the Southern Rocky region in this inventory. The recorded number of glaciers over the whole area is 3575 glaciers. This study estimated the glacier area changes from 1985-2005 as -342.2±98.9 and -235.3±65.2 km2 for the Southern 10 Interior and Southern Rocky areas, respectively, which equates to -15.2±4.4 % and -14.8±4.1 % for these areas, respectively. Estimates of glacier areal loss is about 0.7 % yr-1. 1.4. Outline of chapters This dissertation contains three main chapters. In chapter 2, I answer the first and second research objectives: examining the importance of meteorological forcing data and complexity of the snow models, on winter accumulation in the absence of weather station data. Two different forcing data sets from the Weather Research and Forecasting model (WRF) and the North American Land Data Assimilation System (NLDAS), are used to force two snow models with different complexity (SnowModel and Alpine3D). The results are evaluated based on insitu observations conducted by the UNBC cryosphere group over four glaciers. In addition, the integrated snow depth over the domain is compared with LiDAR measurements. The importance of forcing data, and snow model complexity are investigated by comparison between the results of each simulation for winter 2015-2016. In the third chapter, the effect of wind on the winter accumulation is explored utilizing the second chapter results, in terms of selecting forcing data and the snow model. The results of this chapter answer the third research question: determining the effect of wind on glacier winter mass balance in terms of snow erosion, deposition and sublimation. The wind influence on winter mass balance is demonstrated by implementing two simulations for each glacier: with and without consideration of wind effects. The results show details of snow deposition and sublimation over four glaciers, as well as the magnitude of such effects on winter accumulation. 11 Chapter 4 projects snow cover changes under climate change scenario RCP 8.5 by the end of the century to respond to the fourth research objective: assessing changes in snow depth and SWE over the Upper Columbia Basin by the end of this century using dynamically and statistically downscaled data. Forcing data for the reference period (1979-1994) and future projections for two 15-yr periods (2045-2059, 2085-2099) were derived from downscaled CESM-WRF for four different ensembles. In addition, for the same periods, statistically downscaled data derived from CESM over the region are used to force SnowModel to project future snow cover over the whole Columbia Mountains. The effect of each forcing on the snow cover projections and the related uncertainties are discussed. Chapters 2, 3, and 4 are each based on manuscripts. The paper derived from Chapter 2 is accepted in the journal Hydrological Processes. Chapters 3 and 4 are being prepared as manuscripts for submission to the Journal of Glaciology and International Journal of Climatology, respectively. Each manuscript consists of an introduction, data and methods, results, discussion and conclusions sections. The appropriate permissions have been obtained to replicate the manuscripts herein. As the first author of the papers, I am the primary contributor to the articles including modelling, data processing and analysis, and preparing drafts for submission. Other co-authors include my co-supervisors (Drs. Menounos and Jackson) as leaders for the project and helping to edit the manuscripts. Two other co-authors provided some of the forcing data (Dr. Erler) as well as field data (Dr. Pelto). Chapter 5 provides a summary of the main results obtained in the dissertation, the study limitations, and some suggestions for future research. 12 Chapter 2. The Role of Meteorological Forcing and Snow Model complexity in Winter Glacier Mass Balance Estimation, Columbia River Basin, Canada Mortezapour, M.1, Menounos, B.1, Jackson, P. L.1, Erler A. R.2, and Pelto, B. M.1 1. Natural Resource and Environmental Institute, Prince George, British Columbia, Canada 2. Aquanty, Waterloo, Ontario, Canada Accepted by Hydrological Processes and published Hydrological Processes. 2020;1–19. https://doi.org/10.1002/hyp.13929 13 Abstract Glaciers are commonly located in mountainous terrain subject to highly variable meteorological conditions. High resolution meteorological (HRM) data simulated by atmospheric models can complement meteorological station observations to assess changes in glacier energy fluxes and mass balance. We examine the performance of two snow models, SnowModel and Alpine3D, forced by different meteorological data for winter mass balance simulations at four glaciers in the Canadian portion of the Columbia Basin. The Weather Research and Forecasting model (WRF) with resolution of 1 km and the North American Land Data Assimilation System (NLDAS2) with ~12 km resolution, provide forcing data for the two snow models. Evaluation is based on the ability of the snow models to simulate snow depth at both point locations (automated snow weather stations; ASWSs) and over the entire glacier surface (airborne LiDAR surveys) during the 2015-2016 winter accumulation. When forced with HRM data, both models can reproduce snow depth to within ± 15% of observed values. Both models underestimate winter mass balance when forced by HRM data. When driven with WRF data, SnowModel underestimates winter mass balance integrated over the glacier area by 1 and 10%, while Alpine3D underestimates winter mass balance by 12 and 22% compared with LiDAR and stake measurements, respectively. The overall results show that SnowModel forced by WRF simulated winter mass balance the best. 14 2.1. Introduction Runoff from high mountain areas supplement many river systems in the world. Precipitation and winter temperature patterns and variability can significantly affect spring and summer flows in cold regions (Barnett et al., 2005; Adam et al., 2009). Change in glacier area, for example, influence the hydrology of glacierized basins, which affects water availability for urban consumption, agriculture, and hydroelectricity generation, posing challenges for longterm water management (Khadka et al., 2013). Climate controls the global distribution of alpine glaciers. A glacier is affected by the interaction of weather systems with complex mountain topography (Hofer et al., 2015, Gurgiser et al., 2013; Khadka, 2013); this interaction affects energy and mass fluxes that control the nourishment and depletion of these ice masses. Direct observations provide one approach to evaluate the links between meteorological conditions that favour mass gain or loss (e.g. Brock et al., 2000; Hock and Holmgren, 2005; Andreassen et al., 2008). Prior work examined the relationship between climate and winter accumulation using statistical (e.g. Trachsel and Nesje, 2015) and geodetic (e.g. Sold et al., 2013; Beedle et al., 2015; Pelto et al., 2019) methods. Models that simulate snow accumulation and melt can improve our understanding of physical processes that lead to mass accumulation for alpine glaciers, especially those where sublimation and calving processes are negligible. Many studies model mass loss during the ablation season (e.g. Blӧschl et al., 1991; Arnold et al, 2006; Hock 2003; Lefebre et al., 2003; Hock and Holmgren, 2005; Lott and Lundquist, 2008; DeBeer and Pomeroy, 2009, 2017) whereas fewer simulate winter accumulation (e.g. Winstral and Marks, 2002; Winstral et al., 2013; Reijmer et al., 2008; Franz et al., 2008; Liston et al., 2008; Mernild et al., 2005). Over 15 the last decade, several studies have been conducted to simulate winter accumulation using models of different complexities (Bernhardt et al., 2009; Warscher et al., 2013; Mott et al., 2010; Vionnet et al., 2017). In most studies of winter accumulation on alpine glaciers, physically-based models are forced with observed weather station data. In some studies, however, meteorological data from regional atmospheric models have been used to directly drive surface energy balance models (e.g. Mölg and Kaser, 2011; Gallée et al., 2012; Collier et al., 2013; Vionnet et al., 2014). While coupled models are appropriate tools to investigate how the snowpack interacts with the atmosphere, uncertainties introduced by downscaling atmospheric data directly to the glacier surface may affect spatial and temporal distribution of meteorological variables in the snow models which may affect the snow model results (Machguth et al., 2006a). On the other hand, the complexity of a given modelling approach should be consistent with its application objectives and with temporal and spatial resolution required in the simulation (Kirnbauer et al., 1994; Hock and Jansson, 2005). Our study goal is to address two primary questions: how do differences in the meteorological forcing (HRM data) in modelling snow processes affect the simulated glacier mass balance in regions where direct observations are absent? How does the complexity of snow models affect the simulation of winter mass balance? In this study, we use meteorological fields from mesoscale atmospheric models to simulate winter mass balance. We evaluate simulated temporal and spatially distributed time series of snow depth and SWE using two distributed snow models with different complexities. The paper does not consider the effects of wind transport. We chose four glaciers in British Columbia with different characteristics to evaluate the ability of these snow models to simulate winter accumulation. 16 2.2. Study area and data used in this study 2.2.1. Study area The four glaciers studied lie in the Columbia Mountains of British Columbia, Canada (Fig. 2.1), which are drained by the Fraser, South Thompson, Columbia, and Kootenay rivers (Holland, 1976). These glaciers are: Zillmer (52°41' N, 119°35' W), Conrad (50°49' N, 116°55' W), Nordic (51° 25' 13'' N, 117° 42' 13'' W) and Kokanee (49° 45' N, 117° 8' 30'' W) that have mean elevations of 2094, 2301, 2718 and 2600 m above sea level (m a.s.l.), respectively (Fig. 2.1). The Conrad, Zillmer, Nordic, and Kokanee glaciers cover an area of 11.45, 5.43, 3.39, and 1.8 km2, respectively. In this study, we used a domain larger than each main glacier (the main glaciers plus surrounding areas) to evaluate models in glacier-free areas, especially when comparing with LiDAR data. The topography surrounding each glacier is characterized by steep alpine terrain (elevation range between 420 to 3700 m a.s.l.). The prevailing winds are westerly over this area however, the final wind pattern is strongly affected by the complex local topographic features. Based on the data from the Parameter-Elevation Relationships on Independent Slopes Model (PRISM; 4 km resolution) dataset the annual average temperatures are below 0°C and as cold as -4°C in high elevation areas in the central and northern parts of the Columbia basin for years 1961-1990 (Murdock & Werner, 2011; Daly et al. 1994). The average relative humidity for elevations higher than 1500 m a.s.l. is 74% for the same data (1961-1990 normals). The mean winter precipitation over the entire Columbia basin (19002000) is reported to be 350 mm based on the data derived from ClimateBC for areas below 1500 m a.s.l. (Carver, 2017). 17 Figure 2.1. Top) Study sites: Zillmer, Nordic, Conrad, and Kokanee glaciers and the nearest weather station locations to each glacier (Table 2.1). The blue outline represents the Canadian portion of the Columbia River basin. Bottom) Topography of each individual glacier with LiDAR footprint boundaries (Yellow outlines) and main (ON) glacier boundaries (Red outlines) 18 2.2.2. Input and evaluation data a) Modelled meteorological forcing data Below, we describe the two datasets used to force the snow models. The North American Land Data Assimilation System (NLDAS; V002) is a multimodel land modelling and assimilation system. The non-precipitation land-surface forcing fields for NLDAS2 are derived from the North American Regional Reanalysis (NARR) fields. The precipitation field is a result of a temporal disaggregation of a gauge analysis of Climate Prediction Center (CPC) daily precipitation which is directly performed on the NLDAS2 grid using an orographic adjustment based on the PRISM climatology. The hourly NLDAS2 data have a spatial resolution of 0.125˚ (latitude/longitude) (Cosgrove et al., 2003; Mitchell et al., 2004). Hourly surface meteorological variables from both datasets used to force the snow models include precipitation intensity (mm hr-1), temperature (˚C) and relative humidity (%) at 2 m height, wind speed (m s-1) and wind direction (degrees) at 10 m height, as well as incoming shortwave and longwave radiation (W m-2). The other meteorological driving data used to force the snow models were derived from high resolution regional climate simulations that have been conducted with WRF – the Weather Research and Forecasting model (Skamarock and Klemp, 2008), version 3.6.1. Due to WRF high computation time, a nested domain was used with an outer/inner domain at 7 km and 1 km horizontal resolution (one-way coupling) (Fig S1.2). The outer domain covers most of British Columbia and parts of Alberta, while the inner domain covers only the Canadian Rocky Mountains and the Columbia Mountains, which are the primary regions of interest for this study. The outer domain is forced by meteorological and sea-surface temperature data from the ERA-Interim reanalysis product (Dee et al., 2011) at 6-hour intervals. The WRF simulations 19 were run for an entire model year starting at the beginning of August and integrated continuously until the end of August of the following year. The first month of simulations was used to spin up the model and not analyzed in this study. Furthermore, nudging was applied within nine grid cells of the outer domain boundary and weak spectral nudging was applied above the boundary layer. The model was configured with 42 hybrid terrain-following vertical layers with a model top at 50 hPa (approximately 20 km above sea level) and a weak damping layer starting 8 km below the model top. The following physical parameterization schemes have been employed: the Rapid Radiative Transfer Model for GCMs (RRTMG) radiation scheme (Iacono et al., 2008), the Thompson micro-physics scheme (Thompson et al., 2008), the Noah-MP land surface model (Niu et al., 2011), the MYNN3 planetary boundary and surface layer scheme (Nakanishi and Niino, 2006) and the Grell-3 ensemble cumulus parameterization with subsidence spreading (Grell and Dévényi, 2002); no cumulus scheme was used in the inner domain, since at 1 km convection is explicitly resolved (Liu et al., 2011). The column lake model (Gu et al., 2013), topographic shading and slope effects on radiation, as well as cloud feedbacks on radiation have also been enabled (Wang et al., 2009). The physics configuration is similar to the configuration used by Erler et al. (2015), but a newer version of WRF was employed here. b) Automated Snow Weather Station (ASWS) data We used data from four ASWS (Fig. 2.1 and Table 2.1) to evaluate the HRM datasets and to calibrate the two snow models and optimize their performance by comparing the time series of snow depth and SWE obtained from nearby ASWS sites with simulated data. These weather stations are distributed throughout British Columbia above valley bottoms from the 20 subalpine to treeline at altitudes between 700 to 2200 m a.s.l.; their measurements include snow depth, SWE, cumulative precipitation, and air temperature. We found a few short data gaps and a small proportion of data appeared to be outliers. The four ASWS measurement datasets were preprocessed by examining data (< 5 th percentile and > 95th percentile) to remove outliers and values resulting from sensor malfunction; linear interpolation was used to fill small data gaps (<5 hours). Outliers and gaps were detected by graphical display of data. End of season SWE was evaluated using observed data from ASWSs at the point scale, and by manual probing of snow depths over glaciers. Table 2.1. Automated snow pillow station characteristics. For the Glacier-NMF station, temperature and precipitation are derived from the BC Hydro network (station ID: FDL), while snow depth and SWE are derived from Automated snow pillow network from Ministry of Environment (station ID: 2A17). Station Station Name ID 1E08P Network Elevation name Azure River Parameters1 (m) ENV-ASP 2 ENV-ASP 2 Period of Frequency record 1630 T, P, S, 2003-2020 East Creek 1960 T, P, S, Hourly 2A17 Glacier Np ENV-ASP 1850 S, SWE 2003-2020 Hourly 2D14P Fidelity Redfish Creek 15 km to Conrad 1963-2020 Daily 20 km to Nordic Mt Fidelity FDL 8.7 km to Zillmer SWE 2 to Glacier SWE 2D08P Distance 3 BCH ENV-ASP 2 2080 T, P 1982-2020 Daily T, P, S, 2003-2020 Hourly SWE 1 T: Temperature, P: Precipitation, S: Snow-depth, SWE: Snow water equivalent 2 Ministry of Environment - Automated Snow Pillow Network; 3 BC Hydro 3.4 km to Kokanee c) In situ snow measurements We measured winter surface mass balance at the four glaciers since 2014 (Pelto et al., 2019) using the glaciological method (Cogley et al., 2011). The measurement networks for all four 21 glaciers including snow depth measurement locations, and snow pit/core locations are shown in Fig. 2.2. The glaciological field campaign occurred for the winter in late April to early May and for summer in late August to early September (Table 2.2). We used snow pit measurements with a 100 cm3 box cutter and snow cores using a 7.25 cm diameter Kovacs corer to quantify the mass balance for each glacier. The spring snow density was measured at low, middle, and high elevations, three locations, for each glacier (Fig. 2.2) and above the transient snow line for the fall season. Figure 2.2. Measurement networks for the (a) Zillmer, (b) Nordic, (c) Kokanee, and (d) Conrad glaciers. Snow depth measurement locations (black dots), and snow pit/core locations (blue squares) are identified. Contours represent the elevation bands (m a.s.l.) over glaciers. This figure has been adapted from Pelto et al. (2019). 22 Pelto et al. (2019) observed marginally denser snow in snow pits a than measured in snow cores (0.2±5.7 %). To calculate glacier-wide mass balance, we used the profile method (Escher-Vetter et al., 2009), applied over 100 m hypsometric elevation bins. The area–altitude distribution of a given glacier was obtained using our annual late summer LiDAR measurements. Table 2.2. Glaciological visits (number of locations, dates) and geodetic acquisition information (dates and point density) for each glacier. Field dates are the median date of glacier visit. (Derived from Pelto et al., 2019) Glacier Summer glac. (m/d/y) N Summer ALS (m/d/y) Points (m-2) Winter glac. (m/d/y) N Winter ALS (m/d/y) Points (m-2) Zillmer 8/23/2015 23 10/3/2015 2.75 4/14/2016 46 4/18/2016 3.69 Nordic 8/31/2015 11 9/11/2015 1.99 5/2/2016 28 4/17/2016 3.21 Conrad 9/5/2015 9 9/12/2015 1.35 4/26/2016 44 4/17/2016 2.45 Kokanee 8/27/2015 11 9/12/2015 1.04 4/19/2016 33 4/17/2016 2.77 N: Number of observation locations d) LiDAR measurements The availability of high resolution Light Detection and Ranging (LiDAR) data made it possible to evaluate simulated integrated snow depth over the LiDAR footprint and on each glacier’s area for the accumulation period 2015-2016. We used LiDAR to monitor elevation change over each glacier at a resolution of 1 m and accuracy of better than 15 cm (Arnold et al., 2006; Krabill et al., 1995). LiDAR data were collected twice per year for each glacier (Table 2.2). To estimate mass balance changes during a year, two surface elevation models made from the LiDAR surveys were post-processed, co-registered, bias corrected, and differenced. Full details of the LiDAR surveys can be found in Pelto et al. (2019). The observed snow density was used to convert changes in snow elevation to changes in mass, which was then integrated over each glacier to estimate winter balance. Since the LiDAR maps of snow depth include the effects of ice dynamics (emergence and submergence of the ice surface), they 23 do not represent the spatial distribution of snow depth over the glaciers. These data, however, can be used to assess the performance of snow models for average snow depth over the glaciers. Since the net vertical velocity over the glacier is zero, we can average the LiDAR elevation difference (ΔH) over the spatial domain of the glacier to determine the average snow depth, to compare with the modelled average snow depth. To make this comparison, we downsampled the 1-m LiDAR DEMs to the spatial extent of the model grid cells at 30 m resolution. The model integrated surface mass balance over the glacier surface was estimated for each glacier and compared to the glacier-wide surface mass balance obtained from the LiDAR measurements (Hedrick et al., 2015; Musselman et al., 2015). 2.3. Snow models description and simulation 2.3.1. SnowModel SnowModel (Liston and Elder, 2006) is a physically-based snow model, which can be applied to simulate snow evolution in different climatic settings. This model can be run on grid increments of 1 to 200 m with a temporal resolution of 30 minutes to one day (Liston et al., 2018). SnowModel consists of four main sub-models: MicroMet, EnBal, SnowPack, and SnowTran3D. (1) MicroMet is a submodel for data preprocessing and distributing meteorological forcing data (Liston and Elder, 2005). MicroMet interpolates and distributes forcing data over the domain using Barnes objective analysis (Koch et al., 1983). The interpolation process is started by interpolating the data from spatially irregular stations to a regular grid, using a scheme developed by Barnes (1964, 1973). Then a relationship between meteorological parameters and terrain is used to distribute the data set over the domain (Liston and Elder, 2005). In MicroMet the elevational gradients of temperature and precipitation (i.e. 24 lapse rates) vary monthly through the year for temperature. Precipitation lapse rates are calculated based on a non-linear function of elevation difference. The monthly temperature lapse rate and precipitation adjustment factors can be adapted by the user. (2) EnBal is a submodel to simulate energy fluxes to and from the snowpack (Liston, 1995; Liston et al., 1999b). In addition, ablation by non-blowing snow sublimation is calculated in EnBal and used to adjust snowpack depth. (3) SnowPack is a simple, user defined multi-layer module to simulate snowpack evolution based on density changes (Liston and Hall, 1995). SnowModel uses the surface energy balance equation (Liston, 1995; Liston et al., 1999b; Liston and Elder, 2006), and snowpack evolution occurs based on density changes due to compaction from changes in snow temperature and the weight of overlaying snow. SnowPack simulates changes in the snowpack in response to the energy fluxes and precipitation input given by MicroMet. The energy from MicroMet alters snow density by decreasing snow depth and redistributing meltwater through the snowpack until a maximum snow density is reached. Some internal processes are taken into account by SnowModel such as refreezing, but the model does not consider melting from internal ice deformation, changes in drainage systems, geothermal effects, and sub-glacial frictional melting. In addition, SnowModel is a one-way model forced by meteorological data at each time step without considering exchanges from the surface back to the atmosphere (Mernild et al., 2014, 2017). Within this model, the input liquid equivalent precipitation is considered to be snow if the wet-bulb temperature is less than 1°C. (4) The SnowTran-3D sub-model, which is deactivated in this study, includes state variables related to snow redistribution by wind (Liston and Sturm, 1998, 2002; Liston et al., 2007). SnowModel uses a constant albedo of 0.8 for non-melting snow or user defined values for melting snow. 25 2.3.2. Alpine3D The seasonal and distributed physically-based snow cover model, Alpine3D (Lehning et al., 2006), also was used to evolve winter mass balance for the four glaciers of this study. Alpine3D simulates surface processes and spatial and temporal variability of snow cover in alpine terrain in a multi-layered snowpack. The number of layers is arbitrary, and the model can simulate very thin layers like ice crust and hoar. In the present study, we used the snowpack evolution module (SNOWPACK), and a 3-D radiation energy balance model, which calculates the effects of inclination and terrain exposure on the radiation balance. Alpine3D also manipulates input data using statistical methods through the MeteoIO module using transparent caching, filtering, resampling, and spatial interpolation (Bavay and Egger, 2014). An altitudedependent lapse rate and inverse distance weighting (IDW-LAPSE) was chosen to distribute meteorological data to the model domains. MeteoIO distributes the precipitation and temperature in complex terrain based on a linear regression of these variables with elevation using a constant lapse rate. The albedo is statically calculated by the models based on the snowpack parameters during the simulation. In Alpine3D, the snow cover characteristics and duration are predicted by the finite-element based SNOWPACK model (Bartelt and Lehning, 2002; Lehning et al., 2002a, b). This model solves the equations for non-stationary heat transfer and settlement in a phase-changing snow cover. The heat transfer, vapour diffusion, water transport, and phase changing processes is formulated. In addition, surface hoar and snow metamorphism are described by grain form and size, dendricity and sphericity. In SNOWPACK, the measured snow surface temperature is considered as an upper boundary condition as long as it is below -1.3°C (Lehning et al., 2002b). We did not use Alpine3D’s snow transport and hydrology routines. In Alpine3D, in addition to snow metamorphism, the 26 key difference with SnowModel is the development of multilayer snowpacks based on meteorological conditions (not user-defined layers) and considering the exchange of heat fluxes between layers (Table S1). 2.3.3. Model simulations Topographic data for the four sites were obtained from a digital elevation model (DEM). We used the Shuttle Radar Topography Mission (SRTM; Farr et al., 2007) dataset in 30-m grid cell increments that was merged with high-resolution LiDAR data to construct an updated terrain map over glaciers that accounts for elevation change between acquisition of the SRTM data (February, 2000) and the present day (September, 2015). The vegetation cover was derived from Global Land Cover database (GLC) (Chen et al., 2015) reclassified based on the vegetation cover classification in SnowModel and Alpine3D (Liston and Sturm, 2002; Mernild et al., 2006b; Lehning et al., 2006). The nearest grid points from HRM datasets were used as input for both snow models. We performed the simulations with a 30-m grid increment and one-hourly time step for NLDAS2, and 6-hourly time step for simulations using WRF. We found no substantial changes in snow parameters using interpolated hourly time series derived from WRF data instead of the 6-hourly timestep (less than 2% difference; Fig. S1.4). The simulation spanned two winter seasons, with winter defined as the period including 1 October - 30 April. To optimize the models’ performance and agreement between observed and simulated SWE in winter, we applied corrections to HRM temperature and precipitation to reduce elevation bias between HRM and station data before they were inserted into the snow models. We recognize that the optimized elevation correction function may vary in different locations (Liston and Elder, 2006b). We used monthly climate normals for 1981-2010 from the PRISM dataset (Daly et al., 1994) and applied a linear correction function for temperature 27 (Liston and Elder, 2006b) and precipitation (Pelt et al., 2016) of HRM data to adjust forcing data for local topography. Model performance was evaluated and optimized for one year (2013-2014) before doing the main simulation for winter 2015-2016. We tested different parameters such as lapse rate (Table 2.3), and snow albedo for both models, but found reasonable agreement with observations using default values. Curvature length (curvature length scale is approximately one-half the wavelength of the topographic features within the domain) in MicroMet, and different interpolation methods within MeteoIO, were also examined. The calibration of snow models was done based on minimizing bias, mean absolute error (MAE), and root mean squared error (RMSE) of SWE for the 2013-2014 year (Fig. S1.3). Our analysis shows a warmer and drier winter for 2015-2016 than in 2013-2014 (calibration period). Observations indicate a below normal snowpack in 2015-2016, while the 2013-2014 winter experienced near normal snowpacks in the study area. It is noting that although most days are below freezing in the simulation period, melt can still occur during the autumn and early spring. Table 2.3. The monthly temperature lapse rate and precipitation adjustment factors derived from PRISM, WRF and SnowModel default values (Kunkel, 1989). Temperature lapse rate (˚C km-1) Precipitation adjustment factor (m km -1) Oct Nov Dec Jan Feb Mar Apr May Oct Nov Dec Jan Feb Mar Apr May Default -6.8 -5.5 -4.7 -4.4 -5.9 -7.1 -7.8 -8.1 0.25 0.30 0.35 0.35 0.35 0.30 0.25 0.20 PRISM -3.5 -3.8 -3.2 -2.4 -3.4 -4.3 -4.8 -4.6 0.32 0.40 0.43 0.45 0.35 0.29 0.20 0.19 WRF -6.2 -5.2 -4.8 -4.7 -5.6 -5.7 -6.8 -6.6 0.26 0.31 0.34 0.38 0.39 0.20 0.01 0.06 2.4. Results In the following, we evaluate temperature and precipitation derived from HRM datasets and interpolated by snow model preprocessors, MicroMet and MeteoIO, against observations. We then describe the results of the SnowModel and Alpine3D simulations for the accumulation 28 season (1 October to 30 April) during calibration (2013-2014) and validation (2015-2016) periods. Model evaluation for the main simulations is done at the point and spatial scales using observed data from four ASWSs, stake measurements, and LiDAR data. 2.4.1. Evaluation of temperature and precipitation We compared daily average temperature derived from WRF and NLDAS2 datasets before and after downscaling with MicroMet and MeteoIO with observations for winter 2015-2016 at four ASWS (Fig. 2.3a). We used elevation-corrected data for the comparisons. The relationship between observed and modelled temperature for both WRF and NLDAS2 datasets is strong (R2 > 0.9) with a similar daily variability for both datasets (Fig S1.1). The measured daily mean temperature at all locations that range between -2.4 to -3.3˚C during winter 2015-2016 with biases range between -0.4 and -0.27˚C for Azure-River, +1.5 and +0.6˚C for Glacier-NMF, 0.6 and +0.5˚C for East-Creek, and -1.2 and -0.5˚C for Redfish-Creek with respect to ASWSs for WRF and NLDAS2, respectively. A strong relationship exists between observations and both MicroMet and MeteoIO downscaled temperature for all stations with coefficient of determination higher than 0.81 (p < 0.01), except for slightly lower R2 values at Azure-River. Despite the warm bias (1.5˚C) of WRF at Glacier-NMF, downscaling temperature with MicroMet and MeteoIO reproduce a colder monthly temperature at this station. MicroMet and MeteoIO show a mean cold bias of -1.3 and -1.9˚C, respectively over all stations. In general, both preprocessors, especially MicroMet, produce data with a cold bias for monthly temperature at all stations, with the highest bias at Glacier-NMF. We analyzed the precipitation derived from different datasets against observations to determine the influence of the spatialization with MicroMet and MeteoIO on precipitation distribution in the snow models. Because the daily simulated precipitation does not necessarily 29 match the observations according to the location and time of each individual event, adjusted monthly total precipitation were evaluated (Fig. 2.3b). While monthly trends are reproduced in most months for all datasets, WRF underestimates precipitation at Redfish-Creek (Fig. 2.3b) by an average of 34%. Interpolated precipitation by MicroMet and MeteoIO reduces this bias by 27% and 12%, which shows the difference between using the nearest WRF grid cell and interpolation of WRF data. For the same winter, MicroMet-based precipitation is overestimated in East-creek (23 mm) and Glacier-NMF (84 mm). The difference values vary for MicroMet and MeteoIO due to using different parameterizations and methods of interpolation. The values prepared by the model preprocessors are the final meteorological values to force the snow evolution modules within the models. Figure 2.3. a) Mean monthly temperature (˚C) and b) Mean monthly total precipitation (mm) time series for the accumulation period 2015-2016 at the station locations. Temperature values are derived from stations, WRF outputs, NLDAS2 dataset, and also from pre-processed modules (MeteoIo and MicroMet) of two snow models forced by WRF. Statistical values (R 2, RMSE, and Bias) of WRF and NLDAS monthly mean temperature and monthly mean precipitation against station (refers to Table 2.1) data are shown on the graphs. 30 2.4.2. Evaluation of SWE and Snow depth a) Point scale A pointwise comparison was done to evaluate the simulated daily SWE using ASWS observations (Table 2.4). For the average of the four ASWS, both models appear to perform well with Bias and RMSE scores less than 0.12 and 0.16 m w.e., especially using the WRF dataset, while the largest errors exist for SnowModel and Alpine3D using NLDAS2 data (Table 2.4). Table 2.4. Statistical analysis of daily time series of SWE simulated by SnowModel and Alpine3D forced with WRF fields and NLDAS dataset for the accumulation period 2013-2014 (1 October- 30 April). The correlation coefficients are calculated based on the p < 0.01 Station ID Function SM_WRF SM_NLDAS ALP_WRF ALP_NLDAS 2 0.94 0.88 0.90 0.67 MAE 0.06 0.11 0.06 0.06 Bias 0.06 -0.11 0.06 -0.03 RMSE 0.08 0.13 0.08 0.10 R2 0.98 0.96 0.98 0.96 MAE 0.10 0.12 0.04 0.06 Bias -0.10 -0.12 0.03 -0.02 RMSE R Azure-river East-creek 0.13 0.16 0.05 0.07 2 1.00 1.00 1.00 0.98 MAE 0.11 0.05 0.09 0.09 Bias 0.11 -0.01 -0.07 -0.07 RMSE R Glacier-NMF 0.12 0.06 0.11 0.12 2 1.00 1.00 1.00 0.96 MAE 0.03 0.18 0.12 0.34 Bias 0.02 -0.1 -0.12 -0.34 RMSE 0.04 0.21 0.14 0.40 R2 0.98 0.96 0.98 0.88 MAE 0.08 0.12 0.08 0.14 Bias +0.02 -0.09 -0.03 -0.12 RMSE 0.04 0.06 0.04 0.16 R Redfish-creek Average over all sites 31 A strong agreement between simulated and observed SWE is indicated. Coefficient of determination for all stations are higher than 0.70 (p < 0.01) using WRF-forcing, and are higher than 0.60 (p < 0.01) for NLDAS2 forcing in both calibration and evaluation years. Alpine3D underestimated NLDAS2-forced SWE at all sites with an average of 0.18 m w.e. (25%), while the bias for simulations forced by WRF shows a smaller underestimation of 0.01-0.05 m w.e (~3-9%). The largest percent bias occurred for SnowModel and Alpine3D using NLDAS2 data at Redfish-creek with -0.28 m w.e (~29.7%). The overall RMSE is less than 0.04 m w.e. for all simulations, except NLDAS2-forced SWE by Alpine3D with 0.11 m w.e. Snow accumulation is underestimated at Redfish-creek in March and April (Fig. 2.4). NLDAS2based SWE is well captured by SnowModel, while Alpine3D underestimates SWE at all locations. The largest absolute percent bias in WRF-forced SWE occurs at the East-creek (17%) for SnowModel and Glacier-NMF (14%) for Alpine3D. In the evaluation period (winter 2015-2016), the percent bias in SWE are reduced to 11% at East-creek, but increased to 17% at Redfish-creek. Boxplots of simulated SWE by both snow models with different driving data sets (Fig. 2.5) compared with stake measurements show lower median SWE values for most sites. Conrad Glacier is an exception, which had greater WRF-simulated SWE relative to observations at all elevations. A higher modelled median SWE also is shown by SnowModel forced by WRF against observations at the elevations of 2600-2800 m a.s.l. at Nordic Glacier. Both models estimate lower-than-observed SWE for Kokanee Glacier at all elevations. There is a strong relationship (0.81 < r < 0.99; p < 0.01) between simulated snow depth distributions and terrain elevation at all glaciers for both snow models and forcing data due to distribution of precipitation by elevation using adjustment factors and neglecting redistribution by wind. Both 32 models simulate a higher amount of SWE at higher elevations at the end of winter accumulation. Figure 2.4. Simulated daily SWE time series by SnowModel and Alpine3D forced by NLDAS and WRF data compared with observed snow pillow data at station locations for the accumulation period 2015-2016. Boxplots reveal a higher spatial variability for WRF-forced simulated snow depth by SnowModel compared with Alpine3D and NLDAS2 forcing simulations (Fig. 2.5). However, the simulated spatial variability of snow depth by both snow models is much lower than observations. 33 Figure 2.5. Boxplots comparing snow water equivalent [m] as a function of elevation at the end of the 2016 accumulation season for four glaciers. Black boxes - measurements, Blue boxes - SnowModel forced by WRF, Green boxes - SnowModel forced by NLDAS, Olive boxes - Alpine3D forced by WRF, and Gray boxes Alpine3D forced by NLDAS. The boxes span the interquartile range and the median is marked by a horizontal line inside the boxes. The whiskers are the dashed lines that extend to the highest and lowest values of each dataset. For the number of stake observations at each elevation band, please refer to Fig. 2.2. b) Spatial snow depth and SWE evaluation Overall, both snow models reproduce observed average snow depth from airborne LiDAR campaigns to within 20% for three of four glaciers (Table 2.5). 34 Table 2.5. Comparison of integrated snow depth (m) simulated by models over the whole domain for four glaciers with lidar measurement footprints using WRF and NLDAS forcing data for winter 2016. WRF_Forcing Conrad NLDAS_Forcing Nordic Kokanee Zillmer Average Conrad Nordic Kokanee Zillmer Average Lidar 2.9 2.4 3.02 3.16 2.87 2.9 2.4 3.02 3.16 2.87 SnowModel 3.44 3.19 2.9 3.23 3.19 2.23 2.28 2.99 2.26 2.44 Alpine3D 2.97 2.98 2.79 3.33 3.02 2.52 2.52 2.59 2.77 2.60 %Bias_SM 18.6 32.9 -3.9 2.2 11.15 -23.1 -5.6 -0.88 -28.57 -15.0 %Bias_ALP 2.4 24.1 -7.49 5.3 5.23 -13.1 4.56 -14 -12 -9.41 In terms of forcing data, while distributed snow depth using NLDAS2 data was underestimated for all glaciers, WRF overestimated average snow depth for Nordic and underestimated snow depth in Kokanee when used to force both snow models. WRF produces more variability than NLDAS2 (Fig. 2.6) and appears more similar to the spatial variability in the LiDAR data, but future work is needed to determine whether it represents the spatial complexity more accurately than NLDAS2. But in terms of the model differences, the comparisons suggest that SnowModel captures the snow pattern better than Alpine3D. However, the differences between the models are small compared to the uncertainties. The latter model produces a more uniform snow pattern than SnowModel (Fig. 2.6, Table 2.6). Generally, SnowModel produces deeper snow at higher elevations and shallower snow at low elevations for all glaciers, while there is little difference in snow depth at low and high altitudes in Alpine3D, which is in agreement with boxplots of measured snow depth over terrain at point scales (Fig. 2.5). Mean snow depth of LiDAR measurements and averaged simulated snow depth for all glaciers are shown for both snow models forced with HRM data (Table 2.5). The overall calculated percent differences revealed a good approximation of mean snow depth over domains for both snow models using both driving datasets. 35 Figure 2.6. First column illustrates the LiDAR elevation-change (ΔH) over all glaciers. ΔH represents the difference of two elevation measurements at the beginning and end of accumulation season over each domain. The other columns show simulated snow depth distribution at the end of the accumulation period (end of April 2016) over four glaciers using SnowModel and Alpine3D forced with WRF and NLDAS2 datasets. The black outlines represent the main glacier borders where probing of snow depth was completed. Note: LiDAR measurements are elevation-change measurements including ice dynamics so do not represent snow depth changes at specific points. The wind effects on snow pattern are neglected for the simulations. Although WRF-forced results demonstrate an improvement of snow depth over all domains, an overestimation is clear with the highest percent difference for Nordic Glacier. The overall average simulated snow depth over LiDAR footprints by WRF/NLDAS2 is overestimated/underestimated, respectively by +11% and -15% compared with LiDAR measurements. Alpine3D approximates the average snow depth for both forcing datasets better than SnowModel for all glaciers over the LiDAR footprints, while the results are different on 36 glaciers and in glacier-free areas. To examine model performance over both glacierized and non-glacierized areas separately (due to availability of LiDAR footprints larger than glaciers, but having only stake measurements on the main glaciers), we clipped both WRF and NLDAS2-forced modelled snow depth over each domain using the main glaciers shape files to calculate the averaged snow depth on each glacier and in glacier-free areas. Table 2.6. Comparison of integrated WRF-forced and NLDAS-forced snow depth (m) with lidar measurements on the glacierized and in free glacier areas (same as Table 5 but separated for glacierized and glacier-free areas). WRF_forced_Integrated snow depth (m) Off_glacier NLDAS_forced_Integrated snow depth (m) Off_glacier Conrad Nordic Kokanee Zillmer Average Conrad Nordic Kokanee Zillmer Average Lidar 2.86 2.15 2.89 3.10 2.93 2.86 2.15 2.89 3.10 2.93 Lidar_Sd 1.53 2.03 1.37 1.62 1.54 1.53 2.03 1.37 1.62 1.54 SnowModel 3.42 3.01 2.87 3.21 3.11 3.42 2.25 2.95 2.26 2.42 SM_Sd† 0.53 0.57 0.34 0.65 0.52 0.53 0.40 0.43 0.39 0.39 Alpine3D 2.99 2.97 3.58 3.33 3.21 2.99 2.52 2.60 2.36 2.49 ALP_Sd† 0.29 0.32 0.42 0.45 0.36 0.29 0.24 0.19 0.20 0.19 %Bias_SM 19.60 40.00 -0.70 3.50 13.00 19.50 4.65 2.08 -27.00 -7.20 %Bias_Alp 4.60 38.00 23.00 7.40 16.70 4.55 17.21 -10.03 -23.90 -10.60 WRF_forced_Integrated snow depth (m) ON_glacier NLDAS_Forced_Integrated snow depth (m) ON_glacier Lidar 3.08 3.77 3.55 4.08 3.62 3.08 3.77 3.55 4.08 3.62 SnowModel 3.50 4.13 3.28 3.37 3.57 2.25 2.39 3.49 2.31 2.61 Alpine3D 2.92 3.00 2.94 3.32 3.05 2.46 2.52 2.70 2.34 2.50 %Bias_SM 13.60 9.50 -19.60 -5.07 -1.38 -27.00 -36.60 -1.70 -43.40 -27.90 %Bias_Alp -5.20 -20.40 -27.90 -6.40 -15.70 -20.00 -33.20 -24.00 -42.60 -30.80 †Sd: Standard deviation Averaged over all locations, both models overestimate the snow depth over glacier-free areas, while they underestimate snow depth on the main glaciers. There is a considerable difference between SnowModel and Alpine3D at different locations on the main glaciers. WRF-forced SnowModel underestimated the snow depth overall less than 2%, while Alpine3D 37 underestimated snow depth by about 15% (Table 2.6). The average underestimated snow depth reaches to more than 30% when NLDAS2-forced data is used. Simulated spatial winter mass balance is also evaluated using calculated mass balance from stake and snow pit measurements. Mass balance calculation from LiDAR data was explained previously in section 3.2.d. We compared the simulated, geodetic and glaciological mass balance across all glaciers for winter 2016 (Fig. 2.7). Figure 2.7. Comparison of simulated winter mass balance [m w.e.] by SnowModel and Alpine3D with stake and LiDAR measurements on the main glaciers for the accumulation period of 2015-2016. The notation ‘_B’ denotes the bias (absolute magnitude) of simulated winter mass balance with respect to the observations. Comparisons at individual glaciers revealed that the least and highest percent difference in WRF-forced simulated winter mass balance occurs over Nordic and Kokanee, respectively. Overall winter mass balance over glaciers was calculated to be 1.93 m w.e., while WRF-forced SnowModel and Alpine3D produced 1.74 and 1.51 m w.e, respectively. NLDAS2-forced data reproduced a lower amount of winter mass balance by 0.30 m w.e. averaged over all glaciers. Both snow models underestimate the winter mass balance over glaciers. SnowModel underestimated winter mass balance against stake measurements by about 10% and less than 38 2% compared with LiDAR data, while the difference between observations and Alpine3D was respectively 22% and 13%. Our analysis suggests a better performance for SnowModel by 12% in both comparisons. It is worth noting that simulations using WRF forcing for the accumulation period 2015-2016 decreased the bias by 4% in both models compared to NLDAS2 data. 2.5. Sensitivity analysis and potential uncertainties In this section we investigate the sensitivity of snow amount to several parameters including monthly temperature lapse rates and precipitation adjustment factors, rain/snow threshold for Alpine3D, and multilayer modelling for SnowModel. We tested sensitivity of simulated SWE by Alpine3D using different constant temperature and precipitation adjustment factors (Schlӧgl et al., 2016). Monthly temperature lapse rate and precipitation adjustment factors derived from PRISM and WRF (Table 2.3) were used within the SnowModel and compared with default values to capture the model’s sensitivity to these parameters. Percent difference of simulated daily SWE associated with alternate monthly temperature lapse rates and precipitation adjustment factors derived from WRF and PRISM averaged over all locations is fairly small (< 3%), because forcing the elevation bias corrected data into the snow models, effectively over-riding much of the elevation adjustments made within the snow models. MicroMet specifies precipitation as snow when wet-bulb temperature is less than 1˚C, while snowfall is based on a threshold of 1.2˚C air temperature in Alpine3D (Dai, 2008). The threshold value might affect the snow amount and its frequency. To investigate the impact of different rain/snow thresholds on the simulated SWE by Alpine3D, we ran the model at two locations (Kokanee Glacier and Azure-River station) using a threshold of 2˚C air temperature 39 (instead of the default value of 1.2˚C used in the other simulations) during accumulation season. We found no significant differences (about 0.01 and 0.001 m w.e.) averaged over the glacier area and at station locations, respectively (Fig. S1.5 and Fig. S1.6). We suspect the lack of sensitivity to the rain-snow threshold is because air temperatures are below freezing on the majority of days during the accumulation season for these two locations. Furthermore, we tested SnowModel with a multilayer snowpack simulation (instead of the single layer main simulation) to investigate the effect of increased vertical resolution on simulated snow depth and SWE. We employed six layers for our simulations at Kokanee Glacier (2800 m a.s.l) and Redfish creek station (2000 m a.s.l) locations. We found a mean daily density difference of less than 12 kg.m-3 (~ 2.5%) at both station location and averaged over the glacier which occurs in the late accumulation season, with less difference on the glaciers (Fig. S1.7 and Table S1.5). The difference in SWE was also found to be less than 0.04 m w.e. at station locations with almost no difference (~ 2× 10 -5 m w.e.) on the glacier itself (Fig. S1.8). Uncertainties in forcing data, such as uncorrected bias, and uncertainties in parameterizations in WRF and simplified MicroMet downscaling methods, contribute to bias in simulated snow patterns. Glacier dynamics transport ice down-glacier resulting in higher surface altitudes at the lower elevations and lower surface altitudes at the higher elevations that are not considered within this study. Furthermore, the wind influence on redistributed snow and preferential deposition of snowfall, and avalanche on winter mass balance are not considered in this study and might affect the results, especially on small glaciers. Uncertainties in the measurements also need to be considered. One systematic bias in probing data is the lack of measurements from inaccessible parts of the glaciers. According to some studies, typical 40 errors reported for individual measurements are 0.1-0.3 m w.e. a-1 for snow accumulation (e.g. Gerbaux et al., 2005; Thibert et al., 2008; Huss et al., 2009). Other studies also detected a range of 13-30% error in snow accumulation measurements in different study areas (e.g. Jonas et al., 2009; Sturm and Liston, 2003; Kershaw and McCulloch, 2007). Similarly, errors in mass balance calculated from high resolution LiDAR data with a fixed density is reported in other studies to be 0.15-0.30 m w.e. a-1 (Jóhannesson et al., 2013; Berthier et al., 2014). We found strong agreement between glaciological and geodetic methods for our study glaciers, with glaciological winter balance averaging 1.95±0.08 m w.e., about 4% greater than the geodetic estimate. This agreement suggests that overall, the field sampling scheme reliably represents the glacier-averaged SWE (Pelto et al., 2019). Our findings also revealed a significant agreement between simulated mass balance by SnowModel (about -2%) and Alpine3D (-13%) with geodetic estimates. 2.6. Discussion a) Meteorological data The derived temperature from HRM data accord with observations but show a cold bias for most stations when using WRF data except at one station. We used an identical lapse rate for elevation bias correction at all locations, however the lapse rate can vary between different sites due to topographic complexity in the area (Gao et al., 2012, 2014b; Liston and Elder, 2006). Our results are consistent with other studies, which also found a prominent winter cold bias for both datasets (e.g. Pan et al., 2003; Erler et al., 2015; Lewis and Allen, 2016; Meng et al., 2018), however the NLDAS2 cold bias is lower than the WRF cold bias (refer to Fig. 2.3a). MicroMet and MeteoIO (meteorological preprocessors in SnowModel and Alpine3D) also do 41 not simulate local temperature inversions. MicroMet calculates cloud-cover fraction by converting near surface air temperature and relative humidity to a 700 hPa equivalent elevation, which is not suitable for all conditions, for instance under clear skies at night (Liston and Elder, 2006). The difference in cold bias between the preprocessors is likely due to different lapse rate values and methods utilized for interpolation that cause slightly different meteorological outputs. The winter cold bias might affect the winter mass balance simulation due to the temperature threshold used to divide precipitation into rainfall and snowfall, especially in the beginning and late accumulation season. While this factor could affect the results at station locations at the end of the accumulation season, the glaciers of this study are, on average, 500 m higher than the meteorological stations used in this study and so temperatures there likely remain close to or below freezing. The simulated monthly total precipitation for all stations is relatively consistent with observations. A lower amount of precipitation extracted from NLDAS2 is likely based on the small number of weather stations in these areas of complex terrain. These stations are mostly located at lower elevations and do not capture the full meteorological variability in mountainous areas (Schultz et al., 2002; Pan et al., 2003; Mitchell et al., 2004). An evaluation of WRF performance with the same configuration and parameterization, on climatological monthly mean simulations for Western Canada, can be found in Erler et al. (2015). The downscaled climate simulations were evaluated against several gridded observational datasets: the University of East Anglia Climatic Research Unit (CRU) 0.5° multivariate monthly time series (Harris et al., 2014) and the Global Precipitation Climatology Centre (GPCC) 0.5° monthly precipitation time series (1951-2001), and a 0.25° climatology for 1901 to 2012 (Schneider et al., 2014). The cornerstone of their evaluation used the high-resolution (1/48° ≈ 3.3 km) Parameter-Elevation Regressions on Independent Slopes 42 Model (PRISM) temperature and precipitation climatology, based on the period from 1979 to 2009. Determining precipitation to estimate winter mass balance is challenging. Vionnet et al. (2019) and Quéno et al. (2018) applied a numerical weather prediction system, AROME, to generate high resolution forcing data to force a snow model (e.g. Crocus). Their study demonstrates the need to improve the ability of high-resolution numerical weather systems to produce the precipitation field with higher accuracy for snow and glacier modelling in alpine terrain. Precipitation is affected by local topographic features and wind speed and direction as well as the type of weather systems. Precipitation does not vary as linearly as temperature with elevation (Liston and Elder, 2005; Jarosch et al., 2010), making its extrapolation a more challenging problem in atmospheric and snow simulations. We believe that the accuracy of input data (coarse resolution data, such as NLDAS2) to snow models can be improved by using the power spectral method described by Jarosch et al. (2010) rather than using simple bilinear interpolation, which does not add spatial information, or a fixed lapse rate adjustment, especially for downscaling precipitation data in mountainous area. They downscaled precipitation based on a linear theory for orographic precipitation (Smith and Barstad, 2004) including airflow dynamics, condensed water advection, and downslope evaporation, which was modified for large study regions by including moist air tracking. In addition, the capability and efficiency of snow models can be enhanced by improving solid precipitation prediction at high elevations using data assimilation methods (e.g. Winstral et al., 2019). b) Snow depth and SWE The simulated SWE at the point scale accords with observed data at station locations. In terms of point scale comparison, the results emphasize the applicability of SnowModel and 43 Alpine3D using both datasets and demonstrate improved performance with WRF driving data with respect to SWE compared to observations. Reducing meteorological forcing bias and calibration of models for each station separately could increase the accuracy of the models. Underestimation of SWE by both models in some locations is likely due to underestimated precipitation, while neglect of snow transport by wind also might affect the results. The integrated snow depths derived from the WRF forcing data more closely simulates average snow depth than NLDAS2-forced snow depth over the glaciers of this study. This might be due to underestimated precipitation in NLDAS2, especially in high elevation areas with scarce observations. In addition, NLDAS2 does not correct for gauge undercatch, and therefore probably underestimates actual precipitation (Wright et al., 2017). Gauge undercatch errors in gridded precipitation data sets can be significant, especially for snowfall and extreme events (Adam and Lettenmaier, 2003). The distributed snow depth over each glacier and percent bias show that results are consistent with meteorological forcing data (sections 3.1.1 and 3.1.2). Overestimated snow depth due to overestimation of simulated precipitation by regional climate models is also indicated in other studies (e.g. Bellaire and Jamieson, 2013; Parajca et al., 2010). There is a smaller range of snow depth variation in Alpine3D than in SnowModel, likely due to different meteorological data interpolation methods in MeteoIO compared with MicroMet. Similar to our results, Vogali et al. (2016) reported almost a uniform snow depth (higher snow depth in lower elevations and a lower snow depth at the higher elevations compared with observations) in simulations by Alpine3D. In MeteoIO, the IDW interpolation method used does not include the effects of topography on the meteorological data (Schlögl et al, 2016; Sturm and Wagner, 2010), but vertical changes are considered using a constant lapse rate. Such standard 44 precipitation interpolation methods may not be appropriate to distribute precipitation in complex terrain, and usually result in a relatively homogeneous snow cover (Trujillo et al., 2019). Topography is an important factor in snow accumulation since it regulates the spatial distribution of precipitation both vertically and horizontally (Machguth et al., 2006b; Lehning et al., 2008; Dadic et al., 2010). Overall, winter mass balance is underestimated by both models, compared with the calculated mass balance from stake and snow pit measurements over the glaciers. Underestimation of annual mass balance and SWE by physically-based snow models forced by weather station data has been reported in other studies (e.g. Bavera et al., 2014; Paul and Kotlarski, 2010; Schmucki et al., 2014). In all cases a precipitation correction factor was applied to reduce wind induced undercatch leading to reduction of the mean absolute percent error of 14% and 19% for snow depth and SWE, respectively. They concluded that the most critical part of the modelling process is availability of accurate, gridded precipitation data over Alpine terrain. Averaged WRF-forced simulated snow depth is overestimated in glacier-free areas, while the bias is positive or negative over each individual glacier (Table 2.6). SnowModel neglects the geothermal ground heat flux and in Alpine3D, we used a global average value of 0.6 W m-2. Neglecting heat flux could explain discrepancies over ice-free terrain, however Watson et al. (2006) revealed little sensitivity of Alpine3D to a selected mid-range value (1 W m-2) of geothermal heat. The model error also may be induced by inconsistencies resulting from the model structure and parameterization, especially the preprocessor methods. Along with forcing data, neglecting processes that increase the winter snow mass balance – such as avalanches, and wind-driven processes such as preferential deposition of snow and drift, might be responsible for underestimation of simulated snow depth on the glaciers. Neglecting wind effects in snow simulations also may cause snow density to be underestimated with respect to 45 observations (e.g. Liston and Hiemstra, 2011; Sturm et al., 2010). We found that SnowModel underestimates snow density, by less than 1% and 3% against snow pit measurements and calculated density from ASWS data, respectively. Alpine3D underestimates density by 23% and 24% when compared with observed data from snow pits and ASWS data, respectively (Table 2.7). Table 2.7. Snow density (kg m -3) comparison between simulated density by SnowModel and Alpine3D and observed data from snow pit measurements and density calculated from snow pillow data at the end of accumulation season (April 2016). Density (kg m-3): Snow pillow† data Density (kg m-3): Snow pit measurements RF_C A_R Average over all stations Average over Conrad Nordic Kokanee Zillmer E_C G_NMF all glaciers OBS 480±5 565±10 464±20 413±10 480 500±10 535±5 532±5 549±2 529 SnowModel Alpine3D %B‡-SM %B‡_ALP 543±5 405±10 8 -26 494±20 358±10 -12 -28 446±5 345±2 -4 -23 434±5 352±2 5 -19 479 365 -0.2 -23 509±5 433±8 2 -13 545±5 384±5 2 -28 528±12 356±4 -1 -33 478±5 382 -13 -30 515 389 -2.6 -24 ‡%Bias; †E_C: East_Creek; G_MNF: Glacier-NMF; RF_C: Redfish_Creek; A_R: Azure-River The results show that both snow models lead to similar relative difference values of snow depth when compared to the observations at all locations, with a low variability between the models (~ 4% averaged over all glaciers). However, for SWE and density, the snow models are somewhat less consistent. In particular, the difference in the simulated snow density by the models is considerable (~ 25%), leading to differences in estimating glacier winter mass balance. Underestimation of snow density by Alpine3D, could be due to two reasons: underestimation of new snow density, or slow compaction rate, or both (Saito et al., 2012). Multilayer snowpack in SnowModel does not make a considerable change in the results, given the fact that SnowModel considers only thermal conduction between the layers. Conduction has less influence on ice-covered areas, certainly not as much as forested floors or 46 even bare rock. Our results and use of a 1-layer snowpack could be transferrable to other glaciers but probably not for alpine/forested areas in general. Based on the results obtained, this study suggests that a model with higher order complexity may not be required to simulate winter mass balance over large alpine glaciers. We employed the snow evolution modules of both snow models (e.g. snow metamorphism, layering, and energy exchange between the layers, ground and air), but did not consider effects such as snow avalanches or wind and their effects on snow redistribution. Our analysis reveals a degraded (< 15%) performance for WRF-forced distributed snow depth, possibly because of factors that involve redistribution of snow by gravity and wind, which are probably much more important for small glaciers such as Kokanee and Nordic glaciers. In the absence of wind transport, WRFforced simulations are better than NLDAS2-forced simulations for both snow models at the four glaciers. However, NLDAS2 data still is a valuable and easily accessible dataset for estimation of winter accumulation in areas with scarce observations, recognizing there is uncertainty due to underestimated precipitation. NLDAS2 might be useful if it could be bias corrected, but that would require a much larger domain over which to complete that work. In addition, one needs to ensure that meteorological records used for bias correction are not directly assimilated by NLDAS2. Our future work will assess the effects of wind on winter accumulation in this area. Overall, comparison of SnowModel and Alpine3D revealed that the former model better captured the variability of snow depth against elevation even though gravitation effects, avalanche, and the interaction of wind with terrain were not considered in this study. Different methods of data distribution and interpolation within MicroMet and MeteIO also contribute to these results. 47 2.7. Conclusion The main objective of the study was to evaluate different forcing data and the performance of two forcing spatialization and snowpack evolution models with different complexity in simulating winter mass balance for alpine glaciers in western Canada and also in free glacier areas. Our analysis showed that SnowModel and Alpine3D successfully reproduced short-term variations of the observed snow depth (within ± 15%) and SWE (< 10% for SnowModel and < 20% for Alpine3D) for four glaciers in BC’s Columbia Mountains over an accumulation season with downscaled meteorological forcing data from WRF and NLDAS2. Our results imply that with input data consisting of air temperature, relative humidity, wind speed, incoming shortwave and longwave radiation and precipitation intensity from HRM data, winter mass balance can be estimated with less than 20% bias for glaciers located in remote areas where measured meteorological data are not available. The greatest source of uncertainty relates to precipitation forcing and its spatial variability as a function of altitude and cloudiness (Schmucki et al., 2013). In this study we focused on the accumulation season with limited ablation processes, so error related to uncertainty in radiation fluxes and temperature should be minor. Our results reveal that WRF-forced snow depth and SWE are closer to the observations than snow depth and SWE using NLDAS2 forcing data that has a slightly higher underestimation. SnowModel reproduces winter mass balance with less than 10% bias, compared with a 22% bias for Alpine3D, which is likely due to a greater underestimation of snow density by Alpine3D. The variability of snow depth also is better represented by SnowModel than by Alpine3D due to its different interpolation methods. In this study, snow redistribution by wind is neglected and influence of wind-terrain-interactions on snow pattern should be investigated. 48 Simulations suggest that Alpine3D with a more complex snow evolution module does not necessarily produce a better result for this application. Given the time and computational demands of Alpine3D versus SnowModel (13-fold increase in execution time) the value of a more complex model was not demonstrated for winter mass balance estimation. The accuracy of meteorological forcing data, in particular precipitation, is the most important factor in accurately simulating snow properties, hence the data downscaling and interpolation methods play a key role in simulations when the resolution of input data is not as high as the resolution of model outputs. WRF is a valuable tool for producing meteorological fields as input for snow modelling and can be a good alternative when weather station data are absent or lacking. In addition, the spatial snow pattern might be better represented by implementing snow drift modules using high resolution wind fields from WRF. This study suggests that a one-layer energy balance model is sufficient to estimate winter mass balance, but there may be limitations for simulating snowmelt patterns. Similar studies have demonstrated the SnowModel ability to estimate snow depth and SWE forced by WRF dataset in Spain and Chile (Alonso-González et al., 2018; Réveillet et al., 2020). Therefore, our results suggest that SnowModel can be a better choice for estimating winter mass balance in alpine terrain in terms of simplicity, time and computing demand and data distribution over the domain. 49 Chapter 3. Wind influenced snow redistribution and sublimation over alpine terrain using dynamically downscaled data Mortezapour, M.1, Menounos, B.1, Jackson, P. L.1 1. Natural Resources and Environmental Studies Program, Geography department, Prince George, British Columbia, Canada 2. Natural Resources and Environmental Studies Program, Atmospheric Science department, Prince George, British Columbia, Canada 50 Abstract We model winter mass balance at four glaciers in the Columbia Mountains of British Columbia, Canada using a distributed numerical snow model (SnowModel). The model is forced by dynamically downscaled meteorological data from the Weather Research and Forecasting Model (WRF) and North American Land Data Assimilation System (NLDAS) for winter 2015-2016. The modelling was performed to assess the effects of wind on snow redistribution and winter mass balance. The integrated snow depth and winter mass balance over the surface area of all glaciers were compared with high resolution light detection and ranging (LiDAR) measurements and manual snow probing. Including the snow transport module (SnowTran-3D) into the simulation forced by NLDAS has a negligible effect on the glacier-averaged mass balance. Wind improves the simulated WRF-forced winter mass balance by 2% when compared with glaciological observations. The small-scale snow pattern is not represented by NLDAS data resulting in a 6% greater bias compared with the WRFforced simulation, which is able to represent the snow pattern with realistic small-scale features. Spatial variability of snow drifting sublimation is large and reaches up to 540 mm w.e. summed over the winter accumulation period at the small scales. The strongest contribution of snow reduction by drifting sublimation is about 8% of the snow water equivalent (SWE) with an average of less than 17 mm or 2% over all the sites. 51 3.1. Introduction Snow transport by wind affects the temporal and spatial distribution of snow cover in the winter season. In complex terrain, wind and precipitation interact with topography and results in inhomogeneous snow distribution (e.g. Pomeroy and Gray, 1990; Doorschot et al., 2001; Mott et al., 2010; Mott et al., 2018). A non-uniform snow cover results in large uncertainty in the estimation of snow water storage and, consequently, hydrological resources of alpine catchments (Mernild et al., 2006; Lehning et al., 2008). Snow accumulation is an important part of glacier mass balance in alpine areas that has been less studied than processes affecting summer ablation. Snow distribution in mountainous areas is spatially variable due to spatial differences in factors such as topography, shading, and wind redistribution that cause snow accumulation disparities in alpine terrain resulting in mass balance variations (Winstral and Marks, 2002; Sold et al., 2013; Winstral et al., 2013). There are other processes that affect winter accumulation including avalanching, which has been shown to be highly important for smaller glaciers (e.g. Kuhn, 1995; Mott et al., 2019) or winddriven snow accumulation processes (Dadic et al., 2010; Vionnet et al., 2017; Mott et al., 2018). Snow distribution is also strongly dependent on both large-scale orographic precipitation and local scale processes (e.g. Houze, 2012; Mott et al., 2014). It is not only orographic precipitation that plays an important role but also the effect of the mean wind field on the local precipitation field, as discussed by various studies (i.e. Vionnet et al., 2017; Mott et al., 2018; Gerber et al., 2019). The main physical processes governing snow accumulation include: (a) the transport of deposited snow (snow drift), which includes suspension and saltation processes; and (b) the preferential deposition of precipitation due to topographicinduced wind field perturbations during a snow storm (Lehning et al., 2008; Dadic et al., 2010). 52 The importance of topography on the surface energy balance for a snowpack was investigated in several studies (e.g. Hock and Holmgren, 2005; Arnold et al., 2006). The enhanced snow accumulation on leeward slopes is due to preferential deposition, which is an important source of uncertainty in the surface mass balance of cold regions (e.g. Lehning et al., 2008; Mott and Lehning 2010; Comola et al., 2019) and at smaller scales is affected by local topographic features and wind characteristics (Mott et al., 2010). Snow accumulation is affected by wind redistribution, especially under high wind speed conditions. Such conditions also lead to a large amount of drifting snow sublimation (Pomeroy and Gray, 1990; Liston and Sturm, 1998; Doorschot et al., 2001; Lehning et al., 2006; Pomeroy et al., 2006; Bernhardt et al., 2009; Dadic et al., 2010). Drifting snow also is an important factor to be considered for detailed glacier mass balance assessments (Mott et al., 2008). To accurately simulate snow accumulation in alpine terrain, topography and its interaction with wind need to be considered (Winstral and Marks, 2002; Sold et al., 2013, Mott et al., 2018, Vionnet et al., 2018). Various blowing snow models have been developed during the last three decades, for example, the Prairie Blowing Snow Model (PBSM) (e.g. Pomeroy et al., 1993; Pomeroy and Li, 1997b), PIEKTUK (e.g. Déry and Yau, 1999), SnowTran3D (e.g. Liston et al., 2007), Snow2blow (e.g. Sauter et al., 2013) and iSNOBAL (e.g. Winstral et al., 2013). These models cover a wide range of complexity, but generally consist of two components: 1) a snowpack module to estimate the threshold wind speed for snow transport and erodible snow mass; and 2) an atmospheric component to simulate the spatial and temporal evolution of the wind field and snow transport (Liston et al., 2007). Wind affects snow accumulation through the transport of snow to leeward sides and sheltered areas as well as having an influence on the snowpack due to sublimation and evaporation 53 processes. Some studies estimate a small amount of turbulent sublimation which varies based on the site properties such as topographic configuration and vegetation type (Gordon et al., 2006; Strasser et al., 2007). Sublimation due to snow transported by wind in some cases is not negligible, however (e.g. Strasser et al., 2007; MacDonald et al., 2010; Zwaaftink et al., 2011; Vionnet et al., 2014). These studies demonstrate that drifting snow sublimation has a large spatial variability in complex terrain, especially at crests during a winter snowfall. Therefore, an accurate representation of drifting snow is needed for a reliable estimate of blowing snow sublimation, which can reduce snow accumulation over alpine areas. The main motivation for the present study is to estimate the impact of wind on winter mass balance. The simulation of wind-induced snow transport in alpine terrain requires high-resolution wind fields in regions of complex topography. Consequently, this study applies a distributed snow evolution model: SnowModel (Liston and Sturm, 1998; Liston et al., 2007) forced by high resolution meteorological fields from the Weather Research and Forecasting model (WRF) to represent the effect of wind on the winter mass balance in alpine terrain. We examine snow transport effects and quantify the seasonal pattern of drifting sublimation losses at four individual glaciers over the Columbia Mountains. 3.2. Study sites Four individual glaciers have been chosen in the southeastern interior of British Columbia (Fig. 3.1). The mid elevation of these glaciers ranges are 2094, 2301, 2718 and 2600 m above sea level (m a.s.l.) for Zillmer (52°41' N, 119°35' W), Conrad (50°49' N, 116°55' W), Nordic (51° 25' 13'' N, 117° 42' 13'' W) and Kokanee (49° 45' N, 117° 8' 30'' W), with the areas of ice masses of 5.2, 16.2, 3.4, 1.8 km2, respectively. All glaciers flow east with the exception of the 54 north-facing Nordic Glacier (Fig.3.1). The complex features over the interior ranges of BC has a significant influence on its climate. The combination of changing topographical factors such as elevation, slope, orientation to the mean wind, and topographic exposure alter the local climate on small scales. There is a lower average temperature and increased precipitation at higher elevations, while in the valleys, inversions resulting from cold air drainage increases the occurrence of frost and fog. Moderate temperature, heavy precipitation and snowfall in wintertime are characteristics of the Columbia Mountains (Klock and Mullock, 2002). Figure 3.1. The four glacier locations (stars) over Columbia Mountains, British Columbia, Canada. Red triangles represent the station locations used for model evaluation. The topography of each glacier is also displayed using contours over each simulation domain. The yellow outline shows the LiDAR footprint and main glaciers are represented by red outlines. 3.3. Snow Evolution Model description In this study, SnowModel (Liston and Elder, 2006), a distributed energy balance snow model is used to simulate the time evolution of the winter snowpack. SnowModel consists of a micro-meteorological model (MicroMet) that interpolates and distributes forcing data over 55 the model domain, based on elevation and land cover (Liston and Elder, 2005). The interpolated meteorological fields of precipitation, humidity, pressure, and incoming irradiance are passed to EnBal (Liston, 1995; Liston et al., 1999) and SnowPack (Liston and Hall, 1995) to compute snow evolution. The vertical distribution of the atmospheric under-saturation of water vapour with respect to ice is provided by a modification of equations provided by Pomeroy et al. (1993), and Déry and Taylor (1996) to describe the vertical variation of relative humidity from observations at 2 m. The energy balance is used to calculate the energy available for ablation, while internal snowpack processes, forest-snow interactions, snow sublimation, and wind-snow interactions are computed by SnowPack and SnowTran-3D (Liston and Sturm, 1998, 2002; Liston et al., 2007). A full description of each sub-model, governing equations, approximations, and assumptions can be found in Liston (1995), Liston and Hall (1995), Liston and Sturm (1998), Liston and Elder (2006a, 2006b), Liston et al. (2007), and Liston and Hiemstra (2008). SnowModel is a one-way model forced by meteorological data at each time step (Mernild et al., 2014, 2017). The input water equivalent precipitation is considered to be snow if the wetbulb temperature is less than 1°C. Static surface sublimation, which is the turbulent flux of latent heat from the surface, is simulated by the EnBal module using an energy balance calculation. SnowTran3D (Liston and Sturm, 1998) is a three-dimensional model that calculates the saltation, suspension and sublimation of snow being transported by wind. Blowing snow transport includes three modes of movement: creep - the rolling movement of those particles, which are too heavy for wind to move; saltation - the movement of snow particles along the snow surface through jumping; and suspension - the motion of snow particles in suspended flow with a mean horizontal velocity close to the air speed (Pomeroy 56 and Brun, 2001). SnowModel solves the mass balance equation for the temporal changes in snow depth in the four vertical layers due to melting, precipitation, saltation and suspension (Liston et al., 2006). = [ρ P − + + + +Q ] (1) where ζ is the snow depth (m), ρs and ρw respectively, are snow and water density (kg m-3), P is the precipitation rate (m s-1), and (kg m-1 s-1) are the saltation rate, and amount of turbulent suspended snow, respectively. Qv represents the sublimation of transported snow particles in (kg m-1 s-1). The amount of transported snow at a given time step depends on the shear stress due to wind on the surface. The saltation transport process under equilibrium condition is described according to Pomeroy and Gray (1990). Q = . u∗ (u∗ −u∗ ) ∗ (2) where ρa (kg m -3) is the air density. In general, ∗ will vary in space and time, so Q s-max will also vary. A sufficient fetch length is required for snow transport to reach equilibrium. For the non-equilibrium (and equilibrium) conditions the saltation-transport rate, Qs (kg m-1s-1) is assumed to follow ( ∗) ∗ = ϝ Q_ Q (x ∗ ) = min (x ∗ ) − Q (x ∗ ) Q (x ∗ − Δx ∗ ) Q _ (x ∗ ) for for ∗ ∗ ≥ 0 (2a) ∗ ∗ <0 (2b) where x* (m) is the horizontal coordinate in a reference frame defined by the wind-flow direction, Δx* (m) is horizontal grid resolution; μ is a unitless constant, and f (m) is the equilibrium-fetch distance, which is assumed to be between 300 m to 500 m. However, the 57 actual fetch is variable through time and space, and depends on wind speed, surface roughness, and snow supply. A threshold friction velocity ( ∗ ) must be met to lift the snow for transport. The key determining factor for transport is the shear stress on the surface, produced by the wind which is defined by ∗ (friction velocity) u∗ = u (3) where ur (m s-1) is the wind speed at the reference height zr (m), the surface roughness length is represented by z0 (m), and k is the von Kármán’s constant. To initiate transport, the friction velocity of a given wind speed must exceed the threshold friction velocity value ( threshold friction velocity is lower for fresh and dry snow during snowfall ( s-1) and higher for older, wet and dense snow ( ∗ = 0.25-1.0 m s -1 ∗ ). The ∗ = 0.07-0.25 m ) (Kind, 1981). The threshold velocity can be modified by the roughness length of the surface, which depends on the snow depth fraction and the vegetation type. In this study we used the pre-defined vegetation types with corresponding snow-holding depths. The maximum saltation transport is achieved only if there is sufficient distance to transport the full amount of snow that the wind speed is capable to carry (Liston and Sturm, 1998). The snow concentration at the top of the saltation layer acts as a lower boundary condition for turbulent suspension; in this boundary layer, the snow particles move at a velocity equal to average wind speed. The suspension layer is defined by the concentration of snow within the layer, which is described by Kind (1992) as the interval between the top of saltation layer and top of the suspension layer where snow concentration is zero. Q (x ∗) = ∫ ϕ (z) dz ∗ (4) 58 The integration limits are the top of the saltation layer, h* (m), and the top of the turbulentsuspension layer, zt (m) where the concentration is zero. ϕt (kg m -3) is the mass concentration of the particulate cloud derived from the concentration ϕr at a reference level ztr. ϕ =ϕ ∗ ∗ ∗ +1 − ∗ ∗ (5) where s is the particle-settling velocity (0.3 m s-1 per Schmidt, 1982) and ϕ* is a concentration scaling parameter, defined by: ϕ∗ = β ∗ ϕ , ϕ = (6) ∗ where β is a unitless coefficient (0.5), and up is a constant horizontal velocity in the saltation layer, which is approximated it to be 2.8×u*. Sublimation also is calculated by the following equation: Q (x ∗ ) = ψ ϕ h∗ + ∫ ψ (x∗ , z)ϕ (x ∗, z)dz (7) ∗ where ψ (s-1) is the sublimation-loss rate coefficient, ϕ (kg m-3) is the vertical massconcentration distribution. The sublimation-loss-rate coefficient describes the rate of particle mass-loss as a function of height within the blowing snow profile and it is a function of (a) temperature-dependent humidity gradients between the snow particle and the atmosphere, (b) the energy- and moisture-transfer by conductive and advective mechanisms, (c) particle size, and solar radiation intercepted by the particle. In addition, it is assumed that the mean particle size declines exponentially with height, the relative humidity follows a logarithmic distribution in which it decreases with height, the air temperature in the snow-transport layer is well mixed and constant with height, while the variables defined in the suspension layer vary with height, and the solar radiation absorbed by snow particles is a function of the solar59 elevation angle (latitude, time-of-day, day-of-year) and fractional cloud cover (Liston and Sturm, 1998). Although the sublimation rate from the surface is fairly small, it can be considerable within the snow transportation zone (Hood et al., 1999; MacDonald et al., 2010). In SnowModel, feedbacks between the sublimation rate and the humidity field are neglected that might exaggerate the sublimation. 3.4. Meteorological data SnowModel requires precipitation, air temperature and relative humidity at 2m height, and wind fields at 10 m height above the surface. MicroMet is able to estimate incoming short and longwave radiation. We provided all required data from both WRF (Skamarock et al., 2008), a high resolution, non-hydrostatic model; and NLDAS, which is a product of land surface model forced by observations in combination with an atmospheric model with a spatial resolution of 0.125˚ (Cosgrove et al., 2003; Mitchell et al., 2004). WRF was set up based on a multi-nesting approach with a resolution of 1 km for the inner domain that covers the study area. The outer domain was forced by meteorological data (wind, pressure, temperature, humidity) and sea-surface temperature data from the ERA-Interim reanalysis data (Dee et al., 2011). We used the same land/atmospheric schemes as Erler et al. (2015) and described in Chapter 2. a) Temperature and precipitation The spatiotemporal meteorological data were derived from WRF (1 km) and NLDAS (~ 12 km) datasets with temporal resolution of 6-hourly and 1-hourly, respectively. In our 60 previous study (Mortezapour et al., 2020), we evaluated the snow model performance using observed SWE time series derived from four snow pillow stations within the study area (Table 3.1). The analysis of results from the meteorological models shows a good correlation of the diurnal cycle for air temperature and precipitation with the coefficient of determination (R 2) of 0.87 and 0.60, averaged over all stations. In general, WRF has a cold bias in winter and MicroMet reproduces this colder temperature with a negative bias of about 1.5˚C, while downscaled NLDAS data by MicroMet reduces this cold bias to -0.5˚C. WRF represents a better estimation of precipitation than NLDAS, but both datasets underestimate precipitation in this area (Table 3.1). While monthly trends are reproduced in most months for all datasets, MicroMet-based precipitation is underestimated at Redfish-Creek and overestimated in Eastcreek and Glacier-NMF. Table 3.1. Snow pillow station characteristics. The observed mean temperature and mean monthly precipitation compared with simulated data from WRF and derived data from NLDAS for the accumulation period 2015-2016. Station ID Station Name Elevation Distance [m] Glacier to Mean temperature Mean monthly [˚C] precipitation [mm] OBS WRF NLDAS OBS WRF NLDAS 2D08P East Creek 1960 15 km to Conrad -3.08 -3.59 -3.39 168 124 98 1173191 Glacier NMF 1850 20 km to Nordic -3.32 -1.56 -2.57 121 100 120 2D14P Redfish Creek 2080 3.4 km to Kokanee -2.41 -3.27 -3.49 207 143 120 1E08P Azure River 1630 8.7 km to Zillmer -3.09 -4.07 -3.73 215 140 113 b) Evaluation of wind fields Wind fields are one of the most important meteorological factors in simulating snow transport. To evaluate wind fields generated by WRF and NLDAS, the results were compared with observations derived from the nearest weather station to each glacier, which cover the period 61 of 1 October 2015 to 30 April 2016. The modelled data at the nearest weather station were estimated based on the nearest model grid point (Table 3.2). Table 3.2. Automated weather station characteristics that have been used for wind evaluation for the 2015-2016 accumulation period. Station Name** The Location Elevation (m) Used Parameters* Distance to Glacier Valemount2 52.78˚N, 119.31˚W 1195 WS, WD 20 km to Zillmer Tabernacle 51.75˚N, -117.76˚W 2438 WS, WD 31 km to Nordic Duncan 50.78˚N, -117.18˚W 1376 WS, WD 12 km to Conrad Slocan 49.78˚N, -117.44˚W 1230 WS, WD 20 km to Kokanee *WS: Wind speed, WD: Wind direction ** Data were provided from BC Ministry of Forest, Lands, and Natural Resource Operations. downscaled wind field data by MicroMet (Table 3.3) are the final wind fields are forced into the drifting snow model. MicroMet distributes wind over terrain using a weighing function based on the slope and curvature of each grid cell. We defined the curvature length scale for each individual domain that the curvature calculation is performed on within the wind submodel. Curvature length scale is approximately one-half the wavelength of the topographic features (from a valley bottom to a nearby peak) within the domain, which varies between the sites from 500 to 900 m. All wind speeds less than 1.0 m s-1 are set to 1.0 m s -1 by MicroMet. The simulated mean and 95th percentile wind speed from WRF and downscaled by MicroMet are underestimated with respect to the measured wind speed, especially at the Tabernacle station, which is located at the highest elevation (2438 m a.s.l.) near a ridgetop (Table 3.3). Table 3.3. Observed and simulated mean, 95 th percentile, and standard deviation of wind speed (m s -1) at 10 m height above the surface by WRF and MicroMet at station locations. Winter 2015-2016 Mean wind speed (m s-1) AWS 1 Valemount2 2 Tabernacle 3 Duncan 4 Slocan Average OBS 4.43 10.27 2.93 4.97 5.65 WRF 2.14 1.57 3.07 1.15 1.98 MicroMet 1.33 1.50 1.54 2.65 1.76 95th percentile of wind speed (m s-1) Standard deviation of wind speed (m s-1) OBS 8.95 15.8 7.25 9.69 10.42 OBS 2.66 7.30 2.41 2.82 3.79 WRF 6.03 4.29 5.73 3.25 4.83 MicroMet 2.35 3.38 2.51 6.96 3.80 WRF 1.84 1.35 1.94 0.96 1.52 MicroMet 0.48 0.88 0.54 2.20 1.03 62 Underestimated wind speed might affect the blowing snow transport rate. The WRF simulated wind variability, represented by its standard deviation, is similar to the observed standard deviation, but MicroMet estimates lower variability than WRF. This difference might be due to the imposition of a minimum wind speed in MicroMet. Wind roses of weather station observations and the modelled wind fields (Fig. 3.2) represent the dominant westerly wind flows at Tabernacle reasonably well, but the magnitude is too low. Wind field derived from WRF is affected by its limited representation of the land morphology (hills, valleys, ridges and slopes, etc.) around the stations. In addition, WRF results are influenced by mismatches between actual and model terrain elevation in mountains that might affect the accuracy of the atmospheric conditions in the lower atmosphere there (Santos-Alamillos et al., 2013; Zhang et al., 2013). To summarize, downscaled wind fields represent a reasonable approximation of dominant wind at each station with a lower mean wind speed compared with weather station data. The MicroMet downscaled wind fields reproduce a similar pattern of wind direction compared with WRF wind fields, but with a higher rate of underestimation. The mean percent bias of downscaled mean and 95th percentile of wind speed is about -10% and -20%, respectively. This underestimation may affect the snow pattern as a result of erosion and deposition as well as snow sublimation. 63 Figure 3.2. Wind roses of observations, NLDAS and WRF simulations, as well as downscaled WRF wind data (all at 10 m height) by MicroMet at four station locations for the period of 1 October 2015 to 30 April 2016. Colors indicate the 6-hourly mean wind velocity (m s -1), the frequency of occurrence is represented in percentage circles. 3.5. Model simulation For this study, we used spatial fields of terrain elevation (a Digital Elevation Model – DEM) from the Shuttle Radar Topography Mission (SRTM) and vegetation cover from Global Land Cover (GLC) at 30 m resolution. We updated the DEMs over the glacier surfaces by 64 replacing SRTM surface elevations with elevations obtained from LiDAR surveys in September 2015. Due to a difference between actual elevation and simulated meteorological data, an elevation bias correction was done on both input datasets (i.e. WRF and NLDAS) using a linear regression equation for temperature (Liston and Elder, 2006b) and precipitation (Pelt et al., 2016) with elevation. We used a monthly temperature lapse rate and a precipitation adjustment factor derived from the monthly climate normal from Parameter-Elevation Relationships on Independent Slopes Model (PRISM; Daly et al., 1994) for 1981-2010. More details on the model simulations and evaluation can be found in Mortezapour et al. (2020). The closest meteorological stations to each glacier (Tables 3.1 and 3.2) provide evaluation data for model simulations in winter 2015-2016. Lack of observational data prevents us from evaluating modelled incoming shortwave and longwave radiation. To evaluate simulated wind fields with WRF and derived wind speed and direction from NLDAS, another four stations with wind data were chosen (Table 3.2). The seasonal snow distribution over each glacier was measured by LiDAR and in-situ stake probing (Pelto et al., 2019). The simulations were done twice for each glacier, with and without consideration of wind effects, for two winter seasons 2013-2014 (calibration) and 2015-2016 (validation) from 1 October through 30 April. The flowchart of simulations is illustrated in Fig. 3.3. The simulated snow depths and winter mass balances are evaluated against the probed snow depths and LiDAR measurements over glaciers. The average winter accumulation over each glacier was estimated by differencing coregistered DEMs derived from LiDAR measurements (Pelto et al., 2019) at the end of summer (early September 2015), and the end of winter (end of April 2016). Details concerning the probing data and LiDAR measurement analysis can be found in Pelto et al. (2019). 65 Figure 3.3. SnowModel flowchart with and without snow transport by wind (SnowTran3D module) 3.6. Results 3.6.1. Patterns of snow erosion and deposition First, we analyze the snow depth time series at the station locations to assess how wind might influence the snow depth at specific points. Daily simulated snow depth was compared with observed snow pillow data at each station from October 2015 to April 2016. Those statistics (Table 3.4) show a negligible difference between the results with and without accounting for snow redistribution by wind at the station locations. 66 Table 3.4. Statistical analysis of simulated daily snow depth and SWE time series using SnowModel with and without wind for four snow pillows for the accumulation period 2015-2016 SnowDepth [m] WRF_forced With wind WRF_forced Without wind SWE [m w.e.] WRF_forced With wind WRF_forced Without wind Function/Station ID Correlation (p<0.01) Absolute percent bias Bias Mean absolute error RMSE Correlation (p<0.01) Absolute percent bias Bias Mean absolute error RMSE East-creek 0.76 11.99 -0.06 0.21 0.25 0.76 12.03 -0.07 0.21 0.25 Glacier-NMF 0.78 8.26 0.03 0.16 0.19 0.78 8.71 0.005 0.17 0.20 Redfish 0.71 21.27 0.05 0.44 0.56 0.72 21.28 0.05 0.44 0.56 Azure river 0.76 13.63 0.05 0.19 0.23 0.77 13.43 0.003 0.19 0.22 Average over all Correlation (p<0.01) Absolute percent bias Bias Mean absolute error RMSE Correlation (p<0.01) Absolute percent bias Bias Mean absolute error RMSE 0.72 15.3 -0.09 0.07 0.12 0.71 16.2 -0.10 0.08 0.13 0.86 7.00 0.12 0.11 0.10 0.86 7.90 0.11 0.10 0.12 0.87 12.6 -0.04 -0.03 -0.03 0.85 18.6 -0.02 0.02 0.04 0.75 4.30 0.05 0.04 0.10 0.75 5.20 0.06 0.04 0.08 0.80 9.80 -0.08 0.06 0.09 0.79 11.9 -0.07 0.06 0.09 0.75 13.79 0.02 0.25 0.31 0.76 13.86 -0.002 0.25 0.31 The simulated daily snow depth timeseries forced by WRF are well correlated (r = 0.75; p < 0.05) with observed snow pillow data. The model is able to reproduce the snow depth with the average mean absolute error and RMSE of 0.25 and 0.31 m water equivalent (m w.e.), respectively. Likewise, goodness of fit statistics remained almost the same with or without wind with only small differences in percent bias (Table 3.5). Wind still plays a significant role in spatially redistributing snow within a glacier, however. Using high resolution LiDAR measurements, the snow depth averaged over the area of each glacier around the time of peak accumulation was evaluated and compared with the simulated winter accumulation over each glacier (Table 3.5). Additional evaluation statistics for the accumulation period 2013-2014 are shown in Table S2.1. 67 Table 3.5. Comparison of integrated snow depth forced by WRF and NLDAS over the whole domain of each glacier with LiDAR measurements and the related percent bias for the 2015-2016 accumulation period. Averaged snow depth (m) LiDAR WRF forcing _with wind %Bias WRF forcing _without wind %Bias NLDAS forcing _ with wind %Bias NLDAS forcing _ without wind %Bias Conrad 02.90 03.31 14.10 03.44 18.60 02.24 -32.51 02.23 -32.66 Nordic 02.41 03.04 26.12 03.19 32.23 02.26 -06.31 02.28 -05.60 Kokanee 03.02 02.90 -03.91 02.90 -03.91 03.00 -00.78 02.99 -00.88 Zillmer 03.16 03.19 00.98 03.23 02.22 02.22 -29.93 02.26 -28.57 Averaged over all 02.87 03.11 08.36 03.19 11.15 02.43 -15.33 02.44 -14.98 The averaged snow depth over each domain simulated by NLDAS forcing is underestimated while mean WRF-forced snow is overestimated at all locations, except at the Kokanee Glacier (Table 3.5). On average, the snow accumulation results using WRF forcing improve the simulation results relative to using NLDAS data (Table 3.5). We found negligible differences between glacier-averaged simulated snow depth forced by NLDAS with wind and without wind effects for each glacier. Therefore, we continue our discussion of NLDAS forced data only with considering wind effects. Caution is required when evaluating the spatial pattern of snow accumulation for each glacier (Fig. 3.4) derived from LiDAR since the elevation change between the two DEMs is not corrected for ice dynamics (Pelto et al., 2019). Ice dynamics ( emergence and submergence ice velocities) would need to be accounted for in order to estimate snow depth by differencing the autumn and winter DEMs. Observed winterime submergence and emergence velocities for these glaciers respectively averages about 0.6 to 1.8 m for the accumulation and ablation zones of the glaciers (Pelto, 2020). These velocities, however, strongly vary in space, and correcting for ice dynamics is beyond the scope of this study. 68 Figure 3.4. Comparison of simulated snow depth forced by WRF with and without wind effects, and simulated snow pattern forced by NLDAS with wind, along with LiDAR footprint measurements at four glaciers for the accumulation period 1 October 2015 to 30 April 2016. The last column shows the subtracted LiDAR measurements from WRF derived simulated snow depth with wind effects. Considering wind effects improves the integrated snow depth over the area forced by WRF data at all locations, while the calibration of the model for each individual glacier separately might enhance the results. In this study, the chosen parameters were applied to all sites within SnowModel. Overall, WRF performed 6% better than NLDAS in terms of the mean snow depth bias averaged over all of the glaciers (Table 3.5), and simulations with wind transport 69 reduced the averaged snow depth bias by 3%. To remove the complicating effects of vertical ice dynamics (e.g. submergence or emergence) affecting apparent snow depth over the glaciers, we compared simulated and observed snow depth over ice-free terrain. Boxplots determine the simulated snow pattern and the effect of wind on the snow redistribution (Fig. 3.5). Figure 3.5. Measured vs. simulated WRF-forced snow depth distribution with elevation over glacier free areas, surrounding each glacier for the 2015-2016 accumulation period. Box heights denote 25 th and 75 th percentile. The whiskers extend to the highest and lowest snow depth at each elevation band. The sign (×) shows the medians, and dots represent the outliers. The sample sizes differ for each glacier and elevation band ranges from 20 (e.g. at low-elevations in small glaciers) to more than 8000 samples (e.g. at mid-elevations in large glaciers). For all locations, snow depths without wind increases with elevation. Comparing visualized snow depth distribution with and without wind effects, SnowModel forced by WRF represents the main features of snow depth distribution in glacier free areas (Fig. 3.4). The model simulates erosion at the ridges and deposition on the lee sides of mountains. However, a snow 70 depth difference map of modelled and LiDAR data (Fig. 3.4) illustrates that deposited snow amount on the immediate leeward side of ridges with steep slopes is high in non-glacier areas and may be overestimated. The net transport values over each glacier shows the highest snow reduction at Nordic (Fig. 3.6), while there is a positive net transport for the Zillmer Glacier during the simulation period. The wind graphs (Fig. 3.2) represent the highest westerly monthly wind speed which occurred in November and February with the highest wind speed over Conrad. The model suggests the wind transport of snow over glacier areas to be -30 mm, -20 mm and +7 mm for Nordic, Conrad, and Zillmer, respectively, while the Kokanee glacier experiences negligible snow transport summed over the accumulation period. However, the amount of transported snow within the domains is highly variable. At smaller scales, the summed snow transportation reaches more than 14 m in the lee sides of steep slopes with the lowest snow deposition for Kokanee and largest snow deposition at Conrad Glacier. We also evaluated the simulated SWE distribution with elevation against probing data and pit measurements for about 30-50 points for each glacier (not shown). More details of simulated snow density and pit measurements can be found in Mortezapour et al. (2020) and Pelto et al. (2020). We found that NLDAS-forced SWE was underestimated at all locations with an average of 20% relative to WRF forcing data. The integrated variability of WRF-forced SWE over each glacier area is slightly underestimated by SnowModel, but it is more realistic considering wind transport at all glaciers, than without having wind (Table 3.6). 71 Figure 3.6. Mean monthly drifting snow sublimation, solid precipitation per time step, and mean monthly wind speed at 10 m height above the surface (Left); Relative humidity with respect to ice, accumulated transported snow, and monthly wind direction (Right) averaged over the main glaciers (Fig.3.8) during the 2015 -2016 accumulation period. Table 3.6. Winter mass balance (SWE) comparison between simulations and observations (stake and LiDAR measurements) on the glaciers for the 2015-2016 accumulation period. With-wind Winter Mass Balance (m.w.e) Probing Zillmer Without-wind LiDAR WRF forced with wind WRF forced no wind %Bias -Prob %BiasLiDAR %BiasProb %BiasLiDAR 1.99±0.23 1.68±0.19 1.75±0.20 1.72±0.12 -12.0 +4.1 -15.3 0.0 Nordic 1.79±0.14 1.79±0.22 1.92±0.15 1.96±0.10 +7.2 +7.2 +9.5 +9.5 Conrad 1.88±0.12 1.40±0.18 1.73±0.16 1.62±0.11 -7.9 +17.6 -13.8 +15.7 Kokanee 2.07±0.13 1.98±0.22 1.65±0.17 1.67±0.09 -20.2 -16.6 -19.3 -15.7 1.93±0.16 1.71±0.20 1.76±0.16 1.74±0.11 -8.8 +1.7 -10.0 +1.57 Average all over 72 3.6.2. Surface and drifting snow sublimation Due to minor differences between results forced by NLDAS with and without wind, we only discuss the simulated sublimation forced by WRF. The highest modelled monthly average drifting sublimation rate was 0.12 mm per day at Conrad Glacier (Fig. 3.6). The lowest drifting sublimation and snow transport was at Kokanee (Fig. 3.7). The averaged surface sublimation over Conrad and Zillmer domains are in a similar range as snow drift sublimation, while surface sublimation shows a higher rate at the two other glaciers, especially at Kokanee. Drifting snow sublimation is higher in winter than spring, but surface sublimation is stronger at the end of the accumulation season (Fig. 3.7). Figure 3.7. Cumulative mean reduction of SWE (mm w.e.) averaged over the whole domain due to static sublimation at the surface (black, dashed), and drifting snow sublimation (black, solid) and the total sublimation (red, solid) during the complete simulation for four glacier locations for the 2015-2016 accumulation period. 73 Scatterplots of snow blown sublimation and meteorological data show an expected increase in drifting snow sublimation with wind speed (Fig. 3.8). Sublimation of drifting snow is limited to some areas with a distinctive pattern of maximum values on the windward slopes and minimum values occurring in the valley locations (Fig. 3.9b). The effects of drifting snow sublimation occur on the west and northwest- (i.e. windward) facing slopes for all sites, whereas the surface sublimation happens throughout the whole domain, where snow is present (Fig. 3.9a). Figure 3.8. Scatter plots of daily drifting snow sublimation rate against air temperature, relative humidity with respect to ice, wind speed and wind direction at 10 m height above the surface for the 2015-2016 accumulation period over four glaciers. 74 The surface sublimation also is higher on the north-facing slopes for all glaciers, except at Nordic, which shows the highest values on the west-facing slopes. The highest drifting snow sublimation occurs at the ridges (compare the topography in Fig. 3.1 with Fig. 3.9b), where it can reach 500 mm w.e. for Zillmer and Conrad glaciers. Peak drifting snow sublimation values are lower for Nordic (180 mm w.e.) and the lowest occurs at Kokanee glacier with 45 mm w.e. The average drifting snow sublimation over each model domain is 0.34%, 0.90%, 1.7%, and 2.1% compared with accumulated SWE at the same time for Kokanee, Nordic, Zillmer, and Conrad, respectively during the entire accumulation season. In addition, the average of drifting snow sublimation over all sites on the glaciers is slightly higher than surface sublimation (35 vs. 31 mm w.e.) during winter accumulation. Figure 3.9. Summed static surface (a) and drifting snow sublimation (b) at Conrad, Kokanee, Nordic, and Zillmer glaciers (from left to right) in [mm w.e.]. Fluxes are cumulated over the simulation period (the 2015 -2016 accumulation period). The black and red outlines show the LiDAR footprint and main glaciers boundaries, respectively. 75 Generally, the winter mass balance is underestimated by SnowModel for both WRF and NLDAS, with and without wind effects on the glaciers compared with manual measurements and LiDAR data (Tables 3.5 and 3.6). The percent difference between winter mass balance wind and without wind for the glaciers is less than 4% for all sites, except on the Conrad which is 6% compared with stake measurements. The blowing snow sublimation causes an averaged reduction of less than 2% winter mass over all glaciers. To summarize, the simulations suggest that the high drifting snow sublimation occurs mostly at the ridges and windward slopes and is limited to small areas. Consequently, drifting snow sublimation may not affect glacier winter mass balance significantly in alpine terrains (Strasser et al., 2007). 3.7. Discussion Wind is a key factor in snow distribution over complex terrain (Dadic et al., 2010). The wind field generated by WRF are more representative at alpine sites (e.g. Tabernacle station) compared with the stations located in the valleys, where wind is affected by local topography and thermally forced winds. The wind speed interpolated by MicroMet is lower on average than the forcing data which likely contribute to model error. Underestimation of wind velocity by MicroMet was also reported in other studies (e.g. Gascoin et al., 2013; Musselman et al., 2015; Ravazzani et al., 2016). Therefore, using a more comprehensive method to simulate mechanical effects of terrain on windflow (e.g. Marsh et al., 2020; Vionnet et al., 2020) might lead to more reliable results. 3.7.1. Snow transport We found a negligible difference between simulated snow depth forced by NLDAS with and without wind. This might be due to the coarse resolution of the wind field in NLDAS data 76 which is not able to capture the dominant direction of wind in the complex terrain. Comparisons with LiDAR data display more realistic snow patterns, when the snow transport module is activated with WRF driven winds (Prasad et al., 2001; Gascoin et al., 2013). The WRF forced pattern of snow erosion and deposition is consistent with wind speed and direction over each glacier (Fig. 3.10). Figure 3.10. Downscaled mean wind speed (ms -1) at a height of 10 m above the surface and direction over each domain by MicroMet for the 2015-2016 accumulation period. Snow erosion occurs mostly along west facing slopes that experience stronger wind speeds, while snow is deposited in the leeward sides. We showed a negative net transport of snow by wind, except at Zillmer where the positive net accumulation may be due to snow from 77 surrounding snow-covered slopes. The spatial heterogeneity of simulated snow depth is increased when snow is redistributed by wind. However, the percent bias of integrated mean snow depth over area and winter mass balance of glaciers are not reduced considerably. Wind also may contribute to the release of avalanches affecting winter mass balance, especially in small glaciers. Generally, SnowModel is capable of producing winter mass balance with tolerable uncertainties using WRF forcing data. Our results in terms of amount and patterns of transported snow and drifting sublimation are consistent with Déry et al. (2010) over the Cariboo Mountains, MacDonald et al. (2010) and Musselman et al. (2015) over ridge-tops in the Canadian Rockies and Zwaaftink et al. (2011) over the whole simulation domain in the Swiss Alps using different models and forcing data. 3.7.2. Snow sublimation Wind speed plays a key role in blowing snow sublimation. Higher wind speeds result in a higher rate of mass loss due to drifting sublimation, which is proportional to the square root of wind velocity in SnowModel. High-speed winds create a fully mixed layer above the snow surface that maintains gradients of vapour pressure and facilitates sublimation (Stigter et al., 2018). The strong surface sublimation at the end of the accumulation period is likely related to higher air temperature and more energy available from shortwave radiation. The variability of drifting snow sublimation over the domains is analyzed with respect to the weather conditions, particularly wind direction. The highest rate of drifting snow sublimation occurs where there is a combination of conditions with a warmer temperature, lower humidity environment, and a higher wind speed (Pomeroy and Jones, 1996; Liston and Sturm, 2004; Svoma, 2016) on the windward slopes (N and NW facing slopes) and the top-ridges, where snowfall and wind erosion is significant (Fig. 3.9). All these conditions together lead to large 78 vapour deficits where snow particles are present, resulting in a higher rate of sublimation. The blowing snow sublimation strongly varied over space and time with a range of 1-45% winter mass loss in respect to SWE with the highest rates on the ridges. There is a large range of reported percent mass loss (1-90%) due to drifting snow sublimation in alpine areas in different studies that might be due to different climate conditions, scales, and using the different model setup and resolution (e.g. Pomeroy et al., 1997a; Déry and Yau, 2001; Molotch et al., 2007; Strasser et al., 2007; MacDonald et al., 2010; Bernhardt et al., 2012; Gascoin et al., 2013; Zwaaftink et al., 2013) and less than 6% over an entire domain (e.g. Vionnet et al., 2014). In our study, the greatest magnitude of drifting snow sublimation occurs on the windward slopes, which is consistent with results from Musselman et al. (2015). In contrast, Zwaaftink et al. (2011) demonstrate the strongest blowing snow sublimation occurs on the leeward slopes, due to different topographical features and climate conditions in which adiabatic warming and drying sustain sublimation. Furthermore, Zwaaftink et al. (2011) used a steady state field of the mean wind over the region during a storm event, which could neglect the dominant pattern of wind and drifting sublimation on the region. While underestimated wind speed by MicroMet might affect drifting snow sublimation rate, it can be compensated with neglecting the selflimiting aspect of blowing snow sublimation (thermodynamic feedback) within SnowModel which exaggerates the blowing snow sublimation. 3.7.3. Limitations Uncertainty in the results stems from inadequate information both from model and data limitations. In terms of input data, the simulated results are influenced by systematic errors and uncertainties associated with parameterizations in regional climate models. Limited observations make the performance of the model difficult to evaluate. Specifically, there are 79 few observations of wind in this region, especially over the glaciers we are most interested in. Due to the lack of observations, longwave and shortwave radiation from WRF could not be evaluated. Furthermore, no observations of sublimation rates exist for the glaciers of this study. The model structure and parameterization also might affect the outcomes. Micromet reduces the already too low wind speed from WRF, resulting in underestimation of wind speed in this region. SnowModel does not take into account preferential deposition of snow during precipitation events, and feedbacks between the sublimation rate and the humidity field. Furthermore, in SnowTran3D, the effect of suspended snow due to flow separation along steep ridges, which occurs in clear, windy days on the mountain peaks is neglected by assuming an equilibrium between transport flux with the near-surface winds (Liston et al., 2007). These simplifications might affect the amount of drifting sublimation, with the expectation of a higher rate of actual sublimation in alpine environments. The findings can be also influenced by calibration parameters such as curvature and wind shear in MicroMet. In this study, the model sensitivity to the curvature parameter and slope weight was tested for one glacier. The results show a slight influence of curvature on the snow transported by wind as well as both static and drifting sublimation (<5%), contributing to uncertainty in the model results. In this study, we used a constant wind shear threshold to initiate snow transport; evolution of this parameter during simulation periods should be tested in future studies that may improve the results. The wind shear threshold is a function of wind speed and temperature. The snowpack consolidates over time, consequently rising temperature in the spring will increase the minimum wind shear stress required to initiate snow transport (Gascoin et al., 2013). In addition, the gravitational redistribution (i.e. avalanche and snowslides) and its effect on the winter mass balance calculation, have not been taken into consideration. Avalanching is able to redistribute mass 80 from areas above the glacier onto its surface. Gravitational redistribution may be significant and should be considered in future research. 3.8. Conclusion The main aim of the study was to determine the effect of wind on winter mass balance estimation using NLDAS and WRF forcing dataset. We found negligible differences in winter mass balance and snow depth when forcing SnowModel and its snow transport module with NLDAS data. It is suspected that these minor differences arose from the coarse resolution of wind fields. In contrast, SnowModel forced by WRF data reasonably represents the pattern and magnitude of snow depth over the alpine terrain. Generally, the averaged simulated winter mass balance forced by WRF is underestimated by about 10% and overestimated by 1.7% compared with probing and LiDAR measurements, respectively. The inclusion of wind improved the simulated winter mass balance by 2%. The simulated maps of snow depth reveal windward erosion of snow at the ridges and deposition on the lee sides of steep terrain, which appears to be in good agreement with the general pattern of snow accumulation revealed by LiDAR measurements over ice-free areas. The variability of simulated snow depth is less than stake observations suggest, however. Over the accumulation season, there was a net loss of wind transported snow when averaged over all glaciers, with largest losses (8 cm) at Conrad Glacier. The model results demonstrate strong drifting snow sublimation occurs for those sites with a higher mean wind speed. The greatest reduction of snow is on the windward slopes due to the presence of more snow and exposure to the westerly winds. Averaged over the whole simulation domain within a winter season, drifting snow sublimation is less than 2% of SWE, however at individual locations it can be considerably higher reaching up to 0.5 m (45%) over the accumulation period. We currently lack observations to evaluate simulated snow 81 sublimation, however. We suggest that although drifting snow sublimation is not negligible at smaller scales and varies significantly in time and space, the effect of wind transport on the winter mass balance is not considerable when averaged over the domain and it may be more important for smaller (e.g. cirque) glaciers on the leeward side of mountains. Given the fact that simulation of these processes is computationally expensive, considering wind effects on estimation of glacier-wide average winter mass balance is probably not required for glaciers in the Columbia Basin. However, running the model for a full annual cycle, including the ablation season, may show different results and may need higher complexity models to better represent ablation processes. More accurate modelling of wind fields that captures the local winds, higher resolution meteorological fields, in particular precipitation, and a higher resolution of simulation that can influence snow transport may lead to a more detailed snow pattern and a higher accuracy of winter mass balance. In addition, the effect of gravitational redistribution on the winter mass balance estimation needs to be addressed in future research. 82 Chapter 4. Snow Changes over the Columbia Mountains, Canada under a Future Climate Change Scenario using a Distributed Snow Model Mortezapour, M.1, Menounos, B.1, Jackson, P. L.2, Erler, A. R.3 1. Natural Resources and Environmental Studies Program, Geography department, Prince George, British Columbia, Canada 2. Natural Resources and Environmental Studies Program, Atmospheric Science department, Prince George, British Columbia, Canada 3. Aquanty, Waterloo, Ontario, Canada 83 Abstract In western North America, many communities rely on runoff from mountain snowpacks. Projections of how future climate change will affect the seasonal snowpack are thus of interest to water managers, communities, and policy makers. We investigated projected changes in seasonal snowcover under global warming for the 21 st century and the Canadian portion of the Columbia River basin using a physically-based snow distribution model (SnowModel) at 500 m horizontal resolution. Forcing data for the reference (1979-1994) and future (2045-2059, 2085-2099) periods originate from a 4-member initial condition ensemble of global Community Earth System Model (CESM1) simulations based on the Representative Concentration Pathway (RCP 8.5) scenario. The ensemble has been dynamically downscaled (DD) to 10 km resolution using the Weather Research and Forecasting model (WRF). In addition, we also evaluated the performance of SnowModel using publicly available, statistically downscaled (SD) temperature and precipitation projections. We observe a 38% and 30% decrease in WRF-simulated snow depth and SWE, respectively by the end of this century relative to the reference period. Snow depth and snow water equivalent (SWE) are most affected for elevations below 2000 m asl. The averaged SD-forced snow depth and SWE over the area represents a reduction of 28% and 15% toward the end of the 21st century. Our results indicate that the projected loss of snowpack depends largely on elevation and season, while stronger projected thinning of the snowpack occurs for the dynamically downscaled simulations by the end of the 21st century compared to the statistically downscaled simulation. Nevertheless, both SD- and DD-forced simulated snow depth and SWE show losses in most areas by the end of the century. This has significant implications for water resource availability for various sectors and communities living in these areas 84 4.1. Introduction High latitudes and altitudes are most affected by recent climate change, and changes in seasonal snow is one indicator of this change (Prowse et al., 2006; López-Moreno et al., 2012; Milanes et al., 2018). Mountain snowpacks act to store precipitation during winter, and later, as the weather warms, melting snow provides water supply for many downstream communities in western North America. Therefore, projections of how future climate change will affect water resources is of high interest to policy makers (Erler et al., 2017). Global climate change and its impact on the cryosphere is reported and summarized in a report by the Intergovernmental Panel on Climate Change (2019). Studies indicate that areas of British Columbia (BC) have experienced warming at a rate twice the global average (Gayton, 2008; Warren and Lemmen, 2014); this warming shrank glaciers and thinned seasonal snowpacks during the 20th century (Larsen et al., 2007; Schiefer et al., 2007; Bolch et al., 2010; Menounos et al., 2019). The climate of BC is influenced by maritime and continental air masses and dominated by heavy orographic precipitation (Burn, 2008; Mote and Salathé, 2010). Consequently, climate modelling and downscaling for this area is affected by its complex topography and is therefore sensitive to model resolution (Salathé et al., 2010). Local climate stems from the evolution of weather systems, resulting from the interaction between the atmosphere, land-surface, ocean, snow, and ice at the specific location (Gurgiser et al., 2013; Khadka, 2013). Global Climate Models (GCMs) are too coarse to resolve these high-resolution meteorological conditions and thus, their sole use on understanding how future climate may affect snow and glaciers in mountain environments is limited (Radić and Hock, 2006; Radić and Clarke, 2011). In addition, the projected snow cover changes that have been investigated by the direct output of GCMs has a limited ability to simulate the changes in higher elevations due to their coarse 85 resolution (Mankin and Diffenbaugh, 2015). To bridge the gap between GCM resolution and the spatial resolution needed to resolve mass and energy fluxes in topographically complex terrain, downscaling strategies are commonly employed (Maraun et al., 2010). Statistical downscaling uses statistics derived from empirical relationships between simulated large-scale atmospheric variables by GCMs and observed meteorological variables on a local scale. Statistical downscaling has low computational demands, and, due to the limited availability of observed data, this method is often applied to downscale precipitation and air temperature (Maraun et al., 2010; Teutschbein et al., 2011). Furthermore, statistical downscaling assumes that the relations between coarse resolution variables and those at a local scale that currently exist will also exist in the future. This assumption may not be valid in a changing climate, especially if climate change alters regional climate dynamics. In contrast, dynamical downscaling uses a physically-based model that solves equations for conservation of mass, momentum and energy in the atmosphere and typically includes interactions with the land surface to achieve a finer resolution for climate projections. This method is able to represent detailed information of land-atmosphere interactions, especially orographic precipitation, but is computationally expensive. Given the computational expense, dynamical downscaling has been used in fewer studies of snow regime simulations requiring higher resolution, than statistical downscaling approaches (Christensen et al., 2007; Bavay et al., 2009; Marty et al., 2017; Ishida et al., 2018). Energy-balance models are physically-based, improving simulation by considering the spatial and temporal contributions of energy into the snowpack (e.g. radiation, sensible and latent heat exchange, ground and precipitation heat fluxes) to get the required energy for melting and its response to climate changes and surface conditions. Such studies are important for accurate prediction of snow cover variability in mountainous areas due to the complex topography (Arnold et al., 1996). 86 There are several studies investigating the climate change impacts on the snow cover in North America. A declining trend of snow depth between 1948 and 2006 was reported by Park et al. (2012) in the Arctic with a stronger negative trend in North America than in Eurasia. Using remotely sensed data Choi et al. (2010) revealed that snow area in Western Canada decreased from 1972 to 2008, while the core snow season (consistency of snow cover duration) has increased by approximately 10 days over extreme Western Canada. Kang et al. (2014) report a 19% decline of snow melt contribution to Fraser River runoff over the period 19482006. High-resolution dynamic downscaling using WRF has been used in several studies of current and future climate over North America (e.g. Rasmussen et al., 2011; Liu et al., 2019; Musselman et al., 2018). All studies reported snowfall reduction under a future, warmer climate. This method was also used by Rasmussen et al. (2014) to project the future changes in snowpacks over the Colorado Headwaters under the IPCC A1B scenario by 2050. Their results depicted increased precipitation in future winters by 12% with a decrease in snowfall fraction from 0.83 to 0.74 resulting in a reduction of SWE with runoff occurring 2-3 weeks earlier in the future simulations. More liquid winter precipitation and earlier snow melt in Western Canada was also projected by Erler et al. (2015, 2017) using dynamically downscaled CESM1-WRF, which is also used in the present study. There are many other studies on snow cover projections in the Alps with different Regional Climate Models (RCMs) and emission scenarios (e.g. Rousselot et al., 2012; Steger et al., 2013; Marke et al., 2014; Marty et al., 2017). The projected impact of climate change on the snow hydrology of the Fraser River Basin (British Columbia) and Athabasca river basin (Alberta) has been investigated using different statistical downscaling approaches (e.g. Kerkhoven and Gan, 2011; Shrestha et al., 2017; Leong and Donner, 2015; Islam et al., 2017). Their findings also indicate an increased mean precipitation for the 2050s, while the snowfall is expected to decrease. 87 We investigate how future climate could affect the distribution of snow depth and SWE in the Columbia Basin as simulated using SnowModel (e.g., Liston and Elder, 2006a, 2006b; Liston et al., 2007). SnowModel is forced by dynamically downscaled data (four ensembles) generated by CESM1-WRF and statistically downscaled temperature and precipitation from CESM provided by the Pacific Climate Impacts Consortium (PCIC, 2019). The simulations are done for a historical period and two projected periods based on the RCP 8.5 climate scenario (section 2). The validation of the model and general patterns of climate change are provided in section 3. The results are interpreted in section 4, and the possible uncertainties are discussed in section 5 followed by a conclusion, which summarizes the results of the study. 4.2. Data and Methods 4.2.1. Study region The Canadian part of the Columbia Mountains is bounded by the Rocky Mountain Trench (east), the Columbia River (south), the Interior Plateau (west), and the Fraser River (north) with the peaks higher than 3000 m above sea level (a.s.l.) (Fig. 4.1). The study area makes up 15% of the whole Columbia Basin (668,000 km²), which spans seven US states and the province of British Columbia in Canada including rivers and lakes, along with population centres of different sizes. 88 Figure 4.1. Digital elevation map of the Columbia Mountains. Weather stations used for simulations are marked with points. The snowmelt-driven discharge of the rivers, which plays a key role in water supply in the region, declines rapidly in the late summer. The heavy snowfall in the winter is controlled by mid-latitude cyclones with wet and mild westerly air masses coming from the Pacific Ocean that are intercepted by the Columbia Mountains. High annual precipitation, deep snow accumulation and relatively moderate winter temperatures are characteristic of this region. According to the Parameter-Elevation Relationships on Independent Slopes Model (PRISM) dataset maps (4 km resolution), which is an estimation based on climate station data (19611990), the annual mean temperature in high altitudes (elevation range between 420 to 3700 m a.s.l) in the northern and central Columbia Basin is below freezing and colder than -4˚C. The 89 annual precipitation on the windward side of the ranges has been recorded between 1000 and 2500 mm yr-1, whereas on the leeward side precipitation is reported between 250 and 750 mm yr-1 (Murdock and Werner, 2011). 4.2.2. SnowModel We simulated the evolution of seasonal snow using a spatially distributed modelling system, SnowModel (Liston and Elder, 2006b). The model consists of four modules: SnowPack, a snow depth and SWE evolution model (Liston and Hall, 1995); surface energy exchanges are calculated by EnBal (Liston, 1995; Liston et al., 1999); MicroMet is a sub-model to define the meteorological conditions (Liston and Elder, 2006b); and SnowTran-3D calculates the snow redistributed by wind (Liston and Sturm, 1998, 2002; Liston et al., 2007). SnowTran-3D works for high resolution studies (200 m or finer) but is not used in this study. Temporally distributed meteorological data are required for simulations, which we describe in the next section. MicroMet was used for spatial interpolation of meteorological data including air temperature, incoming long wave radiation, incoming solar radiation, precipitation, relative humidity, surface pressure, wind direction and wind speed. 4.2.3. Input data Spatially distributed land surface data are required for SnowModel. We used a digital elevation model from the Shuttle Radar Topography Mission (SRTM; Farr et al., 2007) at 500 m resolution and land cover data were taken from the Global Land Cover database (GLC; Chen et al., 2015) reclassified based on the SnowModel vegetation cover classification (500 m resolution). a) Dynamically downscaled data 90 Spatio-temporal meteorological fields were provided from dynamically downscaled GCM simulations (Erler et al., 2015, 2017) to force the snow model. Due to high computational demand, only one GCM model was used to downscale climate projections for Western Canada. The Community Earth System Model (CESM1) was used to generate an ensemble of four simulations. CESM1 is a fully coupled global climate model with interactive components of the atmosphere, land, ocean, and sea ice with the nominal resolution of the standard configuration of about 1˚ (Erler et al., 2015). Dynamically downscaled meteorological projections are based on an initial condition ensemble of four independent global forecasts from the same GCM. However, the RCM/WRF ensemble members differ in both initial and boundary conditions (including sea surface temperatures and sea ice, which are prescribed from CESM1). Due to the high computational cost of the regional model, three 15-year segments from the GCM were dynamically downscaled with WRF (V3.4.1) (Erler et al., 2017). Each ensemble member was downscaled separately over Western Canada, using a nested configuration consisting of an outer domain of 30 km resolution and an inner domain at 10 km resolution. The ensemble was integrated for a historical validation period (1979-1994: 1980s) and two projection periods for mid and end of the century (2045-2059: 2050s and 2085-2099: 2090s) based on the representative concentration pathway 8.5 (RCP8.5), a high-emission scenario, for the future trajectory of greenhouse gases. In this scenario the CO2 and CH4 concentrations are 541 ppmv and 2740 ppbv in 2050 and 936 ppmv and 3751 ppbv in 2100, respectively (Erler et al., 2017). The required dynamically downscaled data (DD) for snow simulation derived from WRF include air temperature, precipitation, wind speed and direction, and incoming radiation (short and long wave). The 6-hourly downscaled CESM1-WRF values were averaged to create daily mean values and the simulations were performed based on the daily averaged atmospheric forcing for the whole periods. 91 b) Statistically downscaled data In addition to dynamical downscaling, we also forced SnowModel with statistically downscaled temperature and precipitation data from CESM1 provided by PCIC (PCIC, 2019). The statistically downscaled data are based on historical GCM simulations and gridded historical climate station observational data used to develop the statistical downscaling transfer functions. Statistically downscaled data are based on the Bias Correction/Constructed Analogues with Quantile mapping reordering approach (BCCAQ), which is a combination of the Bias Correction/Constructed Analogues (BCCA; e.g. Maurer et al., 2010), and quantile mapping (QMAP; e.g. Gudmundsson et al., 2012) methods. Daily gridded downscaled precipitation and temperature are accessible at 10 km resolution, which was input to SnowModel with a final resolution of 500 m. For statistically downscaled (SD) data, we used the daily precipitation rate and calculated the average temperature from daily maximum and minimum temperatures downloaded from PCIC. We allowed SnowModel to calculate radiation fluxes for the SD approach and the rest of the variables (wind components and specific humidity) were kept the same as in the DD data. 4.2.4. Station data and forcing data correction The assessment of the changes in temperature and precipitation from DD data is based on four ensemble combinations. For this study, we used meteorological stations over the period 1979-1994, for the entire domain (Table 4.1) to evaluate the performance of downscaled data before using them as input for SnowModel. After initial quality checks, we selected 10 stations based on their altitude and latitude to cover nearly the entire basin. Erler et al. (2015, 2017) provide a validation of the CESM1-WRF simulations for the historical period. 92 Table 4.1. Station information used for meteorological evaluation of forcing projections and evaluation of simulated Snow depth by SnowModel for the historical period (1979-1994). Station Name Bugaboo Creek Lodge Coordinates 50.75˚N, 116.71˚W Cariboo Lodge Duncan Lake Dam East Creek Fish Lake Fred Laing Ridge* Glacier Np Mt Fidelity Nakusp Valemount North 52.71˚N, 119.47˚W 50.23˚N, 116.97W 50.63˚N, 116.93˚W 50.04˚N, 117.17˚W 52.04˚N, 118.57˚W 51.23˚N, 117.70˚W 50.25˚N, 117.8˚W 52.85˚N, 119.25˚W Yoho NP Wapta Lake *No summer records 51.45˚N, 116.33˚W Network Name ECCC3 Used data time period Altitude (m) 1979-1994 1529 ECCC 1980-1994 1979-1994 1980-1994 1979-1994 1989-1994 1979-1994 1979-1994 1979-1989 1095.8 548.6 2030 1070 1080 1890 457 892 1979-1994 1646 ECCC ENV-ASP1 MoTIm2 MoTIe4 ECCC ECCC ECCC ECCC 1 BC Ministry of Environment - Automated Snow Pillow Network BC Ministry of Transportation and Infrastructure (manual) 3 Environment and Climate Change Canada 4 BC Ministry of Transportation and Infrastructure (electronic) 2 The winter precipitation is overestimated, and there is a consistent cold bias at higher elevations in both CESM1 and WRF (Erler et al., 2015). To bias correct the CESM1-WRF output we compared the averaged monthly simulated temperature and precipitation against the monthly observed data over all stations (Fig. 4.2). An elevation difference correction was performed on DD data in the same way as in our previous work (Mortezapour et al., 2020), before comparing them with observations. The average precipitation bias was +30 mm (+27%) and -13 mm (-11%) over the whole accumulation period (September-May) for DD and SD data, respectively. Estimated temperatures from CESM1-WRF and SD data during the accumulation period are about 3˚C and 0.6˚C colder than the observed temperature from stations, respectively. The mean monthly bias between simulated (DD) and observed data was added to the meteorological forcing to correct the variables before using them as input to SnowModel. 93 Figure 4.2. Mean monthly dynamically and statistically downscaled temperature and precipitation against observations for 1979-1994 at station locations (See table 4.1 for station information). The mean monthly offset used for precipitation and temperature adjustment for DD forced data for the three simulation periods from 1979 to 2099. We forced bias corrected data into MicroMet to downscale the resolution of 10 km to 500 m. The evaluation of simulated temperature and precipitation by MicroMet against station data for the historical period is presented in Fig. 4.3. We averaged the mean daily recorded data for the 15-year baseline at all station data and compared with the averaged simulated values using the nearest grid point to a given station (Fig. S3.2). Note that the station data used for SDforced model evaluation are not independent of the data used by PCIC to develop the statistical downscaling as there were no available independent station data for the historical period. There is a strong correlation of 0.93 and 0.98 (p < 0.01) between observations and DD and SD forced mean daily temperature (N = 5400), while the correlation coefficient of mean monthly precipitation (N = 180) are lower values of 0.78 (DD) and 0.89 (SD) (p < 0.01). 94 Figure 4.3. Averaged mean monthly simulated temperature, precipitation, and snow depth against observed data from 10 stations (Table 4.1) for the reference period 1979-1994. The boxes extend from the 25th percentile to the 75th percentile. The whiskers extend to the monthly highest and lowest values for each dataset. The sign (×) shows the medians, and dots represent the outliers. Boxplots of monthly precipitation (Fig. 4.3) reveal a higher median for DD-forced data than observations in cold months. Median monthly SD-forced precipitation values are lower than observed most months. The variability of winter monthly SD precipitation is higher than the 95 variability of DD precipitation. Averaged over all stations for the reference period, there is a strong correlation of 0.83 and 0.73 (p < 0.01) between simulated and observed mean monthly snow depth (N = 180) for DD and SD, respectively. However, snow depth is underestimated in both datasets in some years, especially SD-forced snow depths (Fig. 4.3). Table 4.2 summarizes the performance of SnowModel for snow depth simulations at all ten stations in the Columbia basin for the historical period (1979-1994). The model performance was evaluated using three quantitative statistics including Nash-Sutcliffe efficiency (NSE), Willmott Index of Agreement (WI) and Percent Bias (PBias). Table 4.2. The NSE, WI, percent bias and correlation coefficients r (all significant at the p level of 0.01) of DDforced monthly mean snow depth simulated by SnowModel for the (1979–1994) period for the 10 specified stations. Station Name r* PBias+ NSE‡ WI† Bugaboo Creek Lodge 0.82 16.41 0.62 0.87 Cariboo Lodge 0.75 16.73 0.51 0.77 Duncan Lake Dam 0.98 18.60 0.82 0.99 East Creek 0.61 11.20 0.52 0.83 Fish Lake 0.79 17.10 0.50 0.76 Fred Laing Ridge* 0.65 17.25 0.64 0.86 Glacier Np Mt Fidelity 0.83 19.50 0.73 0.96 Nakusp 0.78 23.61 0.78 0.95 Valemount North 0.67 13.90 0.66 0.88 Yoho NP Wapta Lake 0.63 33.10 0.51 0.78 0.63 0.87 Average 0.75 18.74 * Pearson coefficient +Percent Bias ‡ Nash–Sutcliffe efficiency †Willmott Index of agreement Overall, SnowModel reproduces the seasonal variability of snow depth with the WI of 0.87 and NSE of 0.63, using DD-forced data. SD-forced snow depth shows lower values of the WI (0.78) and the NSE (0.53). Based on the evaluation guidance provided by Moriasi et al. (2007) 96 and Ritter and Muñoz-Carpena (2013), the model simulations can be judged as satisfactory (NSE > 0.50 and PBias < ±25). 4.3. Future projections We examined mid- and end-of-century temperature and precipitation changes calculated by MicroMet at 500 m spatial resolution over the domain for both DD and SD projected datasets. The impact of climate change on the snow depth and SWE was examined by computing the temporally and spatially averaged relative changes of the average SWE of four climate ensemble projections of DD and one projection for SD data. We also calculated the projected solid precipitation difference from the reference period. We then present projected changes in SWE over the Columbia Mountains based on the difference between simulated values of the past period and two projection time periods for mid and end of the century using both forcing data. In addition, interannual variability of snow depth and SWE are presented. Twelve separate simulations were done to cover all ensembles of dynamically downscaled data, in addition to three simulations for statistically downscaled data covering the three periods. All DD-forced values presented and analyzed here represent the average of four ensembles. 4.3.1. Temperature and precipitation We compare seasonally downscaled temperature and precipitation data by MicroMet forced by both dynamically and statistically downscaled data for 15-year periods centred on the 1980s, 2050s, and 2090s (Fig. 4.4). Both datasets show increasing projected temperature in winter. The largest temperature increase occurs in spring and autumn for SD and DD data, 97 respectively. Overall, averaged SD and WRF forced temperature increased by 4˚C and 5˚C from 1979 through 2099 during the accumulation season (September to May), respectively. In contrast, mid-century SD monthly temperatures are colder by ~ 0.2˚C compared to the reference period (Fig. S3.1). DD-forcing data show a higher precipitation rate in the autumn and winter with less temperature changes in spring compared with SD-forcing data (Fig. 4.4). Modelled precipitation rates for both SD and DD forcing data increased in the autumn and spring by 10 and 30%, respectively, while precipitation declined slightly (~5%) throughout the region in winter. Figure 4.4. Downscaled mean daily precipitation rate (left) and temperature (right) by MicroMet over the whole area from September to May for the reference and future periods. The boxes extend from the 25 th percentile to the 75th percentile. The sign (×) represents the medians, and dots show the outliers (< 5% and > 95%). The whiskers extend to the seasonal highest and lowest values for each dataset. 4.3.2. Mean snow depth and SWE The monthly average SWE over the entire domain for the four CESM1-WRF ensembles forced into SnowModel show the variance of the ensemble members for the historical period (Fig. 4.5). The large differences between ensembles arise from internal variability due to different initial and boundary conditions for each ensemble member in downscaling approach 98 within the WRF model. Natural variability (represented by the model-internal variability of a single large-scale model) can play a key role in the variability of future climate forecasts, depending on the component, season, area and time horizon considered (Trentini et al., 2019; Xie et al., 2015, Deser et al., 2012b). Snow water equivalent [m] 0.3 0.25 0.2 0.15 0.1 0.05 0 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 -0.05 Figure 4.5. Monthly averaged SWE [m] over the entire domain from CESM1-WRF forced SnowModel simulations for four ensemble members during historical period (1979-1995). Black line with circle marker shows the average values. We also compare simulated end-of-winter snow accumulation and solid precipitation changes for the 2050s and 2090s (Figs. 4.6 and 4.7). These comparisons are based on the difference between daily averaged projections (2045-2059; 2085-2099) and the reference period (1979-1994). Solid precipitation changes in the 2050s show an absolute decline by 0.11 and 0.19 mm day-1 (~ 7% and 13%) for SD and DD forcing data, respectively over the whole domain. The spatial pattern of solid precipitation (Fig. 4.6a and b) indicates a high rate of loss in mid-elevations across the entire domain using DD data, whereas despite an overall reduction of snow, an increased SD-forced solid precipitation at higher elevations of the northern parts of the Columbia and Rocky Mountains is evident at mid-century. 99 Figure 4.6. The averaged solid precipitation rate change between future simulations and reference period for (a) mid-century and (b) end of the century using statistically downscaled (SD) and dynamically downscaled (DD) forcing. Solid precipitation decreases by an average of 0.41 (mm day -1) (~ 26%) for SD data over the period 2085-2099 with the greatest declines on the eastern side of the Columbia basin, while this difference for DD-forced solid precipitation is -0.36 (mm day-1) (~ 27%) with the highest reduction at mid-elevations (Fig. 4.6). 100 Changes in SWE (Fig.4.7a and b) vary between -75% at lower elevations to +80% at higher elevations by mid-century that indicate an average of 6% and 19% less SWE over the whole area for SD and DD forced data, respectively. Greater decreases for DD-forced snow depth and SWE are consistent with decreased solid precipitation (Fig. 4.6a and b), especially at mid elevations. The snow depth distribution over the Columbia Basin for mid- and end of the century is illustrated in Fig. S3.3. The distribution of changes in SWE by the end of century is in the range between -91% to +64% for SD data, likewise the range of variations for DD-force SWE distribution attains from -79% to +80% over the domain. Over the entire domain, an average of 15% of SD-forced and 30% DD-forced SWE will be lost by end of century. Contrary to the relative decrease, our results show a small absolute decrease in low elevation bands, while the greatest snow reduction occurs at mid elevations, particularly for DD-forcing data (not shown here). 101 Figure 4.7. Averaged snow water equivalent (SWE) changes between future simulations and reference period for (a) mid-century and (b) end of century using statistically downscaled (SD) and dynamically downscaled (DD) forcing data averaged over entire water year. 102 4.3.3. Interannual variability and elevation dependence Global Climate Model projections contain multiple biases such as an interannual variability bias in atmospheric variables. Some studies (e.g. Deser et al., 2012a; Rocheta et al., 2017; Pontoppidan et al., 2019) suggest that there is no clear indication of changes in projected winter temperature and precipitation variability. Nevertheless, other studies indicate a slight decline in temperature variability during winter (Holmes et al., 2016; Rupp et al., 2017; Bathiany et al., 2018) and an increase in precipitation variability in a warmer climate in mid-latitudes (Pendergrass et al., 2017; Rupp et al., 2017; Musselman et al., 2018). We analyzed the difference in snow depth and SWE (average of 15-year simulation) between each projection period and the reference period (Fig. 4.8a) over the whole domain based on the future temperature and precipitation variability. The snow depth maxima at the mid and end of the 21st century is lower than the reference mean snow depth maximum, especially for SD forced snow depth. The pattern of DD-forced snow evolution is the same for all periods, while SDforced snow depth and SWE are deeper during autumn for mid-century and the snow depth peak occurs around two weeks earlier with only one half of the SWE compared with the reference period (Fig. 4.8b). Our analysis shows an increase in SD-forced snow depth and SWE in October and November for mid-century, while an increase in snow depth and SWE occurs in April with DD forcing data. The snow depth and SWE changes in the winter season are less than in other seasons for both datasets. Relative changes in snow depth is time and elevation dependent (Fig. 4.9). 103 Figure 4.8. Mean simulated snow depth and relative difference of snow depth (a) and mean simulated SWE and relative difference of SWE (b) averaged over the entire domain. Analysis were done for the reference period and two projections for mid-century and end of the century forced by statistically downscaled (SD) and dynamically downscaled (DD) forcing data. The sign (_diff) refers to the difference between curves of 2090s and 1980s for both SD and DD dataset. The highest decrease of snow depth happens in lower elevations with a higher percent decrease by end of the 21st century. DD-forced snow depths show a higher percent decrease (~20%) than SD-forced results for both mid and end of the 21st century, in particular, at low and mid elevations (elevations below 1800 m a.s.l.). For low and mid elevations, the percent decrease reaches over 70% and 50% for DD and SD forcing by end of century, respectively. Snow depth above 2000 m a.s.l. changes least in both simulations (Fig. 4.8). 104 Figure 4.9. Relative changes (%) of the snow depth versus elevation by the mid-century (2050s) and the end of century (2090s) averaged over area for the entire water year using two forcing data. The blue represents lower and the red shows the higher frequency of snow depth changes. The black line is the average line. 4.4. Discussion The temperature trend of the region by the end of the 21st century is toward a markedly warmer climate and a higher rate of mean liquid precipitation (Erler et al., 2015; Rupp et al., 2017). It is known that most mesoscale models produce too much precipitation, largely due to frequent, low-magnitude precipitation events (Yhang and Hong, 2008; Gutmann et al., 2011; Yhang et al., 2017). Likewise, the climate projections simulated by GCMs are associated with several uncertainties due to their different formulations and parameterizations, which influence the climate response to the forcing data within the models resulting in different outputs (Stocker et al., 2001), in addition to the uncertainty in future emissions. Some of these uncertainties can be quantified using simulation ensembles (Fig. 4.5). As mentioned before, there is also uncertainty in the RCM downscaling procedure, including WRF (Zubler et al., 2014). On the other hand, a disadvantage of all statistical downscaling approaches is that they depend on the presumption that the historical relationship between predictors and regional 105 climate remains unchanged in future climate conditions (Boé, et al., 2007; Gutmann et al., 2011). Our results show differences between the SnowModel simulated SD and DD forced outcomes, which are due to the difference in forcing data produced by each downscaling method. Our bias correction of modelled temperature and precipitation used the mean monthly offset between the averaged CESM1-WRF modelled outputs of four ensembles and the observed meteorological data. Additionally, we allowed SnowModel to calculate short and long wave fluxes for SD data, while CESM1-WRF downscaled radiation values were used in DD simulations. Our analysis shows that SD incoming short wave and long wave radiation generated by SnowModel are higher than simulated DD-forced radiation by 30 and 100 W m2 , respectively, averaged over the accumulation season (Fig. S3.4). Underestimation of SD- forced snow depth in the reference period can be explained by underestimated SD-forced precipitation, especially at higher elevations, and a warmer temperature during the accumulation season (Fig. 4.3). A higher DD-projected spatio-temporal variation of daily precipitation for each season over the entire domain is because of simulating dynamical changes in the meteorological fields. However statistical downscaling has less spatio-temporal variation because it simply projects changes onto a static dynamical state of the atmosphere, as well as the fact that the SD routines are based on observations that do not span the full range of elevations (Fig. 4.4). The average DD simulated snow depth accords with observed data compared with other studies (Schmucki et al., 2015; Marty et al., 2017), given the fact that the comparisons are based on the nearest grid points to the stations, which do not represent the actual elevation and local topographic effects. Moreover, the measured snow depths at the point scale are for specific locations and elevations that do not generally represent a large area. Our results show a 22% underestimation of DD-forced simulated snow depth, which is in the same range of the results simulated by Schmucki et al. (2015) at the point scale for 11 weather 106 stations (between -15% and 26%) and less than the range of percent differences indicated by Marty et al. (2017) at a regional scale (between -47% and 65%) over the Alps. The snow depth projections using the RCP 8.5 scenario decreases by 38% and 28% over the Columbia River Basin in British Columbia in the 21st century for the averaged four ensemble dynamically downscaled climate projections and for the statistically downscaled climate projections. In general, the inferred magnitude of decrease in snow depth and SWE is plausible in light of strong warming and negligible increases in precipitation. Increased SDforced solid precipitation can be explained by slightly decreased autumn temperature with a higher rate of precipitation (Fig. 4.4). A large temperature increase in spring might cause earlier snow melting in the future compared with the reference period, resulting in a large signal in this season. Pytlak et al. (2018) reported an earlier spring peak runoff for the upper Columbia Basin by 2070s under RCP8.5. Our results are consistent with other studies showing that SWE and snow depth during winter changes least under a warming climate, while snow depth is considerably reduced during autumn and spring (Dye, 2002; Lemke et al., 2007; Ishida et al., 2018). Overall, future snow depth depends primarily on the increasing temperature of the accumulation season as the increase in precipitation is small (Fig. 4.4) (Marty et al., 2017). A colder SD air temperature in autumn and a spring cold bias of DD-forcing data might be responsible for increased snow depth in autumn and spring by mid-century, respectively. A poor representation of snowmelt processes along with excess winter precipitation results in an impact of snow–albedo feedback on temperature within WRF simulations (Fig. 4.8; Erler et al., 2015). The larger loss of snow in DD scenarios may be due to the snow albedo feedback on temperature, which is captured in WRF. The timing of the albedo feedback is not correct in the GCM, thus it does not capture the elevation differences, and melt times are too 107 homogeneous across the domain. In this feedback, WRF warms more in specific regions where snow cover is more prone to receding (Erler et al., 2015). Overall, consistencies can be found in the general decrease of the snow cover depth with solid precipitation patterns for both types of forcing data (Fig. 4.6). Islam et al. (2017) projected a decrease of snowfall (by 40%) in the 2050s with a higher rate of reduction of SWE and snow cover at lower elevations for the Fraser River Basin (FRB) in western Canada. The end of century snow reduction by 60% at low/mid elevations follows other studies over the Alps (69%-75%) (e.g. Rousselot et al., 2012; Marty et al., 2017). Derksen et al. (2019) reported a reduction of 5% to 10% snow cover per decade in seasonal snow accumulation by mid-century for a medium emission scenario across Canada, while an increase of 2% to 5% per decade is reported for British Columbia. We found a small absolute decrease for the end of the century at low elevations, which generally is due to a small amount of snow in this elevation zone. In contrast, the greatest absolute decrease in midelevations occurs because currently snow-covered area during winter months are heavily impacted by higher temperatures in the future. The snow cover at high elevations mainly depends on the precipitation rather than temperature, which cause very small absolute decreases and even some increases in some places using SD-forcing data. As mentioned earlier, the duration and depth of snow is most susceptible to a warming climate during autumn and spring since the snowpack is close to 0˚C. Our findings largely accord with Islam et al. (2017) and Erler et al. (2017), which reported the highest snow cover reduction in the spring over western Canada and USA (e.g. Mote et al., 2018). Our results also show the highest reduction in the autumn season with DD-forcing due to a larger temperature change in the autumn and having a cold bias in the spring. Snow depth and SWE are an integrated response to both temperature and precipitation (Gui et al., 2005; Ding and Qin, 2009; Shi et al., 2011). A warmer 108 climate along with rain-on-snow events are responsible for significant autumn and spring snow changes that elevate the likelihood of rain-on-snow floods (Musselman et al., 2018). Owing that daily mean snow depth is decreasing, the total volume of snow must also diminish (Fyfe et al., 2017; Mote et al., 2018). The large difference between DD and SD projected SWE, is due to differences in forcing data, especially short and long wave radiation (Fig. S3.4). Downscaled radiation from WRF was used in the DD projections, while shortwave and longwave radiations were calculated by SnowModel for SD projections, based on the other meteorological variables. We could not evaluate radiation forcing into the snow model due to lack of observed data. However, comparing simulated radiation over the accumulation period by both forcing datasets reveal higher shortwave and longwave radiation for SD than for DD. 4.4.1. Study limitations We note several limitations in this study: (a) Dynamically downscaled data used in our study represent 15-year epochs and do not reflect century-scale transient simulations of future climate. Various RCM and GCM combinations (Giorgi, 2006; Salathé et al., 2010) and realizations of a single global model (Deser et al., 2012b) can generate different results. In transient simulations, instead of long simulations, the length of the time slices multiple realizations (15-year period × 4 ensembles = 60) have to be performed. At this stage it is not clear if 60 years (four ensembles for 15-year periods, for three periods) of data are adequate to achieve a reliable and stable climatology (Lucas-Picher et al., 2008; Schindler et al., 2015, Erler et al., 2017); (b) In statistically downscaled data, there are important feedbacks that cannot be captured with non-transient simulations. These approaches depend on the presumption that the historical relationship between predictors and regional climate remains unchanged in future climate conditions (Boé et al., 2007; Gutmann et al., 2011); (c) Ground 109 heat flux is not simulated with SnowModel. We expect this to be a minor effect based on the result of Watson et al. (2006) who analyzed the sensitivity of Alpine3D to geothermal heat in the Madison River watershed, USA; (d) MicroMet downscaling methods are simplified in terms of atmospheric conditions such as cloud fraction and temperature inversion layers (Liston and Elder, 2006b). Moreover, effects of wind on snow transport and sublimation of blowing snow is neglected within this study; (e) The SD results used wind components and specific humidity derived from DD data due to lack of these statistically downscaled variables. Therefore, the results might be different if there were available SD data, especially in terms of calculation of surface sublimation and cloudiness within SnowModel; (f) Incoming radiation for SD data were calculated by SnowModel based on the forcing meteorological data. We found a higher rate of incoming shortwave and longwave radiation than the radiation from dynamical downscaling. The lack of observed data prevents us from evaluating simulated radiation; (g) Statistically downscaled data are based on interpolated gridded observations, which may be different from conditions at monitoring stations as biases can exist at high altitudes or other areas with few stations (Eum et al., 2014). The observed data over complex topography in alpine terrain are sparse and limited mostly to the atmospheric conditions and contain uncertainties. Furthermore, another limitation is the use of a single future climate scenario (RCP 8.5) whereas others are also plausible. 4.5. Conclusion We investigated the impacts of climate change (RCP 8.5 scenario) on snow distribution over the Columbia Mountains by a distributed snow model (SnowModel) forced by both regional climate model data (WRF) and statistically downscaled data (PCIC). The model was evaluated by comparing the simulated results with the average of observed data from 10 110 selected stations over the basin. The snow distribution climatology at two future time periods and one historical period was analyzed for CESM1 forcing data both statistically downscaled and dynamically (four different ensembles) downscaled. SWE and snow depth changes are closely related to temperature and precipitation, which differ with latitude and altitude. In addition to the synoptic scale variability, sub-grid processes and the pronounced terrain variation in the Columbia Basin has a strong impact on its climate. At small scales, the changing topographical characteristics such as elevation, orientation to the mean winds, slope, and exposure all combine to change the local climate that affects snow cover (Klock and Mullock, 2002; Boé et al., 2007; Chen et al., 2012). Therefore, application of the highresolution snow model and forcing data (SD/DD dataset) is important to address most of the factors and processes affecting snow cover. The downscaled climate model outputs for the RCP8.5 scenario used in this study demonstrated a clear temperature increase for all time periods using both SD and DD data. There is a slight increase in total precipitation by the end of the century, with the largest changes in spring for SD data. Future seasonal temperature increase during the accumulation season is projected to be highest in spring for SD and in autumn for DD datasets. The snow projections show snowpack reduction through the 21st century with the largest changes for the beginning and end of the accumulation season. The magnitude of the decrease is larger by 10% using DD forcing data. The projected snow depth reduction averaged over the region is about 25% for mid-century and about 40% towards the end of century using DD data. However, under the RCP8.5 scenario the snow depth reduction is projected to reach more than 65% at lower elevations. The mid-elevations are the most vulnerable elevation zones due to higher sensitivity to temperature increase, while higher elevations are less affected by warming. The 111 results indicate that the projected reduction of snow depends largely on elevation and season. Furthermore, other research suggests that there is not a clear preference in using a particular method for downscaling. Both SD and DD have advantages and disadvantages, and within each category there is a wide range of specific methods (Gutmann et al., 2011; Nover et al., 2016). Nevertheless, for both methods our results show a snow depth and SWE reduction in most areas by end of the century that influence water resources for various sectors and communities living in these areas. It is worth noting that our research employed one climate model (CESM1), which is identified as one of the best GCMs representing climate conditions in this region (Mote and Salathé, 2010; Rupp et al., 2017), but is not a comprehensive evaluation of the full range of possibilities. Furthermore, in addition to the uncertainties associated with downscaling approaches in RCMs in terms of configuration, climate change signal and inter-annual variabilities can differ in the various GCM projections. Nonetheless, in such situations, using a multi-model ensemble of regional climate models could be a viable option to simulate precipitation patterns and variability if there is a reliable observed dataset to use for bias correction of RCM outputs, even though they are computationally expensive. Developing statistical downscaling functions for other variables (e.g., wind, pressure, humidity, radiation) to use in land-surface or distributed snow models would lead to a better comparison (Gutmann et al., 2011). 112 Chapter 5: Conclusion 5.1. Dissertation Summary and lessons learned Faced with changing climate conditions, the need for assessment of snowpacks in remote areas provides the motivation for my dissertation research. Western Canada relies on seasonal snow for runoff (Bawden et al., 2015; Bonsal et al., 2019), so quantifying stored freshwater resources is important for adaptation plans. In this dissertation, I investigated glacier winter mass balance and snow distribution using two distributed snow models forced with different input data. The importance of input data and complexity of snow models on glacier mass balance estimation was examined by using high resolution meteorological data from NLDAS and dynamically downscaled data from WRF to force both SnowModel and Alpine3D. I then assessed the effects of wind on snow distribution and winter accumulation by comparing two different simulations of SnowModel, with and without wind transport. I evaluated the results based on in-situ observations and LiDAR measurements. Finally, I simulated future snow depth and SWE under future climate change for much of the Canadian portion of the Columbia Basin using dynamically downscaled CESM1-WRF and statistically downscaled CESM1 from PCIC based on the RCP 8.5 future greenhouse gas concentration scenario. Chapter 2 addressed the first and second research questions that asked how forcing data and snow models complexity affect winter mass balance estimation in the absence of weather station data. We used two snow models (SnowModel and Alpine3D) with different forcing data (WRF and NLDAS) to answer these questions. Both snow models successfully reproduced observed snow depth and SWE when driven with gridded reanalysis and dynamically downscaled data. WRF-forced snow depth and SWE are closer to the observations 113 than simulations forced with NLDAS data. This finding is likely explained by the observed negative precipitation bias observed in NLDAS data over the study area. Our results also revealed that averaged WRF-forced simulated snow depth is overestimated, likely because mesoscale atmospheric models like WRF tend to simulate too much precipitation. In short, these results indicate that accurate downscaling of precipitation data is crucial to correctly estimate snow depth and SWE. In terms of glacier mass balance, SnowModel reproduces winter mass balance with, on average, less than 10% bias, but with a higher variability of snow depth than Alpine3D that had a 22% bias and produced a more uniform snow depth and SWE. When compared to observations, I found that both SnowModel and Alpine3D respectively underestimates snow density by 3% and 24%. This underestimation might arise from underestimation of wind speed (Sturm et al., 2010; Liston and Hiemstra, 2011) by MicroMet in SnowModel and inappropriate parameterization of density calculation in Alpine3D for this region. These results reveal that winter mass balance for glaciers in complex terrain can be adequately simulated with a physically-based model such as SnowModel; higher order-models that require large computational demands may not be required to simulate winter mass balance over large alpine glaciers. The degree to which these models correctly simulate snowmelt over the glaciers remains uncertain, however. In terms of study limitations, Alpine3D is computationally more intensive (~ 13-fold more time) than SnowModel. Therefore, using Alpine3D in event-based simulations might be a better choice than using it in long-term simulations, especially for large areas. In addition, SnowModel and Alpine3D utilize different interpolation methods to prepare data for snow modelling. Employing the same preprocessor for preparing data to force the snow evolution modules makes a better comparison between the models, since the input data have a major 114 influence on simulated snow properties. Since NLDAS data are too coarse to use directly for small glaciers, I recommend making it more site-specific by downscaling in a few steps rather than using its nearest point to the study area. This chapter neglected snow transported by wind, considering wind effects on winter mass balance might lead to a better comparison between the results of the two snow models. Chapter 3 described the effects of wind-topography interactions on snow distribution using WRF and NLDAS meteorological forcing data to answer the third research question (e.g. investigating the effect of wind on glacier winter mass balance in terms of snow erosion, deposition and sublimation). I found that the wind redistribution did not make much difference with the NLDAS data. Comparisons of simulated wind fields to observations showed that MicroMet underestimates wind velocity, this may affect the snow amount results. SnowModel is able to reproduce the erosion on the ridges and deposition on the lee sides of steep terrain. The results indicate negative net transported snow by wind averaged over all glaciers. The model results reveal that strong drifting snow sublimation occurs on the windward slopes likely due to the presence of snow and exposure to the westerly winds. The highest rate of drifting snow sublimation occurs where drifting snow is associated with a higher wind speed, a warmer temperature, and lower humidity conditions (Pomeroy and Jones, 1996; Liston and Sturm, 2004). Although my results show a minor amount of snow loss due to drifting snow sublimation over the accumulation period (2% of SWE averaged over all domains), snow loss can be considerably higher at the smaller scales summed over the accumulation period (reaching up to 500 mm (~100 %) SWE). It is expected that actual sublimation rates would be higher since MicroMet was found to underestimate wind speeds. However, neglecting the 115 thermodynamic feedback within SnowModel could compensate a part of that trend. In terms of wind impacts on winter mass balance, the analysis showed marginal improvement (2%) when averaged for all glaciers with a range of -2.3% to +5.9%. Generally, this study suggests that considering wind effects on the glacier-wide average winter mass balance estimation is not necessary over this study area for large glaciers, given the higher computational demand for simulating such processes in SnowModel. However, blowing snow is a required process to accurately simulate the heterogeneity of winter snow cover in windy environments, especially for small glaciers. I also encountered limitations in this chapter due to the scarcity of suitable wind observations at higher altitudes proximal to each glacier for evaluating simulated wind data. The large number of outliers and significant data gaps made selection of wind data during the winter periods challenging. In addition, MicroMet tends to underestimate the wind speed. Modifying the wind parameterizations in MicroMet may increase the effects of wind redistribution at small scales. Overall, modelling of wind fields with higher accuracy that captures local winds (Vionnet et al., 2020), higher resolution meteorological fields, particularly precipitation (Vionnet et al., 2019), and a higher resolution of simulation that may influence snow transport can lead to a more detailed snow pattern and a higher winter mass balance accuracy. Chapter 4 answered the fourth research question and described projected snow changes in the Columbia Mountains under greenhouse gas emission scenario RCP 8.5. I employed dynamically downscaled CESM1-WRF and statistically downscaled CESM1 provided by PCIC to force SnowModel. The climatology of the snow distribution was 116 compared for three time periods (1979-1994; 2045-2059; 2085-2099). Snow changes differ with both latitude and altitude, which are closely related to temperature and precipitation patterns. Both statistically and dynamically downscaled datasets clearly show increased temperature by the end of the century with a slight precipitation increase that occurs particularly in the spring. The projections show snow depth and SWE decline through the 21st century with the largest changes in the beginning and end of the accumulation season with 10% more reduction using dynamically downscaled data. The highest projected basin average snow reduction is about 25% for the mid-century and slightly less than 40% towards the end of the century. These projections show snow reduction of more than 65% at low elevations (< 2000 m) as these elevation zones are vulnerable to temperature increases and higher rainfall rates. Snow depth at high elevations is projected to remain stable in the winter season. My results also demonstrate the least snow depth changes in winter under a warming climate, while the most snow depth changes happen during autumn and especially spring. While some locations are projected to have increased snow depth through the century. These areas are mostly located at higher altitudes using dynamically downscaled data and at higher latitudes using statistically downscaled data. A bias correction was necessary for the dynamically downscaled data. I used a seasonal average of differences between observed and downscaled data to make the bias correction. Using gridded observation-based dataset such as the climate data for Western North America (ClimateWNA) (Hamann et al., 2013), instead of the 10 stations, may help reduce uncertainty and provide a better representation of interannual variability of snow properties. For SD-forced simulations, there are only two available variables (temperature and precipitation) to force the snow models, while other needed variables are not available. Having statistically downscaled 117 data for wind and humidity may improve snowpack modelling. In addition, I did not evaluate forced and simulated short and long wave radiations due to the lack of observations. Errors or biases in radiation may affect the results, especially at the end of the accumulation season. 5.2. Suggestions for future research I conclude with some recommendations for future work that will help address limitations and improve winter accumulation modelling: 1) Conduct a long-term glacier mass balance simulation along with analysis of each component of the surface energy balance equation to determine the sensitivity of glaciers to climate change, including temperature and cloud coverage for the study area. The climate system is complex, and glaciers respond differently to climatic factors in different locations. To improve our understanding of glacier mass balance and its response to climate change, understanding mass balance sensitivity to subtle changes in the energy balance and accumulation regimes at high altitudes is necessary (Bonekamp et al., 2019). 2) Apply a data assimilation process using sufficient quality and resolution of satellitederived products, that adjusts WRF outputs to better match observations, or use gridded observation-based dataset such as ClimateWNA for bias correction. Uncertainties associated with meteorological forcing are substantial. A cold bias in WRF is often found at high altitudes (e.g., Bonekamp et al., 2018) and would likely impact snowfall and glacier mass balance estimation. In addition, precipitation measurements at high altitudes are sparse and unable to capture spatial variability and the amount of orographic precipitation (Immerzeel et al., 2015). 118 3) Consider the preferential deposition process for winter accumulation simulations. Spatially variable fields of precipitation will increase snow heterogeneity (Lehning et al., 2008). Preferential deposition is an important process in snow patterns (Mott et al., 2014; Vionnet et al., 2018; Gerber et al., 2019) in alpine terrain with steep features where snowfall occurs during strong wind events. This process reduces deposition of snow on windward slopes, while increasing snow deposition on leeward slopes, contributing to the final pattern of snow (Mott et al., 2018). However, the importance of preferential deposition on winter mass balance is unknown and needs to be studied. 4) Take into account snow-avalanching, which may be a significant part of glacier mass input (Dadic et al., 2010), especially for small glaciers. Avalanches redistribute snow to low altitudes, affecting the estimate of mass balance. 5) Use a more comprehensive method to simulate mechanical effects of terrain on windflow. Generating an adequately precise wind field is a key component for accurate blowing snow simulation (Mott et al., 2010; Musselman et al., 2015). Investigating application of methods upon the terrain-curvature wind models for blowing snow (e.g. Liston and Elder, 2006) may improve simulating a detailed snow pattern. 6) Employ a variable horizontal resolution grid (e.g. Marsh et al., 2020) instead of using a fixed resolution spatial discretization method (e.g. Liston and Sturm, 1998). This may lead to a more accurate outcome with less computational demand. To estimate wind redistribution and sublimation of snow at high resolution over large areas, efficient calculation of spatially distributed blowing snow fluxes is required. 119 7) Consider glacier dynamics in the simulations that would improve the results and enable a better comparison with LiDAR measurements of snow patterns. The movement of snow / ice from higher altitudes into the valleys will cause ice and snow to melt faster as the temperature in the low elevations is higher (Johnson and Ohara, 2018) and create higher thinning rates on glaciers (Dehecq et al., 2019). Therefore, a proper coupling between mass balance models and ice flow models (Wijngaard et al., 2019) might make a better estimation of glacier mass balance and final snow pattern. In addition, knowledge of glacier submergence at high elevations and emergence at low elevations would allow a direct comparison with LiDAR measurements at each point on the glacier. 8) Develop statistical downscaling functions for other variables, especially humidity and wind, to make a reliable snow cover modelling in case of using statistically downscaled data. For the use of statistically downscaled data, there are only two available variables (temperature and precipitation) to force the snow models - other variables are needed too. Statistically downscaling wind and humidity would allow us to make snow projections using distributed snow models. 5.3. Research implications The mountain cryosphere is a major source of fresh water in the mountains themselves and in downstream areas. The mountain runoff makes them an important source of freshwater to sustain the ecosystem and supports livelihoods both within and far outside the mountains (Hock et al., 2019). The Canadian portion of Columbia River Basin covers 15% of the total watershed but provides about 40% of the total annual streamflow. The presence of snow, glaciers and permafrost significantly influences the quantity, timing and biogeochemical 120 characteristics of runoff (Kohler et al., 2009). The Columbia River originates in British Columbia and flows into the Pacific Ocean in Oregon, USA. The river is 2000 km long and has over 40 major hydropower stations that generate most of the Pacific Northwest’s electricity (U.S. DIBR, 2016). Adequate and secure water supply is crucial for human health, economic strength and environmental and ecological protection in the Western North America while global climate change poses a major challenge for the sustainability of these resources. In brief, it is evident that decreases in glacier and snow cover in many of the snow-dominated or glacier-fed river basins have influenced and will continue to alter river runoff in amount and seasonality. The average winter runoff is projected to increase, and the maximum spring peak will occur earlier, with resulting decreases in runoff in Western North America no later than the end of the 21st century (Lee et al., 2016; Madani and Lund, 2010). The expected increases in runoff in the near future will affect downstream water management, related hazards and ecosystems (Hock et al., 2019). Increased runoff and sediment transport would increase the risk of overflow (non-productive discharge), especially during winter and spring melt, with the highest impact on run-of-river power plants (Minville et al., 2010; Warren and Lemmen, 2014). Snowpack declines also affect species, particularly those that depend on deep snowpacks. For instance, in British Columbia, the range of mountain caribou has decreased by 40% from historical distribution (Spalding, 2000). Mountain caribou have adapted behaviourally to deep snow conditions by seasonal altitudinal migration. In late winter, deep snowpacks provide a platform for caribou to access food (mainly arboreal lichen) at high elevation ranges. Decreasing snow cover in this region will intensify today's problems for resource managers in British Columbia (integrating the needs for caribou and forestry due to timber harvesting at lower elevations) (Wittmer, 2004). There are many other expected impacts on different sectors 121 such as fishing, tourism, agriculture, ecosystem management, hydropower generation, etc. Therefore, I believe that human society needs to prepare adaptation and mitigation steps for the full range of impacts of snow cover declines and glacier shrinkage in the affected regions. 122 Appendix 1. Supplementary information of Chapter 2 Figure S1.1. Scatter plots of daily mean temperature derived from WRF and NLDAS against observations for 2013-2014 (top) and 2015-2016 (bottom) at four station locations. 123 Figure S1.2. Topography (km) and outline of the outer at 7 km and inner WRF domains at 1 km. 124 Figure S1.3. Simulated time series of daily SWE for the calibration period winter 2013-2014 (top) and snow depth for evaluation period for winter 2015-2016 (bottom) against observed data for four stations. 125 Figure S1.4. Comparison of simulated SWE [mm w.e.] by SnowModel averaged over the entire simulation domain at East-Creek station for winter 2013-2014. Light-blue line represents SWE forced by hourly data and dashed-black line shows simulated SWE forced by 6-hourly data. Figure S1.5. Comparison of simulated SWE (Left) and snow depth (Right) using a rain/snow threshold default value of 1.2˚C (Light blue) and user defined value of 2˚C (Dashed black) by Alpine3D at a point scale for accumulation period winter 2015-2016. 126 Figure S1.6. Spatial scale comparison of simulated SWE [m w.e.] by Alpine3D, averaged over the winter of 20152016 at Kokanee Glacier using different thresholds of 1.2˚C (A) and 2˚C (B). (Same as Fig. S1.4, but at spatial scale) Figure S1.7. Spatial comparison of simulated snow density by SnowModel at Kokanee glacier using multilayer (A) and one-layer (B) methods and their absolute difference (multilayer – 1-layer) (C) for winter accumulation 2015-2016. 127 Figure S1.8. Spatial comparison of simulated SWE [mm w. e.] by SnowModel with multilayer (A) and one-layer (B) averaged over accumulation period 2015-2016. Difference between both methods (A - B) is illustrated in (C). Figure S1.9. Time series of daily snow density [kg m-3] simulated by SnowModel (SM) and Alpine3D (ALP) for winter 2015-2016. 128 Table S1.1. The main differences between snow evolution modules of two snow models, Alpine3D and SnowModel SNOWPACK SnowPack (Alpine3D’s snow evolution core) (SnowModel’s snow evolution core) Soil/Snow/Canopy column – Neglects lateral transfers Snow/Canopy column – Neglects lateral transfers The snow is modelled as three phases (ice/liquid water/water vapour). It can also simulate layers such as ice lenses permafrost, ponding etc. SnowPack is a simple one-layer model; however, it is possible for users to define multiple layers for snow evolution. Number of layers is arbitrary and can simulate very thin layers like ice crust and hoar. State variables include temperature, liquid water content and density that must be known for each layer. From mentioned variables, four microstructure parameters are derived: grain size, bond radius, sphericity and dendricity. SnowPack defines changes in the snowpack in response to the melt fluxes and precipitation input given by MicroMet. Equilibrium growth metamorphism and kinetic metamorphism will drive the temporal evolution of these parameters according to various models. Thermal conductivity and viscosity are computed by microstructure parameters. Furthermore, it includes pressure sintering strain amplification on the necks and its feedback on metamorphism as well as linear and non-linear viscosity ranges. Compaction-based snow density evolutionclosely follows that of Anderson (1976) where density changes with time in response to snow temperature and weight of overlying snow. snow melting alters snow density by decreasing snow depth and redistributing meltwater through the snowpack until a maximum snow density is reached These parameters allow computing mass and energy conservation and settling. Energy balance includes solar radiation, sublimation/deposition of water vapour, melting refreezing and heat conduction. It takes into account refreezing; however, it does not consider melting from internal ice deformation, changes in drainage system, geothermal effects, and sub-glacial friction melting. Non-blowing snow sublimation is calculated in EnBal and used to adjust snowpack depth The heat transfer, vapour diffusion, water transport, and phase changing processes was formulated SNOWPACK has also a detailed description of the interaction with the atmospheric boundary layer is a one-way model without considering from the surface back to the atmosphere 129 Table S1.2. The percent bias and NSE of daily SWE for SnowModel performance using different lapse rates. User lapse rate derived from PRISM data and default lapse rates within the model. Default P&T_LR* Function NSE % Bias East Creek 0.95 11.2 Glacier-NMF 0.96 8.2 Redfish-Creek 0.30 33.6 Azure-Rive 0.93 13.4 User T_LR Default P_LR NSE % Bias 0.96 10.3 0.97 8.6 0.04 40.5 0.93 12.0 User P_LR Default T_LR NSE % Bias 0.92 14.6 0.93 14.4 0.35 32.0 0.93 12.6 User T&P_LR NSE % Bias 0.93 14.4 0.93 14.1 0.10 38.8 0.95 11.0 *T_LR: Temperature lapse-rate, P_LR: Precipitation lapse-rate 130 Table S1.3. Statistics of simulated hourly SWE with respect to the observations at four station locations for the accumulation period from 1 October to 30 April 2015-2016. SWE [m.w.e] Winter (2015-2016) Station ID East-creek Azure-river Redfish-creek Glacier-NMF 1 Function PCC D NSE APB Bias MAE RMSE PCC D NSE APB Bias MAE RMSE PCC D NSE APB Bias MAE RMSE PCC D NSE APB Bias MAE RMSE WRFOUT 0.99 0.99 0.97 6.24 -0.002 0.04 0.05 0.97 0.99 0.94 8.84 -0.01 0.06 0.09 0.97 0.97 0.94 8.84 -0.01 0.06 0.09 0.99 0.95 0.76 24.0 0.11 0.12 0.15 SM1_WRF3 0.98 0.98 0.93 11.57 -0.06 0.07 0.08 0.96 0.97 0.86 12.02 0.05 0.08 0.13 0.96 0.97 0.86 12.02 0.05 0.08 0.13 0.99 0.98 0.93 15.6 0.10 0.10 0.11 SM_NLDAS4 0.98 0.98 0.92 8.6 0.05 0.05 0.09 0.90 0.95 0.8 12.2 0.04 0.08 0.16 0.90 0.95 0.79 12.20 0.04 0.08 0.16 0.99 0.99 0.98 7.20 -0.02 0.05 0.06 ALP2_WRF 0.98 0.99 0.95 6.57 0.02 0.04 0.07 0.98 0.98 0.92 13.47 -0.06 0.09 0.10 0.98 0.98 0.92 13.47 -0.06 0.09 0.10 0.99 0.98 0.93 14.4 -0.07 0.09 0.11 ALP_NLDAS 0.98 0.98 0.94 11.61 -0.05 0.07 0.08 0.83 0.80 0.30 33.7 -0.22 0.22 0.30 0.83 0.80 0.30 33.7 -0.22 0.22 0.30 0.99 0.97 0.92 14.7 -0.07 0.09 0.12 SM: SnowModel; 2ALP: Alpine3D; 3_WRF: WRF forcing; 4_NLDAS: NLDAS forcing Table S1.4. Comparison of simulated density and SWE by SnowModel, averaged over Kokanee Glacier during winter accumulation 2015-2016 using one-layer and multilayer (6 layers) methods. Multilayer One-layer Diff Density [kg m-3] ON OFF SWE [m w.e.] ON OFF 434.7 426.4 8.3 1.4 1.4 0 457.8 445.1 12.6 1.3 1.3 0 131 2. Supplementary information of Chapter 3 Figure S2.1. Flowchart of model simulations with considering wind effects and without wind. 132 Table S2.1. Statistics of daily SWE (m.w.e) forced by WRF, with and without wind, against snow pillow observations for the accumulation period from 1 October to 30 April 2013-2014. SWE [m.w. e] WRF_forced with wind WRF_forced without wind Function/Station ID Correlation Index of agreement Nash-Sutcliffe Absolute percent bias Bias Mean absolute error Root mean square error Correlation Index of agreement Nash-Sutcliffe Absolute percent bias Bias Mean absolute error Root mean square error East Creek 2D08P 0.98 0.99 0.96 10.63 -0.02 0.06 0.07 0.98 0.99 0.96 10.67 -0.02 0.06 0.07 Glacier NP MT 0.99 0.99 0.98 7.06 0.00 0.04 0.05 0.99 0.99 0.98 7.73 -0.01 0.05 0.06 Redfish Creek 2D14P 0.90 0.94 0.78 18.66 0.09 0.15 0.22 0.89 0.92 0.76 24.48 -0.03 0.19 0.23 Azure River 1E08P 0.96 0.98 0.92 14.66 0.00 0.07 0.09 0.96 0.98 0.92 14.26 0.00 0.07 0.09 133 3. Supplementary information of Chapter 4 Figure S3.1. Daily mean temperature of statistically downscaled (SD) and dynamical downscaled (DD) data from GESM1 for three 15-year periods (historical, mid-century and end of century) averaged over the Columbia Basin based on RCP 8.5. at resolution of 500 m. Figure S3.2. Averaged snow depth time series over 10 stations compared with the DD-forced snow depth (mean of 4 ensembles) for the historical period 1979-1995. 134 Figure S3.3. Projected snow depth change (%) averaged over the year using statistically downscaled (SD) and dynamical downscaled (DD) data from CESM1 for mid-century (a) and end of century (b) over the Columbia Basin. 135 Figure S3.4. Simulated incoming short wave (SW) (Right) and long wave (LW) (Left) radiations by SnowModel forced by SD and DD datasets for period 1979-1994. The boxes extend from the 25 th percentile to the 75 th percentile. The whiskers extend to the monthly highest and lowest values for each dataset. The sign (×) shows the medians, and dots represent the outliers. 136 References Adam, J.C., Hamlet, A.F., Lettenmaier, D.P., 2009: Implications of global climate change for snowmelt hydrology in the twenty-first century, Hydrological Processes, 23, 962–972 Andreassen, L.M., Van den Broeke, M.R., Giesen, R.H. and Oerlemans, J., 2008: A 5-year record of surface and mass balance from the ablation zone of Storbreen, Norway, Journal of Glaciology, 54(185), 245–258, doi: 10.3189/002214308784886199 Arnold, N.S., Willis, I.C., Sharp, M.J., Richards K.S., and Lawson, W.J., 1996: A distributed surface energybalance model for a small valley glacier. I. Development and testing for Haut Glacier d’Arolla, Valais, Switzerland. Journal of Glaciology, 42(140), 77-89. Arnold, N.S., Gareth Rees, W., Hodson, A.J., and Kohler, J., 2006: Topographic controls on the surface energy balance of a high Arctic valley glacier, Journal of Geophysical Research, 111, F02011 Barnes, S. L., 1964: A technique for maximizing details in numerical weather map analysis, Journal of Applied Meteorology, 3, 396-409 Barnes, S. L., 1973: Regional objective analysis using weighted timeseries observations. NOAA Tech. Memo. ERL NSSL-62, National Severe Storms Laboratory, Norman, OK, 60 pp Barnett, T.P., Adam, J.C., Lettenmaier, D.P., 2005: Potential impacts of a warming climate on water availability in snow-dominated regions, Nature, 438, 303–309. Bartelt, P., and Lehning, M., 2002: A physical SNOWPACK model for the Swiss avalanche warning, Part I: Numerical model, Cold Region Science and Technology, 35,123-145. Bauder, A., Funk, M., and Huss, M., 2007: Ice volume changes of selected glaciers in the Swiss Alps since the end of the 19th century, Annal of Glaciology, 46, 145–149 Bavay, M., and Egger, T., 2014: MeteoIO 2.4.2: a preprocessing library for meteorological data. Geoscience Model Development, 7, 3135-3151 Bavay, M., Lehning, M., Jonas, T., Löwe, H., 2009: Simulations of future snow cover and discharge in Alpine headwater catchments, Hydrological Processes, 23, 95-108 Bavera, D., Bavay, M., Jonas, T., Lehning, M., De Michele, C., 2014: A comparison between two statistical and a physically-based model in snow water equivalent mapping, Advances in Water Resources, 63, 167– 178. Bawden, A., Burn, D., and Prowse, T., 2015: Recent changes in patterns of western Canadian river flow and association with climatic drivers, Hydrology Research, 46(4), 551-565, doi:10.2166/nh.2014.032. Beedle, M.J., Menounos, B., and Wheate, R., 2015: Glacier change in the Cariboo Mountains, British Columbia, Canada (1952–2005), The Cryosphere, 9, 65-80 Bernhardt, M., Zangl, G., Liston, G.E., Strasser, U., and Mauser, W., 2009: Using wind fields from a high‐ resolution atmospheric model for simulating snow dynamics in mountainous terrain, Hydrological Processes, 23(7), 1064–1075, doi:10.1002/hyp.7208. Bernhardt, M., Schulz, K., Liston, G.E., and Zängl, G., 2012: The influence of lateral snow redistribution processes on snow melt and sublimation in alpine regions, Journal of Hydrology, 424-425, 196-206 137 Berthier, E., Vincent, C., Magnusson, E., Gunnlaugsson, A., Pitte, P., Le Meur, E., Masiokas, M., Ruiz, L., Pálsson, F., Belart, J. M. C., and Wagnon, P., 2014: Glacier topography and elevation changes from Pléiades very high resolution stereo images, The Cryosphere, 8, 4849–4883, doi:10.5194/tcd-8-48492014. Blöschl, G., Kirnbauer, R., & Gutknecht, D., 1991: A spatially distributed snowmelt model for application in alpine terrain. Snow, Hydrology and Forests in high Alpine Areas (Proceedings of the Vienna Symposium, August 1991), IAHS Publ, no. 205 Blӧschl, G., Kirnbauer, R., Gutknecht, D., 1991: Distributed snowmelt simulations in an Alpine catchment. 1. Model evaluation on the basis of snow cover patterns, Water Resources Research, 27, 3171-3179. Boé, J., Terray, L., Habets, L., and Martin, E., 2007: Statistical and dynamical downscaling of the Seine basin climate for hydro-meteorological studies, International Journal of Climatology, 27(12), 1643-1655, doi: 10.1002/joc.1602 Bolch, T., Menounos, B., and Wheate, R., 2010: Landsat-based inventory of glaciers in western Canada, 1985– 2005, Remote Sensing of Environment, 114, 127-137 Bonekamp, P.N.J., Collier, E., and Immerzeel, W.W., 2018: The impact of spatial resolution, landuse and spinup time on resolving spatial precipitation patterns in the Himalayas, Journal of Hydrometeorology, 19: 1565–1591, doi: 10.1175/JHM-D-17-0212.1 Bonekamp, P.N.J., de Kok, R.J., Collier, E., and Immerzeel, W.W., 2019: Contrasting meteorological drivers of the glacier mass balance between the Karakoram and central Himalaya, Frontier in Earth Science, 7, 107, doi: 10.3389/feart.2019.00107 Bonsal, B.R., Peters, D.L., Seglenieks, F., Rivera, A., and Berg, A., 2019: Changes in freshwater availability across Canada; Chapter 6 in Canada’s Changing Climate Report, (ed.) E. Bush and D.S. Lemmen; Government of Canada, Ottawa, Ontario, 261–342. Brahney, J., Weber, F., Foord, V., Janmaat, J., Curtis, P.J., 2017: Evidence for a climate-driven hydrologic regime shift in the Canadian Columbia Basin. Canadian Water Resources Journal, 42 (2), 179-192, doi:10.1080/07011784.2016.1268933. Brock, B.W., Willis, IC., and Sharp, M.J., 2000: Measurement and parameterization of albedo variations at Haut Glacier d’Arolla, Switzerland, Journal of Glaciology, 46(155), 675–688 Brun, E., Martin, E., Simon, V., Gendre C., Coleou, C., 1989: An energy and mass model of snow cover suitable for operational avalanche forecasting, Journal of Glaciology, 35, 333-342. Brun, E., David, P., Sudul M., Brunot G., 1992: A numerical model to simulate snow-cover stratigraphy for operational avalanche forecasting, Journal of Glaciology, 38, 13-22. Burn, D.H., 2008: Climatic influences on streamflow timing in the headwaters of the Mackenzie River Basin. Journal of Hydrology, 352, 225-238 Chen, J., Brissette, F.P., and Leconte, R., 2012: Coupling statistical and dynamical methods for spatial downscaling of precipitation, Climatic Change, 114(3-4), 509-526, doi: 10.1007/s10584-012-0452-2 Chen, J., Chen, J., Liao, A., Cao, X., Chen, L., Chen, X., He, C., Han, G., Peng, S., Lu, M., Zhang, W., Tong, X., Mills. J., 2015: Global land cover mapping at 30 m resolution: A POK-based operational approach, Journal of Photogrammetry and Remote Sensing, 103: 7-27 138 Choi, G., Robinson, D.A., and Kang, S., 2010: Changing Northern Hemisphere snow seasons, Journal of Climate, 23, 5305–5310, doi:10.1175/2010JCLI3644.1 Christensen, J.H., Carter, T.R., Rummukainen, M., Amanatidis, G., 2007: Evaluating the performance and utility of regional climate models: the PRUDENCE project, Climatic Change, 81, 1-6 Chu, Q., Xu, Z., Chen, Y., Han, D., 2017: Evaluation of the WRF model with different domain configurations and spin-up time in reproducing a sub-daily extreme rainfall event in Beijing, China. Hydrology and Earth System Sciences Discussions, doi:10.5194/hess-2017-363 Cogley, J., Hock, R., Rasmussen, L., Arendt, A., Bauder, A., Braithwaite, R., Jansson, P., Kaser, G., Möller, M., and Nicholson, L., 2011: Glossary of glacier mass balance and related terms, IHP-VII technical documents in hydrology No. 86, IACS Contribution No. 2, Int. Hydrol. Program UNESCO Paris. Collier, E., Mӧlg, T., Maussion, F., Scherer, D., Mayer, C., and Bush, A.B.G., 2013: High-resolution interactive modelling of the mountain glacier–atmosphere interface: an application over the Karakoram, The Cryosphere, 7,779-795 Comin, A.N., Schumacher, V., Justino, F., Fernández, A., 2018: Impact of different microphysical parameterizations on extreme snowfall events in the Southern Andes, Weather and Climate Extremes, 21, 65-75 Cosgrove, B.A., Lohmann, D., Mitchell, K.E., Houser, P.R., Wood, E.F., Schaake, J.C., Robock, A., Marshall, C., Sheffield, J., Duan, Q., Luo, L., Higgins, R.W., Pinker, R.T., Tarpley, J.D., Meng, J., 2003: Realtime and retrospective forcing in the North American Land Data Assimilation System (NLDAS) project, Journal of Geophysical Research, 108(D22), 8842 Dadic, R., Mott, R., Lehning, M., and Burlando, P., 2010: Wind influence on snow depth distribution and accumulation over glaciers. Journal of Geophysical Research, 115(F1), F01012, doi: 10.1029/2009JF001261 Dai, A., 2008: Temperature and pressure dependence of the rain-snow phase transition over land and ocean. Geophysical Research Letters, 35, L12802, doi:10.1029/2008gl033295. Daly, C. R., Neilson, R. P., and Phillips, D. L., 1994: A statistical topographic model for mapping climatological precipitation over mountainous terrain. Journal of Applied Meteorology, 33, 140–158, doi: 10.1657/1938-4246-42.1.76 DeBeer, C.M., and Pomeroy, J.W., 2009: Modelling snow melt and snowcover depletion in a small alpine cirque, Canadian Rocky Mountains, Hydrological Processes, 23, 2584 – 2599 DeBeer, C.M., Pomeroy, J.W., 2017: Influence of snowpack and melt energy heterogeneity on snow cover depletion and snowmelt runoff simulation in a cold mountain environment, Journal of Hydrology, 553, 199–213 Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Vitart, F., 2011: The ERAInterim reanalysis: configuration and performance of the data assimilation system, Quarterly Journal of the Royal Meteorological Society, 137(656), 553-597, doi:10.1002/qj.828 Dehecq, A., Gourmelen, N., Gardner, A. S., Brun, F., Goldberg, D., Nienow, P. W., Berthier, E., Vincent, C., Wagnon, P., and Trouvé, E., 2019: Twenty-first century glacier slowdown driven by mass loss in high mountain Asia, Nature Geoscience, 12, 22–27, doi: 10.1038/s41561-018-0271-9 Derksen, C., Burgess, D., Duguay, C., Howell, S., Mudryk, L., Smith, S., Thackeray, C. and Kirchmeier-Young, M., 2019: Canada’s changing climate report: Changes in snow, ice, and permafrost across Canada; 139 Chapter 5 in Canada’s Changing Climate Report, (ed.) E. Bush and D.S. Lemmen; Government of Canada, Ottawa, Ontario, 194–260. Déry, S.J., and Taylor, P.A., 1996: Some aspects of the interaction of blowing snow with the atmospheric boundary layer, Hydrological Processes, 10(10), 1345-1358 Déry, S.J., and Yau, M.K., 1999: A bulk blowing snow model. Boundary Layer Meteorology, 93, 237-251 Déry, S.J., and Yau, M.K., 2001: Simulation of blowing snow in the Canadian Arctic using a double-moment model, Boundary-Layer Meteorology, 99, 297–316 Déry, S.J., Clifton, A., MacLeod, S., and Beedle, M.J., 2010: Blowing snow fluxes in the Cariboo Mountains of British Columbia, Canada. Arctic, Antarctic, and Alpine Research, 42(2), 188-197 Deser, C., Knutti, R., Solomon, S., and Phillips, A. S., 2012a: Communication of the role of natural variability in future North American climate. Supplement, Nature Climate Change, 2, 775–779, doi:10.1038/nclimate1562 Deser, C., Phillips, A., Bourdette, V., Teng, H., 2012b: Uncertainty in climate change projections: the role of internal variability, Climate Dynamics, 38:527–546, doi:10.1007/s00382-010-0977-x Ding, Y. J., and Qin, D.H., 2009: Cryosphere change and global warming: impact and challenges in China. China Basic Science, 34–10 Doorschot, J., Raderschall, N., Lehning, M., 2001: Measurements and one-dimensional model calculations of snow transport over a mountain ridge, Annals of Glaciology, 32(1), 153-158 Dye, D.G., 2002: Variability and trends in the annual snow-cover cycle in Northern Hemisphere land areas. 19722000: Hydrological Processes, 16, 3065-3077 Erler, A.R., Peltier, W.R., and D’orgeville, M., 2015: Dynamically downscaled high-resolution hydroclimate projections for Western Canada, Journal of Climate, 28(2), 423-450, doi: 10.1175/JCLI-D-14-00174.1 Erler, A.R., Peltier, W.R., 2017: Projected hydroclimatic changes in two major river basins at the Canadian west coast based on high-resolution regional climate simulations, Journal of Climate, 30, 8081–8105, doi: 10.1175/JCLI-D-16-0870.1 Etchevers, P., Martin, E., Brown, R., Fierz, C., Lejeune, Y., Bazile, E., Boone, A., Dai, Y.-J., Essery, R., Fernandez, A., Gusev, Y., Jordan, R., Koren, V., Kowalczyk, E., Nasonova, N.O., Pyles, R.D., Schlosser, A., Shmakin, A.B., Smirnova, T.G., Strasser, U., Verseghy, D., Yamazaki, T., and Yang, Z.-L., 2002: SnowMIP, an intercomparison of snow models: first results, Paper presented at the International Snow Science Workshop, Penticton, British Columbia. 29 September - 4 October 2002. Eum, H.I., Dibike, Y., Prowse, T., and Bonsal, B., 2014: Inter-comparison of high-resolution gridded climate data sets and their implication on hydrological model simulation over the Athabasca Watershed, Canada, Hydrological Processes, 28, 4250–4271, doi: 10.1002/hyp.10236 Finger, D., Heinrich, G., Gobiet, A., Bauder, A., 2012: Projections of future water resources and their uncertainty in a glacierized catchment in the Swiss Alps and the subsequent effects on hydropower production during the 21st century, Water Resources Research, 48, W02521, doi:10.1029/2011WR010733. Franz, K.J, Hogue, T.S., Sorooshian, S., 2008: Operational snow modelling: Addressing the challenges of an energy balance model for National Weather Service forecasts, Journal of Hydrology, 360, 48-66 140 Fyfe, J. C., Derksen, C., Mudryk, L., Flato, G.M., Santer, B.D., Swart, N.C., Molotch, N.P., Zhang, X., Wan, H., Arora, V.K., Scinocca, J., and Jiao, Y., 2017: Large near-term projected snowpack loss over the western United States, Nature Communications, 8, 14996, doi: 10.1038/ncomms14996 Gallée, H., Trouvilliez, A., Agosta, C., Genthon, C., Favier, V., Naaim-Bouvet, F., 2012: Transport of snow by the wind: A comparison between observations in Adélie Land, Antarctica, and simulations made with the regional climate model MAR, Boundary Layer Meteorology, 146, 133-147 Gascoin, S., Lhermitte, S., Kinnard, C., Borstel, K., Liston, G.E., 2013: Wind effects on snow cover in PascuaLama, Dry Andes of Chile, Advances in Water Resources, 55, 25-39 Gao, L., Bernhardt. M., Schulz, K., 2012: Elevation correction of ERA‐Interim temperature data in complex terrain, Hydrology and Earth System Sciences, 16(12), 4661–4673, doi: 10.5194/hess-16-4661-2012. Gao, L., Schulz, K., Bernhardt, M., 2014b: Statistical downscaling of ERA‐Interim forecast precipitation data in complex terrain using LASSO algorithm, Advances in Meteorology, 472741, 1–16, doi: 10.1155/2014/472741. Gayton., D., 2008: Impacts of climate change on British Columbia’s Diversity: A literature review. Forrex Forest Research Extension Partnership, Kamloops, BC, Forrex Series 23, url: http://www.forrex.org/publications/forrexseries/fs23.pdf Gerbaux, M., Genthon, C., Etchevers, P., Vincent, C., and Dedieu, J.P., 2005: Surface mass balance of glaciers in the French Alps: distributed modelling and sensitivity to climate change, Journal of Glaciology. 51(175), 561–572 Gerber, F., Mott, R., Lehning, M., 2019: The importance of near-surface winter precipitation processes in complex alpine terrain, Journal of Hydrometeorology, 20, 177-196. Gómez-Navarro, J.J., Raible, C.C., Dierer, S., 2015: Sensitivity of the WRF model to PBL parameterizations and nesting techniques: evaluation of wind storms over complex terrain, Geoscientific Model Development, 8, 3349-3363, doi: 10.5194/gmd-8-3349-2015 Gordon, M., Simon, K., and Taylor, P.A., 2006: On snow depth predictions with the Canadian Land Surface Scheme including a parametrization of blowing snow sublimation, Atmosphere-Ocean, 44(3), 239-255 Grell, G. A., and Dévényi, D., 2002: A generalized approach to parameterizing convection combining ensemble and data assimilation techniques, Geophysical Research Letters, 29(14), 1693. Gu, H., Jin, J., Wu, Y., Ek, M. B., and Subin, Z. M., 2013: Calibration and validation of lake surface temperature simulations with the coupled WRF-lake model, Climatic Change, 129(3-4), 471-483. doi:10.1007/s10584-013-0978-y Gudmundsson, L., Bremnes, J. B., Haugen, J. E., and Engen-Skaugen, T., 2012: Technical Note: Downscaling RCM precipitation to the station scale using statistical transformations – a comparison of methods, Hydrology and Earth System Sciences, 16(9), 3383–3390, doi: 10.5194/hess-16-3383-2012 Gui, C., Yang, Q., Wang, S., 2005: Comparison analysis of the long-term variations of snow cover between mountain and plain areas in Xinjiang region from 1960 to 2003, Journal of Glaciology and Geocryology, 27, 486–90 Gurgiser, W., Marzeion, B., Nicholson, L., Ortner, M., Kaser, G., 2013: Modelling energy and mass balance of Shallap Glacier, Peru, The Cryosphere, 7, 1787-1802 141 Gutmann, E.D., Larson, K.M., Williams, M.W., Nievinski, F.G., and Zavorotny, V., 2011: Snow measurement by GPS interferometric reflectometry: an evaluation at Niwot Ridge, Colorado, Hydrological Processes, 26(19), 2951-2961, doi: 10.1002/hyp.8329 Hamann, A., Wang, T., Spittlehouse, D.L., and Murdock, T.Q., 2013: A comprehensive, high-resolution database of historical and projected climate surfaces for western North America, Bulletin of the American Meteorological Society, 94, 1307–1309 Harris, I., Jones, P.D., Osborn, T.J., and Lister, D.H., 2014: Updated high-resolution grids of monthly climatic observations—The CRU TS3.10 dataset, International Journal of Climatology, 34, 623–642, doi:10.1002/joc.3711. Hedrick, A., Marshall, H.-P., Winstral, A., Elder, k., Yueh, S., and Cline, D., 2015: Independent evaluation of the SNODAS snow depth product using regional-scale lidar-derived measurements, The Cryosphere, 9, 13– 23. Holland, S.S. 1976: Landforms of British Columbia, A Physiographic Outline, Bulletin 48, 64-80. Hock, R., 2003: Temperature index melt modelling in mountain areas, Journal of Hydrology, 282, 104–115. Hock, R., 1999: A distributed temperature-index ice- and snowmelt model including potential direct solar radiation, Journal of Glaciology, 45(149), 101-111 Hock, R., and Holmgren, B., 2005: A distributed surface energy-balance model for complex topography and its application to Storglaciären, Sweden, Journal of Glaciology, 51(172), 25-36 Hock, R. and Jansson, P., 2005: Modelling glacier hydrology, In Anderson MG and McDonnell JJ eds, Encyclopedia of hydrological sciences, 4, 2647–2655 Hock, R., G. Rasul, C. Adler, B. Cáceres, S. Gruber, Y. Hirabayashi, M. Jackson, A. Kääb, S. Kang, S. Kutuzov, A. Milner, U. Molau, S. Morin, B. Orlove, and H. Steltzer, 2019: High Mountain Areas. In: IPCC, Special report on the ocean and cryosphere in a changing climate [H.-O. Pörtner, D.C. Roberts, V. Masson-Delmotte, P. Zhai, M. Tignor, E. Poloczanska, K. Mintenbeck, A. Alegría, M. Nicolai, A. Okem, J. Petzold, B. Rama, N.M. Weyer (eds.)] Hofer, M., Marzeion, B., Mölg, T., 2015: A statistical downscaling method for daily air temperature in datasparse, glaciated mountain environments, Geoscientific Model Development, 8, 579-593 Hood, E., Williams, M., and Cline, D., 1999: Sublimation from a seasonal snowpack at a continental, mid-latitude alpine site, Hydrological Processes, 13, 1781-1797 Houze, R. A. Jr., 2012: Orographic effects on precipitating clouds. Reviews of Geophysics, 50. RG1001, doi:10.1029/2011RG000365 Huss, M., Bauder, A., and Funk, M., 2009: Homogenization of long-term mass-balance time series. Annals of Glaciology, 50, 198–206, doi: 10.3189/172756409787769627 Iacono, M. J., Delamere, J. S., Mlawer, E.J., Shephard, M.W., Clough, S.A., and Collins, W.D., 2008: Radiative forcing by long‐lived greenhouse gases: Calculations with the AER radiative transfer models, Journal of Geophysical Research, 113, D13103. doi:10.1029/2008JD009944. Immerzeel, W.W., Wanders, N., Lutz, A.F., Shea, J.M., and Bierkens, M.F.P., 2015: Reconciling high-altitude precipitation in the upper Indus Basin with glacier mass balances and runoff, Hydrology and Earth System Sciences, 19, 4673–4687, doi: 10.5194/hess-19-4673-2015 142 IPCC, 2019: Climate Change and Land: an IPCC special report on climate change, desertification, land degradation, sustainable land management, food security, and greenhouse gas fluxes in terrestrial ecosystems [P.R. Shukla, J. Skea, E. Calvo Buendia, V. Masson-Delmotte, H.-O. Pörtner, D. C. Roberts, P. Zhai, R. Slade, S. Connors, R. van Diemen, M. Ferrat, E. Haughey, S. Luz, S. Neogi, M. Pathak, J. Petzold, J. Portugal Pereira, P. Vyas, E. Huntley, K. Kissick, M. Belkacemi, J. Malley, (eds.)], In press. Ishida, K., Ercan, A., Trinh, T., Kavvas, M.L., Ohara, N., Carr, K., Anderson, M.L., 2018: Analysis of future climate change impacts on snow distribution over mountainous watersheds in Northern California by means of a physically-based snow distribution model, Science of The Total Environment, 645, 10651082 Islam, S.U., Déry, S.J., Werner, A.T., 2017: Future climate change impacts on snow and water resources of the Fraser River Basin, British Columbia, Journal of Hydrometeorology, 18(2):473-496, doi: 10.1175/JHMD-16-0012.1 Jarosch, A.H., Anslow, F.S., and Clarke, G.K.C., 2010: High-resolution precipitation and temperature downscaling for glacier models, Climate Dynamics, 28(1-2), 391–409, doi: 10.1007/s00382-010-09491 Jóhannesson, T., Björnsson, H., Magnússon, E., Guðmundsson, S.,Pálsson, F., Sigurðsson, O., Thorsteinsson, T., and Berthier, E., 2013: Ice-volume changes, bias estimation of mass-balance measurements and changes in subglacial lakes derived by lidar mapping of the surface Icelandic glaciers, Annals of Glaciology, 54, 63–74, doi:10.3189/2013AoG63A422 Johnson, R., and Ohara, N., 2018: Dynamic equilibrium modelling of snow and inland glaciers under the evolving climate in Wyoming, Journal of Hydrologic Engineering, 23(1), 04017056, doi: 10.1061/(ASCE)HE.1943-5584.0001585 Jonas, T., Marty, C., Magnusson, J., 2009: Estimating the snow water equivalent from snow depth measurements in the Swiss Alps, Journal of Hydrology, 378 (1–2), 161-167. Jordan, R., 1991: A one-dimensional temperature model for a snow cover: Technical documentation for SNTERERM.89. Special Rep. 91-16. Cold Region Research and Engineers Laboratory. U.S. Army Corps of Engineers, Hanover, NH, 61 pp. Kang, D.H., Shi, X., Gao, H., Déry, S.J., 2014: On the changing contribution of snow to the hydrology of the Fraser River Basin, Canada, Journal of Hydrogeometeorology, 15, 1344-1365 Kershaw, G.P., & McCulloch, J., 2007: Midwinter snowpack variation across the Arctic Treeline, Churchill, Manitoba, Canada. Arctic, Antarctic, and Alpine Research, 39(1), 9-15 Khadka, A., 2013: Predicting the effects of different land-use scenarios on water availability using hydrological model, Tropical Resources Bulletin, 32, 32-33. Kind, R. J., 1981: Snow drifting. In Gray, D. M., and Male, D. H. (eds.), Handbook of Snow, Principles, Processes, Management, and Use, Oxford: Pergamon Press, 338–359. Kind, R.J., 1992: One-dimensional aeolian suspension above beds of loose particles – a new concentration-profile equation, Atmospheric Environment, 26A (5), 927-931 Kirnbauer, R., Blӧschl, G. and Gutknecht, D., 1994: Entering the era of distributed snow models, Nordic Hydrology, 25(1–2), 1–24. Klock, R., and Mullock, J., 2002: The weather of British Columbia, Graphic area forecast 31, Nav Canada, 5556. 143 Koch, S. E., M. DesJardins, and P. J. Kocin, 1983: An interactive Barnes objective map analysis scheme for use with satellite and conventional data, Journal of Applied Climatology and Meteorology, 22, 1487-1503 Kohler, T., and Maselli, D., (eds) 2009: Mountains and climate change - From understanding to action. Published by Geographica Bernensia with the support of the Swiss Agency for Development and Cooperation (SDC), and an international team of contributors, Bern. Krabill, W., Thomas, R., Martin, C., Swift, R., and Frederick, E., 1995: Accuracy of airborne laser altimetry over the Greenland ice sheet, International Journal of Remote Sensing, 16(7), 1211–1222 Kuhn, M., 1995: The mass balance of very small glaciers. Zietschrift fur Gletscherkunde und Glazialgeologie, 31, 171–179. Kunkel, K. E., 1989: Simple procedures for extrapolation of humidity variables in the mountainous western United States, Journal of Climate, 2, 656-669 Larsen, C.F., Motyka, R.J., Arendt, A.A., Echelmeyer, K.A., and Geissler, P.E., 2007: Glacier changes in southeast Alaska and northwest British Columbia and contribution to sea level rise, Journal of Geophysical Research, 112, F01007 Lefebre, F., Gallée, H., Ypersele, J.P., and Greuell, W., 2003: Modelling of snow and ice melt at ETH Camp (west Greenland): A study of surface albedo, Journal of Geophysical Research, 108(D8), 4231 Lehning, M., Völksch, I., Gustafsson, D., Nguyen, T.A., Stähli, M., and Zappa, M., 2006: ALPINE3D: a detailed model of mountain surface processes and its application to snow hydrology, Hydrological Processes, 20(10), 2111-2128 Lehning, M., Bartelt, P., Brown, B., Fierz, C., Satyawali, P., 2002a: A physical SNOWPACK model for the Swiss avalanche warning: Part II. Snow microstructure, Cold Regions Science and Technology, 35(3), 147167 Lehning, M., Bartelt, P., Brown, B., Fierz, C., 2002b: A physical SNOWPACK model for the Swiss avalanche warning: Part III: meteorological forcing, thin layer formation and evaluation, Cold Regions Science and Technology, 35(3), 169–184 Lehning, M., Lowe, H., Ryser, M., and Radeschall, N., 2008: Inhomogeneous precipitation distribution and snow transport in steep terrain, Water Resources Research. 44 (W7), W07404. doi:10.1029/2007WR006545 Lemke, P., Ren, J., Alley, R.B., Allison, I., Carrasco, J., Flato, G., Fujii, Y., Kaser, G., Mote, P., Thomas, R.H., and Zhang, T., 2007: Observations: Changes in Snow, Ice and Frozen Ground. In: Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change [Solomon, S., D. Qin, M. Manning, Z. Chen, M. Marquis, K.B. Averyt, M. Tignor and H.L. Miller (eds.)]. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA. Lewis, C.S., Allen, L.N., 2016: Potential crop evapotranspiration and surface evaporation estimates via a gridded weather forcing dataset, Journal of Hydrology, 456, 450-463, doi: 10.1016/j.jhydrol.2016.11.055 Liston, G. E., 1995: Local advection of momentum, heat, and moisture during the melt of patchy snow covers, Journal of Applied Meteorology, 34, 1705-1715 Liston, G. E., and Elder, K., 2004: Representing subgrid snow cover heterogeneities in regional and global models. Journal of Climate, 17, 1381-1397 144 Liston, G.E., Elder, K., 2006a: A meteorological distribution system for high-resolution terrestrial modelling (MicroMet), Journal of Hydrometeorology, 7, 217-234 Liston, G.E., Elder, K., 2006b: A distributed snow-evolution modelling system (SnowModel), Journal of Hydrometeorology, 7(6), 1259-1276 Liston, G. E., and Hall, D.K., 1995: An energy balance model of lake ice evolution, Journal of Glaciology, 41, 373-382 Liston, G.E., Hiemstra, C.A., Elder, K., and Cline, D.W., 2008: Mesocell study area snow distributions for the Cold Land Processes Experiment (CLPX), Journal of Hydrometeorology, 9, 957-976 Liston, G. E., and C. A. Hiemstra., 2011: The changing cryosphere: Pan-Arctic snow trends (1979–2009), Journal of Climate, 24, 5691–5712 Liston, G. E., Haehnel, R. B., Sturn, M., Hiemstra, C. A., Berezovskaya, S. and Tabler, R. D., 2007: Instruments and methods simulating complex snow distributions in windy environments using SnowTran-3D, Journal of Glaciology, 53(181), 241–256 Liston, G. E., and Mernild, S. H., 2012: Greenland freshwater runoff. Part I: A runoff routing model for glaciated and nonglaciated landscapes (HydroFlow). Journal of Climate, 25, 5997-6014, doi: 10.1175/JCLI-D-1100591.1 Liston, G. E., and Sturm, M.,1998: A snow-transport model for complex terrain, Journal of Glaciology, 44, 498516 Liston, G.E., and Sturm, M., 2002: Winter precipitation patterns in arctic Alaska determined from a blowingsnow model and snow-depth observations, Journal of Hydrometeorology, 3, 646–659 Liston, G. E., Winther, J.G., Bruland, O., Elvehøy, H., and Sand, K., 1999b: Below-surface ice melt on the coastal Antarctic ice sheet, Journal of Glaciology, 45, 273-285 Liston, G.E., and Sturm, M., 2004: The role of winter sublimation in the Arctic moisture budget, Nordic Hydrology, 35(4-5), 325–334 Liston, G. E., Winther, J.G., Bruland, O., Elvehøy, H., and Sand, K., 1999: Below-surface ice melt on the coastal Antarctic ice sheet, Journal of Glaciology, 45, 273-285 Liston, G.E., and Hiemstra, C.A., 2008: A simple data assimilation system for complex snow distributions (SnowAssim), Journal of Hydrometeorology, 9, 989-1003, doi: 10.1175/2008JHM871.1 Liu, C., Ikeda, K., Thompson, G., Rasmussen, R., and Dudhia, J., 2011: High-resolution simulations of wintertime precipitation in the Colorado headwaters region: Sensitivity to physics parameterizations, Monthly Weather Review, 139(11), 3533-3553, doi:10.1175/mwr-d-11-00009.1 Liu, L., Ma, Y., Menenti, M., Zhang, X., Ma, W., 2019: Evaluation of WRF modelling in relation to different land surface schemes and initial and boundary conditions: A snow event simulation over the Tibetan Plateau, Journal of Geophysical Research: Atmospheres, 124(1), 209-226, doi.org/10.1029/2018JD029208 Lott, F.C., Lundquist, J.D., 2008: Modelling spatial differences in snowmelt runoff timing. In: Proceedings of the 76th annual Western Snow Conference, Hood River Oregon, 91–97. MacDonald, M. K., Pomeroy, J.W., and Pietroniro, A., 2010: On the importance of sublimation to an alpine snow mass balance in the Canadian Rocky Mountains, Hydrology and Earth System Sciences, 14, 1401-1415. 145 Machguth, H., Paul, F., Hoelzle, M. and Haeberli, W., 2006a: Distributed glacier mass-balance modelling as an important component of modern multi-level glacier monitoring, Annals of Glaciology, 43, 335–343, doi: 10.3189/172756406781812285 Machguth, H., Eisen, O., Paul F., and Hoelzle, M., 2006b: Strong spatial variability of snow accumulation observed with helicopter-borne GPR on two adjacent Alpine glaciers, Geophysical Research Letters, 33 L13503, doi:10.1029/2006GL026576 Mankin, J., and Diffenbaugh, N., 2015: Influence of temperature and precipitation variability on near-term snow trends, Climate Dynamics, 45, 1099–1116, doi:10.1007/s00382-014-2357-4 Maraun, D., Wetterhall, F., Ireson, A.M., Chandler, R.E., Kendon, E.J., Widmann, M., Brienen, S., Rust, H.W., Sauter, T., Themeßl, M., Venema, V.K.C., Chun, K.P., Goodess, C.M., Jones, R.J., Onof, C., Vrac, M., Thiele‐Eich, I., 2010: Precipitation downscaling under climate change: Recent developments to bridge the gap between dynamical models and the end user, Reviews of Geophysics, 48(3), doi:10.1029/2009RG000314 Marke, T., Strasser, U., Hanzer, F., Stötter, J., Wilcke, R. A. I., and Gobiet, A., 2014: Scenarios of future snow conditions in Styria (Austrian Alps), Journal of Hydrometeorology, 16, 261–277, doi:10.1175/jhm-d-140035.1 Marsh, C. B., Pomeroy, J. W., Spiteri, R.J., and Wheater, H. S., 2020: A finite volume blowing snow model for use with variable resolution meshes, Water Resources Research, 56, e2019WR025307. doi.org/10.1029/2019WR025307 Marsh, P., 1999: Snowcover formation and melt: recent advances and future prospects. Hydrological Processes, 13, 14-15, doi: 10.1002/(SICI)1099-1085(199910)13:14/15<2117::AID-HYP869>3.0.CO;2-9 Marty, C., Schlögl, S., Bavay, M., and Lehning, M., 2017: How much can we save? Impact of different emission scenarios on future snow cover in the Alps, The Cryosphere, 11, 517–529 Maurer, E., Hidalgo, H., Das, T., Dettinger, M., Cayan, D., 2010: The utility of daily large-scale climate data in the assessment of climate change impacts on daily streamflow in California, Hydrology and Earth System Sciences, 14, 1125-1138. Mayr, E., Hagg, W., Mayer, C., Braun, L., 2013: Calibrating a spatially distributed conceptual hydrological model using runoff, annual mass balance and winter mass balance, Journal of Hydrology, 487, 40-49 McGrath, D., Sass, L., O'Neel, S., McNeil, C., Candela, S. G., Baker, E. H., and Marshall, H.P., 2018: Interannual snow accumulation variability on glaciers derived from repeat, spatially extensive ground-penetrating radar surveys. The Cryosphere, 12, 3617-3633, doi: 10.5194/tc-12-3617-2018 Meng, X., Lyu, S., Zhang, T., Zhao, L., Li, Z., Han, B., Li, S., Ma, D., Chen, H., Ao, Y., Luo, S., Shen, Y., Guo, J., Wen, L., 2018: Simulated cold bias being improved by using MODIS time-varying albedo in the Tibetan Plateau in WRF model, Environmental Research Letters, 13, 044028. Menounos, B., Hugonnet, R., Shean, D., Gardner, A., Howat, I., Berthier, E., Pelto, B., Tennant, C., Shea, J., Noh, M., Brun, F., Dehecq, A., 2019: Heterogeneous changes in western North American glaciers linked to decadal variability in zonal wind strength, Geophysical Research Letters, 46, 200–209, doi:/10.1029/2018GL080942 Mernild, S.H., Liston, G.E., Hasholt, B., and Knudsen, N.T., 2006: Snow distribution and melt modelling for Mittivakkat Glacier, Ammassalik Island, SE Greenland. Journal of Hydrometeorology, 7, 808-824 146 Mernild, S.H., Hasholt, B., and Liston, G.E., 2006a: Water flow through Mittivakkat Glacier, Ammassalik Island, SE Greenland. Danish Journal of Geography, 106(1), 25-43 Mernild, S.H., Liston, G.E., and Hasholt, B., 2006b: Snow-distribution and melt modelling for glaciers in Zackenberg River Drainage Basin, NE Greenland, Hydrological Processes, 21, 3249-3263, doi: 10.1 002/hyp.6500. Mernild, S.H., Liston, G.E., Hiemstra, C.A., 2014: Northern Hemisphere glaciers and ice caps surface mass balance and contribution to sea-level rise, Journal of Climate, 27(15), 6051–6073, doi: 10.1175/JCLI-D13-00669.1. Mernild, S.H., Liston, G.E., Hiemstra, C.A., 2017: The Andes Cordillera. Part III: glacier surface mass balance and contribution to sea level rise (1979–2014), International Journal of Climatology, 37, 3154–3174 Mitchell, K.E., Lohmann, D., Houser, P.R., Wood, E.F., Schaake, J.C., Robock, A., Cosgrove, B.A., Sheffield, J., Duan, Q., Luo, L., Higgins, R.W., Pinker, R.T., Tarpley, J.D., Lettenmaier, D.P., Marshall, C.H., Entin, J.K., Pan, M., Shi, W., Koren, V., Meng, J., Ramsay, B.H., Bailey, A.A., 2004: The multiinstitution North American Land Data Assimilation System (NLDAS): Utilizing multiple GCIP products and partners in a continental distributed hydrological modelling system, Journal of Geophysical Research, 109, D07S90 Mölg, T., and Kaser, G., 2011: A new approach to resolving climate-cryosphere relations: Downscaling climate dynamics to glacier-scale mass and energy balance without statistical scale linking, Journal of Geophysical Research, 116, D16101 Molotch, N.P., Blanken, P.D., Williams, M.W., Turnipseed, A.A., Monson, R.K., Margulis, S.A., 2007: Estimating sublimation of intercepted and sub-canopy snow using eddy covariance system. Hydrological Processes, 21(12), 1567-1575 Moriasi, D. N., Arnold, J. G., Van Liew, M. W., Bingner, R. L., Harmel, R. D., Veith, T. L., 2007: Model evaluation guidelines for systematic quantification of accuracy in watershed simulations, Transactions of the ASABE, 50 (3), 885–900. Doi:10.13031/2013.23153. Mote, P.W., and Salathé, E.P., 2010: Future climate in the Pacific Northwest, Climatic Change, 102, 29–50, doi: 10.1007/s10584-010-9848-z Mote, P.W., Li, S., Lettenmaier, D.P., Xiao, M., and Engel, R., 2018: Dramatic declines in snowpack in the western US, Climate and Atmospheric Science, 1, 2, doi:10.1038/s41612-018-0012-1 Mott, R., Faure, F., Lehning, M., Lowe, H., Hynek, B., Michlmayer, G., Prokop, A., Schoner, W., 2008: Simulation of seasonal snow-cover distribution for glacierized sites on Sonnblick, Austria, with the Alpine3D model, Annals of Glaciology, 49, 155-160 Mott, R., Schirmer, M., Bavay, M., Grünewald, T., Lehning, M., 2010: Understanding snow-transport processes shaping the mountain snow-cover, The Cryosphere, 4, 545–559 Mott, R., Scipión, D., Schneebeli, M., Dawes, N., Berne, A., Lehning, M., 2014: Orographic effects on snow deposition patterns in mountainous terrain, Journal of Geophysical Research: Atmospheres, 119(3), 1419-1439 Mott, R., and Lehning, M., 2010: Meteorological modelling of very high-resolution wind fields and snow deposition for mountains, Journal of Hydrometeorology, 11, 934-949 Mott, R, Vionnet, V., and Grünewald T., 2018: The seasonal snow cover dynamics: Review on wind-driven coupling processes, Frontiers in Earth Science, 6, 197, doi: 10.3389/feart.2018.00197 147 Mott, R., Wolf, A., Kehl, M., Kunstmann, H., Warscher, M., and Grünewald, T., 2019: Avalanches and micrometeorology driving mass and energy balance of the lowest perennial ice field of the Alps: a case study, The Cryosphere, 13, 1247–1265, doi:10.5194/tc-13-1247-2019 Moya-Álvarez, A.S., Martínez-Castro, D., Flores, J.L., Silva, Y.,: 2018: Sensitivity Study on the Influence of Parameterization Schemes in WRF_ARW Model on Short- and Medium-Range Precipitation Forecasts in the Central Andes of Peru, Advances in Meteorology, vol. 2018, Article ID 1381092, 16 pages, doi:10.1155/2018/1381092 Moyer, A. N., Moore, R.D., and Koppes, M. N., 2016: Streamflow response to the rapid retreat of a lake-calving glacier, Hydrological Processes, 30 (20), 3650-3665, doi:10.1002/hyp.10890. Murdock, T.Q., and Werner, A.T., 2011: Canadian Columbia Basin Climate Trends and Projections: 2007-2010 Update. Pacific Climate Impacts Consortium, University of Victoria, Victoria. BC, 43. Musselman, K.N., Pomeroy, J.W., Essery, R.L.H., and Leroux, N., 2015: Impact of windflow calculations on simulations of alpine snow accumulation, redistribution and ablation, 29(18), 3983-3999, Hydrological Processes, doi: 10.1002/hyp.10595 Musselman, K. N., Lehner, F., Ikeda, K., Clark, M. P., Prein, A. F., Liu, C., Barlage, M., and Rasmussen, R., 2018: Projected increases and shifts in rain‐on‐snow flood risk over western North America, Nature Climate Change, 8(9), 808– 812, doi:10.1038/s41558‐018‐0236‐4 Najafi, M. R., Zwiers, F., and Gillett, N., 2017: Attribution of the observed spring snowpack decline in British Columbia to anthropogenic climate change, Journal of Climate, 30, 4113–4130, doi.org/10.1175/JCLID-16-0189.1. Nakanishi, M., and Niino, H., 2006: An improved Mellor–Yamada level-3 model: Its numerical stability and application to a regional prediction of advection fog, Boundary-Layer Meteorology, 119(2), 397-407. Niu, G. Y., Yang, Z. L., Mitchell, K. E., Chen, F., Ek, M. B., Barlage, M., Kumar, A., Manning, K., Niyogi, D., Rosero, E., Tewari, M., Xia, Y., 2011: The community Noah land surface model with multiparameterization options (Noah‐MP): 1. Model description and evaluation with local‐scale measurements, Journal of Geophysical Research: Atmospheres, 116, D12109. Nover, D.M., Witt, J.W., Butcher, J.B., Johnson, T.A., and Weaver, C.P., 2016: The effects of downscaling method on the variability of simulated watershed response to climate change in five U.S. basins, Earth Interactions, 20(11), 1–27, doi:10.1175/EI-D-15-0024.1 Ohmura, A., 2001: Physical basis for the temperature-based melt-index method, Journal of Applied Meteorology, 40(4), 753-761, doi: 10.1175/1520-0450(2001)040<0753:PBFTTB>2.0.CO;2 Østrem, G. and Brugman, M., 1991: Glacier mass balance measurements: a manual for field and office work, Saskatoon: National Hydrology Research Institute, and Oslo: the Norwegian Water Resources and Electricity Board. Pacific Climate Impacts Consortium (PCIC), University of Victoria, (Feb. 2019). Statistically Downscaled Climate Scenarios, Downloaded from https://data.pacificclimate.org/portal/downscaled_gcms/map/ Pan, M., Sheffield, J., Wood, E. F., Mitchell, K. E., Houser, P. R., Schaake, J. C., Robock, A., Lohmann, D., Cosgrove, B., Duan, Q., Luo, L., Higgins, R. W., Pinker, R. T., Tarpley, J. D., 2003: Snow process modelling in the North American Land Data Assimilation System (NLDAS). Part II: Evaluation of model simulated snow water equivalent, Journal of Geophysical Research, 108(D22), 8850, doi:10.1029/2003JD003994, in press, 2003. 148 Pan, L., Liu, Y., Knievel, J.C., Monache, L.D., and Roux, G., 2018: Evaluations of WRF Sensitivities in Surface Simulations with an Ensemble Prediction System, Atmosphere, 9, 106, doi:10.3390/atmos9030106 Park, H., Yabuki, H., Ohata, T., 2012: Analysis of satellite and model datasets for variability and trends in Arctic snow extent and depth, 1948–2006, Polar Science, 6(1), 23-37 Paul, F., and Kotlarski, S., 2010: Forcing a distributed glacier mass balance model with the regional climate model REMO. Part II: Downscaling strategy and results for two Swiss glaciers, Journal of Climate, 23, 16071620. Pelt, W.J.J., Kohler, J., Liston, G.E., Hagen, J.O., Luks, B., Reijmer, C.H., Pohjola, V.A., 2016: Multidecadal climate and seasonal snow conditions in Svalbard, Journal of Geophysical Research: Earth Surface, 121, 2100–2117, doi:10.1002/2016JF003999. Pelto, B.M., Menounos, B., and Marshall, S.J., 2019: Multi-year evaluation of airborne geodetic surveys to estimate seasonal mass balance, Columbia and Rocky Mountains, Canada, The Cryosphere, 13, 17091727, doi: 10.5194/tc-2019-30. Pelto, B.M., 2020: An Approach to Remotely Monitor Glacier Mass Balance at Seasonal to Annual Time Scales, Columbia and Rocky Mountains, Canada. University of Northern British Columbia. Pohl, S., Marsh, P., and Liston, G.E., 2007: Spatial–temporal variability in turbulent fluxes during spring snowmelt, Arctic Antarctic and Alpine Research, 38, 136–146 Pomeroy, J.W., Marsh, P., Gray, D.M., 1998: Application of a distributed blowing snow model to the Arctic, Hydrological Processes, 11(11), 1451-1464 Pomeroy, J.W., and Gray, D.M., 1990: Saltation of snow, Water Resources Research, 26(7), 1583-1594 Pomeroy, J.W., Gray, D.M., and Landine, P.G., 1993: The Prairie Blowing Snow Model: characteristics, validation, operation, Journal of Hydrology, 144(1-4): 165-192 Pomeroy, J.W., Jones, H.G., 1996: Wind-blown snow: sublimation, transport and changes to polar snow. In Chemical Exchange between the Atmosphere and Polar Snow, Wolff EW, Bales RC (eds); NATO ASI Series I. Springer-Verlag: Berlin, 43, 453-490 Pomeroy, J.W., and Li, L., 1997: Development of the Prairie Blowing Snow Model for application in climatological and hydrological models, Proceedings of the 65 th Annual Western Snow Conference, 1, 86- 197 Pomeroy, J.W., Bewley, D.S., Essery, R.L.H., Hedstrom, N.R., Link, T., Granger, R.J., Sicart, J.E., Ellis, C.R., and Janowicz, J.R., 2006: Shrub tundra snowmelt, Hydrological Processes, 20(4), 923–941, doi:10.1002/hyp.6124. Prasad, R., Tarboton, D., Liston, G.E., Luce, C., Seyfried, M., 2001: Testing a blowing snow model against distributed snow measurements at Upper Sheep Creek, Idaho. United States of America, Water Resources Research, 37, 1341-1356 Pytlak, E., Frans, C., Duffy, K., Johnson, J., Nijssen, B., Chegwidden, O., and Rupp, D., 2018: Climate and Hydrology Datasets for RMJOC Long-Term Planning Studies: Second Edition (RMJOC-II), Part I: Hydroclimate Projections and Analyses, Chapter 6 and 7. Quéno, L., Vionnet, V., Cabot, F., Vrécourt, D., Dombrowski-Etchevers, I., 2018: Forecasting and modelling ice layer formation on the snowpack due to freezing precipitation in the Pyrenees, Cold Regions Science and Technology, 146, 19-31, doi: 10.1016/j.coldregions.2017.11.007 149 Radić, V., and Hock, R., 2006: Modelling future glacier mass balance and volume changes using ERA-40 reanalysis and climate models: A sensitivity study at Storglaciären, Sweden, Journal of Geophysical Research, 111, F03003 Radić, V., and Clarke, G.K.C., 2011: Evaluation of IPCC models performance in simulating late 20th century climatologies and weather patterns over North America, Journal of Climate, 24, 5257-5274 Ravazzani, G., Davolio, S., Ceppi, A., Fiore, A., 2016: Wind speed interpolation for hydrological modelling in complex topography area. Numerical Computations: Theory and Algorithms, The 2nd International Conference and Summer School, NUMTA2016, 19 – 25 June 2016, Club Med Resort “Napitia”, Pizzo Calabro, Calabria, Italy. Reijmer, C.H., Hock, R., 2008: Internal accumulation on Storglaciaren, Sweden, in a multi-layer snow model coupled to a distributed energy- and mass-balance model, Journal of Glaciology, 54(184), 61-72 Ritter, A., and Muñoz-Carpena, R., 2013: Performance evaluation of hydrological models: statistical significance for reducing subjectivity in goodness-of-fit assessments, Journal of Hydrology, 480 (1), 33-45, doi:10.1016/j.jhydrol.2012.12.004 Rogelj, J., Elzen, M., Höhne, N., Fransen, T., Fekete, H., Winkler, H., Schaeffer, R., Sha, F., Riahi, K., and Meinshausen, M., 2016: Paris Agreement climate proposals need a boost to keep warming well below 2°C. Nature, 534 (7609), 631–639, doi:10.1038/nature18307 Rousselot, M., Durand, Y., Giraud, G., Mérindol, L., Dombrowski-Etchevers, I., Déqué, M., and Castebrunet, H., 2012: Statistical adaptation of ALADIN RCM outputs over the French Alps – application to future climate and snow cover, The Cryosphere, 6, 785–805, doi:10.5194/tc-6-785-2012 Rupp, D. E., Abatzoglou, J. T., and Mote, P. W., 2017: Projections of 21st century climate of the Columbia River basin, Climate Dynamics, 49, 1783–1799, doi: 10.1007/s00382-016-3418-7 Salzmann, N., Machguth, H., Linsbauer, A., 2012: The Swiss Alpine glaciers’ response to the global 2˚C air temperature target, Environmental Research Letters, 7, 044001 Santos-Alamillos, F.J., Pozo-Vázquez, D., Ruiz-Arias, J.A., Lara-Fanego, V., and Tovar-Pescador, J., 2013: Analysis of WRF model wind estimate sensitivity to physics parameterization choice and terrain representation in Andalusia (Southern Spain), Journal of Applied Meteorology and Climatology, 52, 1592-1609 Sauter, T., Moller, M., Finkelnburg, R., Grabiec, M., Scherer, D., and Schneider, C., 2013: Snowdrift modelling for the Vestfonna ice cap, north-eastern Svalbard, The Cryosphere, 7, 1287-1301 Sauter, T., and Obleitner, F., 2015: Assessing the uncertainty of glacier mass-balance simulations in the European Arctic based on variance decomposition, Geoscientific Model Development, 8, 3911–3928, doi:10.5194/gmd-8-3911-2015 Schiefer, E., Menounos, B., and Wheate, R., 2007: Recent volume loss of British Columbian glaciers, Canada, Geophysical Research Letters, 34, L16503 Schlögl, S., Marty, C., Bavay, M., and Lehning, M., 2016: Sensitivity of Alpine3D modeled snow cover to modifications in DEM resolution, station coverage and meteorological input quantities, Environmental Modelling and Software, 83, 387–396 Schmucki, E., Marty, C., Fierz, C., Lehning, M., 2014: Evaluation of modelled snow depth and snow water equivalent at three contrasting sites in Switzerland using SNOWPACK simulations driven by different meteorological data input, Cold Regions Science and Technology, 99, 27–37. 150 Schmucki, E., Marty, C., Fierz, C., and Lehning, M., 2015: Simulations of 21st century snow response to climate change in Switzerland from a set of RCMs, International Journal of Climatology, 35, 3262–3273, doi:10.1002/joc.4205 Schneider, U., Becker, A., Finger, P., Meyer-Christoffer, A., Ziese, M., and Rudolf, B., 2014: GPCC’s new land surface precipitation climatology based on quality-controlled in situ data and its role in quantifying the global water cycle, Theoretical and Applied Climatology, 115, 15–40, doi:10.1007/s00704-013-0860-x. Schultz, D.M., Steenburgh, W. J., Trapp, R. J., Horel, J., Kingsmill, D. E., Dunn, L. B., Rust, W. D., Cheng, L., Bansemer, A., Cox, J., Daugherty, J., Jorgensen, D. P., Meitin, J., Showell, L., Smull, B. F., Trap, K., and Trainor, M., 2002: Understanding Utah winter storms: The Intermountain Precipitation Experiment, Bulletin of American Meteorological Society, 83(2), 189-210 Sevruk, B., 1985: Correction of monthly precipitation for wetting losses. TECIMO III. WMO, Instrument and observing methods report, 22, 7-12. Shea, J. M., Marshall, S. J., and Livingston, J. M., 2004: Glacier distributions and climate in the Canadian Rockies, Arctic, Antarctic, and Alpine Research, 36(2), 272-279 Shi, Y., Gao, X., Wu, J., Giorgi, F., 2011: Changes in snow cover over China in the 21st century as simulated by a high-resolution regional climate model, Environmental Research Letters, 6, 045401, doi: 10.1088/1748-9326/6/4/045401 Shrestha, N. K., Du, W., and Wang, J., 2017: Assessing climate change impacts on freshwater resources of the Athabasca River Basin, Canada, Science of the Total Environment, 601-602, 425-440, doi: 10.1016/j.scitotenv.2017.05.013. Shrivastava, R., Dash, S.K., Oza, R.B., Hegde, M.N., 2015: Evaluation of parameterization schemes in the Weather Research and Forecasting (WRF) model: A case study for the Kaiga nuclear power plant site, Annals of Nuclear Energy, 75, 693-702 Skamarock, W. C., and Klemp, J. B., 2008: A time-split nonhydrostatic atmospheric model for weather research and forecasting applications, Journal of Computational Physics, 227(7), 3465-3485. Skamarock, W.C., Klemp, J.B., Dudhia, J., Gill, D.O., Barker, D.M., Duda, M.G., Huang, X.U., Wang, W., Powers, J.G., 2008: A Description of the Advanced Research WRF Version 3. NCAR Technical Note Smith, R.B., and Barstad, I., 2004: A linear theory of orographic precipitation, Journal of Atmospheric Science, 61, 1377–1391 Sold, L., Huss, M., Hoelzle, M., Andereggen, H., Joerg, P.C., Zemp, M., 2013: Methodological approaches to infer end-of-winter snow distribution on alpine glaciers, Journal of Glaciology, 59(218), 1047-1059 Steger, C., Kotlarski, S., Jonas, T., and Schär, C., 2013: Alpine snow cover in a changing climate: a regional climate model perspective, Climate Dynamics, 41, 735–754, doi:10.1007/s00382-012-1545-3 Stergiou, I., Tagaris, E., and Sotiropoulou, R.P., 2017: Sensitivity Assessment of WRF Parameterizations over Europe, Proceedings 2017. 1. 119., doi:10.3390/ecas2017-04138 Stocker, T.F., Clarke, G.K.C., Le Treut, H., Lindzen, R.S., Meleshko, V.P., Mugara, R.K., ... Holtslag, A.A.M., 2001: Physical Climate Processes and Feedbacks. In J. T. Houghton, Y. Ding, D. J. Griggs, M. Noguer, P. J. van der Linden, X. Dai, & K. Maske (Eds.), IPCC, 2001: Climate Change 2001: The Scientific Basis, Contribution of Working Group I to the Third Assessment Report of the Intergovernmental Panel on Climate Change (pp. 417-470), Cambridge and New York: Cambridge University Press. 151 Strasser, U., Bernhardt, M., Weber, M., Liston, G.E., Mauser, W., 2007: Is snow sublimation important in the alpine water balance? The Cryosphere, 2, 53-66 Sturm, M., and Liston, G., 2003: The snow cover on lakes of the Arctic Coastal Plain of Alaska, U.S.A., Journal of Glaciology, 49, 370–380. Sturm, M., and Wagner, A.M., 2010: Using repeated patterns in snow distribution modelling: An Arctic example, Water Resources Research, 46(W12), W12549, doi:10.1029/2010WR009434 Sturm, M., B. Taras, G. E. Liston, C. Derksen, T. Jonas, and J. Lea., 2010: Estimating snow water equivalent using snow depth data and climate classes, Journal of Hydrometeorology, 11, 1380–94 Teutschbein, C., Wetterhall, F., Seibert, J., 2011: Evaluation of different downscaling techniques for hydrological climate-change impact studies at the catchment scale, Climate Dynamics, 37, 2087–2105, doi:10.1007/s00382-010-0979-8 Thibert, E., Baroudi, D., Limam, A., and Berthet-Rambaud., P., 2008: Avalanche impact pressure on an instrumented structure, Cold Regions Science and Technology, 54, 206-215 Thompson, G., Field, P. R., Rasmussen, R. M., & Hall, W. D., 2008: Explicit forecasts of winter precipitation using an improved bulk microphysics scheme. Part II: Implementation of a new snow parameterization, Monthly Weather Review, 136(12), 5095-5115. Trachsel, M., and Nesje, A., 2015: Modelling annual mass balances of eight Scandinavian glaciers using statistical models, The Cryosphere, 9, 1401–1414, doi:10.5194/tc-9-1401-2015 Trentini, F., Leduc, M., and Ludwig, R., 2019: Assessing natural variability in RCM signals: comparison of a multi model EURO-CORDEX ensemble with a 50-member single model large ensemble. Climate Dynamics, 53:1963–1979. doi:10.1007/s00382-019-04755-8 U.S. DIBR: U.S. Department of the Interior Bureau of Reclamation, 2016: West-Wide climate risk assessment, Columbia River Basin climate impact assessment, Final report. Vionnet, V., Martin, E., Masson, V., Guyomarc’h, G., Naaim-Bouvet, F., Prokop, A., Durand, Y., and Lac, C., 2014: Simulation of wind-induced snow transport and sublimation in alpine terrain using a fully coupled snowpack/atmosphere model, The Cryosphere, 8, 395-415 Vionnet, V., Martin, E., Masson, V., Lac, C., Naaim Bouvet, F., and Guyomarc'h, G., 2017: High-resolution large eddy simulation of snow accumulation in alpine terrain, Journal of Geophysical Research: Atmospheres, 122, 11005–11021, doi: 10.1002/2017JD026947. Vionnet, V., Guyomarc, G., Lafaysse, M., Naaim-Bouvet, F., Giraud, G., Deliot, Y., 2018: Operational implementation and evaluation of a blowing snow scheme for avalanche hazard forecasting, Cold Regions Science and Technology, 147, 1-10 Vionnet, V., Six, D., Auger, L., Dumont, M., Lafaysse, M., Quéno, L., Réveillet, M., Dombrowski-Etchevers, I., Thibert, E., and Vincent, C., 2019: Sub-kilometer precipitation datasets for snowpack and glacier modelling in alpine terrain. Frontiers in Earth Science, 7, 182, doi: 10.3389/feart.2019.00182. Wang, W., Bruyere, C., Duda, M., Dudhia, J., Gill, D., Kavulich, M., Keene, K., Lin, H., Michalakes, J., Rizvi, S., Zhang, X., 2009: WRF-ARW Version 3 Modelling System User’s Guide. Mesoscale & Microscale Meteorology Division, National Center for Atmospheric Research, Boulder, USA Warren, F.J., and Lemmen, D.S., editors, 2014: Canada in a Changing Climate: Sector Perspectives on Impacts and Adaptation; Government of Canada, Ottawa, ON. 286p 152 Warscher, M., Strasser, U., Kraller, G., Marke, T., Franz, H., and Kunstmann, H., 2013: Performance of complex snow cover descriptions in a distributed hydrological model system: A case study for the high Alpine terrain of the Berchtesgaden Alps, Water Resources Research, 49, 2619–2637, doi:10.1002/wrcr.20219 Watson, F.G.R., Newman, W.B., Coughlan, J.C., and Garrortt, R.A., 2006: Testing a distributed snowpack simulation model against spatial observations, Journal of Hydrology, 328(3-4), 453-466. Wijngaard, R.R., Steiner, J.F., Kraaijenbrink, P.D.A., Klug, C., Adhikari, S., Banerjee, A., Pellicciotti, F., Beek, L.P.H., Bierkens, M.F.P., Lutz, A.F., and Immerzeel, W.W., 2019: Modelling the response of the Langtang Glacier and the Hintereisferner to a changing climate since the little ice age, Frontiers in Earth Sciences, 7, 143, doi: 10.3389/feart.2019.00143 Winstral, A., and Marks, D., 2002: Simulating Wind Fields and Snow Redistribution Using Terrain-Based Parameters to Model Snow Accumulation and Melt Over a Semi-Arid Mountain Catchment, 59th Eastern Snow Conference Stowe, Vermont. USA Winstral, A., and Marks, D., Gurney, R., 2013: Simulating wind-affected snow accumulations at catchment to basin scales, Advances in Water Resources, 55, 64-79 Winstral, A., Magnusson, J., Schirmer, M., & Jonas, T., 2019: The bias‐detecting ensemble: A new and efficient technique for dynamically incorporating observations into physics‐based, multilayer snow models, Water Resources Research, 55, 613– 631, doi.org/10.1029/2018WR024521 Xie, S.P., Deser, C., Vecchi, G.A., Collins, M., Delworth, T.L., Hall, A., Hawkins, E., Johnson, N.C., Cassou, C., Giannini, A., Watanabe, M., 2015: Towards predictive understanding of regional climate change. Nature Climate Change, 5:921–930, doi:10.1038/nclimate2689 Xue, M., Droegemeier, K.K., and Wong, V., 2000: The Advanced Regional Prediction System (ARPS) ± A multiscale nonhydrostatic atmospheric simulation and prediction model. Part I: Model dynamics and verification, Meteorology and Atmospheric Physics, 75, 161-193. Yhang, Y.B., and Hong, S.Y., 2008: A simulated climatology of the East Asian summer monsoon using a regional spectral model, Asia-Pacific Journal of Atmospheric Sciences, 44(4), 325–339 Yhang, Y.B., Sohn, S.J., and Jung, I.W., 2017: Application of dynamical and statistical downscaling to East Asian summer precipitation for finely resolved datasets, Advances in Meteorology, 2017, 2956373, doi: 10.1155/2017/2956373 Zemp, M., Haeberli, W., 2008: Glaciers and ice caps. Part I: Global overview and outlook. Part II: Glacier changes around the world, In Eamer, J., ed. Global outlook for ice and snow. Nairobi, United Nations Environment Programme (UNEP), Contributors: Bajracharya, S., Chinn, T. J., Fountain, A. G., Hagen, J. O., Huggel, C., Kääb, A., Kaltenborn, B. P., Karki, M., Kaser, G., Kotlyakov, V. M., Lambrechts, C., Li, Z. Q., Molnia, B. F., Mool, P., Nellemann, C., Novikov, V., Osipova, G. B., Rivera, A., Shrestha, B., Svoboda, F., Tsvetkov D. G., Yao, T. D., Zurich Open Repository and Archive, doi.org/10.5167/uzh40427 Zhang, H., Pu, Z., Zhang, X., 2013: Examination of errors in near-surface temperature and wind from WRF numerical simulations in regions of complex terrain, Weather and Forecasting, 28, 893-914 Zubler, E.M., Fischer, A.M., Liniger, M.A., Croci-Maspoli, M., Scherrer, S.C., Appenzeller, C., 2014: Localized climate change scenarios of mean temperature and precipitation over Switzerland, Climatic Change, 125, 237–252, doi: 10.1007/s10584-014-1144-x 153 Zwaaftink, C.D.G., Löwe, H., Mott, R., Bavay, M., Lehning, M., 2011: Drifting snow sublimation: A highresolution 3-D model with temperature and moisture feedbacks, Journal of Geophysical Research, 116, D16107 Zwaaftink, C.D.G., Mott, R., Lehning, M., 2013: Seasonal simulation of drifting snow sublimation in Alpine terrain, Water Resources Research, 49(3), 1581-1590 154