Arts, Faculty of
Geography, Department of
DSpace
UBCV
Chai, Zijie
Master of Science - MSc
University of British Columbia
Glaciers make significant contributions to streamflow in western North America. However, glaciers throughout the region are losing mass and retreating, leading to concerns about summer streamflow, particularly during extended spells of warm, dry weather. An empirical analysis of climatic and streamflow trends in the Canadian portion of the Columbia River basin found a regional climate-driven decline in August streamflow, which was augmented in glacier-fed catchments by glacier retreat. The objective of this study was to use a hydrological model to quantify the relative roles of glacier retreat and climatic trends on streamflow in Illecillewaet catchment. We applied a modified version of HBV-EC within the Raven modelling framework to Illecillewaet River, located in the Columbia River headwaters, from 1977 to 2020. The model was calibrated using daily streamflow, snow water equivalent from snow pillows and geodetic glacier mass change. The model was run for 1977-2020 with four scenarios: (a) declining glacier cover based on historical imagery; (b) constant glacier extent with 1985 glacier cover; (c) constant glacier extent with 2016 glacier cover and (d) no glacier cover. Key conclusions are: (1) in the Illecillewaet catchment, where glacier cover represents about 5% of the total area, glacier ice melt can contribute between 3.5% and 41.1% to the August streamflow; (2) the analysis of glacier ice melt trends showed a small but statistically insignificant increase, attributed to the synchronous rising of the snowline along with the retreat of the glacier terminus; and (3) the decline in August streamflow appears to have been primarily driven by climate influences, particularly decreasing precipitation in August. Currently, the presence of glaciers helps offset this declining trend, but the glacier compensation effect is weakening as the glaciers retreat.
Diagnosing the relative contributions of glacier retreat andclimatic variability to streamflow trends in IllecillewaetcatchmentbyZijie ChaiB.Sc., Hohai University, 2020A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THEREQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Geography)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2023© Zijie Chai, 2023The following individuals certify that they have read, and recommend to the Faculty ofGraduate and Postdoctoral Studies for acceptance, the thesis entitled:Diagnosing the relative contributions of glacier retreat and climatic variability tostreamflow trends in Illecillewaet catchmentsubmitted by Zijie Chai in partial fulfillment of the requirements forthe degree of Master of Sciencein GeographyExamining Committee:Robert Daniel Moore, Professor, Geography, UBCSupervisorGeorg Jost, Adjunct Professor, Geography, UBCSupervisory Committee MemberBrett Eaton, Professor, Geography, UBCAdditional examinerAdditional Supervisory Committee Members:Steven Weijs, Associate Professor, Civil Engineering, UBCSupervisory Committee MemberiiGlaciers make significant contributions to streamflow in western North America. However,glaciers throughout the region are losing mass and retreating, leading to concerns aboutsummer streamflow, particularly during extended spells of warm, dry weather. An empir-ical analysis of climatic and streamflow trends in the Canadian portion of the ColumbiaRiver basin found a regional climate-driven decline in August streamflow, which was aug-mented in glacier-fed catchments by glacier retreat. The objective of this study was to usea hydrological model to quantify the relative roles of glacier retreat and climatic trends onstreamflow in Illecillewaet catchment. We applied a modified version of HBV-EC withinthe Raven modelling framework to Illecillewaet River, located in the Columbia River head-waters, from 1977 to 2020. The model was calibrated using daily streamflow, snow waterequivalent from snow pillows and geodetic glacier mass change. The model was run for1977-2020 with four scenarios: (a) declining glacier cover based on historical imagery; (b)constant glacier extent with 1985 glacier cover; (c) constant glacier extent with 2016 glaciercover and (d) no glacier cover. Key conclusions are: (1) in the Illecillewaet catchment,where glacier cover represents about 5% of the total area, glacier ice melt can contributebetween 3.5% and 41.1% to the August streamflow; (2) the analysis of glacier ice melttrends showed a small but statistically insignificant increase, attributed to the synchronousrising of the snowline along with the retreat of the glacier terminus; and (3) the decline inAugust streamflow appears to have been primarily driven by climate influences, particu-larly decreasing precipitation in August. Currently, the presence of glaciers helps offset thisAbstractdeclining trend, but the glacier compensation effect is weakening as the glaciers retreat.iiiGlaciers in western North America significantly contribute to streamflow, but they arelosing mass and retreating due to climate change. This study used a hydrologic model toanalyze the trend of summer streamflow and roles of glacier melt and changing climate inthe Illecillewaet catchment in the Canadian Columbia River basin. The findings revealedthat over the past 40 years, there has been a decrease in summer rainfall. Fortunately, thepresence of glaciers has helped offset this drier trend by contributing to river flow. But, asthe glaciers continue to retreat, the rivers in the catchment will experience further decreasesLay Summaryin water flow, posing challenges for summer water availability in the future.ivThis MSc thesis represents my original research, with invaluable guidance from my super-visory committee: Dan Moore, Georg Jost, and Steven Weijs. I am also deeply grateful forPrefacethe help provided by Ben Pelto and Matt Chernos.A portion of this research was presented as a poster titled ”Contributions of glaciersto streamflow in Illecillewaet catchment and how they have varied in response to climaticvariability and glacier retreat” at the 2023 Canadian Geophysical Union (CGU).vAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiAcknowledgementsTable of Contents. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation for the study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2.1 Overview of glacier contributions to streamflow . . . . . . . . . . . . 11.2.2 Peak water concept . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2.3 Approaches for quantifying glacier contributions to streamflow . . . 31.2.4 Model-based approaches . . . . . . . . . . . . . . . . . . . . . . . . . 41.3 Research questions and thesis structure . . . . . . . . . . . . . . . . . . . . 62 Study area and methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.1 Study area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.1.1 Physiography and land cover . . . . . . . . . . . . . . . . . . . . . . 72.1.2 Hydroclimate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.2 Data Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2.1 Meteorological data . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2.2 Streamflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.2.3 Snowpack water equivalent . . . . . . . . . . . . . . . . . . . . . . . 112.2.4 Spatial data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.3 Hydrological modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.3.1 Classification of hydrologic response units . . . . . . . . . . . . . . . 132.3.2 Model description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.3.3 Calibration and validation procedure . . . . . . . . . . . . . . . . . . 172.3.4 Extraction of contributions associated with rainfall, snowmelt and icemelt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.3.5 Extraction of glacier mass change . . . . . . . . . . . . . . . . . . . . 182.4 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.1 Split sample model calibration and validation . . . . . . . . . . . . . . . . . 203.1.1 Streamflow simulation . . . . . . . . . . . . . . . . . . . . . . . . . . 20vi3.1.2 Snowpack water equivalent . . . . . . . . . . . . . . . . . . . . . . . 213.1.3 Glacier mass change . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.2 Full period model calibration and performance . . . . . . . . . . . . . . . . 253.2.1 Streamflow simulation . . . . . . . . . . . . . . . . . . . . . . . . . . 253.2.2 Snowpack water equivalent . . . . . . . . . . . . . . . . . . . . . . . 273.2.3 Glacier mass change . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.3 Trends in observed climate and streamflow data . . . . . . . . . . . . . . . . 283.3.1 Climate trend . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.3.2 Streamflow trend . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.4 Ice melt contributions to streamflow . . . . . . . . . . . . . . . . . . . . . . 293.5 Trend analyses of monthly streamflow for four scenarios . . . . . . . . . . . 343.6 Trend analyses of monthly ice melt contributions for three glacier-cover sce-narios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 344 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.1 Utility of using Raven for trend attribution studies . . . . . . . . . . . . . . 384.2 Validity of the model results . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.3 Glacier melt contributions to streamflow in Illecillewaet catchment . . . . . 394.4 Streamflow trends and attribution . . . . . . . . . . . . . . . . . . . . . . . 404.5 Implications for the peak water concept . . . . . . . . . . . . . . . . . . . . 415 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 435.1 Key findings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 435.2 Avenues for future research . . . . . . . . . . . . . . . . . . . . . . . . . . . 43References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50Appendix A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50viiTable 2.1 BEC zones in Illecillewaet catchment. . . . . . . . . . . . . . . . . . . 9Table 2.2 Characteristics of weather stations used in this study. . . . . . . . . . 11Table 2.3 Snow courses in Illecillewaet catchment. . . . . . . . . . . . . . . . . . 12Table 2.4 Selected glacier outlines used for model set-up. . . . . . . . . . . . . . 13Table 2.5 Aspect class.List of Tables. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15Table 3.1 Calibrated parameters for full period calibration and split sample cal-ibration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20Table 3.2 Observed and modelled Illecillewaet catchment glacier mass changefor two time periods. The values are expressed as a water equivalentdepth averaged over the catchment area. . . . . . . . . . . . . . . . . 21Table 3.3 Observed and modelled glacier mass change for two time periods inIllecillewaet catchment, expressed as an average depth of water overthe catchment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28Table 3.4 Theil-Sen slope for temperature and precipitation time series. Thebold and italics values indicate that the Spearman correlation coeffi-cient is significantly different from 0 (p value less than 0.05 and 0.1respectively). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28Table 3.5 Trend analyses for summer streamflow for Illecillewaet River. (bold issignificant in Spearman correlation test, with p value less than 0.05). 28Table 3.6 August ice melt contribution to streamflow, July to September averagedaily temperature, average monthly precipitation, and December toApril average monthly precipitation for 1990, 2003 and 2011. . . . . . 29Table 3.7 Trend analyses for July streamflow for Illecillewaet River. . . . . . . 35Table 3.8 Trend analyses for August streamflow for Illecillewaet River. The boldvalues indicate that the Spearman correlation coefficient is significantlydifferent from 0 (p value less than 0.05). . . . . . . . . . . . . . . . . 35Table 3.9 Trend analyses for September streamflow for Illecillewaet River. . . . 35viiiList of FiguresFigure 2.1 Map of Illecillewaet catchment, including glaciers, the locations ofweather stations (RGR and YRV), and the Water Survey of Canadagauging station. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7Figure 2.2 BEC zones and NDVI class in Illecillewaet catchment. The thresholdindex for bare and vegetated is 0.2. Acquisition and processing BECzone and NDVI data are described in section 2.2.4. . . . . . . . . . . 8Figure 2.3 Climograph for Revelstoke Airport (YRV) and Glacier NP Rogers Pass(RGR) weather stations from 1977 to 2020. The red lines are themonthly average temperature. Blue bars are the monthly total pre-cipitation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10Figure 2.4 Illecillewaet River daily streamflow from 1977 to 2020. Line coloursvary by water year. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11Figure 2.5 The outline to generate HRUs in this study. . . . . . . . . . . . . . . 14Figure 3.1 Observed and simulated hydrographs in randomly selected years in thevalidation period: 2011, 2013, 2016, 2018. . . . . . . . . . . . . . . . . 22Figure 3.2 Scatterplots of simulated versus observed SWE in the calibration andvalidation periods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23Figure 3.3 Snow course locations and the centroids for the corresponding HRUs. 24Figure 3.4 Observed and simulated hydrographs in randomly selected years: 1984,1994, 2006 and 2016. . . . . . . . . . . . . . . . . . . . . . . . . . . . 25Figure 3.5 Monthly model result. . . . . . . . . . . . . . . . . . . . . . . . . . . . 26Figure 3.6 Scatterplots of simulated versus observed SWE. . . . . . . . . . . . . 27Figure 3.7 Temperature and precipitation time series for both RGR and YRVstations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30Figure 3.8 Streamflow time series for July, August and September. . . . . . . . . 31Figure 3.9 Boxplots of ice melt contributions to streamflow by month for May,June, July, August, September and October. The upper panel showsabsolute contributions and the lower panel shows percentage contri-butions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32Figure 3.10Streamflow contributions from glacier ice melt, rainfall and snow meltfor July, August and September. . . . . . . . . . . . . . . . . . . . . . 33Figure 3.11Streamflow time series for four scenarios in July, August and September. 34Figure 3.12Time series of ice-melt contributions to streamflow for three scenariosin July, August and September. . . . . . . . . . . . . . . . . . . . . . 36ixFigure 3.13Time series of modelled area of glacier snow cover (average of Augustfrom 1977 to 2020) and observed glacier areas in 1985, 1998, 2007 and2016. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37Figure A.1 LAI distributions in the Illecillewaet catchment on 20th August 2018(left) and 2019 (right). . . . . . . . . . . . . . . . . . . . . . . . . . . 50Figure A.2 Boxplot of LAI values for each land class on 20th August 2018 (top)and 2019 (bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51xAcknowledgementsI am deeply grateful to all the individuals who made this research possible.First and foremost, I’m really thankful to my supervisor, Dan Moore. Without hisinvaluable guidance, this work would not have happened. He introduced me to the world ofglacial hydrology and guided me in R programming as a powerful research tool. Throughoutthis research journey, he patiently imparted his knowledge in data collection, hydrologicmodelling, analysis of the results, and how to write academically. He also organized a labfield trip that allowed me to see the beauty of my study area, the Illecillewaet catchment.I also acknowledge the opportunity to attend CGU 2023 in Banff, where I presented myposter and connected with other wonderful researchers.I’m also grateful to Georg Jost, who is on my committee. He provided expertise inmodel setup and calibration and gave me helpful feedback on my thesis writing. I extendmy gratitude to Brett Eaton, my external examiner, for his valuable comments during thethesis writing process.The financial support for this research came through a Discovery Grant to ProfessorDan Moore from Canada’s Natural Sciences and Engineering Research Council (NSERC).Additionally, I am supported by a Graduate Support Initiative Award from the UBC Facultyof Arts as well as an International Tuition Award.Special thanks go to Ben Pelto for always being willing to lend a helping hand andprovide valuable advice on resolving technical challenges, especially concerning Raven andOstrich setup, as well as data and methods related to glacier mass balance in this study.I am also grateful to Matt Chernos and Joel Trubilowicz for their helpful assistance andinsightful modeling suggestions.I would also like to express my thanks to the College of Hydrology and Water Resources,Hohai University, where I pursued my bachelor’s degree and found my passion in naturalscience and hydrology.My sincere thanks extend to my lab mates Sierra Collins, Lili Callahan and Billy Brown-ing. Together, we are a supportive and encouraging team during the past two years. Myfriends Albert Hao, Tin Satriawan, Andrea Ku, Tzu-Yi Lu, Tony Zhang, Nastaran Pourne-matollahi, and Ling Han also give me so much encouragement and support and helped memaintain a balanced life amidst the demands of work.Lastly, I thank my family for their love, support, and patience over the years. Theirencouragement has motivated me to follow my passions and interests.xi1 Introduction1.1 Motivation for the studyGlaciers and ice sheets play a significant role in the water cycle, storing water in the form ofsnow and ice during the cold season, and releasing it as meltwater during the warm season.They store about 70% of the world’s freshwater, of which nearly 97% is in Antarctica andGreenland (Lemke et al., 2007). While mountain glaciers only account for less than 0.05% ofthe total volume of ice, they are critical to contribute to and modulate streamflow, whichaffects freshwater availability, hydropower operations, and aquatic ecosystems. Glaciersare experiencing a generalized mass loss worldwide, with an enhanced retreat rate in recentyears, mainly due to climate warming.Climate variability and change can influence glacier-melt contributions to streamflow bydirect effects on energy inputs and glacial water flowpaths over sub-daily to seasonal timescales. Over longer time scales, sustained glacier retreat can limit glacier runoff by reducingthe area available for melt. Reductions in glacier runoff not only produce lower streamflow,but can result in higher water temperatures, with potentially negative implications for cold-and cool-water aquatic species like salmonids. Therefore, quantifying the effects of climatevariability and change on glacier runoff and contributions to downstream river discharge iscritical for managing water supply and aquatic resources.The broad goal of this thesis is to quantify the relative roles of climatic variabilityand glacier retreat on recent streamflow variability for Illecillewaet River, located in theColumbia Mountains. The next section provides a literature review that focuses on identi-fying key knowledge gaps, which will provide the basis for the specific research objectives,which are provided in Section 1.3.1.2 Background1.2.1 Overview of glacier contributions to streamflowIn many mountain regions, large volumes of precipitation are stored on glaciers in winterand released by summer melting, which produces a pronounced seasonality in discharge,in terms of both magnitude and the nature of diel fluctuations. Early in the melt season,glaciers will typically be entirely snow-covered. The relatively high albedo of snow willreduce melt rates. The rate of water flow through the glacier will tend to be low due toa combination of deep snow with a relatively low hydraulic conductivity and the lack of a1well-connected subglacial drainage network (Hock, 2005). As the melt season progresses,snow grain size enlarges, which increases the hydraulic conductivity of the snowpack anddecreases the travel time (Stenborg, 1970; Fountain and Walder, 1998). This effect actsin concert with decreasing snowpack depth to reduce the lag time between melting andthe arrival of the melt water at the base of the snow pack. The melt rate also tends toincrease, which also results in higher percolation rates due to the relation between pore-water velocity and flux rate (Fountain and Walder, 1998). In addition, the englacial andsubglacial drainage systems develop and become more efficient during the melt season(Jansson et al., 2003). All of these influences reduce the water travel time associated withvertical and lateral channels.The glacier snowline plays a key role in controlling the timing, magnitude and diel vari-ability of glacier melt contributions to streamflow. The position of the snowline influencesthe rate of meltwater generation because glacier ice has a lower albedo than seasonal snowor firn, and thus absorbs more incident solar radiation and melts faster once it is exposed.In addition, the glacier runoff process differs between snow-covered and snow-free areas. Insnow-covered areas, melt water and rain experience a lag time due to vertical percolation inunsaturated snow and firn, and then flows downslope in a saturated layer that forms abovethe relatively impermeable glacier ice at the base of the snow and firn. The water in thissaturated layer flows downslope as Darcian flow similar to downslope flow in a permeablesoil underlain by relatively impermeable bedrock. Both processes can cause substantiallag times between the time of melt and the time at which the water exits the glacier tobecome streamflow. By contrast, the flow over the bare ice surface is orders of magnitudefaster than the flow in snow or firn. Consequently, as the snowline rises, the greater areaof exposed bare ice results in both increased rates of melt water generation and a morerapid conveyance of the water to the glacier’s outlet stream (Stenborg, 1970; Fountain andTangborn, 1985; Jansson et al., 2003).At an annual time scale, glacier runoff is dominated by summer melt, and exhibits apositive correlation with summer air temperature and a negative correlation with winterprecipitation (Fountain and Tangborn, 1985; Rothlisberger and Lang, 1987). This charac-teristic can moderate interannual streamflow variability, particularly in highly glacierizedbasins, which is also called the glacier compensation effect (GCE) (Fountain and Tang-born, 1985). This compensation effect means that, during warm and dry periods, the meltfrom glacierized areas could compensate for the lack of precipitation in non-glacierizedarea (Koboltschnig et al., 2009). On the other hand, during cold and wet periods, runofffrom the non-glacierized areas could compensate for the lack of melt from the glacierizedpart (Hopkinson and Young, 1998). Though some previous studies suggested a concave-uprelationship between streamflow variability and glacier cover of different glacierized catch-ments, the real situation can be more complex and not show a specific relationship, asdemonstrated by Tiel et al. (2020).21.2.2 Peak water conceptThe peak water concept posits that glacier runoff response to a shift from a neutral tochronically negative mass balance exhibits a sequence of four phases. The first phaseinvolves enhanced streamflow associated with a more rapid rise of the snowline, which inturn results in greater ice exposure and thus melt water generation. In the second stage,glacier runoff continues to increase until reaching a maximum point, also known as “peakwater,” when the loss of glacier area by terminal retreat counteracts the effect of highersnowlines, longer melt seasons, and higher melt rates. The third phase is marked by adeclining trend in runoff. Less meltwater is generated because the rate of glacier retreat ishigher than the rate of ice melt. Finally, in the fourth phase, streamflow continues to declineas the contribution of glacier wastage to streamflow decreases until the catchment reachesa new glacio-climatic equilibrium or shifts to a non-glacial streamflow regime. These fourphases of glacier runoff response are highly simplified, and the actual response to climatechanges may be more complex (Jansson et al., 2003; Moore et al., 2009; Frans et al.,2018).Casassa et al. (2009) reviewed sixteen data sets from basins with glacier retreat, andfound decreasing streamflow in basins with low elevation and/or low glacier cover, such asin south-central British Columbia, and increasing streamflow tended to be found in higherelevation catchments and/or catchments with high glacierization, such as the TianshanMountains and Tibet in central Asia. The research by Huss and Hock (2018) showedsimilar results: in roughly half of the 56 studied basins from 1980 to 2100, the modelledannual glacier runoff continued to rise until the peak water was reached, after which runoffbegan to decline. In the remaining basins, this maximum point had already been passed.1.2.3 Approaches for quantifying glacier contributions to streamflowDifferent concepts have been used to quantify glacier contributions to streamflow in previousstudies, which can cause ambiguity and lack of comparability of published results. First,in the most general sense, glacier runoff is defined as the runoff from the glacierized area,including the portion of all water inputs to the glacier through melt, rain, or other surfaceor sub-glacier inflow exit at the glacier terminus. The second concept describes glaciercontributions as coming from only the melting of the glacier, including glacier ice, firn andsnow, hence excluding any rainwater or other inflow to the glacier system. Third, glacierrunoff is sometimes defined as the meltwater originating solely from ice and firn melt, andsnow melt on the glacier is calculated separately. This definition is consistent with the viewthat all other components (snow melt, rain, etc.) would exist even if the glacier did notexist. In conclusion, the difference among these three concepts is associated with whetherthe contributions include snow melt and rain (Frenierre and Mark, 2014).Several data-based approaches have been applied in previous studies, in which thestreamflow contributions from glaciers are quantified by methods such as direct dischargemeasurement at glacier outlets (Thayyen and Gergan, 2010; Gascoin et al., 2011), hydro-logical balance approach (Baraer et al., 2009), and hydrochemical tracer approaches (Mark3et al., 2005; Nolin et al., 2010). The most frequently used method is the glaciological ap-proach, which compares the volume of ice melting from glaciers (as water equivalent) tostreamflow measured downstream (Huss, 2011; Pelto, 2011).In catchments where ablation stakes are adequate, mass balance data can be obtaineddirectly from ablation stake networks. Zemp et al. (2009) reviewed the present globalablation monitoring network, and conducted a time series analysis of mass balance for 30glaciers (mainly in Europe and North America) with continuous observation series since1976, which showed an accelerated ice loss from 1976 to 2005. However, this approach canintroduce substantial uncertainty, because the stake networks are often limited to accessibleareas of the glacier surface, which normally cannot represent the accumulation and abla-tion rates of the entire glacier, given the highly heterogeneous mountain terrain. Anothermethod to estimate ice melt is comparison of multi-temporal remote sensing products suchas digital elevation models (DEMs), which can be used to compute changes in ice volumewithin a catchment. For example, Moore et al. (2020) estimated the wastage flux contri-butions from glacier volume changes determined for the periods 1985–1999 and 2000–2018.A limitation is that repeated DEM mapping may only be available over periods of a decadeor more, so the method is not suitable for yearly analysis.Alternatively, empirical analysis of the contrasting responses of glacier-fed and unglacier-ized catchments can provide insight into the effect of glacier runoff on downstream discharge(Stahl et al., 2006). Multiple linear regression is also a commonly used empirical methodto assess the attribution of streamflow trends, in which the air temperature term can beassociated with the magnitude of snow and ice melt and the residuals are related to theeffects of changes within a catchment, such as glacier retreat (Duethmann et al., 2015;Moore et al., 2020).Glacier runoff is sometimes defined as the runoff component that is due to glacier netmass loss, referring only to the water originating from the glacier volume change. Thisquantity is often called the “wastage flux.” Hence, a glacier only affects runoff if there is anet mass loss during the considered time period. In this case, the glacier runoff is equivalentto the glacier’s (negative) mass budget (Radic and Hock, 2014).1.2.4 Model-based approachesHydrological modelling is a frequently used approach. Energy balance or temperature-indexmodels are used to simulate water released by melting snow, firn and ice, firn and ice, whichcan represent the glacier contributions to streamflow (according to different definitions).Energy balance models involve the assessment of energy fluxes towards and from the icesurface, including net radiation, sensible heat flux, latent heat flux and heat conductionbelow the ice or snow surface. While such models are, in principle, more accurate and trans-ferable in time and space, they require data such as solar and longwave radiation that aretypically not available outside of research settings. By contrast, temperature-index models(degree-day models) are based on the empirical relationship between air temperature andsnow or ice melt. Data for these models are much easier to access, mainly air temperature4and stream discharge, allowing their wide use across the world. In many instances, degree-day models have been used for estimating the current and future contribution of glaciermelt water (Braithwaite, 1995; Hock, 2003; Shea et al., 2009). Jost et al. (2012) estimatedglacier contributions to streamflow by running a calibrated temperature-index model withand without glaciers. In the no-glacier runs, all glacier cover was converted to open land,because snowmelt and rainfall runoff in areas covered by glaciers would still be there evenif the glaciers disappeared.The modelling approach has been frequently used and improved over time. Jost et al.(2012) illustrated the significance of accounting for glacier area change in the hydrologymodel, and also using glacier mass changes to assist in parameter estimation. Stahl et al.(2008) simulated the effects of glacier retreat on streamflow in three climate scenarios bycoupling a semi-distributed hydrological model with a volume-area scaling relation model.Frans et al. (2018) evaluated long-term changes in the contribution of glacier melt to theHood River system using a spatially distributed hydrological model coupled with a glacierdynamics model for 184 years. Glacier thickness change is another factor that affects theresults of models, Huss and Hock (2015) applied a ΔH parameterization method on a globalscale to account for the effect of elevation changes through time due to glacier thinning. Inaddition to glacier change, Carnahan et al. (2019) considered vegetation succession whenassessing impacts of glacier mass loss on water resources. They demonstrated that thetime at which basin runoff falls below pre-retreat levels is heavily influenced by the rate ofvegetation establishment following the loss of glacier ice.Though the hydrological modelling approach is relatively mature, numerous sources ofuncertainty remain, including those associated with model structure and parameter uncer-tainty. Equifinality is a common source of model uncertainty, which refers to the potentialfor different combinations of parameter sets to produce comparably good model outputs(Beven, 2006). As the number of calibration parameters increases, the risk of equifinalityincreases (Gao et al., 2012). In addition, previous studies suggested that this uncertaintyis especially obvious when glacier mass balance processes are poorly parameterized, sinceother parameters are likely to be erroneously adjusted to better match measured discharge(Frenierre and Mark, 2014). Using several different data sources rather than single dis-charge to calibrate important model parameters helps to constrain process representationand minimize equifinality, as this approach provides a better assessment of model consis-tency (Frenierre and Mark, 2014; Chernos et al., 2020). Jost et al. (2012) used glacier vol-ume change estimated by a geodetic approach in the calibration procedure, which helped toensure the accuracy of simulated glacier mass balance and reduced parameter uncertainty.Alternatively, application of stepwise calibration (Tiel et al., 2020) and extension of thecalibration period (Wheler et al., 2014) also help to minimize equifinality.Some studies have involved combined approaches to increase reliability and confidence ofthe results. A complementary application of both empirical trend analysis and a conceptual-parametric hydrological model was used by Moyer et al. (2016). The empirical analysisindicated no trend in streamflow, but the model-based approach indicated that streamflowhad declined relative to what it would have been had the glacier not retreated in BridgeRiver catchment. Duethmann et al. (2015) analyzed the attribution of increased streamflow5in Tarim River (Central Asia) by applying multiple linear regression analysis as well as usingsimulation experiments with a constant glacier geometry and detrended temperature andprecipitation time series. They found that in Kakshaal catchment, both temperature andprecipitation are the dominant drivers, while increasing temperature and related glacierretreat played a major role in Sari-Djaz catchment.In the Canadian portion of the Columbia River basin, previous empirical studies basedon a spatial comparison approach suggested that contributions from glacier melt to Augustrunoff have already passed peak water, and decreased through time at a rate that waslinearly related to the loss of glacier area (Moore et al., 2020). However, few studies haveused model-based approaches as a complementary implement to support the attribution oftrends identified through data-based analysis.1.3 Research questions and thesis structureThe review conducted in Section 1.2 identified the utility of combining data-based andmodel-based approaches to improve our understanding of streamflow response to climaticvariability and change in glacierized catchments, but also showed that it has not been widelyused. The main objective of this study is to employ a hydrological model to quantify therespective influences of glacier retreat and climatic trends on streamflow in the Illecillewaetcatchment in Columbia River basin from 1977 to 2020, a period for which a spatial com-parison approach suggested that glacier retreat had exacerbated climate-driven trends tolower August streamflow. To address this objective, the thesis aims to answer the followingresearch questions:1. What are the trends in summer (July to September) streamflow and climatic condi-tions in the Illecillewaet catchment?2. What is the trend in simulated ice melt during the summer months, and at whatstage is the Illecillewaet catchment in terms of the “peak water” cycle?3. What are the relative influences of glacier retreat and climatic variability on stream-flow trends?4. To what extent are model-based results consistent with a data-based spatial compar-ison approach?This thesis is structured as follows: Chapter 2 provides an overview of the study areaand a detailed description of data collection, model setup, and analysis procedures. Chapter3 presents the findings of the study, while Chapter 4 discusses how these results addressthe primary research questions. Finally, Chapter 5 summarizes the key conclusions drawnfrom the study and suggests potential avenues for future research.62 Study area and methodology2.1 Study area2.1.1 Physiography and land coverThe Illecillewaet catchment is located in the Columbia River Basin, which is a 668,000km2 transboundary watershed that spans seven U.S. states and the southeastern portion ofBritish Columbia, Canada. The Canadian portion of the Columbia River basin (top rightplot on Figure 2.1) is dominated by the Columbia Mountains, which are subdivided into theCariboo, Monashee, Selkirk, and Purcell ranges, with the Rocky Mountains defining theeastern slope. The Illecillewaet catchment is located in the Selkirk ranges of the ColumbiaRiver Basin (Figure 2.1) and covers an area of approximately 1151 km2 at the Water Surveyof Canada gauging station.Figure 2.1: Map of Illecillewaet catchment, including glaciers, the locations of weather sta-tions (RGR and YRV), and the Water Survey of Canada gauging station.7The Illecillewaet catchment ranges in elevation from 450 m above sea level (m a.s.l.) inthe valley bottoms to nearly 3300 m a.s.l. at the mountain peaks. Mount Sir Donald (3284m) and Uto Peak (3005 m) stand as prominent features in the landscape. Illecillewaet Riveroriginates from the Illecillewaet Glacier located in Glacier National Park and traverses adistance of 62 km before flowing into the Arrow Lakes. The catchment has been shapedby glaciers over millions of years, resulting in the formation of numerous glacially carvedfeatures such as cirques and U-shaped valleys.The Illecillewaet catchment comprises three Biogeoclimatic Ecosystem Classification(BEC) zones: the Interior Cedar Hemlock Zone (ICH), Engelmann Spruce - Subalpine FirZone (ESSF), and Interior Mountain-heather Alpine Zone (IMA) (Pojar et al., 1987). Table2.1 provides a summary of the characteristics of each zone in the catchment. Figure 2.2shows the distribution of the three BEC zones in Illecillewaet catchment and a NormalizedDifference Vegetation Index (NDVI) distribution in the catchment, with 0.2 as the thresholdto distinguish bare and vegetated areas.Figure 2.2: BEC zones and NDVI class in Illecillewaet catchment. The threshold index forbare and vegetated is 0.2. Acquisition and processing BEC zone and NDVI dataare described in section 2.2.4.8Table 2.1: BEC zones in Illecillewaet catchment.ZoneDominantTreeElevation(m)Precipitation(mm)Temperature(oC) Key FeaturesICH Westernred cedar,Westernhemlock1000 to1800800 to2500-2 to 18 Cool, moist maritimeclimate, higher speciesdiversity than ESSF,critical habitat for manywildlife speciesESSF Engelmannspruce,Subalpinefir1200 to2500500 to2000-5 to 12 Cool, moist subalpineclimate, lower speciesdiversity than ICH,important timberproduction zoneIMA Nodominanttrees, low-growingvegetationAbove1800500 to1000-5 to 5 Cold, dry alpine climate,short growing season,important habitat formountain wildlife andplant speciesThe Illecillewaet River catchment contains more than 83 glaciers, which covered 58.5km2 in Illecillewaet catchment in 2020, or 5% of the total catchment area. These glaciersrange in elevation from 1726 to 2753 m a.s.l., with a mean elevation of 2356 m a.s.l. Theglacier areas range from 0.05 km2 to 7.7 km2, with an average of 0.58 km2. The glaciersin the Illecillewaet catchment experienced significant mass loss and retreat over the pastfew decades. Between 1985 and 2020, the area covered by glaciers in the catchment hasdecreased by 17%, from 70.1 km2 to 58.5 km2. Additionally, the glaciers have thinnedconsiderably in recent years. Between 1985 and 1999, they thinned by 6.0 m w.e., and from2000 to 2018, by 4.5 m w.e. (Schiefer et al., 2007; Menounos et al., 2019). These changeshave resulted in an estimated mass loss of 0.6 gigatonnes of ice from the watershed between1985 and 2020.2.1.2 HydroclimateThe Illecillewaet catchment experiences a transitional climate, combining characteristics ofboth maritime and continental climates. Winters are long and cold, while summers are shortand mild. From 1977 to 2020, monthly average temperatures in YRV ranged from -3.7oC inDecember to 19oC in July and from -9oC to 13.3 oC in RGR (Figure 2.3). During autumnand winter, the catchment is influenced by frontal-cyclonic weather systems originating fromthe northern Pacific, bringing consistent moisture from the ocean. In contrast, summers aredominated by anticyclonic systems, resulting in frequent periods of clear and dry weather(Moore et al., 2020).9Snowfall accounts for approximately 60% of the total annual precipitation at RGRstation (1330 m), averaging 902 mm annually from 1977 to 2020, primarily occurring duringthe snow accumulation period lasting from October to May. In the last few decades,particularly since the 1980s, higher summer temperatures and reduced winter precipitationhave contributed to ongoing glacial wastage in the upper Columbia River basin (Stahl etal., 2006; Bolch et al., 2010; Jost et al., 2012; Schnorbus et al., 2014; Clarke et al., 2015).Figure 2.3: Climograph for Revelstoke Airport (YRV) and Glacier NP Rogers Pass (RGR)weather stations from 1977 to 2020. The red lines are the monthly averagetemperature. Blue bars are the monthly total precipitation.The Illecillewaet River exhibits a proglacial streamflow regime (Eaton and Moore, 2010).Seasonal snowmelt contributes to high flows from April through June, while glacier meltsustains flows from July through September. Streamflow from October through March isdominated by a prolonged recession, with occasional increases associated with mid-winterrain-on-snow events (Figure 2.4).2.2 Data Sources2.2.1 Meteorological dataMeteorological data were compiled from two weather stations operated by Environment andClimate Change Canada: Revelstoke Airport (YRV) and Glacier NP Rogers Pass (RGR)(Figure 2.1; Table 2.2). Both stations have been operating since 1960. The YRV stationis located near the western extent of the catchment at an elevation of 445 m while RGRis located just east of the catchment boundary at Rogers Pass within Glacier National10Figure 2.4: Illecillewaet River daily streamflow from 1977 to 2020. Line colours vary bywater year.Park at an elevation of 1,330 m. Daily minimum and maximum temperatures and dailyprecipitation were acquired for both stations for the period 1960 to 2020.Table 2.2: Characteristics of weather stations used in this study.Name Latitude Longitude ElevationGlacier NP Rogers Pass (RGR) 51.30 -117.51 1330 mRevelstoke A (YRV) 51.00 -118.18 445 m2.2.2 StreamflowLong-term streamflow data were sourced from the Illecillewaet River at Greeley gaugingstation operated by Water Survey of Canada (ID: 08ND013). The station has been opera-tional since 1963. Data were accessed using the hy_daily_flows() function in the tidyhydatpackage in R (Albers, 2017).2.2.3 Snowpack water equivalentSnow water equivalent (SWE) data were obtained from three manual snow survey sites(Table 2.3): Fidelity Mountain (2A17), Glacier (2A02), and Mount Abbot (2A14) (Figure112.1). The SWE data were downloaded from the BC Government’s Water Data Hub(https://www2.gov.bc.ca/gov/content/environment/air-land-water/water/water-science-data/water-data-tools/snow-survey-data). These sites are sampled once a month duringthe snow season, generally between late December and June. After 1998, the SWE dataare mainly from December to April. Station 2A17, at an elevation of 1870 m, has beenactive from 1963 to the present, while stations 2A14 and 2A02, at elevations of 2010 mand 1250 m, have been active since 1959 and 1973, respectively. The SWE data used inthis study cover the period from 1977 to 2020.Table 2.3: Snow courses in Illecillewaet catchment.Name ID Latitude LongitudeFidelity Mountain 2A17 51.23 -117.70Mount Abbot 2A14 51.23 -117.50Glacier 2A02 51.23 -117.502.2.4 Spatial dataThe Digital Elevation Model (DEM) used for the Illecillewaet catchment was obtained fromthe Canadian Digital Elevation Model (CDEM) dataset. The data were downloaded witha resolution of 0.75 arc seconds, and were reprojected to UTM11 with a resolution of 25 m.The CDEM can be downloaded from Geogratis, a free and open data portal managed bythe Government of Canada (https://www.geogratis.gc.ca/).Landsat 8 data were obtained from the USGS (https://earthexplorer.usgs.gov/) for twopurposes. Firstly, Landsat data were used to generate NDVI. Landsat data with cloudcover less than 20% during July to August were selected. Specifically, band 4 (red) andband 5 (near-infrared) images acquired on August 20, 2018, were chosen to calculate NDVIfor the catchment. Secondly, Landsat data were used to differentiate between snow andbare ice on glaciers during the summer months. Data from late July to early Septemberwere selected, as this period corresponds to a period when most of the snow above the icein the ablation area has melted, exposing the bare ice. Additionally, only Landsat imageswith cloud cover below 20% were considered. A total of 13 images from the years 2013 to2020 were selected, including the visible and near-infrared bands (bands 2, 3, 4, and 5).The Biogeoclimatic Ecosystem Classification (BEC) system was developed by the BritishColumbia Ministry of Forests, Lands, and Natural Resource Operations in the late 1970sand early 1980s, with the most recent update in 2013. The BEC zone polygons for theIllecillewaet catchment were obtained using the bec() function in the R package bcmaps.However, the original data did not differentiate between vegetated and open areas andincluded both zone and subzone classes. To address this, the polygon data were processedvia three steps: 1) cropping for the Illecillewaet catchment, 2) merging by zone to re-duce the number of subzone polygons, and 3) intersecting with a vegetated/open polygongenerated from NDVI. The result was four distinct classes for the Illecillewaet catchment:12ICH, ESSF, IMA, and BARE (open areas). Data for the BEC system can be found at:https://www.for.gov.bc.ca/hre/becweb/index.html.Glacier outline data for the Illecillewaet catchment from 1984 to 2020 were obtained fromthe work of Bevington and Menounos (2022). They generated annual glacier outlines usingLandsat mosaics from late summer and applied glacier band ratio thresholds for extraction.The data can be downloaded from the following link: https://zenodo.org/record/5802106.The obtained glacier outline data were first cropped for the Illecillewaet catchment. Then,the glacier outlines were visually inspected using corresponding Landsat imagery in QGISto ensure accuracy. In the model, the glacier outlines are updated for 4 time periods (Table2.4). Glacier outlines were selected for 1985, 1998, 2007, and 2016, which were years withrelatively hot summers, better image quality, and less snow cover. It should be noted thatthe selected 1985 outline was one of the earliest available, though glacier shrinkage mayhave occurred between 1977 and 1985.Table 2.4: Selected glacier outlines used for model set-up.Year Represented time period Glacier area (km2)1985 1977-1985 70.131998 1986-2001 64.472007 2002-2013 60.632016 2014-2020 54.582.3 Hydrological modelling2.3.1 Classification of hydrologic response unitsIn this study, the Raven model framework was used to simulate streamflow using a modifiedHBV-EC formulation, in which a watershed is composed of sub-basins, which are furtherdivided into hydrologic response units (HRUs) (Craig et al., 2020). HRUs are defined asareas that exhibit unique hydrological responses to precipitation events, can be contiguousor non-contiguous. For this study, delineation of sub-basins was not applied. Several reasonsinfluenced this decision. Firstly, modelling the Illecillewaet catchment as a single basinproved to yield satisfactory results in terms of model performance. Secondly, the focusof this study was primarily on the modelled streamflow, icemelt, snowmelt, and similarvariables for the catchment as a whole. Consequently, simulating and comparing smallersub-basins were not within the scope of this study. Thirdly, inclusion of sub-basins in thedefinition of HRUs would have increased their number and thus increased the computationaltime for each run. Lastly, the absence of available observed data for sub-basins meant thatsimulated streamflow for sub-basins could not be used to support model calibration. Theonly available Water Survey of Canada (WSC) station in the study area is situated at theoutlet of the Illecillewaet River, representing the entire catchment. In this application, theHRUs were non-contiguous, and could also be referred to as grouped response units (GRUs)13according to Moore et al. (2013). Each HRU is characterized by a unique combination ofelevation band, land class, aspect class, and slope class (Figure 2.5).Figure 2.5: The outline to generate HRUs in this study.Five land classes were employed to generate the HRUs in the model: ESSF, ICH, IMA,Bare Land, and Glacier. The first four land classes were derived from the BiogeoclimaticEcosystem Classification (BEC) and NDVI from Landsat 8 image. The BEC data pro-vided distinctions among different forest and vegetation types, including zone and subzonelevels. In this study, only the zone level was retained (refer to the study area section forfurther details). However, a limitation of the BEC data was its failure to classify bareland and forest-covered land within each zone. To address this, NDVI was calculated,and a vegetation-bare land threshold of 0.2 was applied to classify forest and bare land.The model incorporates a land class for glaciers, for which the HRUs were updated ap-proximately every 10 years. Because visual inspection of the data revealed a retreat ofapproximately one or two pixels (30-60 m) every decade, a more frequent updating wouldnot have significantly impacted the results. The automatic mapping of glacier outlines wasfollowed by visual inspection and comparison with Landsat data to make minor corrections(e.g., distinguishing shadows from rocks) and select more accurate glacier outlines. Ulti-mately, the glacier outlines from 1985, 1998, 2007, and 2016 were chosen for updating themodel (see Table 2.4). The remaining three classification factors, namely elevation bands,aspect classes, and slope classes, were derived from the Canadian Digital Elevation Model(DEM). Elevation bands were defined at intervals of 100 m, while slope classes were catego-rized into three groups: 0-15o, 15-30o, and 30-90o. Aspect classes followed the classificationprovided in Table 2.5. Finally, the interception of the layers of the factors resulted in atotal of 791 HRUs for the model.14Table 2.5: Aspect class.Aspect class Aspect (degrees)1 0-45 and 315-3602 45-1353 135-2254 225-3152.3.2 Model descriptionHydrological processes within the Illecillewaet catchment were simulated using a modifiedversion of the HBV-EC hydrologic model, which is a Canadian adaptation of the originalHBV hydrological model developed by (Bergström, 1995; Hamilton et al., 2000). The HBV-EC model can be emulated within the Raven Hydrological Modelling Framework (Craig etal., 2020).The HBV-EC model provides a suitable approach for studying hydrological processes indata-sparse regions. This model has been extensively utilized for modelling mountainouscatchments in British Columbia and Alberta, as demonstrated by previous studies (Stahl etal., 2008; Jost et al., 2012). The model uses a combination of empirical and process-basedparameterizations to simulate hydrologic fluxes and changes in storage. Descriptions of themodel algorithms can be found in Hamilton et al. (2000) and Stahl et al. (2008), and inthe Raven User’s Manual. The HBV model operates at a daily time step and uses dailyair temperature (mean of daily minimum and maximum) and precipitation data, whichare spatially distributed across the watershed through nearest neighbor interpolation fromclimate stations. To account for orographic effects, the model incorporates monthly varyingtemperature lapse rates and precipitation gradient calculations based on monthly data fromYRV and RGR stations.The phase of precipitation was determined using the HBV linear transition method,which computes the fraction of precipitation falling as snow based on daily average airtemperature. The model then simulates the routing of precipitation through the forestcanopy, accounting for both snow and rainfall interception. Snow interception follows theapproach described in Hedstrom and Pomeroy (1998), while rainfall interception is deter-mined as a percentage using an exponential LAI-based method, where the interception isrelated to the leaf area index and stem area index (LAI and SAI).Snowmelt occurs when the air temperature exceeds the melting temperature threshold,typically around 0 oC. The HBV-EC algorithm employs a temperature index model with aseasonally varying and spatially adjusted melt factor. This method has demonstrated reli-ability in various regions and environments (Jost et al., 2012), particularly in mountainousareas with limited data availability for energy balance methods, such as solar radiation,wind speed, and humidity observations. Daily snowmelt for a given HRU is computed asfollows:15𝑀(𝑡) = 𝑀𝐹(𝑡) ⋅ [𝑇𝑎(𝑡) − 𝑇𝑚] (2.1)where 𝑀(𝑡) is daily snowmelt (mm day-1) on day-of-year 𝑡, 𝑀𝐹(𝑡) is the melt factor onday-of-year 𝑡 (mm day-1 oC -1), 𝑇𝑎(𝑡) is mean daily air temperature on day-of-year 𝑡, and𝑇𝑚 is a threshold temperature (oC), below which melt ceases.For non-glacier HRUs classed as open sites, the melt factor is computed as𝑀𝐹(𝑡) = 𝐶0(𝑡) ⋅ [1 − 𝐴𝑀 ⋅ 𝑠𝑖𝑛(𝑠) ⋅ 𝑐𝑜𝑠(𝑎)] (2.2)where 𝐶0(𝑡) is a sinusoidal function of day-of-year, with a minimum value 𝐶𝑚𝑖𝑛 on thewinter solstice and a maximum value of 𝐶𝑚𝑖𝑛 + Δ𝐶 on the summer solstice, 𝐴𝑀 is aparameter to be calibrated, 𝑠 and 𝑎 are the mean slope and mean aspect of the cells in theHRU (degrees), respectively.For non-glacier HRUs classed as forest sites, the melt factor is multiplied by a melt-reduction factor 𝑀𝑅𝐹 , which can be calibrated and takes values less than unity. Forglacier HRUs, the open-site melt factor is used when the snow water equivalent in the HRUis positive; when the HRU becomes snow-free, the open-site melt factor is multiplied bythe parameter 𝑀𝑅𝐺, which typically ranges from 1.25 to 1.75, and accounts for the loweralbedo of glacier ice compared to snow.Rainfall and snowmelt infiltrate into the soil using the standard approach of the HBVmodel (Bergström, 1995). The infiltration rate is determined by a power-law functionrelated to the soil water content proportion and the soil’s maximum storage capacity. Oncewater enters the soil layers, it moves upwards through capillary rise and downwards throughpercolation. Capillary rise is linearly proportional to the soil saturation of the upper layersoil, while percolation is driven by the soil moisture gradient. To simplify the model, aconstant percolation rate between adjacent soil layers is used.Soil water returns to the atmosphere and the surface through evapotranspiration andbaseflow, respectively. Soil evapotranspiration equals PET when the soil storage capacityexceeds the tension storage threshold. Otherwise, it is linearly proportional to the soilsaturation. In non-forested areas, evapotranspiration is zero if the snow depth is non-zero.Soil water also contributes to surface water as discharge from the fast and slow reservoirs inthe model. The fast reservoir calculates the discharge rate using a two-parameter power-lawresponse, while the slow reservoir layer employs a discharge rate that is linearly proportionalto the water storage in the reservoir.When the landscape is converted from “Glacier” to “Bare,” the model assumes that allglacial ice has melted, and no additional glacier ice melt or glacier discharge is consideredfor that HRU. The hydrologic balance is applied to HRUs, and the calculated soil water,glacier melt, rainfall runoff are routed to the outlet using a triangular unit hydrographapproach. These flows are then instantly routed to the channel outlet according to themethodology described by Craig et al. (2020).162.3.3 Calibration and validation procedureTwo calibration processes were applied in this study. One calibration was performed from1977 to 2009, followed by a validation period from 2010 to 2020. This calibration andvalidation aimed to demonstrate the model’s capability to accurately simulate streamflowand glacier mass change for an independent time period that is different from the cali-bration period in terms of both climatic conditions and glacier extents to help evaluatethe robustness of the model set-up and calibration. The other calibration focused on theentire time period from 1977 to 2020, without including a separate validation period. Thepurpose of this calibration was to generate parameter values that accurately simulated his-torical streamflow and glacier mass change for comprehensive analysis to provide a basis fordiagnosing the relative influences of climate variability and glacier changes on streamflowtrends over the study period, rather than projecting future trends. For both calibrationapproaches, model parameters were calibrated to streamflow observations from the WSCgauging station Illecillewaet River at Greeley. The software tool Ostrich (Matott, 2017)was employed, utilizing the parallel dynamically dimensioned search (PDDS) algorithmproposed by Tolson and Shoemaker (2007).Selection of sensitive parameters for calibration was based initially on past experiencewith HBV-EC (Hamilton et al., 2000; Stahl et al., 2008; Jost et al., 2012) as well as someexploratory runs with manually adjusted parameter values. A total of 13 parameters wereidentified for calibration. Initial parameter sets were established as guided “first estimates”that adhered to physical realism and literature values. These initial sets were manuallyadjusted to approximate the shape and characteristics of the annual hydrograph.To accurately simulate glacier melt, the aspect correction factor, to which glacier meltis highly sensitive, was manually calibrated using observed glacier mass change data. Forthe second calibration process, this aspect correction factor was initially calibrated usingglacier mass change data from 2000 to 2009, followed by validation using data from 2010to 2019.During both calibration and validation phases, the model’s performance was assessedusing the Nash-Sutcliffe Efficiency (NSE) criterion. Additionally, the simulated SnowWaterEquivalent (SWE) was compared to measurements from snow survey sites to further verifythe model’s accuracy.2.3.4 Extraction of contributions associated with rainfall, snowmelt and icemeltThe Raven modeling framework allows the calculation of streamflow contributions originat-ing from distinct regions within the watershed or various water storage compartments acrossthe landscape (Craig et al., 2020). In this study, the contributions of rainfall, snowmelt,and glacier ice melt to the streamflow were computed on a daily time step utilizing Raven’stracer algorithm. Specifically, the analysis focused on the period from July to September,given that these are the months in which glacier contributions are most important.172.3.5 Extraction of glacier mass changeTo assess the model’s performance in simulating glacier ice melt, a comparison was con-ducted between the observed and modelled glacier mass change for two distinct time pe-riods: 2000 to 2009 and 2010 to 2019. The observed data utilized in this analysis werederived from ASTER height measurements (Hugonnet et al., 2022). The height changeson glaciers within the Illecillewaet catchment were computed for the two time periods, andsubsequently converted into water volume loss by multiplying the glacier areas with an icedensity of 0.85 kg/m3 (Huss, 2013). For the calculation of the simulated mass change, ayearly glacier mass balance was first computed following the approach outlined by Stahl(Stahl et al., 2008). The yearly minimum ice and snow storages (expressed as water equiv-alent) for each glacier HRU were obtained from Raven outputs. By subtracting the currentyear’s minimum ice and snow values from the previous year’s, the annual glacier masschange for each HRU was determined. This process was repeated for each year withinthe two 10-year time periods, and the resulting glacier mass changes were converted intoglacier volume loss by multiplying them by each HRU’s area. Finally, the volumes of theHRUs were summed to obtain the total glacier volume loss. Following the calculation ofthe observed and simulated glacier loss volumes, the results were normalized by dividingthem by the catchment area for easier readability and comparison of the data.2.4 AnalysisIn order to gain a better understanding of the trends in summer streamflow and ice meltcontributions, as well as the role of glacier retreat and climatic variability in shaping thesetrends, the model was executed for four distinct scenarios. Trend analyses were then con-ducted for summer streamflow and ice melt based on the outcomes of these scenarios. Thefour scenarios are as follows:a. Retreating glaciers. This scenario aimed to reflect the most realistic conditions byupdating the glacier outlines approximately every 10 years, as described in Section2.3.2.b. Constant glacier coverage based on the earliest extents. In this scenario, the glacieroutline from 1985, a total glacier cover of 70 km2, was employed throughout the entiretime period. This scenario sought to examine the potential outcomes if the glaciersdid not retreat.c. Constant glacier coverage based on the latest extents. Here, the glacier outline from2016, representing a total glacier cover of 55 km2, was used to portray a scenariowith consistent glacier coverage. This scenario aimed to explore the implications ofa smaller constant glacier coverage compared to the previous scenario. The choiceof 2016 as the representative year was because the glacier outline data for that yearwas more accurate than the data available for 2020, as confirmed through manualverification with Landsat 8 data.18d. No glaciers. In this scenario, all glacier covers within the model were replaced withalpine bare land cover. The objective of this scenario was to assess the influence ofclimatic variability on the trend of streamflow.For each scenario, summer streamflow (July, August, and September) and ice melt wereextracted, and trend analyses were performed by calculating the Spearman correlation coef-ficient between summer streamflow and the corresponding year. The Spearman correlationcoefficient enabled the assessment of trends without assuming specific data distributions orlinearity. The magnitude of the trends was estimated using the Theil–Sen slope estimator,which is robust to outliers and to the lack of normality in the residuals (Helsel and Hirsch,2002).193 Results3.1 Split sample model calibration and validation3.1.1 Streamflow simulationFor the split-sample calibration, the Nash-Sutcliffe efficiency (NSE) for the 1977 to 2009calibration period is 0.893 and even higher at 0.901 for the 2010 to 2020 validation period.As seen in Figure 3.1, four randomly selected years’ hydrographs in the validation period,the model generated reasonably accurate simulations of streamflow associated with summersnow and glacier melt and fall rain events. However, the model underestimated dischargefor some spring snowmelt events and winter base flow. The calibrated parameters are shownin Table 3.1.Table 3.1: Calibrated parameters for full period calibration and split sample calibration.Parameter Parameter descriptionSplitsamplecalibra-tionFullperiodcalibra-tionSnowCorrection in RGR Snow Correction in RGR 1.00 1.00SnowCorrection in YRV Snow Correction in YRV 1.30 1.30AdiabaticLapseRate Adiabatic Lapse Rate 6.07 5.55RAINSNOW_TEMP Rain/Snow halfway transitiontemperature1.00 1.00RAINSNOW_DELTA Range of rain-snow transitionzone0.71 0.44MAX_PERC_RATE Percolation rate 4.96 6.37BASEFLOW_COEFF Linear baseflowstorage/routing coefficient0.02 0.02BASEFLOW_N Baseflow exponent 1.10 1.11MELT_FACTOR MAX-MIN Minimum snow melt factorused in degree day models2.50 2.49MIN_MELT_FACTOR The gap between minimum andmaximum snow melt factors1.20 1.24HBV_MELT_GLACIER_CORR Degree day correction factorfor glacial melt1.75 1.75POROSITY Effective porosity of the soil 0.36 0.4020Parameter Parameter descriptionSplitsamplecalibra-tionFullperiodcalibra-tionHBV_MELT_ASP_CORR HBV snowmelt aspectcorrection0.51 0.533.1.2 Snowpack water equivalentFigure 3.2 presents a comparison between measured snow water equivalent (SWE) fromthree snow courses and the simulated values from the corresponding HRUs during both thecalibration and validation periods. Figure 3.3 shows the locations of the snow courses andthe centroids of the HRUs that represent the snow course locations. Due to limited observeddata from April to June since 1998, fewer data points are available for the validation period.The model generally performed well in simulating SWE, indicating that the calibratedparameters effectively captured SWE dynamics in a separate validation period. Most datapoints align closely to the 1:1 line, particularly for snow accumulation during the wintermonths from December to March in stations A14 and A17. However, the model tendedto overestimate SWE in April, May, and June for all three snow sites, but especially forstation A17.3.1.3 Glacier mass changeAfter manual calibration, specifically adjusting the aspect melt correction parameter 𝐴𝑀to 0.513 to match the glacier mass change between 2000 and 2009 (-0.169 m), the simulatedglacier mass change for the 2010 to 2019 validation period was -0.206 m, which is about7% lower in absolute value than the observed value of -0.219 m (Table 3.2).Table 3.2: Observed and modelled Illecillewaet catchment glacier mass change for two timeperiods. The values are expressed as a water equivalent depth averaged over thecatchment area.Time period Observed (m) Modeled (m)2000-2009 -0.169 -0.2192010-2019 -0.169 -0.20621Figure 3.1: Observed and simulated hydrographs in randomly selected years in the valida-tion period: 2011, 2013, 2016, 2018.22Figure 3.2: Scatterplots of simulated versus observed SWE in the calibration and validationperiods.23Figure 3.3: Snow course locations and the centroids for the corresponding HRUs.243.2 Full period model calibration and performance3.2.1 Streamflow simulationWhen calibrated for the entire study period, the Nash-Sutcliffe index was 0.89. As seen inFigure 3.4, the model generally performed well at capturing streamflow variations in themelt season and during autumn rain events. A scatter plot of predicted versus observedflows indicates that the model provided relatively unbiased streamflow predictions for flowsin excess of about 100 m3s-1 (Figure 3.5). However, the model tended to underestimatelower flows, especially during the winter. The calibrated parameters are shown in Table3.1. As can be seen, the parameters did not show much difference between the split-sampleand full-period calibrations.Figure 3.4: Observed and simulated hydrographs in randomly selected years: 1984, 1994,2006 and 2016.25Figure 3.5: Monthly model result.263.2.2 Snowpack water equivalentAs seen in Figure 3.6, simulated SWE generally agreed well with observed SWE for thehighest-elevation station (A14), although it tended to underestimate SWE from April toJune. For the two lower-elevation stations, the model tended to over-predict SWE in theaccumulation season (January to March) and under-predict in April to June, even more sothan for the split-sample calibration and validation.Figure 3.6: Scatterplots of simulated versus observed SWE.3.2.3 Glacier mass changeAfter manually adjusting the melt aspect correction factor on the glacier mass change data,the model successfully simulated glacier loss for both time periods (2000-2009 and 2010-2019, in Table 3.3). Moreover, it accurately captured the clear increase in glacier mass lossin the latter ten years. Compared to the split sample calibration approach (Section 3.1.3),the full-period calibration slightly improved the model’s performance in simulating glaciermass change during the second time period.27Table 3.3: Observed and modelled glacier mass change for two time periods in Illecillewaetcatchment, expressed as an average depth of water over the catchment.Time period Observed (m) Modeled (m)2000-2009 -0.169 -0.2192010-2019 -0.168 -0.2103.3 Trends in observed climate and streamflow data3.3.1 Climate trendBoth weather stations show increasing temperature trends and declining precipitationtrends in July, August, and September (Figure 3.7 and Table 3.4). The temperature trendsare more significant, with all p-values being less than 0.1, and 4 out of 6 being less than0.05. For precipitation, only the August and July trends at the YRV station are significant(p-values < 0.05 and < 0.1, respectively), but they exhibit higher magnitudes, with a de-crease of approximately 6 to 7 mm per decade. The temperature trend analysis indicatesa warming trend of 0.3 to 0.4 oC per decade during summer, with a slightly higher trendin July.Table 3.4: Theil-Sen slope for temperature and precipitation time series. The bold anditalics values indicate that the Spearman correlation coefficient is significantlydifferent from 0 (p value less than 0.05 and 0.1 respectively).MonthTemperature(RGR)Temperature(YRV)Precipitation(RGR)Precipitation(YRV)July 0.041 0.048 -0.463 -0.644August 0.037 0.032 -0.708 -0.687September 0.034 0.040 0.921 -0.0633.3.2 Streamflow trendFrom 1977 to 2020, there was a notable and significant decline in August streamflow (Figure3.8 and Table 3.5), with a decrease of 0.44 m3s-1 per year. However, no clear trends areobserved for July and September.Table 3.5: Trend analyses for summer streamflow for Illecillewaet River. (bold is significantin Spearman correlation test, with p value less than 0.05).Month Theil-sen slope (m3s-1 per year) Spearman correlation p valueJuly -0.049 0.8328Month Theil-sen slope (m3s-1 per year) Spearman correlation p valueAugust -0.436 0.01September 0.017 0.943.4 Ice melt contributions to streamflowFigure 3.9 presents the monthly modelled ice melt contributions to streamflow. The resultsindicate that glacier ice predominantly begins melting in July, with the melting processpeaking in August and persisting through September. During August, glacier ice meltcontributed up to 21.3 m3s-1 in absolute terms or up to 41% of the total streamflow.Figure 3.10 shows a notable shift in the relative contributions of snowmelt, ice melt, andrainfall to streamflow from July to September. In July, snowmelt remained the primarysource throughout the study period, with relatively minor contributions from ice melt andrainfall. Snowmelt contributions declined in August, while ice melt contributions becamemore substantial. By September, the decline in snowmelt continued, accompanied by aslight reduction in ice melt compared to August. However, the rainfall increased signifi-cantly, constituting nearly half of the total streamflow. Thus, the transition from July toSeptember demonstrates a shift in the relative contributions of snowmelt, ice melt, andrainfall to streamflow dynamics.Table 3.6: August ice melt contribution to streamflow, July to September average dailytemperature, average monthly precipitation, and December to April averagemonthly precipitation for 1990, 2003 and 2011.YearAugust Ice meltcontribution (%)Augusttemperature(oC)Augustprecipitation(mm)WinterPrecipitation(mm)August flow(m3s-1)1990 21.1 14.1 66.5 968.8 66.42017 41.1 14.7 27.3 658.3 51.61999 3.5 13.3 119.4 1117.3 146.4Glacier ice melt contributions in August varied significantly depending on weather condi-tions (Table 3.6). In 1999, with high August precipitation and lower temperatures, coupledwith substantial snow accumulation during the preceding winter (November to April), icemelt only accounted for 3.5% of total streamflow. In contrast, the warmer and drier summerof 2017, accompanied by lower winter precipitation, resulted in glacier melt contributing41% of total streamflow. The year 1990, with average temperature and precipitation, fellin between these extremes.29Figure 3.7: Temperature and precipitation time series for both RGR and YRV stations.30Figure 3.8: Streamflow time series for July, August and September.31Figure 3.9: Boxplots of ice melt contributions to streamflow by month for May, June, July,August, September and October. The upper panel shows absolute contributionsand the lower panel shows percentage contributions.32Figure 3.10: Streamflow contributions from glacier ice melt, rainfall and snow melt for July,August and September.333.5 Trend analyses of monthly streamflow for four scenariosThe three summer months exhibited varying degrees of decline in streamflow (Figure 3.11and Table 3.7, 3.8 and 3.9), with non-significant trends observed in July and September.August streamflow exhibited significant trends for all four scenarios, with three out of fourSpearman correlation p-values being less than 0.05, while the fourth is less than 0.1.Specifically focusing on August, the no-glacier scenario demonstrates the most substan-tial decline in streamflow, indicated by a Theil-Sen slope of -0.51 m3s-1 per year. For theretreating glacier scenario, the slope of the decreasing trend is -0.48 m3s-1 per year, which islower in magnitude than the no-glacier scenario. The two constant glacier-extent scenarios,showed the least declining trends, with -0.27 and -0.31 m3s-1 per year respectively for theearliest and latest extents.Figure 3.11: Streamflow time series for four scenarios in July, August and September.3.6 Trend analyses of monthly ice melt contributions for threeglacier-cover scenariosIce-melt contributions to streamflow exhibited increasing trends during the summer monthsof July, August, and September in association with rising temperatures over the years34Table 3.7: Trend analyses for July streamflow for Illecillewaet River.Streamflow trend Icemelt trendScenario Theil-Sen slope Spearman p value Theil-Sen slope Spearman p valueRetreating glacier -0.22 0.52 0.006 0.64Earliest extents -0.20 0.66 0.058 0.09Latest extents -0.26 0.60 0.035 0.10No glacier -0.37 0.52 NA NATable 3.8: Trend analyses for August streamflow for Illecillewaet River. The bold valuesindicate that the Spearman correlation coefficient is significantly different from0 (p value less than 0.05).Streamflow trend Icemelt trendScenario Theil-Sen slope Spearman p value Theil-Sen slope Spearman p valueRetreating glacier -0.48 0.003 0.08 0.139Earliest extents -0.27 0.093 0.21 0.020Latest extents -0.31 0.044 0.15 0.014No glacier -0.51 0.007 NA NATable 3.9: Trend analyses for September streamflow for Illecillewaet River.Streamflow trend Icemelt trendScenario Theil-Sen slope Spearman p value Theil-Sen slope Spearman p valueRetreating glacier -0.070 0.37 0.02 0.41Earliest extents 0.024 0.85 0.09 0.06Latest extents 0.004 0.95 0.06 0.06No glacier -0.040 0.68 NA NA35(Figure 3.12 and Table 3.7, 3.8 and 3.9). The statistical analysis indicates that the p-valuesfor both the earliest and latest constant glacier outlines are consistently below 0.1. InAugust, these p-values fell below 0.05. However, the ice melt contributions to streamflowdid not exhibit a significant trend in the retreat scenario according to Spearman correlationanalysis. As seen in Figure 3.13, although there was a clear decline in glacier snow cover inAugust from 1977 to 2020, associated with higher snow lines, the exposed ice area (the gapbetween the two lines) has not significantly changed as a consequence of terminal retreat.Figure 3.12: Time series of ice-melt contributions to streamflow for three scenarios in July,August and September.36Figure 3.13: Time series of modelled area of glacier snow cover (average of August from1977 to 2020) and observed glacier areas in 1985, 1998, 2007 and 2016.374 Discussion4.1 Utility of using Raven for trend attribution studiesTwo notable advantages of using Raven in this study are its abilities to update land coverclasses and to trace the streamflow sources. The land class updating feature allows formore realistic simulations of streamflow by accounting for retreating glacier extents overtime. This approach is computationally more straightforward compared to those usedby Jost et al. (2012) and Moyer et al. (2016). Jost et al. (2012) implemented codesthat halted the simulation, read in new HRUs representing updated glacier extents, andthen resumed the simulation. Moyer et al. (2016) ran their model for the entire 31-yearsimulation period 31 times, each time using the glacier extent for a given year throughoutthe entire simulation period, then selected from each run the year associated with theglacier outline, and then combined these into a time series.Duethmann et al. (2015) applied the WASA conceptual hydrological model, which al-lowed the updating of glacier extents within the model on an annual basis using an assumedlinear decrease in glacier area based on two glacier inventories. In contrast to this approach,the advantage of using observed glacier outlines within the Raven framework is that it allowsthe simulation to accommodate nonlinear trends in glacier area (as shown in Figure 3.13),which could have an influence on trends in ice-melt contributions to streamflow. Raven’stracing function enables direct estimation of glacier melt from the model and facilitates theanalysis of its temporal trend. This approach is more direct than the approach used byJost et al. (2012) and Tiel et al. (2020), which involves comparing modelled streamflowfor scenarios with and without glacier covers for the entire simulation period.4.2 Validity of the model resultsThe hydrologic model used in this study, a modified HBV-EC implemented in Raven,demonstrated good performance in simulating streamflow trends. In split-sample calibra-tion, both streamflow and glacier mass loss were reproduced well in the validation period,which supports the validity of the model results for diagnosing the causes of the streamflowtrends. The Nash-Sutcliffe Efficiency (NSE) in the calibration period was 0.89 and evenhigher in the validation period of 0.90. The simulated glacier mass change, expressed as anaverage depth of water over the catchment, was calibrated from 2000 to 2009 at the samevalue as the observed (-0.169 m). In the validation period, the simulated mass change was-0.206 m, which aligns with the observed data and the result from full-period calibration.38However, there was a discrepancy between simulated and observed SWE. The modelledvalues are from the HRUs that correspond to the snow courses, which also include cellsdistributed over the catchment. Consequently, as seen in Figure 3.3, there is a substantialdistance between the actual snow course locations and the centroids of their correspondingHRUs, which means that there would be differences in the temperature and precipitationinputs used to simulate hydrologic fluxes.In the full-period calibration, the NSE is 0.89, and the calibrated parameters did notchange much from the split-sample calibration. The focus of this study was on summerstreamflow, particularly in August, which exhibited a statistically significant decreasingtrend with a slope of -0.44 m3s-1 per year based on observed data. The model successfullycaptured this decreasing trend, with the simulated August streamflow for the retreatingglacier scenario showing a similar slope of -0.48 m3s-1 per year. The simulated glaciermass change was -0.169 m and -0.210 m for the periods 2000 to 2009 and 2010 to 2019,respectively, in full-period calibration. That is, the rate of glacier mass loss increased due tothe warmer and drier conditions in the more recent 10-year period. These findings align withthe results presented by Menounos et al. (2019), which suggest a higher per-unit-area massloss in the Columbia headwaters region during 2009-2018 compared to 2000-2009. However,it was difficult to compare the SWE data between the modelled and the simulated, becauseof the differences in locations of the snow courses and the HRU centroids representing thoselocations (Figure 3.3).Using the temperature-index melt model may introduce potential uncertainty becausethe changing climate conditions might affect melt factors over the years. However, thisapproach has proven robust in hydrological and glaciological simulations especially whencalibrated for daily to seasonal time scales (Hock, 2003). In this study, we used glaciermass loss and streamflow data to constrain the calibration of melt factors. Additionally, ourmodel only focused on historical conditions, not taking any extrapolation into the future,thus ensuring the reliability and validity of the simulated snow and ice melt. Further,the split-sample calibration demonstrated that the calibrated parameters do appear to betransferable within the historic period.4.3 Glacier melt contributions to streamflow in IllecillewaetcatchmentIn the Illecillewaet catchment, where glacier cover represents approximately 5% of the totalarea, glacier ice melt contributed up to 41% of the total streamflow during warm and dryAugust periods. This finding aligns with the fractional contributions reported by Jost et al.(2012) and Chernos et al. (2020). Jost et al. (2012) concluded that ice melt contributed upto 25% of August streamflow in the Mica Basin of the Columbia River headwaters, whichhas a similar glacier cover of about 5%. Similarly, Chernos et al. (2020) found that glacierice melt contributed up to 20-30% of the August streamflow in the Athabasca River Basinwith 2.8% of glacier cover. Conversely, in a wet year in 1999, the contribution from glaciermelt drops to as low as 3.6% of the August streamflow. The findings confirm the general39understanding that warmer and drier summer weather amplifies the importance of glaciermelt in buffering streamflow.4.4 Streamflow trends and attributionSimilar to the study by Duethmann et al. (2015), a major goal of this study was to diagnosethe relative contributions of climatic trends and glacier retreat to observed streamflowtrends. Duethmann et al. (2015) used a scenario driven by detrended precipitation and airtemperature data generated by using the daily pattern for a single year, with precipitationand air temperature scaled so that monthly totals and means conformed in magnitude tolonger-term climatological means. The approach used in this study involved using a no-glacier scenario to represent climate-driven variability, and comparing the trend magnitudefor that scenario to those for the retreating-glacier and constant-glacier scenarios. Theadvantage of the approach used in this study is that it incorporates the effects of interannualand decadal climatic variability in a more realistic manner.From the observed streamflow trends (Figure 3.8 and Table3.5), August streamflowdecreased significantly throughout the years, while no notable trends were observed forJuly and September. Furthermore, the observed climate trends (Figure 3.7 and Table3.4) demonstrate that air temperature exhibited increasing trends for July to September,while only August experienced a significant (negative) trend in precipitation. Consideringall these observations, it can be inferred that the declining streamflow in August may beattributed to reduced precipitation and/or changes in glacier ice melt patterns specific tothat month.According to the analysis of modelled August streamflow trends with the four scenar-ios (Figure 3.11 and Table 3.7, 3.8 and 3.9), the no-glacier scenario exhibited the mostsignificant decreasing trend, with a rate of -0.51 m3s-1 per year. Given that there was noland-cover change over the simulation period, this decline in August streamflow can beattributed to climate variations, specifically the decreasing trend in August precipitation,possibly augmented by an increase in evapotranspiration associated with the air tempera-ture trend. In contrast, the three glacier-covered scenarios show a less pronounced decline,and the magnitude of the decrease diminishes with increasing glacier cover. This resultindicates that the declining streamflow cannot be attributed to a decline in glacier-meltcontributions, and thus can be primarily attributed to less rainfall in August. The increas-ing trend in August ice-melt contributions occurred even despite a substantial decrease inglacier area. As seen in Figure 3.13, the modelled glacier snow cover in August also had adecreasing trend, such that the area of bare ice remained roughly constant. This lack of adecrease in ice area, in combination with increasing air temperature, apparently resultedin an increase in total ice melt.Comparing the streamflow trends among the glacier-covered scenarios, the retreatingscenario exhibits the largest Sen slope of -0.48 m3s-1 per year, slightly less negative than thetrend for the no-glacier scenario, followed by the constant-latest-extent scenario, and theconstant-earliest-extent scenario showed the least decreasing trend with a Sen slope of -0.2740m3s-1 per year. These findings provide three key insights. Firstly, even when the glacierswere not retreating during the study period, the impact of climatic variations still resultedin a decline in August streamflow, albeit to a lesser extent. Secondly, compared to theconstant-earliest-extent scenario, the retreating scenario yielded less glacier melt and had alower ability to offset the decrease in rainfall, leading to less streamflow. This finding alignswith the concept of “realized streamflow” discussed by Moyer et al. (2016); that is, eventhough ice-melt contributions to streamflow increased through time, the effect of glacierretreat was to produce less streamflow than would have been the case had glacier arearemained constant. Thirdly, under similar climate conditions, a catchment with a higherpercentage of glacier cover may have a more stable streamflow trend due to its greatercapacity to counterbalance climatic variations.The trends in August streamflow and ice-melt contributions in the Illecillewaet Rivercatchment are consistent with the simulated trends observed for the Hoh River catchmentin Washington State, as depicted in Figure 9 by Frans et al. (2018). Their analysis revealeddeclining August streamflow from about 1980 to 2000, a period that was associated withan increase in glacier melt contributions to streamflow despite a reduction in glacier area.These findings differ from the conclusions drawn by Moore et al. (2020), who applied aspatial comparison approach that found that catchments with more glacier cover exhibiteda greater declining trend in August streamflow, as shown in Figure 9 of that paper. Theauthors inferred that, except for the most glacierized catchments, the Columbia Riverheadwaters region is already experiencing a post-peak-water decline.A potential explanation for this discrepancy between the current study and that ofMoore et al. (2020) is that the spatial pattern of glacier cover may be correlated with someother variable that influences streamflow trends. For example, the unglacierized catchmentstend to have lower elevations and be located in the south, with glacier cover and elevationsof the mountain peaks generally increasing toward the north (see Figure 1 in that paper).It may be that precipitation trends vary with latitude and/or elevation. Indeed, Luce et al.(2013) reported a declining trend in orographic enhancement in the US Pacific Northwest,immediately south of the Columbia River headwaters. Future research should focus on thespatial distribution of precipitation and glacier melt trends throughout the Columbia Riverbasin to gain a better understanding of these patterns and their implications for streamflowdynamics.4.5 Implications for the peak water conceptThe peak water concept describes the effect of a shift from neutral to chronically negativeglacier mass balance on glacier discharge (as introduced in 1.2.1). A number of studieshave attempted to test this hypothesis using empirical analysis of streamflow data forpartially glacierized catchments (e.g., Stahl et al., 2006; Casassa et al., 2009; Baraer et al.,2012). The results of this analysis, and that of Frans et al. (2018), indicate that streamflowtrends in catchments with relatively low amounts of glacier cover (5% for Illecillewaet Riverand less than 4% for the Hoh River) may be dominated by climate-driven trends in the41unglacierized portions of the catchment. Therefore, future studies should not rely purely ondata-driven analyses, but should include model-based analysis to support the attributionof the trends.425 Conclusions5.1 Key findingsThe main findings from this study are summarized below:a. The analysis of observed meteorological data revealed a significant increase in temper-ature from July to September and a decline in precipitation during July and August.Additionally, the observed August streamflow exhibited a notable reduction from 1977to 2020.b. The hydrological model indicates that glacier ice typically starts melting in July andcontinues until September, with the highest contribution to streamflow occurring inAugust. In the Illecillewaet catchment, where glacier cover represents about 5% ofthe total area, glacier ice melt can contribute between 3.5% and 41.1% to the Auguststreamflow, depending on weather conditions.c. The analysis of glacier ice melt trends in the Illecillewaet catchment showed a smallbut statistically insignificant increase. This is attributed to the synchronous rising ofthe snowline along with the retreat of the glacier terminus.d. By examining August streamflow trends for four simulation scenarios, the declinein August streamflow appears to have been primarily driven by climate influences,particularly decreasing precipitation in August. At the current stage, the presenceof glaciers helps offset this declining trend, but the glacier compensation effect isweakening as the glaciers retreat. However, the combination of retreating glaciers anddeclining precipitation will likely amplify the decline in summer streamflow, leadingto more dramatic decreases in the future in Illecillewaet River catchment.e. The results of this study further support findings from previous studies that the clas-sic glacial peak water cycle is overly simplistic when applied to partially glacierizedcatchments, and does not adequately address the influences of climate-driven stream-flow trends in the unglacierized areas.5.2 Avenues for future researchIn this study, the definition of hydrologic response units did not include a subdivisionby sub-basins. This choice made the model setup easier, allowing for shorter simulationtimes, and also performed well in terms of simulating both streamflow and glacier ice melt.43However, not having sub-basins made it challenging to compare the modelled SWE with theobserved SWE data from snow courses because the spatial centroids of the grouped responseunits that corresponded with the snow courses were located some distance west of the snowcourse locations, such that precipitation interpolated from the weather stations may nothave been representative of that received at the snow courses. Further work should focuson experimenting with the inclusion of sub-basins in the delineation of grouped responseunits to assess whether it can improve the utility of SWE data at snow courses for modelcalibration and testing. Furthermore, adding subbasins near the boundaries of individualglaciers in the catchment would allow a comparison of modelled glacier mass balance withthe observed data for a specific glacier, such as the Illecillewaet glacier.To better understand the discrepancy with the results of Moore et al. (2020) anddetermine the “peak water” phase in the Columbia River Headwaters, future work shouldinvolve applying this model-based method to a wider range of catchments with varyingglacier cover in the Columbia River headwater basins. By doing so, the spatial distributionof summer precipitation trends, late summer streamflow trends, and glacier ice melt trendscan be assessed in this region. This broader analysis would help understand whether thedeclining streamflow can be attributed differently to lower and higher latitudes or if thefinding from the Illecillewaet catchment is an exception in this region.44ReferencesAlbers S. 2017. Tidyhydat: Extract and tidy canadian hydrometric data. The Journal ofOpen Source Software 2 (20): 511 DOI: 10.21105/joss.00511Baraer M, Mark BG, McKenzie JM, Condom T, Bury J, Huh K-I, Portocarrero C, GómezJ, Rathay S. 2012. Glacier recession and water resources in Peru’s Cordillera Blanca.Journal of Glaciology 58 (207): 134–150 DOI: 10.3189/2012jog11j186Baraer M, McKenzie JM, Mark BG, Bury J, Knox S. 2009. Characterizing contributionsof glacier melt and groundwater during the dry season in a poorly gauged catchment ofthe Cordillera Blanca (Peru). Advances in Geosciences 22: 41–49 DOI: 10.5194/adgeo-22-41-2009Bergström S. 1995. The HBV model. Computer models of watershed hydrology.: 443–476Available at: https://www.cabdirect.org/cabdirect/abstract/19961904773Beven K. 2006. A manifesto for the equifinality thesis. Journal of Hydrology 320 (1): 18–36DOI: 10.1016/j.jhydrol.2005.07.007Bolch T, Menounos B, Wheate R. 2010. Landsat-based inventory of glaciers in west-ern Canada, 1985–2005. Remote Sensing of Environment 114 (1): 127–137 DOI:10.1016/j.rse.2009.08.015Braithwaite RJ. 1995. Positive degree-day factors for ablation on the Greenland ice sheetstudied by energy-balance modelling. Journal of Glaciology 41 (137): 153–160 DOI:10.3189/S0022143000017846Carnahan E, Amundson JM, Hood E. 2019. Impact of glacier loss and vegetation successionon annual basin runoff. Hydrol. Earth Syst. Sci.: 15Casassa G, López P, Pouyaud B, Escobar F. 2009. Detection of changes in glacial run-offin alpine basins: examples from North America, the Alps, central Asia and the Andes.Hydrological Processes 23 (1): 31–41 DOI: https://doi.org/10.1002/hyp.7194Chernos M, MacDonald RJ, Nemeth MW, Craig JR. 2020. Current and future projectionsof glacier contribution to streamflow in the upper Athabasca River Basin. CanadianWater Resources Journal / Revue canadienne des ressources hydriques 45 (4): 324–344DOI: 10.1080/07011784.2020.1815587Clarke GKC, Jarosch AH, Anslow FS, Radić V, Menounos B. 2015. Projected deglaciationof western Canada in the twenty-first century. Nature Geoscience 8 (5): 372–377 DOI:10.1038/ngeo2407Craig JR, Brown G, Chlumsky R, Jenkinson RW, Jost G, Lee K, Mai J, Serrer M, SgroN, Shafii M, et al. 2020. Flexible watershed simulation with the Raven hydrolog-ical modelling framework. Environmental Modelling & Software 129: 104728 DOI:10.1016/j.envsoft.2020.104728Duethmann D, Bolch T, Farinotti D, Kriegel D, Vorogushyn S, Merz B, Pieczonka T, JiangT, Su B, Güntner A. 2015. Attribution of streamflow trends in snow and glacier melt-dominated catchments of the Tarim River, Central Asia. Water Resources Research 5145(6): 4727–4750 DOI: https://doi.org/10.1002/2014WR016716Eaton BC, Moore RD. 2010. Regional hydrology., Pike RG, Redding TE, Moore RD, Win-kler RD, Bladon KD (eds).B.C. Ministry of Forest; Range, Forest Science Program,Victoria, B.C.; FORREX Forum for Research; Extension in Natural Resources: Kam-loops; 85–110.Fountain AG, Tangborn WV. 1985. The Effect of Glaciers on Streamflow Variations. WaterResources Research 21 (4): 579–586 DOI: https://doi.org/10.1029/WR021i004p00579Fountain AG, Walder JS. 1998. Water flow through temperate glaciers. Reviews of Geo-physics 36 (3): 299–328 DOI: https://doi.org/10.1029/97RG03579Frans C, Istanbulluoglu E, Lettenmaier DP, Fountain AG, Riedel J. 2018. Glacier Recessionand the Response of Summer Streamflow in the Pacific Northwest United States, 1960-2099. Water Resources Research 54 (9): 6202–6225 DOI: 10.1029/2017WR021764Frenierre JL, Mark BG. 2014. A review of methods for estimating the contribution ofglacial meltwater to total watershed discharge. Progress in Physical Geography: Earthand Environment 38 (2): 173–200 DOI: 10.1177/0309133313516161Gao H, He X, Ye B, Pu J. 2012. Modeling the runoff and glacier mass balance in a smallwatershed on the Central Tibetan Plateau, China, from 1955 to 2008. HydrologicalProcesses 26 (11): 1593–1603 DOI: 10.1002/hyp.8256Gascoin S, Kinnard C, Ponce R, Lhermitte S, MacDonell S, Rabatel A. 2011. Glaciercontribution to streamflow in two headwaters of the Huasco River, Dry Andes of Chile.The Cryosphere 5 (4): 1099–1113 DOI: 10.5194/tc-5-1099-2011Hamilton AS, Hutchinson DG, Moore RD. 2000. Estimating Winter Streamflow UsingConceptual Streamflow Model. Journal of Cold Regions Engineering 14 (4): 158–175DOI: 10.1061/(ASCE)0887-381X(2000)14:4(158)Hedstrom NR, Pomeroy JW. 1998. Measurements and modelling of snow interception in theboreal forest. Hydrological Processes 12 (10-11): 1611–1625 DOI: 10.1002/(SICI)1099-1085(199808/09)12:10/11<1611::AID-HYP684>3.0.CO;2-4Helsel DR, Hirsch RM. 2002. Statistical methods in water resources. Reston, VA. Availableat: http://pubs.er.usgs.gov/publication/twri04A3Hock R. 2003. Temperature index melt modelling in mountain areas. Journal of Hydrology282 (1): 104–115 DOI: 10.1016/S0022-1694(03)00257-9Hock R. 2005. Glacier melt: A review of processes and their modelling. Progress in PhysicalGeography: Earth and Environment 29 (3): 362–391 DOI: 10.1191/0309133305pp453raHopkinson C, Young GJ. 1998. fountain. Hydrological Processes 12 (10-11): 1745–1762DOI: https://doi.org/10.1002/(SICI)1099-1085(199808/09)12:10/11%3C1745::AID-HYP692%3E3.0.CO;2-SHugonnet R, Brun F, Berthier E, Dehecq A, Mannerfelt ES, Eckert N, Farinotti D. 2022.Uncertainty analysis of digital elevation models by spatial inference from stable terrain.IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing15: 6456–6472 DOI: 10.1109/JSTARS.2022.3188922Huss M. 2011. Present and future contribution of glacier storage change to runofffrom macroscale drainage basins in europe. Water Resources Research 47 (7) DOI:10.1029/2010WR010299Huss M. 2013. Density assumptions for converting geodetic glacier volume change to masschange. The Cryosphere 7 (3): 877–887 DOI: 10.5194/tc-7-877-201346Huss M, Hock R. 2015. A new model for global glacier change and sea-level rise. Frontiersin Earth Science 3 DOI: 10.3389/feart.2015.00054Huss M, Hock R. 2018. Global-scale hydrological response to future glacier mass loss.Nature Climate Change 8 (2): 135–140 DOI: 10.1038/s41558-017-0049-xJansson P, Hock R, Schneider T. 2003. The concept of glacier storage: a review. Journalof Hydrology: 14Jost G, Moore RD, Menounos B, Wheate R. 2012. Quantifying the contribution of glacierrunoff to streamflow in the upper Columbia River Basin, Canada. Hydrology and EarthSystem Sciences 16 (3): 849–860 DOI: 10.5194/hess-16-849-2012Koboltschnig GR, Schöner W, Holzmann H, Zappa M. 2009. Glaciermelt of a small basincontributing to runoff under the extreme climate conditions in the summer of 2003.Hydrological Processes 23 (7): 1010–1018 DOI: 10.1002/hyp.7203Lemke P, Ren JF, Alley RB, Allison I, Carrasco J, Flato G, Fujii Y, Kaser G, Mote P,Thomas RH, et al. 2007. Observations: Changes in snow, ice and frozen ground. In:Climate change 2007: The physical science basis. Contribution of working group i tothe fourth assessment report of the intergovernmental panel on climate change. ClimateChange 2007: The Physical Science BasisLuce CH, Abatzoglou JT, Holden ZA. 2013. The Missing Mountain Water: Slower West-erlies Decrease Orographic Enhancement in the Pacific Northwest USA. Science 342(6164): 1360–1364 DOI: 10.1126/science.1242335Mark BG, McKenzie JM, Gómez J. 2005. Hydrochemical evaluation of changing glaciermeltwater contribution to stream discharge: Callejon de Huaylas, Peru / Evalua-tion hydrochimique de la contribution évolutive de la fonte glaciaire à l’écoulementfluvial: Callejon de Huaylas, Pérou. Hydrological Sciences Journal 50 (6) DOI:10.1623/hysj.2005.50.6.975Matott LS. 2017. STRICH: An optimization software tool, documentation and user’s guide,version 17.12.19. 79 pages, university at buffalo center for computational research Avail-able at: http://www.civil.uwaterloo.ca/envmodelling/Ostrich.htmlMenounos B, Hugonnet R, Shean D, Gardner A, Howat I, Berthier E, Pelto B, TennantC, Shea J, Noh M-J, et al. 2019. Heterogeneous Changes in Western North AmericanGlaciers Linked to Decadal Variability in Zonal Wind Strength. Geophysical ResearchLetters 46 (1): 200–209 DOI: 10.1029/2018GL080942Moore RD, Fleming SW, Menounos B, Wheate R, Fountain A, Stahl K, Holm K,Jakob M. 2009. Glacier change in western North America: influences on hydrology,geomorphic hazards and water quality. Hydrological Processes 23 (1): 42–61 DOI:10.1002/hyp.7162Moore RD, Pelto B, Menounos B, Hutchinson D. 2020. Detecting the Effects of SustainedGlacier Wastage on Streamflow in Variably Glacierized Catchments. Frontiers in EarthScience 8: 136 DOI: 10.3389/feart.2020.00136Moore RD, Woods RA, Boyle DP. 2013. Putting PUB into Practice in Mountainous Re-gions. 15 (2)Moyer AN, Moore RD, Koppes MN. 2016. Streamflow response to the rapid retreat of alake-calving glacier. Hydrological Processes 30 (20): 3650–3665 DOI: https://doi.org/10.1002/hyp.10890Nolin AW, Phillippe J, Jefferson A, Lewis SL. 2010. Present-day and future contributions47of glacier runoff to summertime flows in a Pacific Northwest watershed: Implicationsfor water resources. Water Resources Research 46 (12) DOI: 10.1029/2009wr008968Pelto MS. 2011. Skykomish River, Washington: Impact of ongoing glacier retreat onstreamflow. Hydrological Processes 25 (21): 3356–3363 DOI: 10.1002/hyp.8218Pojar J, Klinka K, Meidinger DV. 1987. Biogeoclimatic ecosystem classification in BritishColumbia. Forest Ecology and Management 22 (1-2): 119–154 DOI: 10.1016/0378-1127(87)90100-9Radic V, Hock R. 2014. Glaciers in the earth’s hydrological cycle: Assessments of glaciermass and runoff changes on global and regional scales. Surveys in Geophysics 35 (3):813–837 DOI: http://dx.doi.org/10.1007/s10712-013-9262-yRothlisberger H, Lang H. 1987. In glacio-fluvial sediment transfer: An alpine perspective,wiley-interscience (pp. 207–284). Glacial hydrology.Schiefer E, Menounos B, Wheate R. 2007. Recent volume loss of British Columbian glaciers,Canada. Geophysical Research Letters 34 (16) DOI: 10.1029/2007GL030780Schnorbus M, Werner A, Bennett K. 2014. Impacts of climate change in three hydrologicregimes in British Columbia, Canada. Hydrological Processes 28 (3): 1170–1189 DOI:https://doi.org/10.1002/hyp.9661Shea JM, Moore RD, Stahl K. 2009. Derivation of melt factors from glacier mass-balance records in western Canada. Journal of Glaciology 55 (189): 123–130 DOI:10.3189/002214309788608886Stahl K, Moore RD, Floyer JA, Asplin MG, McKendry IG. 2006. Comparison of approachesfor spatial interpolation of daily air temperature in a large region with complex topogra-phy and highly variable station density. Agricultural and Forest Meteorology 139 (3-4):224–236 DOI: 10.1016/j.agrformet.2006.07.004Stahl K, Moore RD, Shea JM, Hutchinson D, Cannon AJ. 2008. Coupled modellingof glacier and streamflow response to future climate scenarios: MODELLINGOF GLACIER AND STREAMFLOW. Water Resources Research 44 (2) DOI:10.1029/2007WR005956Stenborg T. 1970. Delay of Run-off from a Glacier Basin. Geografiska Annaler: Series A,Physical Geography 52 (1): 1–30 DOI: 10.1080/04353676.1970.11879806Thayyen RJ, Gergan JT. 2010. Role of glaciers in watershed hydrology: a preliminarystudy of a "Himalayan catchment" The Cryosphere 4 (1): 115–128 DOI:10.5194/tc-4-115-2010Tiel M van, Kohn I, Loon AFV, Stahl K. 2020. The compensating effect of glaciers:Characterizing the relation between interannual streamflow variability and glacier cover.Hydrological Processes 34 (3): 553–568 DOI: https://doi.org/10.1002/hyp.13603Tolson BA, Shoemaker CA. 2007. Dynamically dimensioned search algorithm for compu-tationally efficient watershed model calibration. Water Resources Research 43 (1) DOI:10.1029/2005WR004723Wheler BA, MacDougall AH, Flowers GE, Petersen EI, Whitfield PH, Kohfeld KE. 2014.Effects of Temperature Forcing Provenance and Extrapolation on the Performance of anEmpirical Glacier-Melt Model. Arctic, Antarctic, and Alpine Research 46 (2): 379–393DOI: 10.1657/1938-4246-46.2.379Zemp M, Hoelzle M, Haeberli W. 2009. Six decades of glacier mass-balance observations:a review of the worldwide monitoring network. Annals of Glaciology 50 (50): 101–11148DOI: 10.3189/17275640978776959149Appendix AThe Leaf Area Index (LAI) used in the model were obtained from Copernicus(https://land.copernicus.eu/global/products/LAI). Figure A.1 shows LAI distribu-tion in Illecillewaet catchment on 20th August 2018 and 2019. Figure A.2 shows the rangeof LAI for each land class of these two days. In the model, we set the LAI for bare, glacier,ESSF and ICH as 0.1, 0, 3.2 and 4.0 relatively.Figure A.1: LAI distributions in the Illecillewaet catchment on 20th August 2018 (left) and2019 (right).50Figure A.2: Boxplot of LAI values for each land class on 20th August 2018 (top) and 2019(bottom).51
Thesis/Dissertation
eng
Geography
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Moore, R. Daniel
Diagnosing the relative contributions of glacier retreat and climatic variability to streamflow trends in Illecillewaet catchment
Text