Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Evaluation of dynamically downscaled near-surface meteorological variables and energy fluxes at three… Tessema, Mekdes Ayalew 2018

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

Item Metadata

Download

Media
24-ubc_2018_september_tessema_mekdes.pdf [ 15.1MB ]
Metadata
JSON: 24-1.0368955.json
JSON-LD: 24-1.0368955-ld.json
RDF/XML (Pretty): 24-1.0368955-rdf.xml
RDF/JSON: 24-1.0368955-rdf.json
Turtle: 24-1.0368955-turtle.txt
N-Triples: 24-1.0368955-rdf-ntriples.txt
Original Record: 24-1.0368955-source.json
Full Text
24-1.0368955-fulltext.txt
Citation
24-1.0368955.ris

Full Text

Evaluation of dynamically downscaled near-surfacemeteorological variables and energy fluxes at three mountainglaciers in British ColumbiabyMekdes Ayalew TessemaA THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Geophysics)The University of British Columbia(Vancouver)July 2018© Mekdes Ayalew Tessema, 2018The following individuals certify that they have read, and recommend to the Faculty of Grad-uate and Postdoctoral Studies for acceptance, the thesis entitled:Evaluation of dynamically downscaled near-surface meteorological variables and energy fluxesat three mountain glaciers in British Columbiasubmitted by Mekdes Ayalew Tessema in partial fulfilment of the requirements forthe degree of Master of Sciencein GeophysicsExamining Committee Members:Valentina Radic´SupervisorDouw SteynSupervisory Committee MemberGwenn FlowersSupervisory Committee MemberRoland StullAdditional Examining Committee MemberAdditional Supervisory Committee Members:Brian MenounosAdditional Supervisory Committee MemberiiAbstractAll models of glacier melt, regardless of their complexity, must be forced by observed meteoro-logical fields at or in the vicinity of the glacier in question. In the absence of these observations,the forcing is commonly derived from statistical or dynamical downscaling of low resolutionclimate reanalysis models. Here we focus on a dynamical downscaling via Weather Researchand Forecasting (WRF) model, which has previously showed promising results in simulating asurface energy balance (SEB) at several glacierized terrains. Our goal is to evaluate the WRFdownscaling approach at three mountain glaciers in the interior mountains of British Columbiawhere the automatic weather stations (AWSs) recorded data over several summer seasons. TheWRF model, nested within the ERA-Interim global reanalysis produced output fields at 7.5km and 2.5 km spatial resolution, as well as 1 km resolution for one of the sites. We ana-lyze how closely the WRF model output, at sub-daily and daily temporal resolution, resemblesthe observed meteorological fields and SEB fluxes needed to assess seasonal surface melt atthese glaciers. We find that the model at 2.5 km closely simulates the cumulative seasonalmelt (±10% difference) despite large biases in the individual components of the SEB model.Overestimation of the number of clear sky days explains the positive bias in the modeled netshortwave radiation. This positive bias, however, is compensated by a negative bias in the mod-eled net longwave radiation, and by an underestimation of sensible and latent heat fluxes. Theunderestimation in the latter two fluxes, calculated from the bulk aerodynamic method, is dueto underestimated near-surface wind speeds. Radiative fluxes, which are dominant drivers ofseasonal melting, are poorly downscaled with WRF, while successfully simulated by the ERA-Interim at the course spatial resolution. Therefore, we advocate that SEB models be directlyforced with the output from global climate reanalysis. Finally, simulating turbulent heat fluxesat sloped glacier surfaces remains a major challenge, and the 1-km resolution state-of-the-artWRF model is not yet ready to tackle it.iiiLay SummaryMountain glaciers play a key role in the Earth’s hydrological cycle. Melting of mountainglaciers mainly occurs during summer. In order to determine how much melt occurs duringone summer season, a surface energy balance (SEB) model is used. Components of this modelare incoming and reflected sunlight, incoming and emitted infrared radiation, sensible andlatent (evaporative) heat transport, heat from rain and the underlying ground. The sum of theseheat transports is then converted to surface glacier melt. All the components of the SEB canbe directly measured by deploying Automatic Weather Stations (AWS) at the glacier surface,however, these measurements are not available at most glaciers. The goal of this study isto estimate the components of the SEB and surface melting using a computer based weatherforecast model. Results from the model are then compared with observations from AWS atthree mountain glaciers in British Columbia. Our results show that glacier melt estimatedusing the computer model agrees fairly well with the observations despite the poor agreementof the individual components of the SEB model with observations.ivPrefaceThis thesis is based on the output from a Weather Research and Forecasting (WRF) model. Ipreformed all the pre-processing of input data to the WRF model, model set-up and model runs.Post-processing and analysis of the WRF model output and manuscript write-up is performedby myself in collaboration with my supervisor Valentina Radic´ who provided help with thedevelopment of the methods and codes, and with the interpretation of the results.vTable of Contents. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation and research objectives . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Overview of glacier melt models . . . . . . . . . . . . . . . . . . . . . . . . . 31.2.1 Simple positive-degree day (PDD) model . . . . . . . . . . . . . . . . 41.2.2 Enhanced PDD model . . . . . . . . . . . . . . . . . . . . . . . . . . 51.2.3 Semi-empirical surface energy balance (SEB) model . . . . . . . . . . 51.2.4 Surface energy balance (SEB) models . . . . . . . . . . . . . . . . . . 61.3 Climate downscaling methods . . . . . . . . . . . . . . . . . . . . . . . . . . 101.3.1 Application of dynamical downscaling in glacier studies . . . . . . . . 112 Data and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.1 Field measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.2 WRF model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.3 Surface energy balance model . . . . . . . . . . . . . . . . . . . . . . . . . . 20vi2.4 Model evaluation and sensitivity tests . . . . . . . . . . . . . . . . . . . . . . 233 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.1 Simulated versus observed atmospheric conditions and melt rates . . . . . . . . 263.1.1 Towards distributed melt modeling . . . . . . . . . . . . . . . . . . . . 363.2 Cumulative ablation and surface height changes . . . . . . . . . . . . . . . . . 393.3 Sensitivity tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.3.1 Sensitivity to spatial resolution in the WRF model . . . . . . . . . . . . 413.3.2 Sensitivity to land cover . . . . . . . . . . . . . . . . . . . . . . . . . 423.3.3 Sensitivity to parameters in the SEB model . . . . . . . . . . . . . . . 443.3.4 SEB model versus PDD model performance . . . . . . . . . . . . . . . 454 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 475 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52A Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61A.1 Supplementary figures and tables . . . . . . . . . . . . . . . . . . . . . . . . . 61A.1.1 Time series of meteorological variables and SEB fluxes . . . . . . . . . 61A.1.2 Wind rose diagrams . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65A.1.3 Scatter plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67A.2 Guidelines for running the WRF model on Compute Canada computing facilities 81A.2.1 Preparation step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81A.2.2 WRF Pre-processing System (WPS) . . . . . . . . . . . . . . . . . . . 83A.2.3 Model integration Part (real.exe and wrf.exe) . . . . . . . . . . . . . . 88A.2.4 Incorporating new data sets into the WRF model . . . . . . . . . . . . 93A.2.5 Post processing tools . . . . . . . . . . . . . . . . . . . . . . . . . . . 97viiList of TablesTable 2.1 Characteristics of the study glaciers and their automatic weather stations(AWSs) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14Table 2.2 Observed variables and their sensor specifications at Conrad glacier. Onlyvariables used for the purpose of this study are listed. . . . . . . . . . . . . . 15Table 2.3 AWSs elevation (units in meters above sea level): actual vs model for allstudy sites. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19Table 2.4 Specifics of the WRF model setup for the three glaciers . . . . . . . . . . . 19Table 2.5 Roughness lengths (in m) for momentum (z0v), temperature (z0t) and hu-midity (z0q) used in this study, tuned fresh-snow densities (in kg m−3), andcalibrated site-specific melt factors ( fm) in the PDD model (mm w.e. ◦C−1day−1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23Table 3.1 Evaluation results for mod2.5 against the AWS variables: correlation, RMSEand MBE derived from 3-hr timeseries. Bold values are correlations signif-icantly different from zero at a significance level of 0.05. . . . . . . . . . . . 38Table 3.2 Evaluation results for mod2.5 against the AWS variables: correlation, RMSEand MBE derived from daily timeseries. Bold values are correlations sig-nificantly different from zero at a significance level of 0.05. . . . . . . . . . 39Table A.1 Evaluation results for mod1.0 (Castle Creek glacier only), mod2.5, mod7.5and ERA-Interim against the AWS variables: correlation (r), RMSE andMBE derived from daily timeseries. Bold values are correlations signifi-cantly different from zero at a significance level of 0.05. . . . . . . . . . . . 79Table A.2 ESA land category converted to USGS 24-category Land Use Categories . . 81viiiList of FiguresFigure 1.1 (a) Map of the study domain with the glacier sites. Contour topographicalmap is provided for each glacier, with marked site of automatic weatherstation (AWS): (b) Castle Creek glacier in 2010 and 2012 with one on-glacier AWS (pink circle) and two off-glacier AWSs (pink triangles); (c)Nordic glacier in 2014 with one on-glacier AWS (pink circle); (d) Conradglacier with two on-glacier AWSs in 2015 (pink circle and triangle), andtwo on-glacier AWSs in 2016 (blue circle and triangle). Outlines of theglaciers (white lines) are from the Randolph Glacier Inventory. . . . . . . . 2Figure 2.1 WRF domain set up a) d01 = 25 km, d02 = 5 km and d03 = 1 km (for CastleCreek glacier only) and b) d01 = 22.5 km, d02 = 7.5 km and d03 = 2.5 km(for all study sites). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16Figure 2.2 Topography and land category for the domain covering Castle Creek glacierfrom SRTM 30 m DEM and ESA CCI Landcover at 300 m resolution, re-spectively, in comparison to the topography and land category in the WRFnested domains at different spatial resolutions (7.5 km and 2.5 km, and 1km). Categories from the CCI Landcover data are manually converted tothe USGS 24-category Land Cover used in the WRF model (see Appendixfor the land category legend, Table A.2). Land category 24 (white colorin the color bar) is snow/ice. Markers indicate the on-glacier AWS (cir-cle), two off-glacier AWSs (triangles), and the selected grid cell with ice/s-now land category (star) to be used in the sensitivity tests. The outlines ofglaciers (black lines) are from the Randolph Glacier Inventory. . . . . . . . 17Figure 2.3 Topography and land category for the domain covering Nordic glacier fromSRTM 30 m DEM and ESA CCI Landcover at 300 m resolution, respec-tively, in comparison to the topography and land category in the WRFnested domains at different spatial resolutions (7.5 km and 2.5 km) . . . . . 18ixFigure 2.4 Topography and land category for the domain covering Conrad glacierfrom SRTM 30 m DEM and ESA CCI Landcover at 300 m resolution, re-spectively, in comparison to the topography and land category in the WRFnested domains at different spatial resolutions (7.5 km and 2.5 km) . . . . . 18Figure 2.5 Mean daily albedo as observed from AWS (dashed) and as calulated fromthe 3-hr albedo values prescibed for ice, firn and fresh snow to be 0.3, 0.65and 0.8, respectively (solid line). Note that the prescibed albedo valuesdiffer due to differences in observed or modeled snowfall periods. . . . . . 21Figure 3.1 Comparison of observed versus modeled 3-hr meteorological variables atCastle Creek glacier in 2010 and 2012: 2m-air temperature (T2m), 2m-relative humidity (RH2m), 10 m (modeled) and 2 m (observed) wind speed(U), surface pressure (SP) and accumulated precipitation (P). . . . . . . . . 27Figure 3.2 Comparison of observed versus modeled 3-hr meteorological variables atNordic glacier in 2014. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27Figure 3.3 Comparison of observed versus modeled 3-hr meteorological variables atConrad glacier in 2015, station 2. . . . . . . . . . . . . . . . . . . . . . . . 28Figure 3.4 Wind roses of observed and simulated wind speed at Castle Creek glacierin 2012. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29Figure 3.5 Wind roses of observed and simulated wind speed at Nordic glacier in 2014. 29Figure 3.6 Wind roses of observed and simulated wind speed at Conrad glacier in2015, station 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30Figure 3.7 Comparison of observed or AWS-derived versus WRF-derived 3-hr SEBfluxes at Castle Creek glacier in 2010 and 2012: incoming shortwave ra-diation (Kin), incoming longwave radiation (Lin), sensible heat flux (QH),latent heat flux (QE)) and melt energy flux (QM). . . . . . . . . . . . . . . 31Figure 3.8 Comparison of observed or AWS-derived versus WRF-derived 3-hr SEBfluxes at Nordic glacier in 2014. . . . . . . . . . . . . . . . . . . . . . . . 32Figure 3.9 Comparison of observed or AWS-derived versus WRF-derived 3-hr SEBfluxes at Conrad glacier in 2015, station 2. . . . . . . . . . . . . . . . . . . 32xFigure 3.10 Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Castle Creek glacier in 2012 season withcorrelation (r), p-value (p) and mean bias error (MBE) included for each plot. 33Figure 3.11 Comparison of AWS-derived (obs) versus WRF-derived (mod) daily time-series of the selected variables at Castle Creek glacier in 2012 season withcorrelation (r), p-value (p) and mean bias error (MBE) included for each plot. 34Figure 3.12 Comparison of mean daily melt energy flux and its partitioning into netshortwave radiation (Knet), net longwave radiation (Lnet), sensible heat flux(QH)and latent heat fluxes (QE) when observed albedo and roughness lengthsare used for each site. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35Figure 3.13 Comparison of mean daily melt energy flux and its partitioning into netshortwave radiation (Knet), net longwave radiation (Lnet), sensible heat flux(QH)and latent heat fluxes (QE) when modeled albedo and roughness lengths. . . 36Figure 3.14 Left panels: comparison of AWS-derived (obs) versus WRF-derived (mod2.5)daily timeseries of the selected variables from the two stations at Conradglacier in 2016: station 1 in the ablation area, and station 2 in the accumu-lation area. Right panels: the difference in each variable between the twostations (daily timeseries of station 2 minus daily timeseries of station 1),as AWS-derived and WRF-derived. . . . . . . . . . . . . . . . . . . . . . . 37Figure 3.15 Comparison of WRF-derived versus AWS-derived melt, converted to sur-face lowering (zM; set to negative values here), and net surface heightchanges (znet) at each site over the observational period using observedalbedo and roughness length. . . . . . . . . . . . . . . . . . . . . . . . . . 40Figure 3.16 Comparison of WRF-derived versus AWS-derived melt, converted to sur-face lowering (zM; set to negative values here), and net surface heightchanges (znet) at each site over the observational period using prescribedalbedo and roughness length. . . . . . . . . . . . . . . . . . . . . . . . . . 40Figure 3.17 Skill Scores (SSC) calculated for daily timeseries of the selected variablesat each study site. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42xiFigure 3.18 Results of the sensitivity test to varying landcover in WRF: comparison ofmean daily melt energy flux and its partitioning into net shortwave radiation(Knet), net longwave radiation (Lnet), sensible heat flux(QH) and latent heatfluxes (QE) when modeled albedo and roughness lengths are used in theSEB model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43Figure 3.19 Results of the sensitivity test to the choice of four parameterization schemesfor albedo (α) and roughness length (z0) in the SEB model forced bymod2.5: comparison of mean daily melt energy flux and its partitioninginto net shortwave radiation (Knet), net longwave radiation (Lnet), sensibleheat flux(QH) and latent heat fluxes (QE). . . . . . . . . . . . . . . . . . . 44Figure 3.20 Melt converted to surface lowering (zM) as simulated by the simple PDDmodel and by the SEB model (with observed albedo and roughness lengths)for each site. The PDD model is forced by 2m-temperature from AWSand from WRF at different spatial resolution. Results of the sensitivity testwhere the calibrated site-specific melt factor in the PDD model is perturbedby ± 1 mm w.e. day−1 ◦C−1 are shown for mod2.5 (PDD− and PDD+). . . 45Figure A.1 Comparison of observed versus modeled 3-hr meteorological variables atConrad glacier in 2015, station 1: 2m-air temperature (T2m), 2m-relativehumidity (RH2m), 10 m (modeled) and 2 m (observed) wind speed (U),surface pressure (SP) and accumulated precipitation (P). . . . . . . . . . . 61Figure A.2 Comparison of observed or AWS-derived versus WRF-derived 3-hr SEBfluxes at Conrad glacier in 2015, station 1: incoming shortwave radiation(Kin), incoming longwave radiation (Lin), sensible heat flux (QH), latentheat flux (QE) and melt energy flux (QM). . . . . . . . . . . . . . . . . . . 62Figure A.3 Comparison of observed versus modeled 3-hr meteorological variables atConrad glacier in 2016, station 1. . . . . . . . . . . . . . . . . . . . . . . . 63Figure A.4 Comparison of observed or AWS-derived versus WRF-derived 3-hr SEBfluxes at at Conrad glacier in 2016, station 1. . . . . . . . . . . . . . . . . 63xiiFigure A.5 Comparison of observed versus modeled 3-hr meteorological variables atConrad glacier in 2016, station 2. . . . . . . . . . . . . . . . . . . . . . . . 64Figure A.6 Comparison of observed or AWS-derived versus WRF-derived 3-hr SEBfluxes at at Conrad glacier in 2016, station 2. . . . . . . . . . . . . . . . . 64Figure A.7 Wind roses of observed and simulated wind speed at Conrad glacier in2015, station 1. Colorbar shows the wind speed (m s−1). . . . . . . . . . . 65Figure A.8 Wind roses of observed and simulated wind speed at Castle Creek glacierin 2010. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65Figure A.9 Wind roses of observed and simulated wind speed at Conrad glacier in2016, station 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66Figure A.10 Wind roses of observed and simulated wind speed at Conrad glacier in2016, station 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66Figure A.11 Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Nordic glacier in 2014 with correlation(r), p-value (p) and mean bias error (MBE) included for each plot. Thevariables are: T2m (◦C), RH2m (%), U (m s−1), SP (hPa), Kin (W m−2), Lin(W m−2), QH (W m−2), QE (W m−2) and QM (W m−2). . . . . . . . . . . 67Figure A.12 Comparison of AWS-derived (obs) versus WRF-derived (mod) daily time-series of the selected variables at Nordic glacier in 2014. . . . . . . . . . . 68Figure A.13 Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Conrad glacier in 2015, station 1. . . . . . 69Figure A.14 Comparison of AWS-derived (obs) versus WRF-derived (mod) daily time-series of the selected variables at Conrad glacier in 2015, station 1. . . . . . 70Figure A.15 Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Conrad glacier in 2015, station 2. . . . . . 71Figure A.16 Comparison of AWS-derived (obs) versus WRF-derived (mod) daily time-series of the selected variables at Conrad glacier in 2015, station 2. . . . . . 72Figure A.17 Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Conrad glacier in 2016, station 1. . . . . . 73xiiiFigure A.18 Comparison of AWS-derived (obs) versus WRF-derived (mod) daily time-series of the selected variables at Conrad glacier in 2016, station 1. . . . . . 74Figure A.19 Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Conrad glacier in 2016, station 2. . . . . . 75Figure A.20 Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Conrad glacier in 2016, station 2. . . . . . 76Figure A.21 Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Castle Creek glacier in 2010. . . . . . . . 77Figure A.22 Comparison of AWS-derived (obs) versus WRF-derived (mod) daily time-series of the selected variables at Conrad glacier in 2010. . . . . . . . . . . 78xivAcknowledgmentsFunding support for this study was provided through the Natural Sciences and EngineeringResearch Council (NSERC) of Canada. Data for model evaluation is from the field campaignsperformed by Noel Fitzpatrick, Valentina Radic´ and Brian Menounos.My heartfelt appreciation and thanks goes to my supervisor Valentina Radic´ who has beensupportive and encouraging throughout my years of study. Without her tremendous supportand proper guidance, this work would not have been possible.My great thanks goes to Douw Steyn for helping me get this opportunity and his continuoussupport. Gwenn Flowers and Brian Menounos, thank you for your guidance and constructivefeedback.My deepest gratitude goes to my parents, my husband Antonio and son Nahom for provid-ing me unconditional support and encouragement.Special thanks to all my families and friends for checking on me all the time and giving methe spiritual and emotional support that I needed.xv1. Introduction1.1. Motivation and research objectivesMountain glaciers and ice caps play a key role in the Earth’s hydrological cycle. On a globalscale, retreat of glaciers is one of the key indicators of climate change while the worldwideglacial wastage significantly contributes to ongoing sea level rise (Garner et al., 2015). Thecurrent rate of glacier contribution to sea level rise is about 1 mm/yr, which is 30% of theobserved current rate (the other contributors being ocean thermal expansion and the mass lossfrom the Antarctic and Greenland ice sheets). Projections of glacier contribution to sea levelrise until the end of the century range from ∼50 to ∼200 mm (e.g., Marzeion et al., 2012;Radic´ et al., 2014; Huss and Hock, 2015) depending on different climate scenarios and differentmodels used. On regional and local scales, mountain glaciers and ice caps are recognized assignificant contributors to seasonal stream flow, as they store fresh water in the form of snowand ice during winter, that is then released during summer. As an example, contribution ofglacier run-off to water supply for a country like Bolivia can reach up to 14% in the wet seasonand 27% in the dry season (Soruco et al., 2015). Modeling of the magnitude and timing ofglacier run-off is therefore crucial to successfully predict the behavior of melt driven streamflow in order to improve water resource management (Hock, 2005).Since melting at a glacier surface is a dominant process in glacier ablation for all non-calving glaciers, most of the models of glacier ablation focus on simulating glacier surfacemelt. The melt models can be applied to glacier surfaces on local (e.g., Johannesson et al.,1995; Braithwaite and Yu, 2000; Hock, 2003, Woul and Hock, 2005; Radic´ and Hock, 2006;MacDougall et al., 2011), regional and global scales (e.g., Gregory and Oerlemans, 1998;Radic´ et al., 2014; Clarke et al., 2015). Any model of glacier melt, whether empirical orphysics based, requires meteorological variables as input data. A conventional approach toobtain these data is by direct measurements, i.e. by the deployment of an automated weatherstation (AWS) at or in a close vicinity to a glacier. However, meteorological observations thatspan at least one ablation season are available on fewer than 100 glaciers worldwide, mostof which are in the European Alps and Scandinavia (Zemp et al., 2009; Radic and Hock,2014). The remoteness of glaciers, terrain complexity and harsh weather conditions are themain reasons for the lack of long term (several decades) meteorological observations at theglacier-atmosphere interface. A common alternative approach to obtain input data for meltmodels is to use global and/or regional climate reanalysis products. The meteorological fields1from these data-sets are available at relatively coarse spatial resolution (> 30 km) and cannotdirectly be used to force glacier melt models that simulate processes taking place at much finerscales (< 0.1 km). This issue of ‘scale mismatch’ is commonly circumvented by the use ofstatistical and/or dynamical downscaling methods (Lawrence, 1985; Fowler and Wilby, 2007).Figure 1.1: (a) Map of the study domain with the glacier sites. Contour topographicalmap is provided for each glacier, with marked site of automatic weather station(AWS): (b) Castle Creek glacier in 2010 and 2012 with one on-glacier AWS (pinkcircle) and two off-glacier AWSs (pink triangles); (c) Nordic glacier in 2014 withone on-glacier AWS (pink circle); (d) Conrad glacier with two on-glacier AWSsin 2015 (pink circle and triangle), and two on-glacier AWSs in 2016 (blue circleand triangle). Outlines of the glaciers (white lines) are from the Randolph GlacierInventory.The main goal of this study is to investigate the performance of a dynamical downscal-ing approach to derive meteorological fields needed to force a physics-based model of glacier2melt. This downscaling implements the Weather Research Forecasting (WRF) model, whichis a high-resolution regional climate model, forced by a global climate reanalysis data. Themodel is run for a domain covering three mountain glaciers in the interior of British Columbia:Castle Creek glacier, Nordic glacier and Conrad glacier (Figure 1.1). At these glaciers, auto-matic weather stations (AWSs) intermittently recorded data over several summer seasons from2010-2016. Our main objective is to evaluate the performance of the WRF model against theobservations from the AWSs. Specifically, we ask the following questions:1. How do downscaled sub-daily and daily timeseries of meteorological variables comparewith the observations at each site?2. How sensitive are the downscaled fields to the choice of spatial resolution and landcoverdata in the WRF model?3. If a surface energy balance (SEB) model is used to calculate point-scale daily and sea-sonal surface melt at these glaciers, how successfully can the SEB model be forced bythe WRF data instead of the AWS data?The four sections of the study are organized as following: Section 1 provides (i) an overviewof the most widely used melt models for glacier surfaces, (ii) brief review of downscaling meth-ods, focusing on the dynamical downscaling, and (iii) an overview of recent applications of theWRF model at glacierized terrains. Section 2 introduces the data obtained from automaticweather stations (AWS) and explains the setup for the WRF model including the boundaryconditions, topography and landcover input data, and parameterization schemes. This sectionalso introduces the SEB model and a set of sensitivity tests applied to the WRF model and theSEB model. Results and discussion are presented in section 3 and 4, respectively. Conclusionsaddressing the listed research questions are given in section 5.1.2. Overview of glacier melt modelsModels of glacier mass balance which quantify mass fluxes (accumulation and ablation), varyin their complexity, scope, and input data requirements. For example, models for glacier sur-face melt (as part of the glacier mass balance) range from low complexity ones (e.g., tem-perature index model) to high complexity ones (e.g., physics-based surface energy balancemodels). The performance of these melt models is particularly important in glacierized wa-tersheds where glacier melt water production can directly affect the surrounding environment,for example, by triggering seasonal flooding events (Hock, 2005). In this section we provide abrief overview of most commonly used glacier melt models.31.2.1. Simple positive-degree day (PDD) modelThe PDD model is based on an empirical relationship between surface melt and the sum ofpositive daily near-surface air temperature (Braithwaite and Olesen, 1989). Mathematically,the model is described by Cuffey and Paterson (2010) as:D(∆t) =∫∆tTCH(TC)dt with H(TC) =1 if TC > 00 if TC < 0 (1.1)where D is the positive degree day index for time period ∆t measured in days, and TC is near-surface air temperature in degrees Celsius. The D index, often symbolized as PDD in theliterature, approximates the integrated heat energy that drives melt over a time interval of in-terest (∆t). The total surface melt, ms, over time ∆t is parametrized asms = fmD (1.2)where fm (mm day−1 ◦C−1) denotes the degree-day factor, which is an empirically determinedcoefficient, often assumed constant in time. Generally, fm values are larger for an ice surfacethan a snow surface, since the lower albedo of ice can result in greater absorbed energy andmore energy available for melt (Braithwaite and Yu, 2000). Therefore, a common approach isto calculate D separately for the periods of snow cover (Ds) and exposed ice (Di), modifyingthe equation (1.2) toms = fmsDs + fmiDi (1.3)where fms and fmi are melt factors for snow and ice, respectively. The melt factors implicitlyaccount for all components of the surface energy balance. Since these components vary in theirrelative contributions to the melt energy, the melt factors are expected to vary in space and time(Hock, 2005). The melt factors are often found to be glacier specific and can substantially varyfrom point to point on a single glacier, thus their calibration with observed melt is required (e.g.,Johannesson et al., 1995; Hock, 1999; Hock, 2003; Radic´ and Hock, 2006; MacDougall et al.,2011; Gillett and Cullen, 2011; Jonsell et al., 2012; Parajuli et al., 2015). When the modelis calibrated, however, its performance can be equal or better than the performance of morecomplex models, especially at a catchment scale (Rango and Martinec, 1995). Nevertheless,the PDD models are shown not to be suitable for shorter time intervals (hours) and for small4spatial scales, i.e. pixels in satellite images with different exposures to solar radiation (Rangoand Martinec, 1995).1.2.2. Enhanced PDD modelMany attempts have been made to strengthen the physical foundation of the PDD method byincorporating more variables, such as wind speed, vapor pressure or radiation (e.g., Hock,1999). One of these approaches, developed by Hock (1999), introduced direct solar radiationinto the PDD model to capture melt dependency on surface slope, aspect and topographicshading. This ’enhanced’ PDD model is expressed asM =(1n fm +aI)TC if TC > 00 if TC ≤ 0 (1.4)where a is a radiation coefficient different for snow and ice surfaces, n is the number of timesteps per day and I is potential clear sky direct solar radiation at the ice or snow surfaces. Themelt factor ( fm) and the radiation index (a) are empirical coefficients, I is calculated as a func-tion of top of atmosphere solar radiation, assumed atmospheric transmissivity, solar geometryand topographic characteristics. This model is mainly driven by the observed near-surface airtemperature, while the direct solar radiation can be computed from solar geometry and topog-raphy. The model has been applied on several glaciers worldwide (eg., Hock, 1999; Pellicciottiet al., 2005; Shea et al., 2009; Carenzo et al., 2009; MacDougall et al., 2011; Jonsell et al.,2012; Gabbi et al., 2014) as well as on a regional scale i.e. for a full suite of glaciers within agiven region (Clarke et al., 2015).1.2.3. Semi-empirical surface energy balance (SEB) modelAn alternative to the PDD models is to consider a more physics-based approach, where energyavailable for melting is derived from the most relevant fluxes in the surface energy balance(SEB). One example of such an approach is a semi-empirical SEB model developed by Oer-lemans (2001), considers net solar radiation while the longwave radiation and turbulent fluxesare parametrized as a linear function of near-surface air temperature:M = (1−α)G+ c0 + c1Ta (1.5)5where G is global radiation, α is albedo for snow and ice separately, c0 and c1 are calibratedconstants and Ta is near surface air temperature. A study by Giesen and Oerlemans (2012)showed that this model successfully simulates inter-annual variability in seasonal ablationwhen hourly observations of glacier surfaces are used. Similar to the PDD models, however,calibration of model parameters is necessary.1.2.4. Surface energy balance (SEB) modelsWhile the temperature index models, when calibrated, can successfully simulate glacier melt,they belong to a category of ’black box’ models, and thus cannot provide much insight intophysical processes driving glacier melt. Since empirical models require a glacier specific cali-bration with observations, the transferability of model parameters is questionable for sites withno such observations. Unlike empirical models, SEB models when driven by observed fluxesat a glacier surface do not require any calibration and are therefore more transferable from siteto site and potentially applicable on regional (and global) scales. A main drawback, however,is their requirement for more than just one input variable, i.e. air temperature. Some of theseadditional variables (e.g., radiative fluxes) are not necessarily part of standard meteorologicalmeasurements, i.e., the simplest automated weather station installed at a glacier. The SEBquantifies different energy fluxes at the glacier-atmosphere interface: incoming and outgoinglongwave radiation fluxes, incoming and reflected shortwave radiation fluxes, sensible and la-tent heat fluxes, heat flux due to rainfall, heat conduction into the snow/ice, and energy flux putinto melting. To successfully quantify the energy available for melting, each component of theSEB needs to be carefully assessed. On a temperate glacier in mid-latitudes, the net radiationflux is the dominant source of melt energy on annual and seasonal scales while turbulent heatfluxes can dominate over short time intervals of hours and days (Oerlemans and Klok, 2002;Hock, 2005; Hock and Holmgren, 2005; Andreassen et al., 2008; Munro and Marosz-Wantuch,2009; MacDougall and Flowers, 2010; Gillett and Cullen, 2011).The SEB models range in their complexity from those that consider only one surface layer(Ontiveros-Gonzalez et al., 2015; Radic´ et al., 2017) to those that consider multiple (surface andsubsurface) layers (e.g., Klok and Oerlemans, 2002; Oerlemans et al., 2009; MacDougall andFlowers, 2010; Sun et al., 2012). Also, the SEB models can be driven by directly measuredenergy fluxes at a glacier surface (e.g., Conway and Cullen, 2013) or as a combination ofmeasured and calculated (parametrized) fluxes (e.g., Munro, 1989; Oerlemans and Klok, 2002;Klok and Oerlemans, 2002; Hock and Holmgren, 2005; Andreassen et al., 2008; Munro andMarosz-Wantuch, 2009; MacDougall and Flowers, 2010; Gillett and Cullen, 2011).A general SEB model can be expressed in the following form:6QM = Kin +Kout +Lin +Lout +QH +QE +QP +QG (1.6)where QM is the energy available for melting, Kin and Kout are incoming and reflected short-wave radiation, respectively, Lin is downward longwave radiation, Lout is emitted longwaveradiation, QH is sensible heat flux, QE is latent heat flux, QG is subsurface energy flux andQP is heat flux from rain. By convention, the fluxes into the glacier surface are positive andthe fluxes out of the surface are negative and all the fluxes have units of W m−2. In glacierstudies, by convention, QE refers only to latent heat associated with water vapor (i.e. the heatconsumed by evaporation and sublimation and released by condensation and deposition) whilethe latent heat associated with the melting and refreezing is part of QM.If not directly measured, a common approach to derive Lout is by the Stefan-Boltzmannlaw:Lout = εσT 4s (1.7)where ε is emissivity of the surface (assumed to be unity), σ is the Stefan Boltzmann constant(5.67× 10−8 Wm−2K−4 ) and Ts is the glacier surface temperature which, for the meltingsurface, is assumed to be 0◦C.The turbulent heat fluxes are often described by gradient-flux relations (Stull, 2003; Cuf-fey and Paterson, 2010) in which heat transfer by turbulent eddies is viewed as analogous tomolecular diffusion, with eddies playing the role of molecules. Therefore, the vertical flux ofheat is assumed to be proportional to the vertical gradient of mean (the average over a certaintime window) air temperature and to a coefficient specifying the effectiveness of the transferprocess, known as the eddy diffusivity for heat KH :QH = ρacpKH∂T∂ z(1.8)where ρa is air density at standard pressure (ρa = 1.29 kg m−3) and cp is the specific heatcapacity at constant pressure (cp = 1005 J kg−1K−1).The vertical flux of vapor is proportional to the vertical gradient of mean humidity q and acoefficient known as the eddy diffusivity for water vapour KE :7QE = ρaLvKE∂q∂ z(1.9)where Lv = 2.514×106 Jkg−1 is the latent heat of vaporization of snow or ice. Both KH and KEvary with height z and can be parameterized as functions of wind shear, and atmospheric sta-bility. Integrating equations (1.8) and (1.9) between the surface and height z above the surface(e.g., the height of the measurements), and using the mixing length theory to parameterize KHand KE yieldsQH = ρa cp CH uz(Tz−Ts)(1.10)QE =0.622P0ρa LvCE uz (ez− es)where uz is mean wind speed at height z, Tz is mean temperature at height z, Ts is mean sur-face temperature, ez is vapor pressure at z, es is surface vapor pressure, P0 is the air pressureat standard sea level (1013 hPa), P is the actual air pressure and CH and CE are the exchangecoefficients for temperature and humidity, respectively. In equation (1.10), we converted thespecific humidity into vapor pressure using ez = P00.622qz. If relative humidity is measured in-stead of specific humidity, vapor pressure (ez) can be directly calculated from the ClausiusClapeyron equation (WMO, 2008) as:ez =RHz100(6.112)exp(17.62Tz243.12+Tz)(1.11)where RHz is the relative humidity at height z.There is a range of different parametrizations for the exchange coefficients CH and CE (e.g.,Hock, 2005; Conway and Cullen, 2013; Radic´ et al., 2017; Fitzpatrick et al., 2017) but here wesummarize the three most commonly used for simulating turbulent fluxes at glacier surfaces(e.g., Conway and Cullen, 2013; Radic´ et al., 2017; Fitzpatrick et al., 2017):1. Parametrization that assumes a vertical logarithmic profile of wind speed and tempera-ture (valid for statically neutral atmospheric boundary layer):CH log =k2ln( zuz0v ) ln(ztz0t)(1.12)8CE log =k2ln( zuz0v ) ln(zqz0q)(1.13)where k is the von Ka´rma´n constant, zu, zt and zq are measurement heights for windspeed, temperature and water vapor, respectively and z0v, z0t and z0q are roughnesslengths for momentum, temperature and water vapor, respectively.2. Parameterization that includes a stability correction based on the bulk Richardson num-ber (Rib)CH Rib =k2ln( zuz0v ) ln(ztz0t)(1−5Rib)2, for stable conditions (Rib > 0)k2ln( zuz0v ) ln(ztz0t)(1−16Rib)0.75, for unstable conditions (Rib < 0)(1.14)whereRib =g(Tz−Ts)ztTzu2(1.15)and g is acceleration due to gravity. CE Rib has the same form as (1.14) except the zt andz0t are replaced by zq and z0q, respectively.3. Parametrization that includes integrated stability functions (Ψ) based on the Monin-Obukhov similarity theory (Monin and Obukhov, 1954) with stability parameter zLCH MO =k2[ln( zuz0v )−Ψm(zuL )][ln( ztz0t )−Ψh(ztL )] (1.16)where Ψm( zuL ) and Ψh(zL) = Ψq(zL) are the vertically integrated stability functions formomentum and heat respectively. L is the Monin-Obukhov length given asL =ρcpu3∗TzkgQH(1.17)9where u∗ is the friction velocity.Inter-comparison of the above methods on several glaciers showed that all three parameteri-zations simulate daily values of QH and QE reasonably well when roughness lengths, derivedfrom eddy-covariance measurements, are used in the expressions (e.g., Conway and Cullen,2013; Fitzpatrick et al., 2017; Radic´ et al., 2017).However, on sub-daily scales, especially dur-ing prevailing katabatic flows, all methods overestimate QH despite the stability corrections.The rain heat flux (Qp) is commonly modeled as Hock (2005)QP = ρwcwR(Tr−Ts) (1.18)where ρw is the density of water, cw is the specific heat of water, R is the rainfall rate and Tr isthe temperature of rain (assumed to be equal to near surface air temperature).The heat flux within the ice layer of depth h is given byQG =∫ h0ρhcpi∂T∂ tdz (1.19)where∂T∂ tis the rate of change of ice temperature, ρh is the ice density and cpi is the specificheat of ice (2009 J kg−1 K−1). Due to their small contribution to total melt over a melt seasonat temperate glaciers, it is common to ignore both QP and QG (Sicart et al., 2005; Andreassenet al., 2008; Gillett and Cullen, 2011).1.3. Climate downscaling methodsIn the absence of measured meteorological variables at glacier surfaces, a common alterna-tive approach to obtain the input data for melt models is to use global and/or regional climatereanalysis (e.g., Radic´ et al., 2014; Huss and Hock, 2015). Nevertheless, the output fieldsfrom climate reanalysis models have too coarse spatial resolution to adequately represent at-mospheric conditions at a local (glacier) scale. The performance of climate reanalysis modelsis also error prone over complex mountainous terrain mainly because the processes in the atmo-spheric boundary layer are too complex to be adequately parameterized and the observationsneeded for the data assimilation in the reanalysis models are often absent (e.g., Jarosch et al.,2010; Blacutt et al., 2015). The scale mismatch and relatively poor performance of climatereanalysis over a complex terrain have triggered a development of correction schemes to be10applied prior to the implementation of reanalysis data into glacier melt and mass balance mod-els (e.g., Machguth et al., 2009; Kotlarski et al., 2010; Mo¨lg and Kaser, 2011; Jarosch et al.,2010). These correction schemes, more commonly known as climate downscaling methods,can be split into two categories: (1) statistical downscaling - which relies on empirical (statisti-cal) relationships between large scale predictors (e.g., output of a global climate model (GCM)or climate reanalysis model) and the local scale variables observed at a given site, and (2) dy-namical downscaling - which applies a high-resolution regional climate model (RCM) nestedwithin a GCM to directly simulate the local scale variables.Generally, the statistical downscaling is computationally inexpensive and, therefore, usedto project a long-term (e.g., centuries) transient runoff and mass balance changes in responseto climate scenarios from GCMs (e.g., Radic´ and Hock, 2006; Matulla et al., 2009; Springeret al., 2013; Clarke et al., 2015). A critical disadvantage of statistical downscaling, however,is its assumption that the established empirical relationships between large-scale and local-scale variables in the present climate remain unchanged in the future climate, regardless of thechanging climate conditions (e.g., Zorita and Storch, 1999). Also, the method depends on theavailability of local-scale observations which is often a serious limitation for glacierized com-plex terrain. Dynamical downscaling, on the other hand, can explicitly model processes takingplace in the atmosphere and at the surface-atmosphere interface, utilizing physics-based equa-tions that should equally hold in any climate conditions. While observations are still neededfor model evaluation, they are not the key ingredient as is the case for statistical downscaling.To run a state-of-the-art RCM on a sub-kilometer spatial grid, however, is computationallyexpensive (e.g., it takes several days to generate a week-long simulation output on a super-computer) which makes this approach unsuitable for long-term projections. Furthermore, theperformance of dynamical downscaling is sensitive to the choice of RCMs and their parame-terizations (Mearns et al., 2014).1.3.1. Application of dynamical downscaling in glacier studiesApplication of dynamical downscaling in studies of glacier melt and mass balance has found itsuse relatively recently (Machguth et al., 2009; Paul and Kotlarski, 2010; Kotlarski et al., 2010).In one of the first applications (e.g., Paul and Kotlarski, 2010), the dynamical downscaling wasaccompanied by statistical corrections (interpolation methods). More recent studies circum-vented the need for statistical correction by the use of multiple grid nesting (e.g., Mo¨lg andKaser, 2011; Mo¨lg et al., 2012; Claremar et al., 2012; Collier et al., 2013; Collier et al., 2015;Schanke et al., 2015; Wrzesien et al., 2015; Aas et al., 2016; Wu et al., 2016). In this approach,a smaller RCM domain with finer spatial resolution is placed or ‘nested’ inside the original11RCM domain, and the model is set to run simultaneously over original and nested domains.The nesting can be iteratively applied within each domain until a desired spatial resolution isreached over a location of interest. This approach is considered to resolve the mountain glacierscale reasonably well (Mo¨lg and Kaser, 2011; Mo¨lg et al., 2012) and to improve the represen-tation of processes (e.g., orographic precipitation, katabatic winds) taking place over complexmountainous and glacierized terrain (e.g., Maussion et al., 2011).A study of glacier mass balance at Kilimanjaro (Mo¨lg and Kaser, 2011) is the first todemonstrate a successful application of dynamical downscaling at the atmosphere-glacier in-terface. The study applied the Weather Research and Forecast (WRF) model (Powers et al.,2017) to force a SEB model for one month in a dry and a wet season. The WRF model wasrun with a multiple nesting scheme (from 39 km to 0.812 km resolution). while the ERA-Interim climate reanalysis represented boundary conditions. The output model variables (e.g.,near surface air temperature, relative humidity, wind speed, surface air pressure, cloud fractionand precipitation) showed strong correlation with corresponding variables measured at a AWSdeployed at the glacier surface. This successful application prompted the use of WRF modelfor other different glacierized domains. For example, Claremar et al., (2012) applied the WRFmodel, with a nesting scheme of 24 km, 8 km, and 2.7 km to simulate a 2-year surface massbalance for three glaciers in Svalbard. Strong correlations with AWS data were obtained formost downscaled variables except the near-surface wind speed. The increased horizontal reso-lution of the inner-most domain, from 8 km to 2.7 km spatial grid did not improve the modelperformance because the model failed to resolve the key topographic effects (e.g., valley circu-lation, clouds, shading, and albedo). On the other-hand, increased vertical resolution improvedthe model simulations probably due to more correct vertical profiles of wind speed and mois-ture transport from the open sea. Further advancement in the WRF model application in glacierstudies came from Collier et al., (2013) and Collier et al., 2015) who developed a high reso-lution interactive model of the glacier-atmosphere interface and applied it over the Karakoramregion in the Northwestern Himalaya. The model is run for the months of June - August 2004(where the initial 25 days were discarded as a model spin-up time) using three nested domains(33 km, 11 km and 2.2 km spatial grid). The modeled near-surface air temperature and windspeed, over two summer months, agreed well with observations from an AWS in the glaciervicinity, however, the incoming short-wave radiation was overestimated due to poorly simu-lated cloud cover, humidity and topographic shading. Precipitation was also poorly simulated,showing that the resolution of the inner-most nest of 2.2 km failed to resolve processes dictat-ing snow accumulation over a complex terrain (e.g., orographic uplift and microscale complexflow features).122. Data and MethodsThe main goal of our study is to investigate how precisely the WRF model can capture theobserved point- scale meteorological fields and energy fluxes on a glacier surface. To thatend, we aim to compare the downscaled timeseries of these variables and fluxes with theirobserved equivalents. We then simulate point-scale surface melting and ablation using a SEBmodel forced by, separately, the observations and the WRF model output. We also compare theperformance of the SEB model relative to a simple positive degree day (PDD) model when bothmodels are forced by the downscaled variables. In this section we first introduce the data usedfor the WRF model evaluation, i.e. the field measurements from the glacier sites. Secondly,we describe the WRF model setting and briefly introduce the SEB model used to derive themelting and ablation. We also describe the set of evaluation metrics. Finally, we introduce aset of sensitivity tests, categorized as: (1) sensitivity of the downscaled variables to varyingspatial resolution and landcover in the WRF model, and (2) sensitivity of the modeled melt todifferent settings for albedo and surface roughness lengths in the SEB model.2.1. Field measurementsMeteorological and glaciological measurements at three mountain glaciers (Castle Creek, Nordicand Conrad) in the interior of British Columbia (Fig. 1.1) serve to force the SEB model (sec-tion 2.3) and to evaluate the WRF model (section 2.2). Castle Creek Glacier is located in theCariboo Mountains, BC (Figure 1.1) and contributes meltwater to Castle Creek, a tributary ofthe Fraser River. Monitoring of the glaciers annual mass balance has been performed since2009 ((Beedle et al., 2015)). In addition to the on-glacier AWS that operated during the ab-lation seasons of 2010 and 2012, two AWSs in the glacier vicinity (Fig. 1.1) have been inoperation since 2007/2008 (Dery et al., 2010) measuring meteorological variables at two levels(3 and 5 m above surface). The stations were situated near the glaciers transient snowline andin the glacier forefield. At the upper off-glacier AWS, the mean annual air temperature over2007-2010 was -2.6 ◦ C while summer (June-August) mean air temperature was 6.6 ◦ C. Overthe same period, total precipitation during summer (June-August) was 94 mm at the upperoff-glacier AWS, while mean monthly wind speeds often exceed 5 m s−1 (Dery et al., 2010).The lower part of the glacier where the on-glacier AWS is located, is gently sloping with anapproximate mean gradient of 7 ◦. Nordic and Conrad glaciers lie in the Purcell Mountainsin eastern BC and are located within the Columbia river basin. Their seasonal mass balancehas been observed since 2012, but not yet published. Surface area and elevation range of the13Table 2.1: Characteristics of the study glaciers and their automatic weather stations(AWSs)Glacier Area (km2) Elevation range (m) AWS coordinates and elevation (m a.s.l) Observational periodCastle Creek 9.5 1900 - 2800 AWS: 53◦ 3’ 2.99” N, 120◦ 26’ 39.57” W; 1967 01 Aug - 12 Aug 2010AWS: 53◦ 3’ 2.99” N, 120◦ 26’ 39.57” W; 1967 21 Aug - 16 Sep 2012Nordic 5 2000 - 2900 AWS: 51◦ 26’ 3.62” N, 117◦ 41’ 59.02” W; 2280 12 Jul - 28 Aug 2014Conrad 15 1800 - 3200 AWS1: 50◦ 49’ 29.5” N, 116◦ 55’ 20.9” W; 2138 15 Jul - 05 Sep 2015AWS2: 50◦ 49’ 23.0” N, 116◦ 55’ 16.6” W; 2164 16 Jul - 07 Sep 2015AWS1: 50◦ 49’ 22.9” N, 116◦ 55’ 11.7” W; 2163 19 Jun - 28 Aug 2016AWS2: 50◦ 46’ 55.9” N, 116◦ 54’ 43.1” W; 2909 16 Jun - 22 Aug 2016three glaciers are listed in Table 2.1, as well as the location coordinates and operation dates ofthe automated weather stations (AWSs). Observations from all three glaciers consist of incom-ing and outgoing components of shortwave and longwave radiation fluxes, air temperature andhumidity, wind speed and direction, turbulent fluxes, and surface height changes as indicatorof solid precipitation and ablation. At each site, eddy-covariance method was used to derivethe estimates of roughness lengths for momentum, temperature and humidity. At Castle Creekglacier, the measured surface air pressure at the lower and upper off-glacier stations is interpo-lated to derive the pressure at the on-glacier AWS, while rainfall rate is taken from the loweroff-glacier AWS. A melting ice surface was present during all the observations at Castle Creekglacier, with some intermittent fresh snowfall events. At Nordic glacier, a transitional snowsurface was present for the first four days, with partial snow cover diminishing to a fully bareice surface.The details about the data (sensor type, sampling frequency, accuracy control) and atmo-spheric conditions at Castle Creek glacier are given in Radic´ et al., (2017), while the detailsfrom Nordic glacier are in Fitzpatrick et al., (2017). At Conrad glacier, a total of four AWS de-ployments were executed during 2015 and 2016: two stations in the ablation zone from July toSeptember 2015, and one in each of the ablation and accumulation zones from June to August2016 (Fig. 1.1). A melting ice surface was present during observations at both stations in 2015,and for most of the 2016 observation period at the AWS in the ablation zone. At the AWS in theaccumulation zone a snow surface was present throughout the 2016 observation period, and forthe first 10 days at the AWS in the ablation zone. The meteorological sensors were housed ona four-legged quadpod, which provided a stable platform (as monitored by inclinometer) thatlowered as the ice melted, and maintained a nearly constant height of the sensors above the14surface. Near-surface meteorological variables and fluxes were measured at a nearly constantheight above the surface (wind speed at 2.7 m, temperature and relative humidity at 2 m, andradiation variables at 1.6 m). All variables were saved as 1 min averages except for the rainfallwhich was saved as 1 min totals. Table 2.2 provides details on the sensor type and accuracycontrol at Conrad glacier. A time lapse camera in close proximity (10-30 m) to each AWSwas used to observe the surface and atmospheric conditions over a season, and to monitor thestation’s behavior.Table 2.2: Observed variables and their sensor specifications at Conrad glacier. Onlyvariables used for the purpose of this study are listed.Variable Sensor AccuracyWind speed/direction Young 05103ap Wind Monitor ± 0.3 ms−1Air temperature/humidity Rotronic HC2 Probe and Shield (2015) ± 0.1◦C / 0.8 %Aspirated Rotronic HC2 Probe and Shield (2016) ± 0.1◦C / 0.8 %Atmospheric Pressure Vaisala PTB110 Barometer ± 0.3 hPaPrecipitation Texas Electronics Tipping Bucket Gauge ± 1 % (up to 10 mm hr−1)Radiation fluxes Kipp and Zonen CNR4 Radiometer 10 - 20 Wm−2 (pyranometer)5 - 15 Wm−2 (pyrgeometer)Surface height CSI SR50A Sonic Ranger (three sensors per AWS) ± 0.01 m2.2. WRF modelWe run the advanced research WRF model, version 3.8.1 (Powers et al., 2017), configuredwith three nested domains, with the inner-most domain covering the sites of all three studyglaciers. The parent domain (d01) of 22.5 km horizontal grid spacing covers the bulk of NorthAmerica, while the nested domains (d02 and d03) have a horizontal grid spacing of 7.5 kmand 2.5 km, respectively (Figure 2.1 b). Only model output from the two nested domains isused for analysis and henceforth referred to as mod2.5 and mod7.5 output, respectively. ForCastle Creek glacier only, we also run the WRF model configured with three nested domainson 25, 5 and 1 km horizontal grid spacing (Figure 2.1 a). Only model output on 1 km gridspacing is used (referred as mod1.0) in order to test the sensitivity of results to an increase inspatial resolution from the 2.5 km to 1 km. All model runs use the same model specificationsand input data (Table 2.4) following those from (Mo¨lg and Kaser, 2011; Collier et al., 2013;15Collier et al., 2015). As boundary conditions for the parent domains we use ERA-Interimreanalysis (Dee et al., 2011) with 6-hr temporal resolution. ERA-Interm has been commonlyused in other studies (e.g., Mo¨lg and Kaser, 2011; Collier et al., 2013; Collier et al., 2015) andis also shown to perform better than other reanalysis for Western Canada ((?)). Air temperature,wind speed and relative humidity values are taken from 37 pressure levels, in the range from1 hPa (top level) to 1000 hPa (bottom level). Lower initial conditions include soil temperatureand moisture at four levels, surface and sea-level pressure, and near-surface wind, temperatureand humidity. Sea surface temperature (SST) is set to a constant value, i.e., the time varyingoption is turned off. Similarly to other studies that used WRF on glacierized terrains, the one-way nesting approach is applied rather than the two-way nesting. In the one-way nesting, theinformation flows only from the outer (parent) domain to the nested domain, while there is nofeedback from the nest domain solution to the parent. We employ the default lateral boundarycontrol that includes a four-grid-point relaxation zone with the outermost grid-point specified.No nudging or re-initialization is used within the parent and nested domains. All three modeldomains use the same vertical resolution (60 layers up to a model top of 50 hPa) and physicalparameterizations (Table 2.4), with the exception of the cumulusa) b)Figure 2.1: WRF domain set up a) d01 = 25 km, d02 = 5 km and d03 = 1 km (for CastleCreek glacier only) and b) d01 = 22.5 km, d02 = 7.5 km and d03 = 2.5 km (for allstudy sites).convection scheme that is omitted in the inner most domain (2.5 and 1 km). All simulationswere performed using a restart option to prevent the interruption of model runs due to memory16insufficiency. Due to inconsistencies between the physics in the RCMs and those given by theinitial conditions, a model spinup period is required to make sure that all RCM componentsreach a physical equilibrium (Montavez et al., 2017). We, therefore, apply a model spinup of≈3 days for each model run, and the WRF output is saved on 3-hr time step covering the periodof observations with AWSs (Table 2.1). Figures 2.2 - 2.4 show the WRF model topography andlandcover category data for each of the three glaciers in the nested domains of different spatialresolution. To compare the WRF output with the AWS data for each glacier, an output from thegrid cell with the minimum horizontal distance to the AWS geographical coordinates is used.Despite the overlap in the location, the elevation of the selected grid cell differs from the actualAWS elevation (see Table 2.3). Furthermore, none of the selected grid cells that represent theAWS within the ablation area correctly capture the snow/ice land category (Figures 2.2 - 2.4).Instead, the representative grid cells are classified either as the evergreen needle-leaf forest orbare ground tundra. In addition to the selected AWS grid cell for each glacier, we identify agrid cell that is the closest to AWS site and has the snow/ice land category. The WRF outputfrom these grid cells, henceforth labeled as mod∗1.0 and mod∗2.5, is used in the sensitivitytests.Figure 2.2: Topography and land category for the domain covering Castle Creek glacierfrom SRTM 30 m DEM and ESA CCI Landcover at 300 m resolution, respectively,in comparison to the topography and land category in the WRF nested domains atdifferent spatial resolutions (7.5 km and 2.5 km, and 1 km). Categories from theCCI Landcover data are manually converted to the USGS 24-category Land Coverused in the WRF model (see Appendix for the land category legend, Table A.2).Land category 24 (white color in the color bar) is snow/ice. Markers indicate theon-glacier AWS (circle), two off-glacier AWSs (triangles), and the selected grid cellwith ice/snow land category (star) to be used in the sensitivity tests. The outlines ofglaciers (black lines) are from the Randolph Glacier Inventory.17Figure 2.3: Topography and land category for the domain covering Nordic glacier fromSRTM 30 m DEM and ESA CCI Landcover at 300 m resolution, respectively, incomparison to the topography and land category in the WRF nested domains atdifferent spatial resolutions (7.5 km and 2.5 km)Figure 2.4: Topography and land category for the domain covering Conrad glacier fromSRTM 30 m DEM and ESA CCI Landcover at 300 m resolution, respectively, incomparison to the topography and land category in the WRF nested domains atdifferent spatial resolutions (7.5 km and 2.5 km)18Table 2.3: AWSs elevation (units in meters above sea level): actual vs model for all studysites.Glacier Actual mod2.5 mod7.5 mod1.0 mod∗1.0 mod∗2.5Castle Creek 2010/2012 1967 2185 2044 1973 2371 2288Nordic 2014 2208 2255 2021 - - 2325Conrad 2015 stn1 2138 2160 2250 - - 2488Conrad 2015 stn2 2164 2160 2250 - - 2488Conrad 2016 stn1 2163 2160 2245 - - 2488Conrad 2016 stn2 2909 2634 2194 - - 2488Table 2.4: Specifics of the WRF model setup for the three glaciersDomain setupHorizontal grid spacing 22.5, 7.5, 2.5 km (domain 1-3)25, 5, 1 km (domain 1-3) (only for Castle Creek glacier)Time step 90, 30, 10s (domain 1-3)75, 15, 3s (only for Castle Creek glacier)Vertical levels 60Model top pressure 50 hPaModel physicsRadiation RRTMG (Iacono et al., 2008)Cumulus Grell 3D Ensemble (GRELL, 1993; Grell and Devenyi, 2002)Microphysics Thompson 2-moment (Thompson et al., 2008)Atmospheric Boundary Layer Mellor-Yamada (MYNN) Level 3 (Nakanishi and Niino, 2006)Surface layer Mellor-Yamada (MYNN) Level 2.5 (Nakanishi and Niino, 2006)Land surface Noah-MP (Niu et al., 2011; Yang et al., 2011)DynamicsTop Boundary conditions Upper level diffusive layerHorizontal advection 2nd order diffusion on model levelsLateral boundariesSpecified boundary width 5Forcing ERA Interim, 0.75◦×0.75◦ (Dee et al., 2011)Land Cover data ESA CCI (300 m) (ESA, 2018)Topography GMTED2010 (Danielson and Gesch, 2011)192.3. Surface energy balance modelThe selection of a SEB model is motivated by our main goal, i.e. to evaluate the WRF modelin simulating the most relevant components of surface energy balance at the on-glacier AWSlocations. To that end, we chose a SEB model of low complexity and force it first with the AWSdata and then with the WRF output of the selected grid cells. In addition to 3-hr melt rates wealso assess the surface height changes, as a difference between ablation and accumulation. Themodeled surface height changes are compared to those measured by sonic rangers at each studysite. In our study, the available energy flux for melting (QM) at a given point is estimated as:QM = (1−α)Kin +Lin +Lout +QH +QE (2.1)where the terms correspond to, from left to right: broadband albedo, incoming shortwave ra-diation, incoming and outgoing longwave radiation, and turbulent fluxes of sensible and latentheat. We neglect the ground heat flux and heat flux from rainfall, since both gave negligiblecontributions to the total melt over a melting season on Nordic glacier (Fitzpatrick et al., 2017).Since the selected grid cells in the WRF model do not represent the ice/snow land category, αis computed using observed daily Kout and Kin at each site (Figure 2.5).If a glacier surface is at a melting point, a positive QM given in equation (2.1) drives meltand a negative QM drives freezing (if water is available at the surface). The total melt M (inunits of water equivalent thickness) integrated over a given time period is then calculated as:M =QMρwL f, QM > 00, QM ≤ 0(2.2)where ρw is density of water, L f = 3.34 × 105 J kg−1 is the latent heat of melting (or fusion).The turbulent heat fluxes QH and QE in equation (2.1) are parametrized using the bulkmethod given in equation (1.10). Among the three parametrizations for the exchange coef-ficient C, the approach given in 1.16 is chosen since this parametrization is found to be thebest performing for Castle Creek (2010 and 2012) and Nordic glacier (2014) (Fitzpatrick et al.,2017; Radic´ et al., 2017). The roughness lengths (z0v, z0t , z0q) at each AWS site are estimatedfrom the eddy-covariance method (see Radic´ et al., 2017; Fitzpatrick et al., 2017 for details)and assumed constant over the observational period (Table 2.5).20220 230 240 250α00.51Castle 2010 & 2012mod 7.5mod 2.5mod 1.0AWSobserved200 210 220 230α00.51Nordic 2014200 210 220 230 240α00.51Conrad st1 2015200 210 220 230 240α00.51Conrad st2 2015Time (day)180 200 220 240α00.51Conrad st1 2016Time (day)180 200 220α00.51Conrad st2 2016Figure 2.5: Mean daily albedo as observed from AWS (dashed) and as calulated fromthe 3-hr albedo values prescibed for ice, firn and fresh snow to be 0.3, 0.65 and0.8, respectively (solid line). Note that the prescibed albedo values differ due todifferences in observed or modeled snowfall periods.The observed albedo and roughness lengths estimates are our default setting in the SEBmodel forced either by the AWS data or the WRF model output (Fig. 2.5. In addition, weintroduce an alternative setting for albedo where the 3-hr mean albedo for ice, firn and freshsnow is prescribed to 0.3, 0.65 and 0.8, respectively (Fig. 2.5). These values have been usedin the previous studies with WRF in SEB on glaciers (Mo¨lg et al., 2012; Collier et al., 2013;Collier et al., 2015). An alternative setting for surface roughness is z0v = 10−3 m, while thescalar roughness lengths are determined from the surface renewal model of Andreas (1987).Combinations of the two settings for albedo and the two settings for roughness lengths yieldfour different choices of parameter settings to run the SEB model and derive QM: (1) observedalbedo and observed roughness, (2) prescribed albedo and observed roughness, (3) observedalbedo and prescribed roughness, and (4) prescribed albedo and prescribed roughness.21We convert the modeled melt (M) into surface lowering, zM, as:zM = Mρ f/iρw, (2.3)where ρ f/i is ice density (ρi=900 kg m−3) or the density of the snow from the previous winter,where the latter only applies at the AWS in the accumulation zone of Conrad glacier. In theabsence of measured snow density in the accumulation area of Conrad glacier we tune its valueto achieve the best match between measured (via SR50) and modeled zM, using the AWS dataas input to the SEB model, over the whole observational period. This tuning is performed byminimizing the root-mean-square-error (RMSE) between the measured and modeled zM usinga range of rho f from 600 to 900 kg m−3 with 50 m−3 increment. The tuning resulted in ρ f =800kg m−3.To directly compare the net surface height change (znet) with SR50 data, in addition tothe modeled melt, we take into account the modeled sublimation and fresh snow accumula-tion. Prior to the comparison, we applied filtering on the SR50 raw data to remove the high-frequency noise. For the sites with more than one SR50 sensor (all except Castle sites), wecarefully inter-compared the filtered data among the sensors and removed any erroneous dataassociated with the infrastructure titling (usually occurring towards the end of the observationalperiod). As the final SR50 data we use the mean filtered signal across the sensors at the site.Due to the applied filtering, the final SR50 timeseries only intermittently captures diurnal cyclein the surface lowering. The sublimation (zS) is calculated only when the latent heat flux (QE)is negative:zS = |QE |ρ f/iρw Ls, (2.4)where Ls is latent heat of sublimation. Finally, we introduce a simple accumulation model toaccount for any surface rising due to a fresh snowfall, zA (m w.e.), using the precipitation totals(P):zA =Pρsρw, Tz ≤Φ0, Tz >Φ,(2.5)where ρs is a fresh-snow density and Φ, a threshold temperature for differentiating betweenliquid and solid precipitation, is set to 0 ◦C. The net surface height change (znet) is then derivedas znet = zA− zM − zS, where negative (positive) znet stands for net surface lowering (rising).22Table 2.5: Roughness lengths (in m) for momentum (z0v), temperature (z0t) and humidity(z0q) used in this study, tuned fresh-snow densities (in kg m−3), and calibrated site-specific melt factors ( fm) in the PDD model (mm w.e. ◦C−1 day−1)Glacier and year z0v z0t z0q ρs fmCastle Creek 2010 10−2.5 10−4.5 10−4.0 100 7.42Castle Creek 2012 10−2.5 10−4.5 10−4.0 200 5.99Nordic 2014 10−2.5 10−5.0 10−6.0 200 5.29Conrad 2015 AWS1 10−2.5 10−4.0 10−3.5 200 6.80Conrad 2015 AWS2 10−3.0 10−4.5 10−3.5 200 5.69Conrad 2016 AWS1 10−3.0 10−4.5 10−4.5 200 7.63Conrad 2016 AWS2 10−2.5 10−5.0 10−5.0 50 4.09Similarly to the determination of ρ f , the fresh-snow density (ρs) is tuned within a range ofexpected values (i.e. from 10 to 250 kg m−3; Cuffey and Paterson, 2010) to achieve roughlythe best match between the measured (via SR50) and modeled zA during the observed snowfallevents. The tuning is performed by minimizing RMSE between measured and modeled znetusing a range of ρs from 10 to 250 kg m−3 with 10 kg m−3 increment. Final values for thesnow densities at each AWS are listed in Table 2.5. Note that the model can only give a crudeestimate of surface height changes since the mass balance model is oversimplified and we lackany actual measurements of fresh and old snow densities. The comparison of modeled versusmeasured surface height changes, nevertheless, serves as a bulk quality control for the SEBmodel at each AWS site.2.4. Model evaluation and sensitivity testsOur first step in the analysis of model performance is to compare the timeseries of downscaledvariables with the observed variables at each study site. The model evaluation is performedwith the use of standard evaluation metrics: root-mean-square-error (RMSE), mean-bias-error(MBE) and Pearson correlation coefficient (r). These metrics are quantified for the 3-hr anddaily means of the following variables: near-surface air temperature (T2m), relative humidity(RH2m), wind speed (U), surface air pressure (SP), precipitation (P), incoming shortwave ra-diation flux (Kin), incoming longwave radiation flux (Lin), sensible (QH) and latent heat (QE)fluxes as assessed by the bulk method, and melt energy flux (QM) as derived from the SEBmodel (Eq.(2.1)). Note that the WRF model output gives 10-m wind speed, while the AWSwind measurements are from ≈ 2 m height above the surface. While one could apply a wind-speed correction for the height difference based on the logarithmic-wind profile assumption(e.g., Claremar et al., 2012), we chose not to since our goal is to evaluate the model skill with-23out the use of any data correction. Specifically for the wind speed, the logarithmic-wind profilelikely misrepresents the actual wind profiles during the prevailing katabatic conditions, withlow-level wind speed maxima, which are commonly observed at our sites (e.g., Fitzpatricket al., 2017; Radic´ et al., 2017).In addition to the aforementioned evaluation metrics, we analyze the downscaling per-formance relative to the reference data that is dynamically downscaled, i.e the ERA-Interimreanalysis. To that end, all the variables needed for forcing of the SEB model are taken fromthe ERA-Interim output at 0.75◦ x 0.75◦ spatial resolution, i.e. the dataset representing theboundary and initial conditions to the WRF model. Following Wilks 2006 we calculate a skillscore (SSC):SSC = 1− MSEMSEref, (2.6)where MSE is the mean squared error between measurements (AWS) and model (WRF modelat 7.5, 2.5 and 1 km spatial grid) and MSEref between measurements and a reference model(ERA-Interim). If SSC > 0 for a given variable, the model performs better than the referencemodel. Inter-comparison of SSCs for different spatial resolution in WRF can explain whetherthe finer horizontal spacing adds skill to the simulation, i.e. sensitivity of the downscaledvariables to the choice of spatial resolution in the WRF model.Another sensitivity test consists of analyzing the impact of landcover in WRF on the down-scaled fluxes. We therefore look into WRF output from the closest snow/ice grid cell (labeledas mod∗2.5 mod∗1.0) to the grid cell representing the AWS at each site, which is classified as forestor bare tundra (Figures 2.2 to 2.4). This analysis, however, only partially addresses the effectof landcover on model performance since model output is also affected by grid cell elevation.The comparison of model output with measurements also depends on how close the selectedgrid cell is to the AWS location. To isolate the effect of landcover, we manually alter the landcategory of the AWS grid cell to snow/ice category in the model’s innermost domain and thenrerun the WRF model. In this way we can directly compare the output for the same grid cellbut with a different land category. This sensitivity test is only performed for Nordic glacierusing the run with the innermost nest of 2.5 km grid spacing, and the modified AWS grid celloutput is labeled as mod2.5−modified.The next step is to analyze the downscaling performance in the melt modeling over theobserved period for each site. We compare the cumulative melt (zM), as well as the cumulativenet surface changes (znet), derived when the SEB model is forced by the AWS measurementsand then by the WRF model output. We investigate the energy partitioning, where the daily24melt energy flux (QM), averaged over the observational period, is partitioned into the main SEBcomponents: net shortwave radiation flux, net longwave radiation flux, and sensible and latentheat fluxes. We then test the sensitivity of the results to the SEB model setup, using the foursettings for albedo and roughness lengths, i.e. observed versus modeled (set) values.Finally, we compare the performance of the SEB model with a simple positive degree-day(PDD) model when downscaled fields are forcing both models. This analysis is motivatedby the fact that the current models of glacier mass balance on regional and global scales usetemperature-index models (e.g., Marzeion et al., 2012; Radic´ et al., 2014; Huss and Hock,2015) or semi-empirical SEB models (e.g., Giesen and Oerlemans, 2012), all requiring fewerinput variables than the SEB model used in our study. On the other hand, these empirical mod-els require glacier-specific parameters (e.g., melt factors) or a calibration with available massbalance observations. The melt factors in the temperature-index models are most commonlycalibrated using observations of seasonal mass balance or surface height changes (Hock, 2003).While we do have measurements of surface height changes at our sites, these measurementsare post-processed and as mentioned above, do not consistently resolve diurnal fluctuationsin the surface lowering. Furthermore, in order to represent only the measurements of surfacelowering due to ice surface melting, the SR50 measurements would need to be corrected forthe fresh-snow accumulation and melting episodes. Instead of using the manipulated SR50data for the calibration of the PDD model, we use the daily melt rates derived from the SEBmodel when forced with the AWS observations. By linearly regressing the daily SEB melt ratesagainst the daily sum of positive 3-hr observed temperatures, we derive a melt factor value foreach site and observational period (Table 2.5). These melt factors are then used in the PDDmodel forced with the WRF model temperatures at 2.5 km resolution to derive the cumulativemelting over the observational period at each site.253. Results3.1. Simulated versus observed atmospheric conditions andmelt ratesHere we present the results of the inter-comparison among the observed (AWS) variables, thesimulated (WRF) variables from the three nested domains (7.5 km and 2.5 km for all the sites,and 1 km for Castle Creek glacier only), and the simulated variables from ERA-Interim (≈ 80km spatial resolution). The simulated variables are given for the grid cell representing the AWSlocation at each study site. In Figures 3.1 - 3.6, we show the 3-hr timeseries of the selectedvariables for Castle 2010 and 2012 season, Nordic 2014 season and Conrad 2015 season forAWS2 while the figures of the remaining sites are available in Appendix A.1. We also showscatter plots (modeled versus observed values) of 3-hr and daily selected variables for CastleCreek 2012 season (Figures 3.10 and 3.11), while scatter plots for the remaining sites areavailable in the Appendix A.1. Tables 3.1 and 3.2 show the results of this inter-comparison for3-hr and daily timeseries (RMSE, MBE and Pearson’s correlation). The following results areconsistent among all the sites and observational periods:Air temperature and relative humidity: Both ERA-Interim and the WRF model cap-ture well the observed inter-diurnal variability in 2-m air temperature (T2m). Amplitude of thediurnal cycle in temperature is, however, more pronounced in the model runs than in the obser-vations because the grid cell representing the AWS has a different category from snow/ice. Thelocal cooling effect due to the presence of snow/ice at the surface is, therefore, not captured inthe modeled timeseries. In the sensitivity tests, when the neighboring snow/ice covered gridcell is used (mod∗2.5, mod∗1.0) instead of the AWS grid cell, the diurnal cycle of temperatureresembles more closely the observed one (not shown). Most sites display a warm temperaturebias in the WRF model (maximum MBE of 1.6 K for mod2.5 across the site; Tables 3.1 and3.2). Modeled 3-hr and daily timeseries of 2-m air temperature are strongly correlated withobservations (r > 0.7 at all sites; Tables 3.1 and 3.2). Timeseries of relative humidity (RH2m),due to its dependence on temperature, shows similar comparison of modeled versus observedvalues as the timeseries of air temperature. A statistically significant correlation of r > 0.5(p-value < 0.05) is achieved for all sites except for Castle Creek 2012 (Tables 3.1 and 3.2).26T 2m (o C)01020Castle 2010 & 2012RH2m (%)406080100U (m/s)0510SP (hPa)780800820840860Time (day)215 220 225 230 235 240 245 250 255 260P (mm)051015ERAmod7.5mod2.5mod1AWSFigure 3.1: Comparison of observed versus modeled 3-hr meteorological variables atCastle Creek glacier in 2010 and 2012: 2m-air temperature (T2m), 2m-relative hu-midity (RH2m), 10 m (modeled) and 2 m (observed) wind speed (U), surface pres-sure (SP) and accumulated precipitation (P).T 2m (o C)01020RH2m (%)20406080100U (m/s)246810SP (hPa)780800820Time (day)195 200 205 210 215 220 225 230 235 240P (mm)01020ERAmod7.5mod2.5AWSFigure 3.2: Comparison of observed versus modeled 3-hr meteorological variables atNordic glacier in 2014.Surface air pressure: The highest correlations in 3-hr and daily timeseries among allvariables are achieved for the surface air pressure (r > 0.9 at all sites; Tables 3.1 and 3.2),27while the constant bias in air pressure reflects the elevation difference and can, for example, becorrected by using a hydrostatic balance equation. Since the variability in surface air pressure(SP) is associated with the atmospheric circulation at synoptic scales, it is not surprising thatERA-Interim and downscaled timeseries are highly correlated with the observed ones.T 2m (o C)01020RH2m (%)406080100U (m/s)2468SP (hPa)780800820Time (day)200 205 210 215 220 225 230 235 240 245P (mm)01020 ERAmod7.5mod2.5AWSFigure 3.3: Comparison of observed versus modeled 3-hr meteorological variables atConrad glacier in 2015, station 2.Precipitation: The timing and magnitude of precipitation events are generally well cap-tured by both WRF and ERA-Interim. ERA-Interim and WRF (7.5 km) generally overestimatethe frequency of small precipitation events, while the WRF (2.5 km) underestimates precip-itation frequency (not shown). Observed rainfall rates are to be taken with caution because:(1) the rain gauge at the sites might suffer from under-sampling (Fitzpatrick et al., 2017); and(2) the rain gauge does not adequately capture the snowfall rate that contributes to the totalprecipitation rate in the model.Wind speed and direction: The prevailing wind direction for the AWS stations is downglacier as expected for alpine glaciers that experience katabatic flow (e.g., Radic´ et al., 2017).Wind directions from WRF do not resemble the observations (see Figures 3.4, 3.5, 3.6 and alsofigures in the appendix A.1) leading us to a conclusion that the downscaled wind fields can notcorrectly simulate this downslope wind direction.288.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)024681012141618202224Castle 2012: obs Castle 2012: mod2.5Castle 2012: mod*2.5Castle 2012: mod1.0 Castle 2012: mod*1.0Castle 2012: mod7.5windspeed(m/s)Figure 3.4: Wind roses of observed and simulated wind speed at Castle Creek glacier in2012.8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)0123456789101112Nordic 2014: mod2.5Nordic 2014: mod7.5Nordic 2014: mod*2.5Nordic 2014: obswindspeed(m/s)Figure 3.5: Wind roses of observed and simulated wind speed at Nordic glacier in 2014.298.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)01234567891011128.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)Castle 2012: mod7.5Castle 2012:  mod2.5Castle 2012: mod*2.5Castle  2012: obswindspeed(m/s)Figure 3.6: Wind roses of observed and simulated wind speed at Conrad glacier in 2015,station 2.The WRF model, regardless of the spatial resolution, underestimates the wind speed (MBE onaverage being -2 m s−1 for mod2.5 daily wind speed; Table 3.2), with the largest underestima-tion during clear sky conditions, i.e. during the prevailing katabatic flow. This underestimationhappens despite the height differences (10 m height in the model versus 2 m height at AWS)which theoretically should lead to an overestimation in the model rather than the underesti-mation. The correlations between modeled and observed wind speeds at 3-hr and daily scalesare low (r < 0.3) and for most sites not statistically significant at 5% confidence level (Tables3.1 and 3.2). On the other hand, a closer resemblance between modeled and observed dailywind speed occurs during storm events associated with the surface air pressure drop and pre-cipitation. Among all the sites, the one in the accumulation area at Conrad glacier in 2016and the one at Castle Creek glacier in 2012 display statistically significant correlation betweenobserved and modeled (at 2.5 km) timeseries of daily wind speeds (r > 0.5; p-value < 0.05; Ta-bles 3.1 and 3.2). Note that, for Castle 2012, the finest resolution model run (mod1.0) simulateshigher wind speeds than other model runs and, in part, overestimates the observed wind speed.These higher wind speeds, however, are not associated with the downslope wind direction (seeFigure 3.4).Incoming radiation fluxes: For the observed clear sky conditions, the modeled incoming30shortwave (Kin) and longwave (Lin) radiation fluxes closely resemble the observed fluxes. TheWRF model, however, consistently overestimates the frequency of clear sky days in the ob-servational period, which results in the overestimation of seasonal Kin, and underestimation ofseasonal Lin. High correlation for 3-hr timeseries in Kin (r > 0.75; Table 3.1) is mainly due tothe well captured daily cycle. For the daily timeseries, however, the correlation drops (0.3 <r < 0.68; Table 3.2) as the overestimation of occurrences of clear sky days becomes apparent(positive MBE across the sites). If the storm events are well captured in the model, modeledKin and Lin closely resemble the observed fluxes. Generation of local convective clouds, on theother hand, seems to represent the main challenge for the model to successfully capture. Thispattern of model performance yields an overestimated mean Kin (MBE for mod2.5 in the rangefrom 14 to 68 W m−2 across the sites; Tables 3.1 and 3.2) and an underestimated mean Linover the observational period (MBE for mod2.5 in the range from -29 to -10 W m−2 across thesites; Tables 3.1 and 3.2). The biases are the most pronounced in the model run with the finestresolution (mod2.5 and mod1.0).Kin (W/m2 )0200400600800Castle 2010 & 2012L in (W/m2 )200250300350QH (W/m2 )-1000100QE (W/m2 )-100-50050Time (day)215 220 225 230 235 240 245 250 255 260QM (W/m2 )0200400600800ERAmod7.5mod2.5mod1AWSFigure 3.7: Comparison of observed or AWS-derived versus WRF-derived 3-hr SEBfluxes at Castle Creek glacier in 2010 and 2012: incoming shortwave radiation(Kin), incoming longwave radiation (Lin), sensible heat flux (QH), latent heat flux(QE)) and melt energy flux (QM).31Kin (W/m2 )0200400600800L in (W/m2 )250300350QH (W/m2 )0100200QE (W/m2 )050100150Time (day)195 200 205 210 215 220 225 230 235 240QM (W/m2 )0200400600800ERAmod7.5mod2.5AWSFigure 3.8: Comparison of observed or AWS-derived versus WRF-derived 3-hr SEBfluxes at Nordic glacier in 2014.Kin (W/m2 )0200400600800L in (W/m2 )250300350QH (W/m2 )050100150QE (W/m2 )050100150Time (day)200 205 210 215 220 225 230 235 240 245QM (W/m2 )0200400600800ERAmod7.5mod2.5AWSFigure 3.9: Comparison of observed or AWS-derived versus WRF-derived 3-hr SEBfluxes at Conrad glacier in 2015, station 2.Turbulent fluxes: The simulation of turbulent fluxes, calculated by the bulk method, de-pends on how the model simulates wind speed, 2-m air temperature (for QH) and 2-m specifichumidity (for QE). Since the wind speed is mainly underestimated (MBE for mod2.5 is from32-1.1 to -2.3 m s−1; Tables 3.1 and 3.2) and poorly correlated with observed values, the turbu-lent fluxes are also underestimated and poorly correlated with AWS-derived fluxes. For mod2.5run, there is a consistent mean negative bias in QH for all the AWS sites (MBE from -12 to 42W m−2; Tables 3.1 and 3.2), while the correlation for daily timeseries spans from statisticallyinsignificant negative correlations to statistically significant correlations of r = 0.77 and r = 0.76(Table 3.2) for Castle 2012 and Conrad 2016 stations in the accumulation area, respectively.The underestimation of wind speed is likely the key reason for the underestimation of QH sincethe mean bias in 2-m air temperatures is close to zero or slightly positive (within 2◦C) acrossthe study sites.T2m (° C) obs-5 0 5 10T 2m (°  C) mod2.5-50510r=0.73 p=0.00 MBE=-0.8RH2m (%) obs60 80 100RH2m (%) mod2.56080100r=0.07 p=0.29 MBE=4.2U (m/s) obs2 4 6 8 10U (m/s) mod2.5246810r=0.27 p=0.00 MBE=-1.1SP (hPa) obs780 790 800SP (hPa) mod2.5780790800r=0.92 p=0.00 MBE=-21.3Kin (W/m2) obs0 200 400 600 800K in (W/m2 ) mod 2.50200400600800r=0.79 p=0.00 MBE=20.4Lin (W/m2) obs250 300 350L in (W/m2 ) mod 2.5250300350r=0.27 p=0.00 MBE=-28.8QH (W/m2) obs-50 0 50QH (W/m2 ) mod 2.5-50050r=0.59 p=0.00 MBE=-24.1QE (W/m2) obs-50 0 50QE (W/m2 ) mod 2.5-50050r=0.31 p=0.00 MBE=-4.5QMelt (W/m2) obs0 200 400QMelt (W/m2 ) mod 2.50200400r=0.85 p=0.00 MBE=-22.1Figure 3.10: Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Castle Creek glacier in 2012 season with corre-lation (r), p-value (p) and mean bias error (MBE) included for each plot.33T2m (° C) daily obs-4 0 4 8T2m (°  C) daily mod2.5-4048r=0.77 p=0.00 MBE=-0.8RH2m (%) daily obs70 80 90RH2m (%) daily mod2.5708090r=0.03 p=0.88 MBE=4.2U (m/s) daily obs2 4 6U (m/s) daily mod2.5246r=0.56 p=0.00 MBE=-1.1SP (hPa) daily obs780 790 800SP (hPa) daily mod2.5780790800r=0.92 p=0.00 MBE=-21.3Kin (W/m2) daily obs50 100 150 200 250Kin (W/m2 ) daily mod2.550100150200250r=0.38 p=0.06 MBE=20.4Lin (W/m2) daily obs240 260 280 300 320Lin (W/m2 ) daily mod2.5240260280300320r=0.37 p=0.07 MBE=-28.8QH (W/m2) daily obs-40 -20 0 20 40QH (W/m2 ) daily mod2.5-40-2002040r=0.77 p=0.00 MBE=-24.0QE (W/m2) daily obs-40 -20 0 20QE (W/m2 ) daily mod2.5-40-20020r=0.45 p=0.02 MBE=-4.4QMelt (W/m2) daily obs0 100 200QMelt (W/m2 ) daily mod2.5050100150200r=0.81 p=0.00 MBE=-22.0Figure 3.11: Comparison of AWS-derived (obs) versus WRF-derived (mod) daily time-series of the selected variables at Castle Creek glacier in 2012 season with corre-lation (r), p-value (p) and mean bias error (MBE) included for each plot.Melt rates: All sites show statistically significant correlations, at 5% confidence level,between AWS-derived and WRF-derived daily timeseries of energy available for melt (QM;0.60 < r < 0.89 for mod2.5; Tables 3.1 and 3.2). The results here are for the observed albedoand observed roughness length schemes in the SEB model, but the findings are the same forthe other three schemes where the albedo and z0 are prescribed. There is no consistent patternin the MBE across the sites, i.e. for some sites the model slightly overestimates QM, while forother sites there is an underestimation (see Tables 3.1 and 3.2). To further analyze the modelperformance for each site we plot the partitioning of QM, averaged over the observationalperiod, into the net shortwave (Knet) and net longwave (Lnet) radiative fluxes, and the turbulentheat fluxes (QH and QE) (Figure 3.12). At all sites, observed Knet is the dominant contributorto the observed melt energy, followed by QH . Thus the model skill in representing the net34shortwave fluxes dominates the model performance in simulating the melt rates.QM Knet Lnet QH QEenergy flux (W/m2 )-50050100150200250 Castle 2010ERAmod 7.5mod 2.5mod 1.0AWSQM Knet Lnet QH QEenergy flux (W/m2 )-50050100150 Castle 2012QM Knet Lnet QH QE-50050100150200 Nordic 2014QM Knet Lnet QH QE-50050100150200250 Conrad st1 2015QM Knet Lnet QH QE-50050100150200 Conrad st2 2015QM Knet Lnet QH QE-50050100150200250 Conrad st1 2016QM Knet Lnet QH QE-40-20020406080100 Conrad st2 2016Figure 3.12: Comparison of mean daily melt energy flux and its partitioning into netshortwave radiation (Knet), net longwave radiation (Lnet), sensible heat flux(QH)and latent heat fluxes (QE) when observed albedo and roughness lengths are usedfor each site.The WRF model at 2.5 km resolution, as already mentioned, overestimates the clear sky con-ditions, and therefore yields overestimated mean daily Knet (up to 30% of observed Knet) andunderestimated Lnet (down to 60% of observed Lnet). Together with the underestimation ofQH (down to 80% of AWS-derived QH), these biases in the most relevant SEB componentsact to compensate each other, resulting in a relatively successful simulation of the mean dailymelt rates (difference within 10% of AWS-derived QM) over the observational period. Thus,despite the poorest performance of mod2.5 run in simulating Kin and Lin relative to mod7.5 andERA-Interim, mod2.5 gives the closest QM to the AWS-derived values due to the compensationof biases in the SEB model. Overall, ERA-Interim displays the best skill in simulating Knet andQH , and therefore successfully simulates QM. Only for the AWS at the accumulation site ofConrad glacier in 2016, ERA-Interim substantially overestimates QM (up to 40% of the AWS-derived QM) due to overestimating the 2-m air temperature which leads to an overestimationof QH . For this site, mod2.5 best simulates QM, again due to the cancellation of biases in theindividual components of the SEB. The same overall pattern in the SEB partitioning is foundfor other schemes of albedo and surface roughness lengths (Figure 3.13).35QM Knet Lnet QH QEenergy flux (W/m2 )-50050100150200 Castle 2010ERAmod 7.5mod 2.5mod 1.0AWSQM Knet Lnet QH QEenergy flux (W/m2 )-50050100150 Castle 2012QM Knet Lnet QH QE-50050100150200250 Nordic 2014QM Knet Lnet QH QE-50050100150200 Conrad st1 2015QM Knet Lnet QH QE-50050100150200 Conrad st2 2015QM Knet Lnet QH QE-50050100150200250 Conrad st1 2016QM Knet Lnet QH QE-50050100150 Conrad st2 2016Figure 3.13: Comparison of mean daily melt energy flux and its partitioning into netshortwave radiation (Knet), net longwave radiation (Lnet), sensible heat flux(QH)and latent heat fluxes (QE) when modeled albedo and roughness lengths.3.1.1. Towards distributed melt modelingConsidering that surface melting is spatially variable, a natural extension of our analysis isto investigate how well the SEB model, driven by the WRF model output, can resolve thisspatial variability. At Conrad glacier, two stations operated simultaneously during 2016, onestation in the ablation (station 1) and other in the accumulation area (station 2), each onerepresented by a different grid cell in mod2.5 (Figure 2.4). For each site we plot the timeseriesof daily mean values of the selected variables, observed versus mod2.5 (Figure 3.14). The SEBmodel here is used with the observed albedo and observed roughness length scheme. For eachvariable we also plot the difference in the observed (modeled) timeseries between the upperand lower station (Figure 3.14). Note that for both the 2-m temperature and QM, the differenceis converted to a gradient, i.e. difference between the stations is divided by their elevationdifference. The results reveal smaller difference in Kin between upper and lower stations inthe observations than the model simulations. Observed time series of daily mean Lin show astronger correlation between the two stations (r = 0.93), than is the case in mod2.5 (r = 0.60).36Figure 3.14: Left panels: comparison of AWS-derived (obs) versus WRF-derived(mod2.5) daily timeseries of the selected variables from the two stations at Conradglacier in 2016: station 1 in the ablation area, and station 2 in the accumulationarea. Right panels: the difference in each variable between the two stations (dailytimeseries of station 2 minus daily timeseries of station 1), as AWS-derived andWRF-derived.As measured at the stations throughout the observational period, the upper station consistentlyreceives less longwave radiation than the lower station, while this difference is only intermit-tently captured in mod2.5. The temperature gradient in observations is consistently smaller andmore variable (mean ± standard deviation: -0.43 ± 0.18 (100 m)−1) than in mod2.5 (-0.89 ±0.08 (100 m)−1). The observed temperature gradient co-varies with the wind speed measuredat the lower station (r = 0.79; not shown), a pattern indicative of advective cooling during akatabatic flow. The pattern is not captured in mod2.5. About 90% of the time, the upper stationreceives higher precipitation rates than the lower station, however, the precipitation increasewith elevation is replicated for less than 30% of the time in mod2.5. In addition, modeled andobserved precipitation are more correlated at the lower station (r = 0.64) than at the upperstation (r = 0.24). Results for the melt rates show that QH is the main contributor to the inter-diurnal variability of QM. At the lower station, the correlation between the AWS-derived QMand QH (T2m) is 0.68 (0.74), and at the upper station: 0.77 (0.80). The equivalent correlation inmod2.5 at the lower station is 0.22 (and 0.63), while at the upper station is r = 0.52 (r = 0.85).37Table 3.1: Evaluation results for mod2.5 against the AWS variables: correlation, RMSEand MBE derived from 3-hr timeseries. Bold values are correlations significantlydifferent from zero at a significance level of 0.05.Castle Creek Nordic Conrad 2015 Conrad 20162010 2012 2014 stn1 stn2 stn1 stn2Variables Units Correlation (r)T2m 0.79 0.73 0.85 0.79 0.80 0.69 0.84RH2m 0.63 0.07 0.55 0.51 0.51 0.55 0.63U -0.21 0.27 0.00 0.05 0.03 -0.05 0.27SP 0.94 0.92 0.97 0.96 0.96 0.93 0.94Kin 0.83 0.79 0.86 0.86 0.86 0.88 0.86Lin 0.32 0.27 0.41 0.43 0.42 0.27 0.29QH -0.15 0.59 0.03 0.21 0.20 0.02 0.48QE -0.06 0.31 0.36 0.41 0.41 0.23 0.58QM 0.86 0.85 0.86 0.88 0.89 0.89 0.93RMSET2m ◦C 2.6 2.9 2.5 2.6 2.6 2.6 1.9RH2m % 14.1 22.3 15.2 16.7 17.8 17.8 12.5U m s−1 2.3 2.9 2.1 2.0 2.2 2.2 2.4SP hPa 4.9 19.1 21.4 2.4 1.5 1.5 22.8Kin W m−2 162.0 161.0 159.0 169.6 171.5 171.5 174.3Lin W m−2 43.6 46.8 35.5 37.5 35.2 35.2 38.9QH W m−2 54.5 38.3 55.8 51.0 41.0 41.0 23.2QE W m−2 18.7 20.6 15.8 18.1 15.0 15.0 11.6QM W m−2 120.1 86.1 103.0 116.1 100.9 100.9 31.7MBET2m ◦C 0.1 -0.8 1.4 -0.2 -0.4 1.6 0.4RH2m % 8.8 4.2 5.9 9.4 11.7 -0.6 2.3U m s−1 -1.3 -1.1 -1.3 -1.2 -1.4 -1.8 -1.9SP hPa -19.1 -21.3 -4.8 -2.1 -1.0 0.7 22.8Kin W m−2 27.6 20.4 37.9 65.8 67.3 53.8 38.8Lin W m−2 -23.1 -28.8 -18.9 -26.2 -20.6 -24.1 -9.9QH W m−2 -42.9 -24.1 -35.4 -31.9 -27.1 -26.4 -12.5QE W m−2 -12.8 -4.4 -1.2 1.0 2.6 -5.5 0.8QM W m−2 -45.1 -22.1 -22.4 1.6 9.2 -8.0 -7.5The better model skill in downscaling wind speeds and temperature at the upper station yields ahigher correlation between WRF-derived and AWS-derived QM (r = 0.88) relative to the lowerstation (r = 0.73). The variability in the AWS-derived melt gradient (∆QM/∆z) is mainly drivenby temperature variability observed at the upper station (r = 0.64).The variability in modeled ∆QM/∆z, however, is mainly driven by the variability in modeledKin at the lower station. Modeled mean ∆QM/∆z over the observational period is -0.27 Wm−2 m−1 which is larger than the AWS-derived value of -0.18 W m−2 m−1. Larger gradientin melting is explained by overestimated Kin in mod2.5: the difference between WRF-derivedand AWS-derived ∆QM/∆z is strongly correlated with the difference between modeled andobserved Kin at the lower station (r = -0.89).38Table 3.2: Evaluation results for mod2.5 against the AWS variables: correlation, RMSEand MBE derived from daily timeseries. Bold values are correlations significantlydifferent from zero at a significance level of 0.05.Castle Creek Nordic Conrad 2015 Conrad 20162010 2012 2014 stn1 stn2 stn1 stn2Variables Units Correlation (r)T2m 0.92 0.77 0.94 0.91 0.90 0.81 0.92RH2m 0.89 0.03 0.69 0.62 0.63 0.70 0.76U -0.38 0.56 -0.06 0.10 0.05 -0.21 0.59SP 0.95 0.92 0.98 0.98 0.98 0.94 0.95Kin 0.60 0.38 0.68 0.63 0.61 0.46 0.32Lin 0.49 0.37 0.50 0.57 0.53 0.27 0.36QH -0.12 0.77 -0.02 0.31 0.31 -0.05 0.76QE -0.01 0.45 0.48 0.51 0.51 0.26 0.73QM 0.61 0.81 0.66 0.75 0.76 0.73 0.89RMSET2m ◦C 1.8 2.5 1.8 1.5 1.6 2.4 1.3RH2m % 10.7 19.6 12.5 13.9 15.3 6.5 8.3U m s−1 2.6 1.5 1.9 1.7 1.9 2.2 2.1SP hPa 19.2 21.4 4.9 2.3 1.4 1.6 22.8Kin W m−2 80.7 77.6 75.3 93.4 95.3 88.2 94.6Lin W m−2 31.4 38.9 27.5 31.5 28.2 29.8 27.7QH W m−2 49.0 32.7 49.7 44.8 36.7 34.4 18.9QE W m−2 15.8 14.6 11.2 12.9 10.8 9.7 9.4QM W m−2 81.0 45.5 56.1 59.4 50.8 47.5 17.6MBET2m ◦C -0.1 -0.8 1.4 -0.2 -0.4 1.6 0.4RH2m % 9.4 4.2 5.8 9.4 11.7 -0.6 2.3U m s−1 -2.3 -1.1 -1.3 -1.2 -1.4 -1.8 -1.9SP hPa -19.1 -21.3 -4.8 -2.1 -1.0 0.7 22.8Kin W m−2 13.7 20.4 37.9 65.8 67.3 53.8 38.8Lin W m−2 -24.5 -28.8 -18.9 -26.2 -20.6 -24.1 -9.9QH W m−2 -42.3 -23.9 -35.4 -31.9 -27.1 -26.4 -12.4QE W m−2 -12.2 -4.4 -1.2 1.0 2.6 -5.5 0.7QM W m−2 -55.3 -22.0 -22.4 1.6 9.2 -8.0 -7.53.2. Cumulative ablation and surface height changesFigure 3.15 shows the modeled cumulative melt, converted to surface lowering (zM), over theobservational period at each site. For all the sites except Castle, the net melt derived whenthe SEB model is forced with mod2.5 is within 10% of the net melt derived when the SEBmodel is forced with the AWS observations. We place less confidence in the evaluation resultsat the Castle site, since the sampling period is shorter than at other sites, and the sensorsmeasuring radiative fluxes are older and less precise than at other sites. As shown earlier inthe evaluation of QM, ERA-Interim yields the closest simulation of the cumulative melt (with10% difference with observed melt) at each site except the accumulation site at Conrad 2016(up 40% difference). Similar results are found when the modeled albedo and modeled surface39roughness lengths are used instead of the observed values (see Figure 3.16).zM (m)-0.8-0.6-0.4-0.200.2Castle 2010 & 2012Time (day)220 240 260znet (m)-0.8-0.6-0.4-0.200.2ERAmod 7.5mod 2.5mod 1.0AWSSR50-2.5-2-1.5-1-0.50Nordic 2014Time (day)200 220 240-2.5-2-1.5-1-0.50-2.5-2-1.5-1-0.50Conrad st1 2015Time (day)200 220 240-2.5-2-1.5-1-0.50-2.5-2-1.5-1-0.50Conrad st2 2015Time (day)200 220 240-2.5-2-1.5-1-0.50-4-3-2-10Conrad st1 2016Time (day)180 200 220 240-4-3-2-10-2-1.5-1-0.50Conrad st2 2016Time (day)180 200 220-2-1.5-1-0.50Figure 3.15: Comparison of WRF-derived versus AWS-derived melt, converted to sur-face lowering (zM; set to negative values here), and net surface height changes(znet) at each site over the observational period using observed albedo and rough-ness length.zM (m)-0.8-0.6-0.4-0.200.2Castle 2010 & 2012Time (day)220 240 260znet (m)-0.8-0.6-0.4-0.200.2ERAmod 7.5mod 2.5mod 1.0AWSSR50-2.5-2-1.5-1-0.50Nordic 2014Time (day)200 220 240-2.5-2-1.5-1-0.50-3-2.5-2-1.5-1-0.50Conrad st1 2015Time (day)200 220 240-3-2.5-2-1.5-1-0.50-2.5-2-1.5-1-0.50Conrad st2 2015Time (day)200 220 240-2.5-2-1.5-1-0.50-4-3.5-3-2.5-2-1.5-1-0.50Conrad st1 2016Time (day)180 200 220 240-4-3.5-3-2.5-2-1.5-1-0.50-2-1.5-1-0.50Conrad st2 2016Time (day)180 200 220-2-1.5-1-0.50Figure 3.16: Comparison of WRF-derived versus AWS-derived melt, converted to sur-face lowering (zM; set to negative values here), and net surface height changes(znet) at each site over the observational period using prescribed albedo and rough-ness length.Next we look at the net surface height changes (znet) that, in addition to melting, take into40account the sublimation and fresh-snow accumulation. We compared the modeled znet to thesurface height changes measured by sonic rangers (SR50) at each site (Figure 3.15). Sincethe sublimation contributes less than 1% to the net surface changes, the largest difference inthe modeled versus measured znet comes from the fresh-snow accumulation, as is displayed atsites on Castle 2012 and at Conrad 2016 in the accumulation zone. Surface height changesfrom mod2.5 yield a close resemblance (difference within 10%) to those modeled by the AWSobservations and those measured by SR50. In fact, at the accumulation site of Conrad 2016,the WRF model (mod2.5) gives closer resemblance to measured znet than the AWS observationsdo. Using the SEB model settings with prescribed albedo (prescribed values for fresh and oldsnow, and ice) the albedo is altered from its initial value to the albedo value for fresh snow(0.8) only during the 3-hr period with the snowfall detected at AWS or modeled in WRF. ForCastle 2012, using the prescribed albedo in mod2.5 improves the skill in simulating znet, whileat Conrad 2016 it degrades the skill. At both sites, however, mod2.5 performs better than ERA-Interim and mod7.5 which both fail to correctly capture the snowfall events throughout theobservational periods. Since the variance in zM drives the bulk of the variance in znet, thereis a close resemblance between WRF-derived and measured (SR50) znet for mod2.5 regardlessof how well the model simulates the accumulation. Nevertheless, we show here that correctlydownscaling temperature and precipitation plays an important role in zM and znet modelingthrough the albedo feedback.3.3. Sensitivity tests3.3.1. Sensitivity to spatial resolution in the WRF modelHere we analyze the downscaling performance, in the simulations of daily variables, relativeto the reference ERA-Interim simulations. In particular, we analyze whether the increase ofspatial resolution (from 7.5 to 2.5 and 1 km) adds skill to the WRF model simulation. Theskill scores (SSC), calculated for the variables used in the SEB model, show that the down-scaling only improves (SSC > 0) in the simulation of 2-m air temperature (Figure 3.17). Inter-comparison of SSCs for model runs of different spatial resolution reveals a relatively complexpattern: for the radiative fluxes (Kin and Lin) the model run with a finer resolution consistentlyyields poorer skill across all the sites.41SSC00.20.40.60.81T 2m-10-50RH-2-101USSC-5-4-3-2-101Prec-3-2-101K in-8-6-4-20LinSitesCa10 Ca12 N14 C15 1 C15 2 C16 1 C16 2SSC-1.5-1-0.500.51QHSitesCa10 Ca12 N14 C15 1 C15 2 C16 1 C16 2-1.5-1-0.500.51QESitesCa10 Ca12 N14 C15 1 C15 2 C16 1 C16 2-1.5-1-0.500.51QMmod 7.5mod 2.5mod 1.0Figure 3.17: Skill Scores (SSC) calculated for daily timeseries of the selected variablesat each study site.For example, precipitation is better simulated in mod2.5 than in mod7.5 for all the sites exceptCastle 2012, while for wind speed only half of the sites are better simulated in mod2.5. Theeffect of refining the model resolution from 2.5 to 1 km, as tested on the Castle site in 2010and 2012, does not improve the skill for any variable. While one might expect a katabaticflow to be better resolved at 1 km resolution, this is certainly not the case for our sites wherethe SSC shows more negative values for mod1.0 than for mod7.5 run. Finally, at most sites,mod2.5 has a better skill than mod7.5 in simulating QM. As previously shown, the success insimulating QM depends on the bias compensations in the components of the SEB. The effectof bias cancellation works the best in the mod2.5 run.3.3.2. Sensitivity to land coverTo analyze the effect of landcover on the model performance, we look into the model outputof the grid cell with snow/ice landcover, labeled as mod*2.5 and mod*1.0, that is spatially theclosest to the grid cell representing the AWS. In addition, for Nordic site in 2014, we extractthe model output at 2.5 km resolution where the landcover category of the AWS grid cell wasmanually altered from the bare tundra to snow/ice (see figure 3.18). By doing so we altered the42grid cell’s albedo that the WRF model ‘sees’ from 0.18 (bare tundra) to 0.67 (snow/ice). Thealbedo in the SEB model, however, remains unaltered, i.e. set to the observed daily albedo atthe AWS. The results of these runs, in terms of the mean daily melt energy and its partitioning,reveal that the sensitivity to the landcover alteration, as tested for the Nordic site, is negligible,i.e. the new run produced a small increase (up to 1%) in each energy flux (QM, Kin, Lin andturbulent heat fluxes).QM Knet Lnet QH QEenergy flux (W/m2 )-50050100150200Castle 2010QM Knet Lnet QH QEenergy flux (W/m2 )-50050100150Castle 2012AWSmod2.5mod*2.5mod1.0mod*1.0mod2.5-modifiedQM Knet Lnet QH QE-50050100150200250Nordic 2014QM Knet Lnet QH QE-50050100150200Conrad st1 2015QM Knet Lnet QH QE-50050100150200Conrad st2 2015QM Knet Lnet QH QE-50050100150200250Conrad st1 2016QM Knet Lnet QH QE-40-20020406080100120Conrad st2 2016Figure 3.18: Results of the sensitivity test to varying landcover in WRF: comparison ofmean daily melt energy flux and its partitioning into net shortwave radiation (Knet),net longwave radiation (Lnet), sensible heat flux(QH) and latent heat fluxes (QE)when modeled albedo and roughness lengths are used in the SEB model.The largest difference in the new run relative to the original run is in the 2-m air tempera-ture (not shown), where the diurnal cycle is decreased relative to the original run and betterresembles the one observed at the AWS. The same effect on the diurnal cycle is observed inthe mod*2.5 and mod*1.0 at all the sites. On the other hand, mod*2.5 and mod*1.0 yield largerdifferences in the energy fluxes from the original run at each site (on average 10% differenceacross the sites) than is the case in the sensitivity test on Nordic. This finding is not surpris-ing considering that the model output is being compared between the grid cells of differentelevations and proximity to the AWS site, in addition to the landcover difference.433.3.3. Sensitivity to parameters in the SEB modelHere we analyze the sensitivity of modeled QM and its partitioning to the four SEB model’sschemes for albedo and roughness lengths (see figure 3.19): (1) observed albedo and observedroughness, (2) set albedo and observed roughness, (3) observed albedo and set roughness, and(4) set albedo and set roughness. Results reveal that switching the albedo schemes has moreeffect on altering QM than switching between the roughness length schemes, as expected dueto the dominant effect of Knet on available energy for melting. As shown earlier, the largestsensitivity of daily mean QM to the albedo setting is at the sites with observed fresh snowfallduring the observational period.QM Knet Lnet QH QEenergy flux (W/m2 )-50050100150200Castle 2010αobs-z0obsαmod-z0modαmod-z0obsαobs-z0modQM Knet Lnet QH QEenergy flux (W/m2 )-50050100150Castle 2012QM Knet Lnet QH QE-50050100150200Nordic 2014QM Knet Lnet QH QE-50050100150200250Conrad st1 2015QM Knet Lnet QH QE-50050100150200Conrad st2 2015QM Knet Lnet QH QE-50050100150200250Conrad st1 2016QM Knet Lnet QH QE-40-20020406080100Conrad st2 2016Figure 3.19: Results of the sensitivity test to the choice of four parameterization schemesfor albedo (α) and roughness length (z0) in the SEB model forced by mod2.5:comparison of mean daily melt energy flux and its partitioning into net shortwaveradiation (Knet), net longwave radiation (Lnet), sensible heat flux(QH) and latentheat fluxes (QE).At Castle 2012, the difference between set vs observed albedo yields a difference of 20% inQM, while for Conrad 2016 (in the accumulation zone) the difference in QM is 30%. On theother hand, the difference in QM at other sites is within 10%. As expected, the larger thedifference between the set and observed albedo, the larger the difference (bias) in QM. Sincethe WRF model underestimates the occurrence of snowfall over the observational period, theoverall albedo is underestimated, leading to overestimation of Knet relative to the one derivedwith the observed albedo in the SEB model (Figure 3.19). The effect of different roughness44length schemes on altering QH, and consequently on altering QM, is small because (1) z0v isof the same order of magnitude (10−3 m) in both schemes (Table 2.5), and (2) the Andreas(1987) model is known to simulate well the order of magnitude difference between scalar andmomentum roughness lengths for these glacier surfaces (Fitzpatrick et al., 2017; Radic´ et al.,2017). As long as the correct order of magnitude is used, the choice of roughness length formomentum has a small effect on mean daily QH (< 5% difference) and consequently negligibleeffect on mean daily QM (< 1% difference).3.3.4. SEB model versus PDD model performanceThe simple PDD model (equation (1.2)) forced by mod2.5 2m-temperature, using the calibratedmelt factors ( fm, Table 2.5), simulates the cumulative melt over the whole observation periodwithin 10% difference from the AWS-derived melt. Despite the good performance of the PDDmodel, its major limitation is the required calibration of the melt factors. When altering thecalibrated value of fm by ± 1 mm w.e. day−1 ◦C−1 at each site, for example, daily mean meltrate is changed by up to 40% from its original value (Figure 3.20).Time (day)220 240 260zM (m)-0.6-0.4-0.200.20.40.6Castle 2010 & 2012SEB (AWS)PDD (mod 2.5 )PDD (mod 1.0 )PDD (AWS)PDD-  (mod 2.5 )PDD+  (mod 2.5 )Time (day)200 220 240-3.5-3-2.5-2-1.5-1-0.50Nordic 2014Time (day)200 220 240-3-2.5-2-1.5-1-0.50Conrad st1 2015Time (day)200 220 240-2.5-2-1.5-1-0.50Conrad st2 2015Time (day)180 200 220 240-5-4-3-2-10Conrad st1 2016Time (day)180 200 220-1.6-1.4-1.2-1-0.8-0.6-0.4-0.200.2Conrad st2 2016Figure 3.20: Melt converted to surface lowering (zM) as simulated by the simple PDDmodel and by the SEB model (with observed albedo and roughness lengths) foreach site. The PDD model is forced by 2m-temperature from AWS and from WRFat different spatial resolution. Results of the sensitivity test where the calibratedsite-specific melt factor in the PDD model is perturbed by ± 1 mm w.e. day−1◦C−1 are shown for mod2.5 (PDD− and PDD+).Note that the perturbation in the melt factor of ± 1 mm w.e. day−1 ◦C−1 is relatively smallconsidering that, for a set of neighboring glaciers in the region, reported glacier-specific fm forice and for snow are in the range 4.0 - 9.7 mm w.e. day−1 ◦C−1 and 2.4 - 6.6 mm w.e. day−1◦C−1, respectively (Radic´ and Hock, 2011). In addition, as has been already noted before (seeHock, 2005), it is questionable whether the melt factor for the same location can be treated as45a constant in time. In our study, the melt factor for ice differs by more than 1 mm w.e. day−1◦C−1 when calibrated for the approximately same site in two different seasons (Castle site for2010 and 2012, and Conrad site for 2015 and 2016; Table 2.5).464. DiscussionSeveral recent studies showed a successful application of the WRF model in dynamicallydownscaling meteorological fields needed to force a melt model for glacier surfaces. Herewe focused on evaluating the WRF model in simulating meteorological variables and fluxesneeded to force a single-point SEB model at three mountain glaciers in the interior of BC. Asreference data for model evaluation we used observations from the automated weather stationsoperating intermittently at these glaciers in summers of 2010-2016. The ERA-Interim reanal-ysis was dynamically downscaled with the WRF model in three domains with resolutions of7.5 km, 2.5 km, and 1 km. The inner nest with 1 km resolution was only used for one site.We also analyzed the sensitivity of modeled daily melt rates to the choice of spatial resolutionand landcover in the WRF model, and to the choice of parameterization schemes for albedoand surface roughness in the SEB model. With the exception of near-surface air temperature,all meteorological variables measured by AWS in the ablation areas of the glaciers are bettersimulated by ERA-Interim than by WRF. In contrast, WRF performs best for simulating dailymelt for our single accumulation area site. The better performance of ERA-Interim over theWRF model is, in part, expected since the data assimilation is incorporated in the reanalysismodel, while the WRF model does not include any nudging to available in-situ and/or remotesensing observations for the domain.WRF at 2.5 km resolution overestimates the frequency of clear sky days over the obser-vation periods, overestimating the seasonal incoming shortwave radiation (Kin), and underes-timating the seasonal incoming longwave radiation (Lin). These results corroborate previousstudies that showed the overestimation of Kin in a domain with finest resolution indicates aproblem with convective cloud simulation over complex terrain (e.g., Claremar et al., 2012;Collier et al., 2013; Schanke et al., 2015). An increase in the complexity of the topographyaffects the cloudiness. If this effect is not sufficiently well described by the model resolution itmay induce more errors in the cloud cover relative to a coarser resolution. Furthermore, we donot use a cumulus parameterization in the innermost domain (2.5 and 1 km) assuming that cu-mulus convection is explicitly resolved. However, previous studies indicate that a grid spacingon the order of 100 m (Bryan et al., 2003; Petch, 2006) or even 10 m (Craig and Doernbrack,2008) is needed to capture the dominant length scales of moist cumulus convection.The finer model resolution (2.5 km versus 7.5 km) has better skill in simulating melt rates,but this is mainly due to the cancellation of biases in the SEB. The positive bias in Kin is com-pensated by the negative biases in Lin and sensible heat fluxes (QH). The underestimation ofQH down to -80% is due to underestimated near-surface wind speeds used in the bulk aerody-47namic method. As shown for our sites in the ablation area, WRF fails to simulate the prevailingdownslope katabatic flow and, therefore, fails to capture the adiabatic cooling in the ablationarea. The latter effect is less important for the simulation of QH , thus, the poor simulation ofwind speeds explains the poor performance in estimating QH . On the other hand, better WRFperformance in both temperature and wind speed is found for the site in the accumulation area.While WRF at few-kilometer resolution is able to simulate katabatic or other slope winds atlarge outlet glaciers (> 100 km2) and ice caps (Claremar et al., 2012), it fails to do so at oursites even at 1 km spatial resolution. We conclude that, for relatively small glaciers (5-15 km2)in complex terrain, the 1 km grid spacing might not be fine enough to simulate the local topo-graphic wind effect which highly depends on local gradients in the topography that correctlycaptures ice/snow cover. As indicated in Claremar et al., (2012), increasing the vertical reso-lution in the atmospheric boundary layer (e.g., 2 m grid in the lowest 10s of meters; Soderbergand Parmhed, 2006) might improve the simulation of near-surface wind speed and direction.It is likely that, for our sites, an increase in both horizontal and vertical resolution would beneeded for better simulation of katabatics, but this would also increase computational time ofthe simulations. Another problem in simulating QH is the use of inadequate bulk methodsrooted in the Monin-Obukhov theory, which in inadequate in its application to sloped glaciersurfaces (e.g., Radic´ et al., 2017). While we used the downscaled temperature and wind speedto force the bulk method, the method is based on the same theory used for parameterizing theboundary layer turbulence in the WRF model (e.g., Janjic, 1996; Nakanishi and Niino, 2009).To potentially improve the simulation of katabatics and QH at our sites in the ablation area,without increasing the resolution in WRF or running an eddy-simulator, one could make useof the successfully downscaled variables at the accumulation area. Downscaled temperatureand wind speed in the accumulation area could serve, for example, as boundary conditions toa katabatic flow model (Prandtl, 1942) solved with an elevation-varying eddy viscosity in theflux-gradient method (e.g., Grisogono and Oerlemans, 2002).A set of sensitivity tests showed that modifying the landcover category in WRF at 2.5km for the grid cell representative of the AWS site has negligible effect on the downscaledvariables. This sensitivity analysis, however, consisted of altering the landcover of only onegrid cell in the domain, which is likely too small a perturbation in the input conditions tosubstantially impact the model output. Furthermore, since the albedo in the SEB model istreated independently of WRF, the modeled melt rates are also not affected by this alterationin the landcover. It remains to be seen, however, how important the landcover alteration wouldbe in coupled WRF-SEB modeling approach (e.g., Mo¨lg and Kaser, 2011). The modeledmelt rates are shown to be sensitive to the albedo scheme in the SEB model, indicating the48importance of (i) accurate representation of albedo, and (ii) a correct simulation of snowfallevents throughout the ablation season. While the first point can be addressed by incorporatinga better albedo model (e.g., Schmidt et al., 2017), simulating snowfall, and precipitation ingeneral, remains a major challenge for WRF. Underestimation of frequency and intensity ofprecipitation is consistent with an overestimation of Kin and, as pointed out by Collier (2013),may reflect errors in the forcing data at the lateral boundaries and/or in the WRF resolutionto fully resolve orographic uplift or microscale complex flow features that affect precipitationat the sites. The modeled melt rates are shown to be insensitive to the choice of roughnesslengths as long as the roughness length for momentum (z0) is within one order of magnitudeof its observed value. The value used here (z0=10−3 m) is commonly assigned for glacier icesurfaces (e.g., Fitzpatrick et al., 2017). On the other hand, the underestimation of the near-surface wind speed in WRF, leading to a substantial underestimation of QH , makes the correctchoice for z0 almost irrelevant. Due to erroneous input wind speeds to the bulk method, a betterfirst order approximation for QH would be to parametrize it solely as a function of downscaledtemperature, assuming a constant drag coefficient and setting the mean wind speed to somereasonable constant value (e.g. 3 m s−1) throughout the observation period.495. ConclusionsThe main goal of this study was to investigate the performance of a dynamical downscalingapproach to derive meteorological fields needed to force a physics-based model of glaciermelt. Our main conclusions, as answers to the questions from the introduction section, arethe following:1. How do downscaled sub-daily and daily timeseries of meteorological variables comparewith the observations at each site? ERA-Interim climate reanalysis at 80 km resolutioncan successfully simulate the glacier-scale radiative fluxes that are the dominant driversof daily melt rates during the ablation season. Dynamical downscaling of ERA-Interim,with the state-of-the-art WRF model, to 2.5 km and 1 km resolution was successful fornear-surface air temperature, somewhat successful for radiative fluxes, and unsuccess-ful for near-surface wind speed and turbulent heat fluxes. WRF at 2.5 km resolutionoverestimates the frequency of clear sky days over the observation periods, overestimat-ing seasonal incoming shortwave radiation (Kin), and underestimating seasonal incominglongwave radiation (Lin). The sensible heat flux (QH) is underestimated down to -80%due to underestimated near-surface wind speeds.2. How sensitive are the downscaled fields to the choice of spatial resolution and landcoverdata in the WRF model? A set of sensitivity tests showed that an increase of spatialresolution in WRF from 7.5 km to 2.5 km, and 1 km for one tested site, did not improvethe simulation of downscaled variables other than near-surface air temperature. For rel-atively small glaciers (5-15 km2) located in a complex terrain, the 2.5 and 1 km gridspacing might not be fine enough to simulate the local cloud convection and topographicwind effects (katabatics). If this complex terrain is inadequately captured in the WRFtopography, a model at higher resolution may induce more errors in the cloud cover rel-ative to a coarser resolution. While the elevation of the WRF grid cell that covers theAWS site agrees well with the actual station elevation, the grid-cell land category is oftennot classified correctly as ice or snow. The sensitivity of model results to the wronglyassigned land category is nevertheless small.3. If a surface energy balance (SEB) model is used to calculate point-scale daily and sea-sonal surface melt at these glaciers, how successfully can the SEB model be forced bythe WRF data instead of the AWS data? The best performance in simulating seasonalmelt rates (within± 10% difference from the measured melt) in ablation area is obtainedwhen the SEB model is forced by the ERA-Interim fluxes, and not by the dynamically50downscaled fluxes. For the site in the accumulation area, on the other hand, the WRFmodel at 2.5 km resolution yielded best match between simulated and AWS-derived sea-sonal melt rates. The finer WRF resolution (2.5 km versus 7.5 km) is shown to havebetter skill in simulating the seasonal melt rates, but this is mainly due to the cancella-tion of biases in the SEB model. The positive bias in seasonal Kin is compensated by thenegative biases in seasonal Lin and sensible heat fluxes (QH). The modeled melt ratesare shown to be sensitive to the albedo scheme in the SEB model, indicating the impor-tance of an accurate representation of albedo, and a correct simulation of snowfall eventsthroughout the ablation season.Considering the limitations with the WRF downscaling, we advocate for the SEB modelingdirectly forced with the output from global climate reanalysis for the past and present climate,and from global climate models (GCMs) for future climate scenarios. Those GCMs, however,need to be evaluated against the climate reanalysis prior to their applications in the projections.Temperature downscaling remains important for the simulation of turbulent fluxes, especiallyif one is interested in capturing the variability in melt rates throughout the ablation season.More importantly, however, we found that the WRF failure to resolve the katabatic flow in theablation area of the relatively small glaciers (5-15 km2) strongly deteriorates the simulationof turbulent heat fluxes even when near-surface air temperature is successfully downscaled.A way forward, as we see it, is to use the air temperature and wind fields, downscaled to theaccumulation area of a glacier, as boundary conditions to a locally nested katabatic flow model.51BibliographyAas, K. S., T. Dunse, E. Collier, T. V. Schuler, T. K. Berntsen, J. Kohler, and B. Luks2016. The climatic mass balance of Svalbard glaciers: a 10-year simulation with a coupledatmosphere-glacier mass balance model. Cryosphere, 10(3):1089–1104. → page 11Andreas, E.1987. A Theory for the scalar rougness and the scalar transfer-coefficients over snow andsea ice. Boundary-Layer Meteorology, 38(1-2):159–184. → pages 21, 45Andreassen, L. M., M. R. Van Den Broeke, R. H. Giesen, and J. Oerlemans2008. A 5 year record of surface energy and mass balance from the ablation zone ofStorbreen, Norway. Journal of Glaciology, 54(185):245–258. → pages 6, 10Beedle, M. J., B. Menounos, and R. Wheate2015. Glacier change in the Cariboo Mountains, British Columbia, Canada (1952-2005).Cryosphere, 9(1):65–80. → page 13Blacutt, L. A., D. L. Herdies, L. G. G. de Goncalves, D. A. Vila, and M. Andrade2015. Precipitation comparison for the CFSR, MERRA, TRMM3B42 and CombinedScheme datasets in Bolivia. Atmospheric Research, 163(SI):117–131. → page 10Braithwaite, R. J. and O. B. Olesen1989. Calculation of glacier ablation from air temperature, west greenland. In GlacierFluctuations and Climatic Change, J. Oerlemans, ed., Pp. 219–233, Dordrecht. SpringerNetherlands. → page 4Braithwaite, R. J. and Z. Yu2000. Sensitivity of mass balance of five Swiss glaciers to temperature changes assessed bytunning a degree-day model. Journal of Glaciology, 46(152):7–14. → pages 1, 4Bryan, G., J. Wyngaard, and J. Fritsch2003. Resolution requirements for the simulation of deep moist convection. MonthlyWeather Review, 131(10):2394–2416. → page 47Carenzo, M., F. Pellicciotti, S. Rimkus, and P. Burlando2009. Assessing the Transferability and robustness of an enhanced temperature-indexglacier-melt model. Journal of Glaciology, 55(190):258–274. → page 5Claremar, B. o¨., F. Obleitner, C. Reijmer, V. Pohjola, A. Waxeg a˚rd, F. Karner, andA. Rutgersson2012. Applying a Mesoscale Atmospheric Model to Svalbard Glaciers. Advances inMeteorology, 2012:1–22. → pages 11, 12, 23, 47, 48Clarke, G. K., A. H. Jarosch, F. S. Anslow, V. Radic´, and B. Menounos2015. Projected deglaciation of western Canada in the twenty-first century. NatureGeoscience, 8:372–377. → pages 1, 5, 1152Collier, E., F. Maussion, L. Nicholson, T. Mo¨lg, W. Immerzeel, and A. Bush2015. Impact of debris cover on glacier ablation and atmospheric-glacier feedbacks in theKarakoram. The Cryosphere, 9:1617–1632. → pages 11, 12, 16, 21Collier, E., T. Mo¨lg, F. Maussion, D. Scherer, C. Mayer, and A. Bush2013. High-resolution interactive modelling of the mountain glacier-atmospheric interface:an application over the Karakoram. The Cryosphere, 7:779–795. → pages11, 12, 15, 16, 21, 47, 49Conway, J. and N. Cullen2013. Constraining turbulent heat flux parameterization over a temperate maritime glacierin New Zealand. Annals of Glaciology, 54(63):41–51. → pages 6, 8, 10Craig, G. C. and A. Doernbrack2008. Entrainment in Cumulus Clouds: What Resolution is Cloud-Resolving? Journal ofThe Atmospheric Sciences , 65(12):3978–3988. → page 47Cuffey, K. M. and W. Paterson2010. The Physics of Glaciers, 4 edition. Academic Press. → pages 4, 7, 23Danielson, J. and D. Gesch2011. Global multi-resolution terrain elevation data 2010 (gmted2010). Technical report,U.S. Geological Survey. → page 19Dee, D. P., S. M. Uppala, A. J. Simmons, P. Berrisford, P. Poli, S. Kobayashi, U. Andrae,M. A. Balmaseda, G. Balsamo, P. Bauer, P. Bechtold, A. C. M. Beljaars, L. van de Berg,J. Bidlot, N. Bormann, C. Delsol, R. Dragani, M. Fuentes, A. J. Geer, L. Haimberger, S. B.Healy, H. Hersbach, E. V. Holm, L. Isaksen, P. Kallberg, M. Koehler, M. Matricardi, A. P.McNally, B. M. Monge-Sanz, J. J. Morcrette, B. K. Park, C. Peubey, P. de Rosnay,C. Tavolato, J. N. Thepaut, and F. Vitart2011. The ERA-Interim reanalysis: configuration and performance of the data assimilationsystem. Quarterly Journal of the Royal Meterological Society, 137:553–597. → pages16, 19Dery, S. J., A. Clifton, S. MacLeod, and M. J. Beedle2010. Blowing Snow Fluxes in the Cariboo Mountains of British Columbia, Canada. ArcticAntarctic and Alpine Research, 42(2):188–197. → page 13ESA2018. Esa climate change initiative (cci) data. Technical report, Center for EnvironmentalData Analysis. → page 19Fitzpatrick, N., V. Radic, and B. Menounos2017. Surface Energy Balance Closure and Turbulent Flux Parameterization on aMid-Latitude Mountain Glacier Purcell Mountains, Canada. Frontiers in Earth Science, 5.→ pages 8, 10, 14, 20, 24, 28, 45, 4953Fowler, H. J. and R. L. Wilby2007. Beyond the downscaling comparison study. International Journal of Climatology,27(12):1543–1545. → page 2Gabbi, J., M. Carenzo, F. Pellicciotti, A. Bauder, and M. Funk2014. A comparison of empirical and physically based glacier surface melt models forlong-term simulations of glacier response. Journal of Glaciology, 60(224):1140–1154. →page 5Garner, K. L., M. Y. Chang, M. T. Fulda, J. A. Berlin, R. E. Freed, M. M. Soo-Hoo, D. L.Revell, M. Ikegami, L. E. Flint, A. L. Flint, and B. E. Kendall2015. Impacts of sea level rise and climate change on coastal plant species in the centralCalifornia coast. PEERJ, 3. → page 1Giesen, R. H. and J. Oerlemans2012. Calibration of a surface mass balance model for global-scale applications.Cryosphere, 6(6):1463–1481. → pages 6, 25Gillett, S. and N. J. Cullen2011. Atmospheric controls on summer ablation over Brewster Glacier, New Zealand.International Journal of Climatology, 31:2033–2048. → pages 4, 6, 10Gregory, J. and J. Oerlemans1998. Simulated future sea-level rise due to glacier melt based on regionally and seasonallyresolved temperature changes. Nature, 391(6666):474–476. → page 1GRELL, G.1993. Prognostic Evaluation of Assumptions used by Cumulus Parameterizations.Monthly Weather Review, 121(3):764–787. → page 19Grell, G. and D. Devenyi2002. A generalized approach to parameterizing convection combining ensemble and dataassimilation techniques. Geophysical Research Letters, 29(14). → page 19Grisogono, B. and J. Oerlemans2002. Justifying the WKB approximation in pure katabatic flows. Tellus Series A-DynamicMeteorology and Oceanography, 54(5):453–462. Climate Conference 2001, Utrecht Univ,Utrecht, Netherlands, Aug, 2001. → page 48Hock, R.1999. A distributed temperature-index ice-and snowmelt model including potential directsolar radiation. Journal of Glaciology, 45(149):101–111. → pages 4, 5Hock, R.2003. Temperature index melt modelling in mountain areas. Journal of Hydrology,282:104–115. → pages 1, 4, 2554Hock, R.2005. Glacier melt: a review of processes and their modelling. Progress in PhysicalGeography, 29(3):362–391. → pages 1, 3, 4, 6, 8, 10, 45Hock, R. and B. Holmgren2005. A distributed surface energy-balance model for complex topography and itsapplication to Storglaciaren, Sweden. Journal of Glaciology, 51(172):25–36. → page 6Huss, M. and R. Hock2015. A new model for global glacier change and sea-level rise. Frontiers in Earth Science,3(54):1–22. → pages 1, 10, 25Iacono, M. J., J. S. Delamere, E. J. Mlawer, M. W. Shephard, S. A. Clough, and W. D. Collins2008. Radiative forcing by long-lived greenhouse gases: Calculations with the AERradiative transfer models. Journal of Geophysical Research-Atmosphere, 113(D13). →page 19Janjic, Z.1996. The mellor-yamada level 2.5 turbulence closure scheme in the ncep eta model. WorldMeteorological Organization-Publications-WMO TD, Pp. 4–14. → page 48Jarosch, A. H., F. S. Anslow, and G. K. Clarke2010. High-resolution precipitation and temperature downscaling for glacier models. ClimDyn, 38:391–409. → pages 10, 11Johannesson, T., O. Sigurdsson, T. Laumann, and M. Kennett1995. Degree-day glacier mass-balance modellling with applications to glaciers in Iceland,Norway and Greenland. Journal of Glaciology, 41(138):345–358. → pages 1, 4Jonsell, U., F. Navarro, M. Ban˜o´n, J. Lapazaran, and J. Otero2012. Sensitivity of a distributed temperature-radiation index melt model based on AWSobservations and surface energy balance fluxes, Hurd Peninsula glaciers, Liningston Island,Antarctica. The Cryosphere, 6:539–552. → pages 4, 5Klok, E. L. and J. Oerlemans2002. Model study of the spatial distribution of the energy and mass balance ofMorteratschgletscher, Switzerland. Journal of Glaciology, 48(163):505–518. → page 6Kotlarski, S., F. Paul, and D. Jacob2010. Forcing a Distributed Glacier Mass Balance Model with the Regional Climate ModelREMO. Part I: Climate Model Evaluation. Journal of Climate, 23:1589–1606. → page 11Lawrence, G. W.1985. The use of general circulation models in the analysis of the ecosystem impacts ofclimatic change. Climatic change, 7:267–284. → page 255MacDougall, A., B. Wheler, and G. Flowers2011. A preliminary assessment of glacier melt-model parameter sensitivity andtransferability in a dry subarctic environment. The Cryosphere, 5:1011–2011. → pages1, 4, 5MacDougall, A. H. and G. Flowers2010. Spatial and Temporal Transferability of a Distributed Energy-Balance Glacier MeltModel. Journal of Climate, 24:1480–1498. → page 6Machguth, H., F. Paul, S. Kotlarski, and M. Hoelzle2009. Calculating distributed glacier mass balance for the Swiss Alps from regional climatemodel output: A methodical description and interpretation of the results. Journal ofGeophysical Research, 114:1–19. → page 11Marzeion, B., A. Jarosch, and M. Hofer2012. Past and future sea-level change from the surface mass balance of glaciers. TheCryosphere, 6:1295–1322. → pages 1, 25Matulla, C., E. Watson, S. Wagner, and W. Schoener2009. Downscaled GCM projections of winter and summer mass balance for Peyto Glacier,Alberta, Canada (2000-2100) from ensemble simulations with ECHAM5-MPIOM.International Journal of Climatology, 29(11):1550–1559. → page 11Maussion, F., S. D., R. Finkelnburg, J. Richters, W. Yang, and T. Yao2011. WRF simulation of a precipitation event over the Tibetan Plateau, China - anassessment using remote sensing and ground observations. Hydrology and Earth SystemSciences, 15:1795–1817. → page 12Mearns, L., M. Bukovsky, S. C. Pryor, and V. Magaa2014. Downscaling of Climate Information. → page 11Mo¨lg, T. and G. Kaser2011. A new approach to resolving climate-cryosphere relations: Downscaling slimatedynamics to glacier-scale mass and energy balance without statistical scale linking. JournalofGeophysical research, 116:779–795. → pages 11, 12, 15, 16, 48Mo¨lg, T., F. Maussion, W. Yang, and S. D.2012. The footprint of Asian monsoon dynamics in the mass and energy balance of aTibetan glacier. The Cryosphere, 6:1445–1461. → pages 11, 12, 21Monin, A. and A. Obukhov1954. Basic laws of turbulent mixing in the surface layer of the atmosphere. Tr. Akad. NaukSSSR Geophiz, 24(151):163–187. → page 9Montavez, J. P., J. M. Lopez-Romero, S. Jerez, J. J. Gomez-Navarro, and P. Jimenez-Guerrero2017. How much spin-up period is really necessary in regional climate simulations?Geophysical Research Abstracts. → page 1756Munro, D.1989. Surface rougness and bulk heat transfer on a glacier: comparison with eddycorrelation. Journal of Glaciology, 35(121):343–348. → page 6Munro, D. and M. Marosz-Wantuch2009. Modeling Ablation on Place Glacier, British Columbia, from Glacier and Off glacierData Sets. Arctic, Antarctic, and Alpine Research, 41(2):246–256. → page 6Nakanishi, M. and H. Niino2006. An improved mellor-yamada level-3 model: Its numerical stability and application toa regional prediction of advection fog. Bounday-Layer Meteorology, 119(2):397–407. →page 19Nakanishi, M. and H. Niino2009. Development of an Improved Turbulence Closure Model for the AtmosphericBoundary Layer. Journal of the Meteorological Society of Japan, 87(5):895–912. → page48Niu, G.-Y., Z.-L. Yang, K. E. Mitchell, F. Chen, M. B. Ek, M. Barlage, A. Kumar,K. Manning, D. Niyogi, E. Rosero, M. Tewari, and Y. Xia2011. The community Noah land surface model with multiparameterization options(Noah-MP): 1. Model description and evaluation with local-scale measurements. Journal ofGeophysical Research-Atmosphere, 116. → page 19Oerlemans, J.2001. Glaciers and Climate Change. A.A.Balkema, Lisse, Abingdon, Exton (PA), Tokyo,145 pp. → page 5Oerlemans, J., R. H. Giesen, and M. R. Van den Broeke2009. Retreating alpine glaciers: increased melt rates due to accumulation of dust (Vadretda Morteratsch, Switzerland). Journal of Glaciology, 55(192):729–736. → page 6Oerlemans, J. and E. Klok2002. Energy Balance of Glacier Surface: Analysis of Automatic Weather Station Datafrom the Morteratschgletscher, Switzerland. Arctic, Antarctic, and Alpine Research,34(4):477–485. → page 6Ontiveros-Gonzalez, G., H. Delgado-Granados, and J. Cortes-Ramos2015. Surface Energy Balance model for high-altitude glacial system at 19 degrees N onGlaciar Norte, Mexico. Geofisica Internacional, 54(4):299–314. → page 6Parajuli, A., M. B. Chand, R. B. Kayastha, J. M. Shea, and P. K. Mool2015. Modified temperature index model for estimating the melt water discharge fromderbis-covered Lirung Glacier, Nepal. Remote Sensing and GIS for Hydrology and WaterResources, 368:409–414. → page 457Paul, F. and S. Kotlarski2010. Forcing a Distributed Glacier Mass Balance Model with the Regional Climate ModelREMO. Part II: Downscaling Strategy and Results for Two Swiss Glaciers. Journal ofClimate, 23:1607–1620. → page 11Pellicciotti, F., B. Brock, U. Strasser, P. Burlado, M. Funk, and J. Corripio2005. An enhanced temperature-index glacier melt model including the shortwaveradiation balance: development and testing for Haut Glacier dA´rolla, Switzerland. Journalof Glaciology, 51(175):573–587. → page 5Petch, J. C.2006. Sensitivity studies of developing convection in a cloudresolving model. QuarterlyJournal of the Royal Meteorological Society, 132(615):345–358. → page 47Powers, J. G., J. B. Klemp, W. C. Skamarock, C. A. Davis, J. Dudhia, D. O. Gill, J. L. Coen,D. J. Gochis, R. Ahmadov, S. E. Peckham, G. A. Grell, J. Michalakes, S. Trahan, S. G.Benjamin, C. R. Alexander, G. J. Dimego, W. Wang, C. S. Schwartz, G. S. Romine, Z. Liu,C. Snyder, F. Chen, M. J. Barlage, W. Yu, and M. G. Duda2017. The Weather Research and Forcasting Model: Overview, System Efforts, and FutureDirections. Bulletin of the American Meteorological Society, 98(8):1717–1737. → pages12, 15Prandtl, L.1942. Bemerkungen zur Theorie der freien Turbulenz . Zeitschrift Angewandte Mathematikund Mechanik, 22:241–243. → page 48Radic´, V., A. Bliss, A. C. Beedlow, R. Hock, E. Miles, and J. G. Cogley2014. Regional and global projections of twenty-first century glacier mass changes inresponse to climate scenarios from global climate models. Clim Dyn, 42:37–58. → pages1, 10, 25Radic´, V. and R. Hock2006. Modelling future glacier mass balance and volume changes using ERA-40 reanalysisand climate models: A sensitivity study at Storglaciaren, Sweden. Journal of GeophysicalResearch, 111:1–12. → pages 1, 4, 11Radic´, V. and R. Hock2011. Regionallly differentiated contribution of mountain glaciers and ice caps to futuresea-level rise. Nature Geoscience, 4:91–94. → page 45Radic, V. and R. Hock2014. Glaciers in the Earth’s Hydrological Cycle: Assessments of Glacier Mass and RunoffChanges on Global and Regional Scales. Surveys in Geophysics, 35(3):813–837. → page 1Radic´, V., B. Menounos, J. Shea, N. Fitzpatrick, M. A. Tessema, and S. J. Dery2017. Evaluation of different methods to model near-surface turbulent fluxes for a58mountain glacier in the Cariboo Mountains, BC, Canada. Cryosphere, 11(6):2897–2918.→ pages 6, 8, 10, 14, 20, 24, 28, 45, 48Rango, A. and J. Martinec1995. Revisiting the Degree-Day Method For Snowmelt Computation. Water ResoucesBulletin, 31(4):657–669. → pages 4, 5Schanke, K., T. K. Berntsen, J. Boike, B. Etzelmuller, J. E. Kristjansson, M. Maturilli, T. V. Schuler, F. Stordal, and S. Westermann2015. A Comparison between Simulated and Observed Surface Energy Balance at theSvalbard Archipelago. Journal of Applied Meteorology and Climatology, 54:1102–1119.→ pages 11, 47Schmidt, L. S., G. Adalgeirsdottir, S. Gudmundsson, P. L. Langen, F. Palsson, R. Mottram,S. Gascoin, and H. Bjornsson2017. The importance of accurate glacier albedo for estimates of surface mass balance onVatnajokull: evaluating the surface energy budget in a regional climate model withautomatic weather station observations. Cryosphere, 11(4). → page 49Shea, J. M., R. Moore, and K. Stahl2009. Derivation of melt factors from glacier mass-balance records in western Canada.Journal of Glaciology, 55(189):123–130. → page 5Sicart, J., P. Wagnon, and P. Ribstein2005. Atmospheric controls of the heat balance of Zongo Glacier (16 degrees S, Bolivia).Journal of Geophysical Research-Atmosphere, 110(D12). → page 10Soderberg, S. and O. Parmhed2006. Numerical modelling of katabatic flow over a melting outflow glacier.Boundary-Layer Meteorology, 120(3):509–534. → page 48Soruco, A., C. Vincent, A. Rabatel, B. Francou, E. Thibert, E. J. Sicart, and T. Condom2015. Contribution of glacier runoff to water resouces of La Paz city Bolivia (16oS).Annals of Glaciology, 56(70):147–154. → page 1Springer, C., C. Matulla, W. Sch o¨ner, R. Steinacker, and S. Wagner2013. Downscaled GCM projections of winter and summer massbalance for CentralEuropean glaciers (2000-2100) fromensemble simulations with ECHAM5-MPIOM.International Journal of Climatology, 33:1270–1279. → page 11Stull, B.2003. An introduction to Boundary Layer Meteorology. Kluwer Academic, Dordrecht,Boston, London, 308 pp. → page 7Sun, W., X. Qin, J. Ren, X. Yang, T. Zhang, Y. Liu, X. Cui, and W. Du2012. The Surface Energy Budget in the Accumulation Zone of the Laohugou Glacier No.12 in the Western Qilian Mountains, China, in Summer 2009. Arctic Antarctic and AlpineResearch, 44(3):296–305. → page 659Thompson, G., P. R. Field, R. M. Rasmussen, and W. D. Hall2008. Explicit Forecasts of Winter Precipitation Using an Improved Bulk MicrophysicsScheme. Part II: Implementation of a New Snow Parameterization. Monthly WeatherReview, 136(12):5095–5115. → page 19Wilks, D. S.2006. Statistical Methods in the Atmospheric Sciences, International Geophysics Series,2nd edition. Elsevier Academic Press. → page 24WMO2008. Guide to Meteorological Instruments and Methods of Observation. WMO, WorldMeteorological Organization, wmo-no.8 edition. → page 8Woul, M. D. and R. Hock2005. Static mass-balance sensitivity of Arctic glaciers and ice caps using a degree-dayapproch. Annals of Glaciology, 42:217–224. → page 1Wrzesien, M. L., T. M. Pavelsky, S. B. Kapnick, M. T. Durand, and T. H. Painter2015. Evaluation of snow cover fraction for regional climate simulations in the SierraNevada. International Journal of Climatology, 35(9):2472–2484. → page 11Wu, X., Y. Shen, N. Wang, X. Pan, W. Zhang, J. He, and G. Wang2016. Coupling the WRF model with a temperature index model based on remote sensingfor snowmelt simulations in a river basin in the Altay Mountains, north-west China.Hydrological Processes, 30(21):3967–3977. → page 11Yang, Z.-L., G.-Y. Niu, K. E. Mitchell, F. Chen, M. B. Ek, M. Barlage, L. Longuevergne,K. Manning, D. Niyogi, M. Tewari, and Y. Xia2011. The community Noah land surface model with multiparameterization options(Noah-MP): 2. Evaluation over global river basins. Journal of GeophysicalResearch-Atmosphere, 116. → page 19Zemp, M., M. Hoelzle, and W. Haeberli2009. Six decades of glacier mass-balance observations: a review of the worldwidemonitoring network. Annals of Glaciology, 50(50):101–111. → page 1Zorita, E. and H. v. Storch1999. The Analog Method as a Simple Statistical Downscaling Technique: Comparisonwith More Complicated Methods. Journal of Climate, 12:2474–2489. → page 1160A. AppendixA.1. Supplementary figures and tablesA.1.1. Time series of meteorological variables and SEB fluxesT 2m (o C)01020RH2m (%)406080100U (m/s)2468SP (hPa)780800820Time (day)200 205 210 215 220 225 230 235 240 245P (mm)01020ERAmod7.5mod2.5AWSFigure A.1: Comparison of observed versus modeled 3-hr meteorological variables atConrad glacier in 2015, station 1: 2m-air temperature (T2m), 2m-relative humid-ity (RH2m), 10 m (modeled) and 2 m (observed) wind speed (U), surface pressure(SP) and accumulated precipitation (P).61Kin (W/m2 )0200400600800L in (W/m2 )250300350QH (W/m2 )0100200QE (W/m2 )0100200Time (day)200 205 210 215 220 225 230 235 240 245QM (W/m2 )0200400600800ERAmod7.5mod2.5AWSFigure A.2: Comparison of observed or AWS-derived versus WRF-derived 3-hr SEBfluxes at Conrad glacier in 2015, station 1: incoming shortwave radiation (Kin),incoming longwave radiation (Lin), sensible heat flux (QH), latent heat flux (QE)and melt energy flux (QM).62T 2m (o C)01020RH2m (%)406080100U (m/s)246SP (hPa)780800820Time (day)180 190 200 210 220 230 240P (mm)0510ERAmod7.5mod2.5AWSFigure A.3: Comparison of observed versus modeled 3-hr meteorological variables atConrad glacier in 2016, station 1.Kin (W/m2 )05001000L in (W/m2 )250300350QH (W/m2 )050100150QE (W/m2 )-2002040Time (day)180 190 200 210 220 230 240QM (W/m2 )0200400600800ERAmod7.5mod2.5AWSFigure A.4: Comparison of observed or AWS-derived versus WRF-derived 3-hr SEBfluxes at at Conrad glacier in 2016, station 1.63T 2m (o C)01020RH2m (%)406080100U (m/s)2468SP (hPa)720760800Time (day)170 180 190 200 210 220 230P (mm)051015 ERAmod7.5mod2.5AWSFigure A.5: Comparison of observed versus modeled 3-hr meteorological variables atConrad glacier in 2016, station 2.Kin (W/m2 )0200400600800L in (W/m2 )200250300350QH (W/m2 )0100200QE (W/m2 )-50050Time (day)170 180 190 200 210 220 230QM (W/m2 )0200400ERAmod7.5mod2.5AWSFigure A.6: Comparison of observed or AWS-derived versus WRF-derived 3-hr SEBfluxes at at Conrad glacier in 2016, station 2.64A.1.2. Wind rose diagrams8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)0123456789101112windspeed(m/s)Conrad st1 2015:obsConrad st1 2015:mod7.5 Conrad st1 2015:mod*2.5Conrad st1 2015:mod2.5Figure A.7: Wind roses of observed and simulated wind speed at Conrad glacier in 2015,station 1. Colorbar shows the wind speed (m s−1).8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)0123456789101112Castle 2010: obs Castle 2010: mod2.5 Castle 2010: mod7.5Castle 2010: mod*2.5 Castle 2010: mod1.0 Castle 2010: mod*1.0windspeed(m/s)Figure A.8: Wind roses of observed and simulated wind speed at Castle Creek glacier in2010.658.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)0123456789101112Conrad st1 2016:obs Conrad st1 2016:mod2.5Conrad st1 2016:mod7.5 Conrad st1 2016:mod*2.5windspeed(m/s)Figure A.9: Wind roses of observed and simulated wind speed at Conrad glacier in 2016,station 18.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)8.3333%16.6667%25%33.3333%0% E(90o)W(270o) N(0o)S(180o)0123456789101112Conrad st2 2016:obsConrad st2 2016:mod7.5 Conrad st2 2016:mod*2.5Conrad st2 2016:mod2.5windspeed(m/s)Figure A.10: Wind roses of observed and simulated wind speed at Conrad glacier in 2016,station 2.66A.1.3. Scatter plotsT2m obs5 10 15T 2m mod 2.551015r=0.85 p=0.00 MBE=1.4RH2m obs40 60 80 100RH2m mod 2.5406080100r=0.55 p=0.00 MBE=5.9U obs2 4 6 8 10U mod2.5246810r=0.00 p=0.98 MBE=-1.3SP obs770 780 790SP mod2.5770775780785790r=0.97 p=0.00 MBE=-4.8Kin obs0 200 400 600 800K in mod 2.50200400600800r=0.86 p=0.00 MBE=37.9Lin obs250 300 350L in mod 2.5250300350r=0.41 p=0.00 MBE=-18.9QH obs50 100 150 200Q H mod 2.550100150200r=0.03 p=0.62 MBE=-35.4QE obs0 50 100Q E mod 2.5050100r=0.36 p=0.00 MBE=-1.2QMelt obs0 200 400 600 800Q Melt mod 2.50200400600800r=0.86 p=0.00 MBE=-22.4Figure A.11: Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Nordic glacier in 2014 with correlation (r), p-value (p) and mean bias error (MBE) included for each plot. The variables are:T2m (◦C), RH2m (%), U (m s−1), SP (hPa), Kin (W m−2), Lin (W m−2), QH (Wm−2), QE (W m−2) and QM (W m−2).67T2m daily obs5 10 15T 2m daily mod2.5468101214r=0.94 p=0.00 MBE=1.4RH2m daily obs40 60 80RH2m daily mod2.5405060708090r=0.69 p=0.00 MBE=5.9U daily obs1 2 3 4 5U daily mod2.512345r=-0.06 p=0.71 MBE=-1.3SP daily obs775 780 785SP daily mod2.5775780785r=0.98 p=0.00 MBE=-4.8Kin daily obs100 200 300K in daily mod2.5100200300r=0.68 p=0.00 MBE=37.9Lin daily obs260 280 300 320 340L in daily mod2.5260280300320340r=0.50 p=0.00 MBE=-18.9QH daily obs20 60 100Q H daily mod2.520406080100120r=-0.02 p=0.90 MBE=-35.4QE daily obs0 20 40 60Q E daily mod2.50204060r=0.48 p=0.00 MBE=-1.2QMelt daily obs100 200 300Q Melt daily mod2.5100150200250300350r=0.66 p=0.00 MBE=-22.4Figure A.12: Comparison of AWS-derived (obs) versus WRF-derived (mod) daily time-series of the selected variables at Nordic glacier in 2014.68T2m obs0 5 10 15T 2m mod 2.5051015r=0.80 p=0.00 MBE=-0.2RH2m obs40 60 80 100RH2m mod 2.5406080100r=0.52 p=0.00 MBE=9.4U obs2 4 6 8U mod2.52468r=0.06 p=0.25 MBE=-1.2SP obs775 780 785 790SP mod2.5775780785790r=0.97 p=0.00 MBE=-2.1Kin obs0 200 400 600 800K in mod 2.50200400600800r=0.87 p=0.00 MBE=65.8Lin obs250 300 350L in mod 2.5250300350r=0.44 p=0.00 MBE=-26.2QH obs0 50 100QH mod 2.5050100r=0.17 p=0.00 MBE=-25.3QE obs0 40 80QE mod 2.5-20020406080r=0.41 p=0.00 MBE=0.6QMelt obs0 200 400 600QMelt mod 2.50200400600r=0.87 p=0.00 MBE=5.9Figure A.13: Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Conrad glacier in 2015, station 1.69T2m daily obs0 5 10T 2m daily mod2.50510r=0.91 p=0.00 MBE=-0.2RH2m daily obs50 60 70 80 90RH2m daily mod2.55060708090r=0.62 p=0.00 MBE=9.4U daily obs1 2 3 4 5U daily mod2.512345r=0.10 p=0.48 MBE=-1.2SP daily obs775 780 785 790SP daily mod2.5775780785790r=0.98 p=0.00 MBE=-2.1Kin daily obs100 200 300K in daily mod2.5100200300r=0.63 p=0.00 MBE=65.8Lin daily obs240 260 280 300 320L in daily mod2.5240260280300320r=0.57 p=0.00 MBE=-26.2QH daily obs0 20 40 60 80QH daily mod2.5020406080r=0.27 p=0.06 MBE=-25.3QE daily obs-10 0 10 20 30QE daily mod2.5-100102030r=0.50 p=0.00 MBE=0.6QMelt daily obs100 200 300QMelt daily mod2.550100150200250300r=0.62 p=0.00 MBE=5.9Figure A.14: Comparison of AWS-derived (obs) versus WRF-derived (mod) daily time-series of the selected variables at Conrad glacier in 2015, station 1.70T2m obs0 5 10 15T 2m mod 2.5051015r=0.80 p=0.00 MBE=-0.4RH2m obs40 60 80 100RH2m mod 2.5406080100r=0.52 p=0.00 MBE=11.7U obs2 4 6 8U mod2.52468r=0.03 p=0.52 MBE=-1.4SP obs775 780 785 790SP mod2.5775780785790r=0.97 p=0.00 MBE=-1.0Kin obs0 200 400 600 800K in mod 2.50200400600800r=0.86 p=0.00 MBE=67.3Lin obs250 300 350L in mod 2.5250300350r=0.42 p=0.00 MBE=-20.6QH obs0 50 100Q H mod 2.5050100r=0.20 p=0.00 MBE=-26.1QE obs0 50 100Q E mod 2.5050100r=0.41 p=0.00 MBE=2.4QMelt obs0 200 400 600Q Melt mod 2.50200400600r=0.89 p=0.00 MBE=10.1Figure A.15: Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Conrad glacier in 2015, station 2.71T2m daily obs0 5 10T 2m daily mod2.50510r=0.90 p=0.00 MBE=-0.4RH2m daily obs50 60 70 80 90RH2m daily mod2.55060708090r=0.63 p=0.00 MBE=11.7U daily obs1 2 3 4 5U daily mod2.512345r=0.05 p=0.73 MBE=-1.4SP daily obs775 780 785 790SP daily mod2.5775780785790r=0.98 p=0.00 MBE=-1.0Kin daily obs100 200 300K in daily mod2.5100200300r=0.61 p=0.00 MBE=67.3Lin daily obs240 260 280 300 320L in daily mod2.5240260280300320r=0.53 p=0.00 MBE=-20.6QH daily obs0 20 40 60 80QH daily mod2.5020406080r=0.31 p=0.03 MBE=-28.6QE daily obs-10 0 10 20 30QE daily mod2.5-100102030r=0.54 p=0.00 MBE=2.2QMelt daily obs100 200 300QMelt daily mod2.550100150200250300r=0.60 p=0.00 MBE=9.8Figure A.16: Comparison of AWS-derived (obs) versus WRF-derived (mod) daily time-series of the selected variables at Conrad glacier in 2015, station 2.72T2m obs0 5 10 15T 2m mod 2.5051015r=0.70 p=0.00 MBE=1.6RH2m obs50 100RH2m mod 2.5406080100r=0.56 p=0.00 MBE=-0.6U obs2 4 6U mod2.5246r=-0.06 p=0.17 MBE=-1.8SP obs775 785 795SP mod2.5775780785790795r=0.93 p=0.00 MBE=0.7Kin obs0 500 1000K in mod 2.505001000r=0.88 p=0.00 MBE=53.8Lin obs250 300 350L in mod 2.5250300350r=0.27 p=0.00 MBE=-24.1QH obs0 50 100QH mod 2.5050100r=0.02 p=0.71 MBE=-27.9QE obs0 40QE mod 2.5-200204060r=0.22 p=0.00 MBE=-6.0QMelt obs0 200 400 600QMelt mod 2.50200400600r=0.89 p=0.00 MBE=-9.4Figure A.17: Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Conrad glacier in 2016, station 1.73T2m daily obs2 6 10T 2m daily mod2.524681012r=0.81 p=0.00 MBE=1.6RH2m daily obs50 60 70 80RH2m daily mod2.550607080r=0.70 p=0.00 MBE=-0.6U daily obs1 2 3 4 5U daily mod2.512345r=-0.21 p=0.08 MBE=-1.8SP daily obs780 785 790SP daily mod2.5780785790r=0.94 p=0.00 MBE=0.7Kin daily obs100 200 300K in daily mod2.5100200300r=0.46 p=0.00 MBE=53.8Lin daily obs260 280 300 320L in daily mod2.5260280300320r=0.27 p=0.02 MBE=-24.1QH daily obs20 40 60 80QH daily mod2.520406080r=-0.06 p=0.60 MBE=-27.9QE daily obs0 10 20 30QE daily mod2.50102030r=0.25 p=0.03 MBE=-6.0QMelt daily obs100 200 300QMelt daily mod2.5100150200250300r=0.49 p=0.00 MBE=-9.4Figure A.18: Comparison of AWS-derived (obs) versus WRF-derived (mod) daily time-series of the selected variables at Conrad glacier in 2016, station 1.74T2m obs-5 0 5 10T 2m mod 2.5-50510r=0.84 p=0.00 MBE=0.4RH2m obs40 60 80 100RH2m mod 2.5406080100r=0.63 p=0.00 MBE=2.3U obs2 4 6 8U mod2.52468r=0.27 p=0.00 MBE=-1.9SP obs710 720 730 740SP mod2.5710720730740r=0.94 p=0.00 MBE=22.8Kin obs0 200 400 600 800K in mod 2.50200400600800r=0.86 p=0.00 MBE=38.8Lin obs200 250 300L in mod 2.5200250300r=0.29 p=0.00 MBE=-10.0QH obs-50 0 50Q H mod 2.5-50050r=0.48 p=0.00 MBE=-12.5QE obs-40 0 40Q E mod 2.5-60-40-2002040r=0.58 p=0.00 MBE=0.8QMelt obs0 100 200 300Q Melt mod 2.50100200300r=0.93 p=0.00 MBE=-7.5Figure A.19: Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Conrad glacier in 2016, station 2.75T2m daily obs0 4 8T 2m daily mod2.5-202468r=0.92 p=0.00 MBE=0.4RH2m daily obs60 80 100RH2m daily mod2.560708090100r=0.76 p=0.00 MBE=2.3U daily obs2 4 6U daily mod2.5246r=0.59 p=0.00 MBE=-1.9SP daily obs710 720 730 740SP daily mod2.5710720730740r=0.95 p=0.00 MBE=22.8Kin daily obs100 200 300K in daily mod2.5100200300r=0.32 p=0.01 MBE=38.8Lin daily obs240 260 280 300 320L in daily mod2.5240260280300320r=0.36 p=0.00 MBE=-10.0QH daily obs-40 0 40Q H daily mod2.5-40-200204060r=0.76 p=0.00 MBE=-12.5QE daily obs-40 -20 0 20Q E daily mod2.5-40-20020r=0.73 p=0.00 MBE=0.8QMelt daily obs20 60 100 140Q Melt daily mod2.52060100140r=0.89 p=0.00 MBE=-7.5Figure A.20: Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Conrad glacier in 2016, station 2.76T2m (° C) obs0 5 10T 2m (°  C) mod2.50510r=0.79 p=0.00 MBE=0.1RH2m (%) obs60 80 100RH2m (%) mod2.560708090100r=0.63 p=0.00 MBE=8.8U (m/s) obs1 2 3U (m/s) mod2.5123r=-0.20 p=0.06 MBE=-2.3SP (hpa) obs780 790 800SP (hpa) mod2.5780785790795800r=0.94 p=0.00 MBE=-19.1Kin (W/m2) obs0 200 400 600 800K in (W/m2 ) mod 2.50200400600800r=0.83 p=0.00 MBE=27.6Lin (W/m2) obs250 300 350L in (W/m2 ) mod 2.5250300350r=0.32 p=0.00 MBE=-23.1QH (W/m2) obs0 10 20 30Q H (W/m2 ) mod 2.50102030r=-0.16 p=0.16 MBE=-42.9QE (W/m2) obs-10 0 10 20 30Q E (W/m2 ) mod 2.5-100102030r=-0.07 p=0.54 MBE=-12.8QMelt (W/m2) obs0 200 400 600Q Melt (W/m2 ) mod 2.50200400600r=0.86 p=0.00 MBE=-45.1Figure A.21: Comparison of AWS-derived (obs) versus WRF-derived (mod) 3-hr time-series of the selected variables at Castle Creek glacier in 2010.77T2m (° C) daily obs2 4 6 8 10T2m (°  C) daily mod2.5246810r=0.92 p=0.00 MBE=-0.1RH2m (%) daily obs70 80 90 100RH2m (%) daily mod2.5708090100r=0.89 p=0.00 MBE=9.4U (m/s) daily obs1 1.5 2 2.5U (m/s) daily mod2.511.522.5r=-0.38 p=0.24 MBE=-2.3SP (hpa) daily obs780 790 800SP (hpa) daily mod2.5780785790795800r=0.95 p=0.00 MBE=-19.1Kin (W/m2) daily obs100 200 300Kin (W/m2 ) daily mod2.5100200300r=0.60 p=0.05 MBE=13.7Lin (W/m2) daily obs250 300 350Lin (W/m2 ) daily mod2.5250300350r=0.49 p=0.13 MBE=-24.5QH (W/m2) daily obs5 10 15QH (W/m2 ) daily mod2.551015r=-0.12 p=0.73 MBE=-42.3QE (W/m2) daily obs0 5 10 15QE (W/m2 ) daily mod2.5051015r=-0.01 p=0.97 MBE=-12.2QMelt (W/m2) daily obs50 100 150 200Q Melt (W/m2 ) daily mod2.550100150200r=0.61 p=0.05 MBE=-55.3Figure A.22: Comparison of AWS-derived (obs) versus WRF-derived (mod) daily time-series of the selected variables at Conrad glacier in 2010.78Table A.1: Evaluation results for mod1.0 (Castle Creek glacier only), mod2.5, mod7.5 and ERA-Interim against the AWS vari-ables: correlation (r), RMSE and MBE derived from daily timeseries. Bold values are correlations significantly differentfrom zero at a significance level of 0.05.7980Table A.2: ESA land category converted to USGS 24-category Land Use CategoriesA.2. Guidelines for running the WRF model on ComputeCanada computing facilitiesA.2.1. Preparation step1. Apply for a [WestGrid account].2. Choose among the available computing facilities such as Arbutus, Cedar, Orcinus etcbased on your requirement of available software (such as WRF), number of processors,memory and storage.3. The list of different systems with their specification can be found at [Westgrid systems].4. Once you choose the system, then use the ssh command to connect to that specific systemssh [username]@[systemname eg. orcinus].westgrid.ca or81ssh [username]@[systemname eg. cedar].computecanada.ca5. Check the availability of WRF model as well as other post processing tools such asNCL and NCO using the command module avail.. This command lists all the availablesoftware.6. In some cases, there is more than one version of WRF available, so choose the one thatfits your interest (recent versions are better due to the improvements and bug fixing thatis happening in each of the versions).7. Once you choose the WRF version to use, then create a directory where you want tocopy all the WRF files (e.g., WRF) related to the chosen version. The files includepreprocessing files such as geogrid, ungrib and metgrid; also model integration part suchas real.exe, wrf.exe.8. Change the directory to WRF using cd WRF9. Find the path for the wrf version that you want to copy usingwhich [wrf version eg. wrf3.8.1]which wrf (eg. For orcinus)10. The above command gives the path (eg. for orcinus in westgrid, the path is /global/software/wrf-3.8.1/bin/wrf), then copy the directory from that path using the command below. Eventhough this is the full path for the directory, the wrf file exists on the wrf-[version direc-tory] (e.g., /global/software/wrf-3.8.1/ for orcinus)cp -r [path] .cp -r /global/software/wrf-3.8.1/ . (for orcinus)11. The wrf directory copied from the above step should contain files such as bin, data,geogrid, metgrid and ungrib (both preprocessing and model integration files). The di-rectory bin contains all the executables and also some visualization tools such as ge-ogrid.exe, metgrid.exe, ungrib.exe, real.exe and wrf.exe12. Create two main directories, one for preprocessing (Eg. WPS-RUN) and the other forthe WRF model integration part (e.g., WRF-RUN).13. Move (copy) the geogrid, metgrid and ungrib directories to the WPS-RUN directory.Also copy everything from the bin directory except real.exe and wrf.exe.8214. Move (copy) the data directory except for WPS-geog-data (if there is any) also the exe-cutables real.exe and wrf.exe from the bin directory to the WRF-RUN directory.15. The important files in the WPS-RUN directory are the executables (geogrid.exe, un-grib.exe, metgrid.exe), WPS geographical input data and boundary input data.16. Create a directory for the wps geography input data (ex. geog) and boundary input data(e.g., DATA) under the WPS-RUN directory.17. In case of orcinus, wps geography input data is provided and it is located in the datadirectory and is named as WPS-geog-data. If this data is not provided, it can be down-loaded from [WPS-geog-data].18. The boundary input can be found from different sources but the commonly used data isEra-interim which can be downloaded from [ERA-Interim data].A.2.2. WRF Pre-processing System (WPS)1. Change the directory to WPS-RUN where all the preprocessing of the input data will beperformed.2. Load the specific wrf model that you want to use using the commandmodule load wrf3. In case of Cedar and Graham, wps and wrf should be loaded separately (when runningwps, load wps and vice versa).4. Specification for the model run such as time period, domain location, size and resolutionare specified through a namelist file called namelist.wps. The components of this file are&share (common for the three executables), &geogrid, &ungrib and &metgrid. Samplenamelist.wps file can be found from [sample namelist.wps files]. Example of basicnamelist.wps file is given below.&sharewrf_core = ’ARW’,max_dom = 3,start_date = ’2012-08-19_00:00:00’,’2012-08-19_00:00:00’,end_date = ’2012-09-19_00:00:00’,’2012-09-19_00:00:00’,83interval_seconds = 21600,io_form_geogrid = 2,debug_level = 0,/&geogridparent_id = 1, 1,parent_grid_ratio= 1, 5,i_parent_start=1, 50,j_parent_start=1, 50,e_we= 140, 176,e_sn= 140, 176,geog_data_res= ’10m’,’2m’,dx=25000,dy=25000,map_proj=’lambert’,ref_lat=53.0433,ref_lon=-120.4443,stand_lon=-120.4443,truelat1 = 60,truelat2 = 30,geog_data_path = ’./geog’/&ungribout_format = ’WPS’,prefix = ’ERA’,/&metgridfg_name = ’ERA’,io_form_metgrid = 2,/5. Modify the above namelist.wps file based on your domain information. Detailed expla-84nation about namelist. wps variables can be obtained from [Description of Namelist.wpsvariables].6. Before you start running the WPS programs, the domain setup can be viewed and ad-justed using the NCAR Command Language (NCL) visualization and scientific dataprocessing tool. Check if it is available on the remote server using the commandmodule avail.7. If it is available, then use the following command to load and start using itmodule load ncl_ncarg/6.3.0 (example is for orcinus)8. There are different NCL scripts that can be used for different purposes. The specificscript used to check the model set up (domain size and location) isplotgrids_new.ncl9. This script is usually provided with the wrf download file also can be can be downloadedfrom [NCL scripts].10. Use the following command to run the scriptncl plotgrids_new.ncl11. Run the three executables such as geogrid.exe, ungrib.exe and metgrid.exe. Order doesn’tmatter for geogrid.exe and ungrib.exe but metgrid.exe should be run after geogrid.exe andungrib.exe.A.2.2.1. Running ungribThe main function of the ungrib program is simply to unpack the required fields for wrfmodel run and writes them out into a format that the metgrid program can read. It doesn’thave dependency on model domain size. It uses Varable tables (Vtables) to identify thefields that should be unpacked from the given grib files. Each and every boundary inputdata has its own Vtable and it is usually provided along with the wrf download files. TheVtables are found in the ungrib directory. Example of a Vtable for Era-Interim pressurelevel data is named as Vtable.ERA-interim.pl.85(a) First link the correct Vtable using the following command (example is for Era-Interim data)ln -sf ungrib/Variable_Tables/Vtable.ERA-interim.plVtable(b) The file named Vtable should be created in the WPS-RUN directory.(c) The raw input data (DATA) can be viewed using the following programs (g1print.exefor grib1 and g2print.exe for grib2 data)g2print.exe \DATA\.. or g1print.exe \DATA\..(d) Second link the boundary input data (Ex. Era-Interim) to be processed by ungribusing the command./link_grib.sch ./DATA/* (use * when there are morethan one file in the DATA directory)(e) After running the above script, file of the form GRIBFILE.AAA, GRIBFILE.AAB,...should be created.(f) Run the program ungrib.exe directly from the command line asungrib.exe(g) The result from the above run should produce files of the form, for example WPS:2014-08-23 00, WPS:2014-08-23 06,...(Same number of files as the 6 hourly boundaryinput data). There is also a log file named ungrib.log that can be used to check ifthere are problems with the run.(h) The intermediate files produced by ungrib.ex can be viewed using the programrd_intermediate.exe [e.g., WPS:2014-08-23$\_$06](i) The above method works only for basic runs (e.g., for few days and coarse res-olution). When perfoming higher resolution run for longer periods a batch scriptshould be used. The first option is to use a script of the form [program name].pbswhen using Westgrid facilities such as orcinus.westgrid.ca. In this script, you canspecify the number of processors and also memory required for the specific run.See the following sample script for orcinus86#PBS -S /bin/bash#PBS -l pmem=2000mb#PBS -l nodes=2:ppn=4#PBS -l walltime=02:00:00#PBS -M [email address]#PBS -m beacd $PBS_O_WORKDIRmodule load wrfmpirun ungrib.exe(j) A slightly different script ([program name].sh) is used when using cedar and gra-ham with login node cedar.computecanada.ca and graham.computecanada.ca re-spectively. See the following sample script#!/bin/bash#TBATCH--account=def-[user name for the primaryaccount holder]#SBATCH --nodes=2#SBATCH --ntasks-per-node=32#SBATCH --mem=128000M#SBATCH --time=00-10:00#SBATCH --mail-user=[email address]#SBATCH --mail-type=BEGIN#SBATCH --mail-type=END#SBATCH --mail-type=FAIL#SBATCH --mail-type=REQUEUE#SBATCH --mail-type=ALLsrun ungrib.exe(k) Submit the script using the command qsub [Name of the script].pbs in case of(A.2.2.1 i) or sbatch [Name of the script].sh for (A.2.2.1 j).(l) The progress of the job can be checked using the command showq -u [user name]in case of (A.2.2.1 i) or [scc user name] for (A.2.2.1 j).(m) The output from this runs includes the files mentioned in (h) but the error files arenamed rsl.error.0000, rsl.error.0001,...(same number of error files as the number ofprocessors used).87A.2.2.2. Running geogridMain functions of geogrid are: defining map projection, geographic location and dimen-sion of domains, horizontally interpolating various terrestrial data sets to the model grids(e.g. Topography height, Land use category, Soil type, Vegetation fraction).(a) Run the program geogrid.exe from the command line or using the script from(A.2.2.1 i or A.2.2.1 j).(b) Output from the above run includes geo em.d01.nc, geo em.d02.nc,... (same num-ber of files as the number of domains). Also, an error file named geogrid.log andrsl.error.0000, rsl.error.0001,... if batch script is used.A.2.2.3. Running metgridMain function of metgrid is to horizontally interpolate meteorological data (extracted byungrib) to simulation domains (defined by geogrid).(a) Run the program metgrid.exe from the command line or using the script from(A.2.2.1 i or A.2.2.1 j).(b) Output from the above run includes met em.d01.2014-07-08 00:00:00.nc,...(samenumber of files as the 6 hourly boundary input file for each domain). Also, an errorfile named metgrid.log or rsl.error.0000, rsl.error.0001,... if batch script is used.A.2.3. Model integration Part (real.exe and wrf.exe)A.2.3.1. running real.exeMain functions of real.exe are reading boundary and geographical input data from WPS andinterpolate data to the model’s vertical coordinate.1. Change the directory to the WRF-RUN2. Similar to the namelist.wps file in section A.2.2, the wrf model integration part usesa namelist file named namelist.input. This file contains information such as time step,options of physics and dynamics schemes. Sample namelist.input file can be obtainedfrom [Sample namelist.input files]. Example of basic namelist.input file is shown below88&time_controlrun_days = 0,run_hours = 12,run_minutes = 0,run_seconds = 0,start_year = 2010, 2010, 2010,start_month = 07, 07, 07,start_day = 31, 31, 31,start_hour = 00, 00, 00,start_minute = 00, 00, 00,start_second = 00, 00, 00,end_year = 2010, 2010, 2010,end_month = 08, 08, 08,end_day = 05, 05, 05,end_hour = 00, 00, 00,end_minute = 00, 00, 00,end_second = 00, 00, 00,interval_seconds = 21600,input_from_file = .true.,.true.,.true.,history_interval = 360, 360, 180,frames_per_outfile = 1000, 1000, 1000,restart = .true.,restart_interval = 720,io_form_history = 2,io_form_restart = 2,io_form_input = 2,io_form_boundary = 2,debug_level = 0,auxinput4_inname = "wrflowinp_d<domain>",auxinput4_interval = 360,io_form_auxinput4 = 2,/&domainstime_step = 180,89time_step_fract_num = 0,time_step_fract_den = 1,max_dom = 3,s_we = 1, 1, 1,e_we = 140, 181, 211,s_sn = 1, 1, 1,e_sn = 140, 181, 211,s_vert = 1, 1, 1,e_vert = 30, 30, 30,p_top_requested = 5000,num_metgrid_levels = 38,num_metgrid_soil_levels = 4,dx = 33000, 11000, 3666.6667,dy = 33000, 11000, 3666.6667,grid_id = 1, 2, 3,parent_id = 1, 1, 2,i_parent_start = 1, 40, 60,j_parent_start = 1, 40, 60,parent_grid_ratio = 1, 3, 3,parent_time_step_ratio = 1, 3, 3,feedback = 0,smooth_option = 0sfcp_to_sfcp = .false.,/&physicsmp_physics = 10, 10, 10,ra_lw_physics = 1, 1, 1,ra_sw_physics = 2, 2, 2,radt = 33, 33, 33 ,sf_sfclay_physics = 2, 2, 2,sf_surface_physics = 2, 2, 2,bl_pbl_physics = 2, 2, 2,bldt = 0, 0, 0,90cu_physics = 3, 3, 3,cudt = 5, 5, 5,isfflx = 1,ifsnow = 1,icloud = 1,surface_input_source = 1,num_soil_layers = 4,num_land_cat = 24,sf_urban_physics = 0, 0, 0, 0,sst_update = 1,/&dynamicsw_damping = 1,diff_opt = 1, 1, 1,km_opt = 4, 4, 4,diff_6th_opt = 2, 2, 2,diff_6th_factor = 0.12, 0.12, 0.12,base_temp = 290.damp_opt = 3,zdamp = 5000., 5000., 5000.,dampcoef = 0.2, 0.2, 0.2,khdif = 0, 0, 0,kvdif = 0, 0, 0,non_hydrostatic = .true., .true., .true.,moist_adv_opt = 1, 1, 1,scalar_adv_opt = 1, 1, 1,/&bdy_controlspec_bdy_width = 5,spec_zone = 1,relax_zone = 4,specified = .true., .false.,.false.,nested = .false., .true., .true.,91/&namelist_quiltnio_tasks_per_group = 0,nio_groups = 1,/3. Description of all the namelist variables can be found at [Description of namelist.inputvariables].4. Modify the existing namelist.input file or create one based on the domain information.5. First, link output files from section A.2.2 i.e the met em.d01.2014-07-08 00:00:00.ncfiles to the WRF-RUN directory using the commandln sf ../WPS-RUN/met_em. *.6. Use the batch script from (A.2.2.1 i or A.2.2.1 j) but using more memories and numberof processors. See the following sample script.#PBS -S /bin/bash#PBS -l pmem=4000mb#PBS -l nodes=2:ppn=12#PBS -l walltime=12:00:00#PBS -M [email address]#PBS -m bea# cd PBS_O_WORKDIRmodule load wrfmpirun real.exe ( mpirun wrf.exe when running wrf.exe)7. The above run should produce files such as wrfbdy d01 (only for domain 1) and wrfinput d01,wrfinput d02,... (same number as the number of domains). Also error files such asrsl.error.0000, rsl.error.0001,....8. For longer runs, it is recommended to use [restart option] which can be specified atnamelist.input file before running real.exe.92A.2.3.2. Running wrf.exeMain function of wrf.exe is to read initial and boundary data from real.exe., read data formnamelist.input and generates the model forecast.1. Use the batch script from (A.2.2.1 i or A.2.2.1 j) to submit the job.2. Output from the model run looks like wrfout d01 2012-08-19 00:00:00, wrfout d02 2012-08-19 00:00:00,...(same as number of specified output interval over the given period foreach domain).3. There are a number of programs that can be used to visualize the wrf model output andalso produce different plots. The list of these programs and their features can be foundat [Post processing tools].4. Since the direct wrf model output takes lots of storage, it is possible to extract only vari-ables of interest using the NCL script called wrfout to cf.ncl which can be downloadedfrom [Re-sampling wrf output].A.2.4. Incorporating new data sets into the WRF modelIt is a common practice to use ones own land cover and topography data to run the WRF modelsince in most cases, the default data sets that come with the WRF model may not have the rightspatial resolution for the specific study of interest. The steps on how to incorporate land coveror topography data to the WRF model run is explained below.1. Download the Landcover [e.g.,ESA-300m] and Topography data e.g., SRTM-90m)2. Data usually comes as tif file so convert the data in to binary format (geogrid format)used by WRF for example using [Convert GeoTIFF to geogrid format used by WRF].3. Create an index file to define projection and dimension of data.4. Description of index files can be found at [Description of index Options].5. Sample index file for ESA land cover data is given below:projection = regular_lltype=categoricalcategory_min=1category_max=2493known_x = 1known_y = 9967known_lat = 69.500343known_lon = -142.727829dx = 2.775223e-03dy = 2.777125e-03signed = yesunits="category"description="24-category ESA"wordsize = 2tile_x = 100tile_y = 100tile_z = 1tile_bdr = 3missing_value = 0.000000scale_factor = 1.000000row_order = bottom_topendian = little6. Add entry for the data in the GEOGRID.TBL file. Example on how to add entry for theheight data (*) and land use data (**) to the GEOGRID.TBL file is shown below.===============================name = HGT_Mpriority = 1dest_type = continuoussmooth_option = smth-desmth_special; smooth_passes=1fill_missing=0.interp_option = gmted2010_30s:average_gcell(4.0)interp_option = gtopo_30s:average_gcell(4.0)interp_option = gtopo_2m:four_ptinterp_option = gtopo_5m:four_ptinterp_option = gtopo_10m:four_ptinterp_option = default:average_gcell(4.0)* interp_option = SRTM:average_gcell(4.0)rel_path = gmted2010_30s:topo_gmted2010_30s/94rel_path = gtopo_30s:topo_30s/rel_path = gtopo_2m:topo_2m/rel_path = gtopo_5m:topo_5m/rel_path = gtopo_10m:topo_10m/rel_path = default:topo_gmted2010_30s/* rel_path = SRTM:SRTM-90m/===============================name=LANDUSEFpriority=1dest_type=categoricalz_dim_name=land_catdominant = LU_INDEXlandmask_water = nlcd2006_9s:17landmask_water =nlcd2006_30s:17landmask_water = nlcd2011_9s:17landmask_water = nlcd2006:17landmask_water = ssib_10m:16landmask_water = ssib_5m:16landmask_water = modis_15s:17landmask_water = modis_30s:17landmask_water = usgs_30s:16landmask_water = usgs_lakes:16,28landmask_water = modis_lakes:17,21landmask_water = usgs_2m:16landmask_water = usgs_5m:16landmask_water = usgs_10m:16landmask_water = default:17,21interp_option = nlcd2006_9s:average_gcell(0.0)interp_option =nlcd2006_30s:average_gcell(0.0)interp_option = nlcd2011_9s:average_gcell(0.0)interp_option = nlcd2006:nearest_neighborinterp_option = ssib_10m:four_ptinterp_option = ssib_5m:four_ptinterp_option = modis_15s:nearest_neighborinterp_option = modis_30s:nearest_neighbor95interp_option = usgs_30s:nearest_neighborinterp_option = usgs_lakes:nearest_neighborinterp_option = modis_lakes:nearest_neighborinterp_option = usgs_2m:four_ptinterp_option = usgs_5m:four_ptinterp_option = usgs_10m:four_ptinterp_option = default:nearest_neighbor** interp_option = ESA:nearest_neighborrel_path = nlcd2006_9s:nlcd2006_ll_9s/rel_path =nlcd2006_30s:nlcd2006_ll_30s/rel_path = nlcd2011_9s:nlcd2011_ll_9s/rel_path = nlcd2006:nlcd2006_ll_30s/rel_path = ssib_10m:ssib_landuse_10m/rel_path = ssib_5m:ssib_landuse_5m/rel_path = modis_15s:modis_landuse_20class_15s/rel_path = modis_30s:modis_landuse_20class_30s/rel_path = usgs_30s:landuse_30s/rel_path = usgs_lakes:landuse_30s_with_lakes/rel_path = modis_lakes:modis_landuse_21class_30s/rel_path = usgs_2m:landuse_2m/rel_path = usgs_5m:landuse_5m/rel_path = usgs_10m:landuse_10m/rel_path = default:modis_landuse_21class_30s/** rel_path = ESA:ESA-300m/===============================7. Specify the the two data sets in the namelist.wps file geogrid section (see *** below)otherwisethe default topography and land use data will be used.&geogridparent_id = 1, 1, 2,parent_grid_ratio= 1, 3, 3,i_parent_start=1, 40, 57,j_parent_start=1, 40, 53,e_we= 140, 211, 241,e_sn= 140, 211, 241,96geog_data_res= ’SRTM+ESA’,’SRTM+ESA’,’SRTM+ESA’, (***)dx=22500,dy=22500,map_proj=’lambert’,ref_lat=53.0433,ref_lon=-120.4443,stand_lon=-120.4443,truelat1 = 60,truelat2 = 30,geog_data_path = ’/home/mekdes12/scratch/WRF/my_wps_run/geog’/&ungribout_format = ’WPS’,prefix = ’ERA’,/&metgridfg_name = ’ERA’,io_form_metgrid = 2,/8. Then run WRF following the same step as section A.2.A.2.5. Post processing tools• For post processing of the WRF model output you can use IDV, NCL scripts as well asMatlab. Other post processing tools can be found at [Post processing utilities].• IDV is one of the tool which enable visualization of different variables and can be down-loaded from [IDV download].• To re sample the output variables from the original wrf output, the following ncl scriptcan be used. [Re-sampling wrf model output].97

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0368955/manifest

Comment

Related Items