ESTIMATION OF SNOW WATER EQUIVALENT IN MOUNTAINOUS TERRAIN
USING AIRBORNE LASER ALTIMETRY AND EMPIRICAL MODELLING, COAST
MOUNTAINS, B.C., CANADA
by
Sergio C. Vazquez Tagle Gallegos
B. Eng., Autonomous University of the State of Hidalgo, 2018

THESIS SUBMITTED IN PARTIAL FULFILLMENT OF
THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE
IN
NATURAL RESOURCES AND ENVIRONMENTAL STUDIES

UNIVERSITY OF NORTHERN BRITISH COLUMBIA
April 2024
© Sergio Vazquez Tagle Gallegos, 2024

Abstract
This research aims to enhance scholarly understanding of snow dynamics, the remotely
sensed snowpacks, and the calculation of basin-wide snow water equivalent in mountainous
terrain.
Mountain snow is a critical source of meltwater. However, forecasting snow
distributions and total water equivalence in mountain basins is limited due to complex terrain,
challenging environmental conditions, and lack of observations. Laser altimetry can provide
detailed observations of snow depth, but an estimate of snow density is required to evaluate the
total basin water equivalence. This study uses laser altimetry surveys and empirically modelled
snow densities to estimate mountainous basin-wide snow-water equivalent (SWE).
Between 2017 and 2020, seven laser altimeter surveys during late winter and spring
were conducted in the LaJoie Basin, Coast Mountains, British Columbia (B.C.), a strategic
hydroelectric power reservoir. The laser-derived snow depths averaged between 1.4 and 2.1 m
for non-glacierized terrain, while glacierized terrain weighted averages ranged between 2.2 and
5.4 m.
The laser-derived depths were combined with empirical snow density models to derive
distributed SWE for the Lajoie Basin. Ten linear and three non-linear snow density empirical
models were tested and developed, from which (a) a snow course multi-parameter, non-linear
relation and (b) snow pillow robust (Huber loss) linear regressions yielded this study’s lowest
root mean squared errors (51.65 and 74.12 kg m-3, respectively).
For non-glacierized terrain, the multi-parameter, non-linear model produced basin-wide
SWE averages between 0.56 and 1.06 m.w.e. and propagated uncertainties from ± 0.1 to ± 0.14
ii

m.w.e. Conversely, glacierized terrain exhibited weighted SWE averages between 0.82 and 2.93
m.w.e., with estimates of uncertainty ranging from ± 0.31 to ± 0.52 m.w.e.
The robust linear regressions yielded non-glacierized SWE averages from 0.49 to 0.92 m.w.e.,
with uncertainties between ± 0.14 and ± 0.2 m.w.e. The weighted SWE averages in glacierized
terrain ranged between 0.77 and 1.92 m.w.e., with estimates of uncertainty between ± 0.29 and
± 0.45 m.w.e. The SWE estimates from lidar and modelled density are comparable to snow
pillow observations at the watershed, which demonstrates the effectiveness of these coupled
techniques and improves our forecasting capabilities.

iii

Table of Contents

Abstract ................................................................................................................................... ii
Table of Contents .................................................................................................................. iv
List of Tables ......................................................................................................................... vi
List of Figures ....................................................................................................................... vii
Acknowledgments .................................................................................................................. x
1.

2.

3.

Introduction .................................................................................................................... 1
1.1.

Dissertation motivation ............................................................................................. 1

1.2.

Thesis objectives and structure ................................................................................. 3

Literature Review .......................................................................................................... 4
2.1.

Snow nature and metamorphism ............................................................................... 4

2.2.

SWE in mountainous regions .................................................................................... 7

2.3.

Laser altimetry for snow depth.................................................................................. 8

2.4.

Empirical snow density modelling .......................................................................... 11

Study area and methods .............................................................................................. 13
3.1.

Study area ................................................................................................................ 13

3.2.

Snow depth: Lidar altimetry data ............................................................................ 17

3.2.1.

Data acquisition ............................................................................................... 17

3.2.2.

Trajectory processing ....................................................................................... 19

3.2.3.

Point cloud projection and filtering ................................................................. 20

3.2.4.

Bayesian alignment algorithms ........................................................................ 20
iv

3.2.5.

Ice dynamics .................................................................................................... 21

3.2.6.

Vertical uncertainty assessment ....................................................................... 22

3.3.

3.3.1.

Observational data preprocessing .................................................................... 25

3.3.2.

Existing empirical snow density models.......................................................... 30

3.3.3.

Least squared regressions ................................................................................ 33

3.3.4.

Multi-parameter, non-linear equation .............................................................. 35

3.3.5.

Model uncertainty assessment ......................................................................... 36

3.4.
4.

Basin-wide SWE estimations and uncertainty propagation .................................... 37

Results ........................................................................................................................... 40
4.1.

Lidar data summaries .............................................................................................. 40

4.1.1.

Snow depth over glaciers ................................................................................. 43

4.1.2.

Vertical snow accumulation gradients ............................................................. 45

4.1.3.

Orientation of snow depth ................................................................................ 46

4.1.4.

Uncertainty assessment .................................................................................... 47

4.2.

Snow density modelling .......................................................................................... 48

4.2.1.

Previous empirical snow density models ......................................................... 48

4.2.2.

Linear regression models from snow pillow data ............................................ 51

4.2.3.

Snow course-derived linear regressions........................................................... 53

4.2.4.

Multi-parameter, non-linear equation .............................................................. 55

4.2.5.

Model evaluation and selection ....................................................................... 56

4.3.

Snow water equivalent estimations ......................................................................... 59

4.3.1.

Basin-wide summary ....................................................................................... 59

4.3.2.

Point-based validation ...................................................................................... 66

4.4.
5.

Snow density: empirical modelling ......................................................................... 24

Uncertainty propagation .......................................................................................... 69

Discussion...................................................................................................................... 77
5.1.

Lidar-derived snow depth........................................................................................ 77

5.1.1.

General observations ........................................................................................ 77

5.1.2.

Comparative analysis of snow depth uncertainties .......................................... 78

5.2.

Empirical snow density modelling .......................................................................... 79
v

5.2.1.

Major findings .................................................................................................. 79

5.2.2.

Comparative analysis of density modelling errors........................................... 80

5.3.

Propagated basin-wide SWE and error ................................................................... 80

5.3.1.

Major findings .................................................................................................. 80

5.3.2.

Comparative analysis ....................................................................................... 81

5.4.

Study limitations ..................................................................................................... 82

6.

Conclusions ................................................................................................................... 86

7.

References ..................................................................................................................... 90

8.

Appendices .................................................................................................................... 98

List of Tables
3.1: Biogeoclimatic zones in the LaJoie, from MacKenzie & Meidinger (2018). ..................16
3.2: Lidar acquisition date, covered area, and point density. ..................................................18
3.3: Regression coefficients calibrated by Jonas et al. (2009) .................................................32
4.1: Minimum, maximum, median, and weighted average depth over glaciers. .....................43
4.2: Bias, standard deviation, and standard error identified through 24 pseudo-GCP. ...........47
4.3: Snow density model performance ....................................................................................57
4.4: Comparison of RMSE values for different sample sets ...................................................58
4.5: Non-glaciated terrain average snow depth (m.) and SWE (m.w.e) estimates ..................62
4.6: Average snow depth and SWE estimates weighted by glacier size .................................62
4.7: Integrated basin-wide SWE averages ...............................................................................62
4.8: Temporal distribution of the snow course observations by depth ....................................64
4.9: Snow pillow observations against SWE estimates (m.w.e.) with the snow course
derived multi-parameter, non-linear snow density model. ......................................................67
vi

4.10: Pillow observations compared with SWE estimates (m.w.e) with the snow pillowderived robust linear regressions model. .................................................................................68
4.11: SWE uncertainty ranges from the snow course derived non-linear equation.................75
4.12: Modelled SWE uncertainty ranges with the snow pillow derived linear regression. .....76

List of Figures
2.1: The snow crystal morphology ............................................................................................. 5
2.2: Phase and wavelength differences lights ............................................................................. 9
2.3: Lidar parameters required for range measurements .......................................................... 10
2.4: Lidar-derived cloud point and underlying hillshaded terrain model. ................................ 11
3.1: Elevation, the Downton Lake, and the glacier coverage at the LaJoie Watershed, Coast
Mountains, B.C. ........................................................................................................................ 13
3.2: Topographic slope within the LaJoie basin with the Downton Lake ................................ 14
3.3: Pixel frequency by terrain orientation in the LaJoie watershed ........................................ 15
3.4: Terrain orientation (i.e., aspect) in the LaJoie watershed. ................................................. 15
3.5: Biogeoclimatic zones in LaJoie ......................................................................................... 16
3.6: LaJoie watershed and lidar flight paths ............................................................................ 19
3.7: Reflectance colour-coded point cloud ............................................................................... 20
3.8: Same survey snow depth grids .......................................................................................... 21
3.9: Pseudo GCP and snow pillow locations in the study area ................................................. 23
3.10: Site selection for the pseudo ground control points. ....................................................... 23
3.11: Snow pillow and course locations ................................................................................... 26
3.12: Temporal distribution of provincial snow course samples .............................................. 27

vii

3.13: Provincial snow course density data ................................................................................ 28
3.14: Provincial snow course data spanning seven decades ..................................................... 29
3.15: Snow depth and density subsets ...................................................................................... 33
3.16: Linear regression performance after data filters .............................................................. 34
4.1: Lidar-derived snow depth over ice-free terrain. ................................................................ 41
4.2: Vertically corrected snow depth estimations ..................................................................... 42
4.3: Mean snow depth by survey, colour coded by survey. ...................................................... 43
4.4: Glacier polygons by their averaged snow depth; colour-coded by survey. ....................... 44
4.5: Vertical gradients of the observed snow depths at the LaJoie basin ................................. 45
4.6: Azimuthal distribution of snow depths in non-glacial portions of the basin ..................... 46
4.7: Observed snow densities compared to existing empirical snow density models .............. 49
4.8: Observed densities against model residuals; empirical snow density models................... 50
4.9: Observed snow densities compared to the snow pillow linear regression model ............. 51
4.10: Observed densities against model residuals; snow pillow linear regressions. ................ 52
4.11: Observed snow densities compared to the snow course linear regression ...................... 53
4.12: Observed densities against model residuals; snow course linear regressions. ................ 54
4.13: Observed snow densities compared to the non-linear equation model. .......................... 55
4.14: Observed densities compared to the model residuals from the non-linear equation. ...... 55
4.15: RMSE between the observed and the modelled snow densities ...................................... 57
4.16: Basin-wide SWE obtained with the non-linear equation density model ......................... 60
4.17: Basin-wide SWE from with the pillow-derived robust linear regressions ...................... 61
4.18: Mean snow water equivalent; multi-parameter, non-linear equation .............................. 63
4.19: Mean snow water equivalent; snow pillow derived robust linear regressions ................ 63

viii

4.20:Vertical gradient observed in the SWE estimates ............................................................ 65
4.21: Snow pillows at the LaJoie Basin. ................................................................................... 66
4.22: SWE observations compared to the non-linear model .................................................... 68
4.23: SWE observations compared to pillow linear regressions .............................................. 69
4.24: Snow course derived non-linear model – SWE uncertainty ............................................ 71
4.25: Snow course derived non-linear model – relative SWE uncertainty ............................... 72
4.26: Snow pillow derived robust linear regression model – SWE uncertainty ....................... 73
4.27: Snow pillow derived robust linear regression model – relative SWE uncertainty .......... 74
4.28: Histogram of SWE uncertainty in non-glacierized terrain; non-linear equation ............. 75
4.29: Histogram of SWE uncertainty in non-glacierized terrain; pillow-derived linear
regressions................................................................................................................................. 76
5.1: Densification pattern recorded by two snow pillows in one snow season ........................ 84

ix

Acknowledgments

I have been no one, and I have worn many faces. I have been naïve, but I conquered
many places. I felt looped in hell, disappeared like a small fish in a tank, and acted like less than
another brick in the wall.
The perfect song is framed with silence means you need darkness to appreciate the light. Things,
however, look brighter when you put God in the center of your life. If you forget where the Sun
is you can always ask the plants- our Mother watches over us all.
I danced through the night, but my roots kept my feet on the ground. My parents Pichur
and Pichon, my sisters Piscis and Peach, and my grandparents Mamita, Tito, Roi, Licha y el
Gordito became my beacon of light. My aunties, cousins, and the Gallegos fam would never let
me go astray: some say we have eyes everywhere.
Birthright could not bring me here, but many opportunities created this path. The
financial support from Mitacs, my employers, and my loved ones gave me the tools, but my
mentors – Joseph Shea and Brian Menounos- taught me how to fight back.
I’ll never know if they knew I was in hell, but they taught me to debate and respectfully disagree
with my demons. They kept it up even in my most uncertain times and looked at my many rabbit
holes with a smile and a potential way out. Three different North American instruments and just
look at the symphony we made.
I extend heartfelt gratitude to my esteemed committee members, Stephen Dery and Rita
Winkler, whose wisdom linked my contributions to a sky full of stars.

x

I thank those who stayed, and I bless those who left; Life is an ever-changing wave. I
thank my friends, and many many cats, for the reminders of who I am and for fueling my fire.
I thank Adam, Anna, Georg, and Frank for lending me a hand with their countless insights. I
thank Sara, Brayden, Fernanda, and Alejandra for offering me a roof, a meal, and a muchneeded tribe. I acknowledge the support from the Finance Department, while I thank my school
for my students, the lessons, and the time.
I am the third of my name, but the only one not in Mictlān. I long for my grandpa El
Gallo or my most beloved Dr. Emmett Brown, but I dreamed of future times. I forgot what these
stories look like, but I have my people, will, and hands. I will take any opportunities I have to
enrich this tapestry of life. Here and now.
Luckily for me, I live right inside the middle of the A of the Wonderland sign, and many peers
have led me to fly. Bottled gratitude to Mrs. Yola, Josefina Casillas, Alberto Blanco, Edgar
Uribe, Pam MacQuarrie, Patricia Montiel, and Rosa Maria Prol. Thanks for encouraging the
most talkative and scattered student you’ve had.
Hope and truth kept me through and made me finish this out of genuine love and lust
for life. This is no trick of the light; I will try to keep on making us proud. I will learn, but I
will, above all things, create.

“The magic, the wonder, the mystery, and the innocence of a child’s heart are the seeds
of creativity that will heal the world” (Jackson, 1993). I trust that water -from the oceans to the
ice- will lure the fire of those willing to fight.
xi

1.

Introduction
1.1.

Dissertation motivation

Snow is a sensitive indicator of climatic fluctuations, affecting local, regional, and
global climate. Ground snow affects the environment due to its high albedo, emissivity, and
low thermal conductivity (Brun et al., 2013; Huss et al., 2017).
The snow-covered area between Canada, Alaska, and Siberia reflects around two
zettajoules of shortwave energy per year (Flanner et al., 2011), an equivalent to five thousand
times the global energy consumption of 9 555 million tonnes of oil equivalent (International
Energy Agency, 2018). Recent history trends suggest smaller snow cover duration and extent
(Huss et al., 2017), progressively restricting the snow cover to mountainous terrain.
Mapping mountainous snowpacks is essential for hydrological modelling, water
availability forecasting, and efficient reservoir management. The estimated global mountain
runoff (12,184 km3/year) contributes approximately 31% of the net freshwater discharge
(38,136 km3/year; Meybeck et al., 2001). This resource is used to meet the operational
requirements of hydroelectric reservoirs (International Energy Agency, 2018) and the needs of
society.
Furthermore, mountain snowpacks are critical to regulating water, temperature,
nutrients, and sediment supply within mountainous watersheds. These ecosystems, in turn,
support vital industries such as fishing, farming, and forestry (Huss et al., 2017).
Accurate mountainous snow mapping becomes even more critical for climate
modelling, forecasting, and risk management (Pigeon & Jiskoot, 2008; Viviroli et al., 2007).
The escalating risks posed by the progressive tropospheric warming shown by ice cores,
1

thinner snowpacks and in-situ observations (Intergovernmental Panel on Climate Change,
2014) lead to higher hazard incidence, losses over the cryosphere and irreversible ecological
and economic damage (Lutz & Howarth, 2015).
The decline in Arctic sea ice, snow cover, and the associated rise in methane emissions
have imposed a cumulative cost on society of between $7.5 and $93 trillion (Euskirchen et al.,
2013). These impacts manifest themselves through global sea levels, hemispheric circulation
patterns, and the global efficiency to reflect solar radiation.
Accurate assessments of snowpack volumes and fluctuations enhance water allocation,
climate modelling, forecasting, and risk management efforts. Understanding the spatial
distribution of mountain snowpack dynamics is essential for effective climate change
adaptation strategies.
Laser altimetry, paired with estimates of snow density, provides viable alternatives to
fulfill snow mapping needs, particularly in areas lacking a snow observation network (Deems
et al., 2013). These advanced techniques can characterize mountainous snowpacks and
estimate basin-wide snow water content (SWE).
Leveraging these innovative approaches helps us better understand hydrologic
dynamics, enhance current resource management practices, and adapt to the challenges posed
by changing environmental conditions (Pigeon & Jiskoot, 2008). A comprehensive
mountainous snow mapping opens possibilities for research and facilitates informed decisionmaking processes.

2

1.2.

Thesis objectives and structure

The motivation of this study is to develop and assess methods to estimate basin-wide
SWE and its uncertainty in a glacierized watershed in southern British Columbia. The steps
required to complete this task include:
1. Acquiring and processing airborne lidar data.
2. Developing and testing empirical models of snow density.
3. Calculating distributed SWE and its associated uncertainties.

The first Chapter of this dissertation establishes its relevance and structure. Chapter 2
provides a concise review of relevant literature. Chapter 3 describes the study area, data, and
methods employed. Chapter 4 details the results from the lidar analysis, the empirical density
modelling, and the basin-wide SWE propagation, while Chapter 5 discusses and interprets
these findings.
Finally, Chapter 6 offers extensive conclusions, acknowledges limitations, and
suggests further research directions. The appendices contain specific ground control points,
glacier extents, snow pillows, and course data locations used in this study.

3

2.

Literature Review
2.1.

Snow nature and metamorphism

Snow consists of ice crystals formed in the atmosphere that precipitate and accumulate
in snow columns to form a snowpack. A snowpack is constituted of snow and pores that contain
air and water vapour, reaching a three-phase system once melt conditions are reached.
Atmospheric conditions and suspended particles determine the nature of the initial
snow crystals. The snow requires the collision of thermally contrasting moist masses of air.
Warm air masses get colder as they are lifted upward, allowing water droplets to condense.
The water droplets may reach a supercooled state by remaining liquid below 0 °C to nearly 40 °C.
The supercooled droplets eventually evaporate, but any suspended particles in the
atmosphere serve as a solid surface from which one droplet can freeze. The freezing
temperature depends on the suspended dust particle, but it starts generally with temperatures
colder than -6 °C.
A frozen droplet starts to attract water vapour, leading to surface deposition. The
deposition of the surrounding vapour forms a crystal-lattice structure fed by water vapour, i.e.,
a snowflake (Figure 2.1; Libbrecht, 2005).

4

Figure 2.1: The snow crystal morphology against water content and temperature. From Libbrecht
(2005)

The snow crystal array and the water, air, vapour, and ice content determine the
snowpack density. Snow stability is also influenced by its crystal nature, with weakly bonded
crystals contributing to unstable layers and the risk of avalanches. Additionally, complexshaped crystals (such as dendrites) tend to have more surface area and capture more air and
moisture (Dingman, 2002).
Typical densities of freshly fallen snow range from 30 to 300 kg m -3. The wind may
break snowflakes and pack them into dense layers, but it can otherwise be steady enough to
allow for interstitial pores. Generally, higher temperatures and wind speeds lead to denser
snowpacks than colder and more stable conditions (McKay & Findlay, 1971; Verschuren,
1972).

5

On a seasonal scale, snow metamorphism occurs from the snowflake deposition until
the snow completely melts. There are four different mechanisms for snow metamorphism
(Dingman, 2002):
•

Gravitational settling. Whenever the snowpack accumulates more snow, the

overlying stress of the snow column progressively increases. Gravitational settling can increase
density from 2 to 50 kg m-3 daily (Anderson & States, 1976).
•

Destructive metamorphism. Variations in vapour pressure over convex and

concave snow surfaces lead to destructive metamorphism. In convex surfaces, higher vapour
pressures lead to sublimation. The moisture is deposited over concave surfaces, leading to
progressively denser and more spherical grains over time (Dingman, 2002).
This metamorphism can increase the snow density by 1% per hour, but it becomes
negligible whenever densities exceed 250 kg m-3 (Anderson & States, 1976). Often, this
metamorphism can only be found in fresh snow (Dingman, 2002).
•

Constructive metamorphism. If the snow grains connect through sintering, a

snowpack’s effective porosity will be reduced. For greater distances, the temperature gradients
within the snowpack may cause sublimation at warmer portions and deposition in colder parts,
leading to changes in the snowpack density.
•

Melt metamorphism. Melt metamorphism is caused by liquid water added to

the snowpack. The thaw at the surface level or any percolated rain can cause layers of solid ice
across the snowpack once frozen. Additionally, changes in the mean grain size may result from
adding water, either by dissolving them or increasing their size.

6

All metamorphism processes lead to a progressive increase in density throughout the
season. Snow metamorphism largely depends on atmospheric conditions; the same area may
witness different snowpack properties over the years, but other parts of a watershed will
simultaneously undergo comparable metamorphism (Dingman, 2002).

2.2.

Snow in mountainous regions

Snow accumulates in a watershed through precipitation and is redistributed by wind or
gravity. The total snowfall a watershed may receive largely depends on the regional climate,
latitude, and continentality. In mountain watersheds, wind redistribution and avalanching
cause much of the spatial heterogeneity in snow depths and SWE (Elder et al., 1991; Vionnet
et al., 2021).
Snow ablation refers to the loss of snow through melting, sublimation, and erosion and
is controlled by radiative and energy exchanges between the snow and the atmosphere.
Shortwave and longwave radiation are the principal radiative inputs that must be considered
for SWE losses, followed by turbulent energy exchanges (Dozier et al., 1988).
Within a mountainous watershed, the forest cover significantly influences snow
accumulation and ablation rates (Varhola et al., 2010). In forested areas, snow accumulation
may be reduced by up to 40%, and the snowmelt rates are up to 70% lower than nearby
reference sites. The reduced accumulation can be attributed to the dense canopy, which
intercepts up to 60% of the snowfall, from which up to 40% gets sublimated (Hedstrom &
Pomeroy, 1998; Pomeroy et al., 1998).

7

The ablation rates are largely dependent on a snowpack’s energy budget. While canopy cover
prevents shortwave radiation from reaching the surface and reduces wind speed, effectively
minimizing the sensible and latent heat fluxes (Boon, 2009; Teti, 2008).
In both forested and open alpine areas, diurnal melt can percolate the snowpack and
refreeze, directly altering snow density and depth but not the SWE. In contrast, if an ice lens
is created at the surface level, the snow beneath can melt or displace, keeping a seemingly
constant snow depth (from ground level up to the recrystallized layer) but affecting the density
and SWE (Elder et al., 1991). Therefore, these magnitudes must be accounted for
independently.

2.3.

Laser altimetry for snow depth

The Light Detection and Ranging (Lidar) technique is a remote-sensing technology that
emits a photon and measures the time difference between the time of emission and its return
to the sensor after being reflected from a surface.
Laser light differs from white light because of its (a) temporal coherence – the light is
composed of a single colour; and because of (b) its spatial coherence - all light waves are sent
in the same phase (Figure 2.2; Csele, 2004).

8

Figure 2.2: Phase and wavelength differences between white (a), monochromatic (b) and laser (c)
lights. From Csele (2004).

The polarized light produces intense, focused beams of light that are powerful enough
to reach the surface and be reflected. The intensity of the reflection may be used to characterize
the nature of the surface (NOAA, 2012), while the distance to the surface is estimated through
the time of flight.
A photon’s time of flight changes with different atmospheric densities, making it
necessary to account for local atmospheric conditions (air humidity or pressure). The return
times also depend on the sensor-specific force, angular rate, and orientation. Thus, the scan
angle, platform altitude, beam divergence, and swath width must be known to transform the
laser returns into altitudinal data (Figure 2.3; Deems et al., 2013).

9

Figure 2.3: Lidar parameters required for range measurements (R); platform height (h); beam
divergence (γ); laser spot footprint (AL); swath width (SW) and amplitude (q). From Deems et al.
(2013).

The return times are processed with an inertial measurement unit (IMU) to correct positional
and atmospheric biases. The data is processed through a differential global positioning system
(GPS) triangulation to assign each measured distance into a fixed elevation on a coordinate
system (Deems et al., 2013), creating a cloud of points representing every reflected surface
(Figure 2.4; Painter et al., 2016).

10

Figure 2.4: Lidar-derived cloud point and underlying hillshaded terrain model. Modified from the
National Oceanic and Atmospheric Administration (NOAA) Coastal Services (2012).

Once classified, a cloud point is interpolated to create a terrain surface or digital
elevation model (DEM). Estimating the snow depth requires subtracting a bare-earth DEM
from a snow-covered DEM, as the snowpacks can be interpreted as changes in the surface
elevation above sea level (Deems et al., 2013).
To difference a snow lidar dataset from a bare-earth model, both datasets must be coregistered to remove bias (Nuth & Kääb, 2011) since horizontal offsets may compromise the
snow depth values.

2.4.

Empirical snow density modelling

Lidar snow depth grids must be paired with an estimate of snow density to estimate
snow water equivalence. In situ observations of snow density are sparse and can be challenging
11

to estimate over large areas (Frei et al., 2012). As an alternative, physical or empirical models
can be used for snow density estimation (Pistocchi, 2016).
While physically based models rely on fundamental physics laws and could be applied
to any watershed, the models often require specific variables such as wind speed, air
temperature or relative humidity. Most snow-covered areas lack a dense enough network of
sensors to retrieve these atmospheric parameters (McCreight & Small, 2014; Pistocchi, 2016).
Empirical snow density models, in contrast, rely on observations such as historical
snow depth, SWE information, and even simple meteorological factors. Sturm et al. (2010)
and Pelto (2020) contrast historic SWE against their snow depth counterpart, thus inferring the
snow density. Shea (2010) developed a model for snow density based on accumulated positivedegree days, accounting for progressive snow metamorphism. Empirical snow density models
may not be transferrable between regions, so empirical models should be tested and calibrated
with region-specific data.

12

3.

Study area and methods
Repeat lidar surveys from the LaJoie watershed, Gold Bridge, British Columbia, were

processed into vertically and horizontally corrected digital elevation models and snow depth
grids. The basin-wide snow depth grids are paired with provincially calibrated empirical snow
density models to produce SWE estimations. The uncertainty in the SWE estimations was
propagated from the snow depth and snow density-independent and random errors.

3.1.

Study area

Located in the southern Coast Mountains of British Columbia, the Lajoie watershed
spans 962.5 km2 (Figure 3.1). The watershed consists of a partially glacierized landscape (183
km2 or 19.01%; Appendix B) between 800 and 2400 m.a.s.l.

Figure 3.1: Elevation, the Downton Lake (black central polygon), and the glacier coverage (hatched
blue polygons) at the LaJoie Watershed, Coast Mountains, B.C. (UTM 10N; ESPG:32610). These
glacier polygons were used in this study and were extracted from the Randolph Glacier Inventory (RGI
Consortium, 2017).

13

The watershed is partially glacierized alpine terrain marked by pronounced elevation
gradients, with altitudes ranging from 1000 to 2000 m within less than 1 km. (Figure 3.2). The
seasonal runoff from this watershed ultimately flows to Downton Reservoir, which is part of
the provincial hydroelectric power generation network.

Figure 3.2: Topographic slope within the LaJoie Basin with the Downton Lake (black)

The predominant terrain orientation is towards the NE-E direction, with the fewest
slopes facing toward the W-NW direction (Figures 3.3 and 3.4).

14

Figure 3.3: Pixel frequency (Millions of 5 m pixels) by terrain orientation in the LaJoie watershed
from a bare earth terrain model generated from lidar data.

Figure 3.4: Terrain orientation (i.e., aspect) in the LaJoie watershed.

Beneath the treeline, the forests consist of montane trees, primarily Engelmann spruce
and subalpine fir. The watershed’s alpine terrain is considered Interior Mountain heather,
15

reflecting cold continental conditions with abundant moisture and snowfalls (Table 3.1, Figure
3.5; MacKenzie & Meidinger, 2018).

Figure 3.5: Biogeoclimatic zones in LaJoie, from MacKenzie & Meidinger (2018).
Table 3.1: Biogeoclimatic zones in the LaJoie, from MacKenzie & Meidinger (2018).
BGC ID

Zone name

ESSF

Engelmann Spruce-Subalpine Fir

MS

Montane Spruce

Subzone

Subzone description

Area
(km2)

dv

Dry Very Cold

119.75

xvp

Very Dry, Very Cold Parkland

0.18

xv

Very Dry Very Cold

0.02

dvw

Dry, Very Cold Parkland

130.58

mw

Moist Warm

17.12

dvp

Dry, Very Cold Parkland

258.03

mwp

Moist Warm Parkland

0.02

mw

Moist Warm

40.42

dc

Dry Cold

62.01

CMA

Coastal Mountain-heather Alpine

unp

Undifferentiated and Parkland

2.48

IMA

Interior Mountain-heather Alpine

unp

Undifferentiated and Parkland

0.09

un

Undifferentiated

283.58

IDF

Interior Douglas-Fir

xc

Very Dry Cold

2.22

dc

Dry Cold

44.91

un

Undifferentiated

1.10

BAFA

Boreal Altai Fescue Alpine

16

The Downton Lake is the uppermost hydroelectric reservoir on the Bridge River. The
lake has a surface of nearly 2400 ha with a maximum and mean depth of 80 and 30 m,
respectively. The Downton Lake acts as a seasonal reservoir since a) the mean water residence
time is 6.3 months, and b) the lake acts as the leading supplier for the adjacent La Joie
hydroelectric dam (Hirst, 1991).
The hydrologic resources of the basin are regulated by BC Hydro for hydroelectric
production. The Lajoie Dam was finished by 1957 with an installed capacity of 22 MW,
generating an average of 170 GWh/year (Hirst, 1991). This installed capacity satisfies 0.22%
of the energy consumed in the province as of 2017 (National Energy Board, 2019).

3.2.

Snow depth: Lidar altimetry data

Laser-altimetry surveys of a mountainous watershed in British Columbia were
processed into digital elevation models (DEMs) to reflect acquisition dates between 2017 and
2020.
3.2.1. Data acquisition
Eleven surveys (Table 3.2) were conducted with a Riegl 780 infrared (1064 nm) laserscanning system mounted on a fixed-wing aircraft (Nolan et al., 2015; Painter et al., 2016).

17

Table 3.2: Lidar acquisition date, covered area, and point density.
Survey date
(year/month/day)
2017-09-22

Watershed coverage
(%)
47

Point density
(Per m2)
1.50

2017-09-27

60

1.47

2017-05-04

52

0.87

2017-05-07

56

0.60

2018-02-19

100

1.48

2018-04-22

100

1.02

2019-03-18

100

1.05

2019-04-30

49

1.03

2019-05-04

55

1.11

2020-03-09

100

1.02

2020-05-10

100

1.03

The acquisition dates were chosen to quantify the snow depth during snow
accumulation and melt seasons, aiming to measure the maximum snow depth every year.
Acquisition dates additionally depended on suitable flying conditions, optimizing safety, and
minimizing atmospheric humidity, a common cause of lidar noise. Logistical considerations
(equipment malfunction or aircraft maintenance) sometimes required two-day surveys (Table
3.2).
The atmospheric conditions at the time of the survey were obtained from the nearest
weather stations (e.g., Pemberton, Lillooet). The aircraft’s position and roll, pitch and yaw
were recorded with an Applanix POS AV Global Navigation Satellite (GNSS) device and with
an inertial navigation (IMU) system, coincident with the laser scans (Deems et al., 2013). The
IMU trajectory was recorded in WGS84 ellipsoidal height.

18

3.2.2. Trajectory processing
The aircraft’s trajectory was post-processed with RiPROCESS 1.8.2 software (RIEGL
Laser Measurement Systems GmbH, 2017) and the Applanix PP-RTX software. Neither
ground GPS stations nor field measurements were required as these devices triangulate their
positions through a satellite system.
After the aircraft's position triangulation, the aircraft's position achieved error ranges
below 0.05 m. The corrected positions were projected into a continuous trajectory line (Figure
3.6).

Figure 3.6: LaJoie watershed and lidar flight paths, April 2018 survey. Darker line: Vancouver to
Pemberton, surveying half the watershed. Brighter line: return flight following the same path.

19

The aircraft’s post-processed flight trajectory and roll, pitch and yaw are stored in a
series of Position & Orientation (POF) files and later used to produce the laser point clouds.
3.2.3. Point cloud projection and filtering
At this stage, and before creating a gridded array, the point clouds contained artifacts
where the laser pulse interacted with features not associated with a surface, such as water
vapour, smoke, birds, or insects (Figure 3.7).

Figure 3.7: Reflectance colour-coded point cloud after applying a -12 dB reflectance filter (Yellow
indicates noise). Manual cleaning was required for the heavily noised 2017-May 7th survey.

Noise was identified and removed using a threshold based on the intensity of a given laser
return. Trial and error identified this lower threshold between -12 and -15 dB depending on a
given survey. Noise can also be a by-product of inexistent “erroneously reflected” beams
(Figure 3.7). As a result, manual noise filtering was also necessary.
3.2.4. Bayesian alignment algorithms
The lidar point clouds were paired with their corrected flight trajectories and
orthogonally corrected positions (POF files) to perform a ‘Bayes alignment.’ This process

20

involves a series of Bayesian adaptive sequence alignment algorithms that identify and
minimize the vertical differences between overlapping LAS flight lines (e.g., Figure 3.6;
Jalobeanu et al., 2014; BayesMap Solutions LLC., 2020).
The Bayesian alignment eliminated systematic bias caused by vertical or horizontal
offsets between subsequent flight lines. These offsets can arise, for example, due to
acquisitions that occur over multiple days (Figure 3.8).

Figure 3.8: Same survey snow depth grids: a) Unaligned and b) Bayesian aligned. Misalignment
disrupts continuity in (a). These images show vertically uncorrected values (section 3.2.6).

After the flight line alignment, each dataset was exported into a Digital Elevation
Model. For uniformity, all surveys were exported into a Universal Transverse Mercator
projection (UTM 10N; ESPG:32610) with a spatial resolution of 5 m.
3.2.5. Ice dynamics
Snow depths cannot be calculated on individual glacierized pixels due to glacier
dynamics, which redistribute mass from high elevations to low elevations. However, the effect
21

of submergence or emergence is zero when integrated over a glacier surface (Mortezapour et
al., 2020; Pelto et al., 2019).
Thus, the average snow depth of an individual glacier can be calculated considering
negligible ice loss. Based on the Randolph Glacier Inventory 6.0 (Appendix A; RGI
Consortium, 2017), there are 104 unique glaciers for the Lajoie. These polygons delimit the
average snow depth for each survey over the glacierized terrain.
3.2.6. Vertical uncertainty assessment
Standard co-registration (Nuth & Kääb, 2011) uses bare-earth datasets and compares
their aspect and slope to remove elevation biases. However, the slope of a snow-covered model
differs from the underlying bare earth because of unevenly distributed snow (by gravity, wind,
or ablation).
Thus, snow-covered scenes cannot be co-registered with bare-earth scenes. Instead, a
vertical bias can be estimated by finding snow-free pixels in all scenes and using these as
pseudo-ground control points to correct the elevation bias. Once vertically corrected, one
summer survey (bare-earth DEM) was subtracted from all the winter surveys to yield basinwide snow depth grids.
Elevation models should reveal no vertical differences in stable, snow-free terrains
such as ridges, mountain peaks or wind-swept locations. Multispectral scenes from Landsat-7
OLI and Sentinel2 were used to identify 24 (Figure 3.9) snow-free sites in the watershed.
The supervised selection of points (pers. comm, B. Menounos) prioritized a distributed
number of pseudo-GCP within the watershed (Appendix A and Figure 3.10).

22

Figure 3.9: Pseudo GCP and snow pillow locations in the study area. Sentinel-2 (ESA) image courtesy
of the U.S. Geological Survey.

Figure 3.10: Site selection for the pseudo ground control points.

For each of the 24 pseudo ground control points, the median elevation among all
surveys was used as the reference elevation. The elevation bias of each survey was then
23

calculated as the median difference between the reference elevations and the measurements of
the individual study. The standard error of the elevation biases is used here as a conservative
estimate of the snow depth uncertainty.
A negative bias indicates that the pseudo-GCP heights are below the reference median
elevation.

3.3.

Snow density: empirical modelling

The densification of a snowpack is influenced by several energy and precipitation
factors including, but not limited to, soil and landcover properties, precipitation mass, air
temperature, atmospheric pressure, wind speed, and incoming long and shortwave radiation
(Hedrick et al., 2018). The main downside inherent to physically based models comes from
the availability of the required inputs since most snow-covered areas often lack a densely
enough network of sensors to retrieve them (McCreight & Small, 2014; Pistocchi, 2016).
Practical observations show that the range of observed snow depths is greater than the bulk
density counterpart (Steppuhn, 1976). Consequently, modelling the conservative parameter
(snow density) and measuring the dynamic parameter (snow depth), is perhaps the most
practical and potentially accurate way to calculate SWE without enough energy, vegetative or
precipitation observations (Sturm et al., 2010).
Three different approaches were used to empirically estimate snow densities, which
are then used to calculate SWE for each non-glacierized grid point, average snow depths and
SWE for each glacier, and a basin-wide SWE value. The first snow density approach consists
of previously published empirical models (Hill et al., 2019; Jonas et al., 2009; Pistocchi, 2016).

24

The second approach estimates snow density on the date of each survey through a least
square regression from synchronous provincial observations, as in Pelto et al. (2019) and Sturm
et al. (2010).
The third approach uses local data to fit five parameters in a non-linear equation created
to represent physical processes of snow densification from the start of the season to the date of
maximum snowmelt.
All density estimation approaches were validated with provincial snow course and
pillow data. The root mean squared error (RMSE) between the modelled and observed snow
density was used to determine the error in snow density.
3.3.1. Observational data preprocessing
The empirical snow density models require calibration with local data to maintain
regional representability. The observational data relevant to this study comprises the snow
pillow network and the historic snow course measurements available from the Province of
British Columbia (Figure 3.11).

25

Figure 3.11: Snow pillow and course locations (BC, 2021b, 2021a). Yellow boxes contain coordinates
in WGS 1984 (EPSG: 4326). Axis ticks represent UTM 10N zone limits (EPSG: 32610).

Provincial snow course data are manual observations of SWE and snow depth
(Province of British Columbia, 2021b) collected at specific locations since 1933.
Snow course observations are sampled monthly or bi-monthly between late December
to June (Figure 3.12)

26

Figure 3.12: Temporal distribution of provincial snow course samples (BC, 2021b). The x-axis
shows 30-day intervals.

For data consistency, the snow depth (cm) and the water equivalent (mm) were
converted into meters to estimate snow density (Equation 3.1).

ߩ‫= ݓ݋݊ݏ‬

(ߩ‫)ܧܹܵ ∗ ݎ݁ݐܽݓ‬
‫ݐ݌݁݀ ݓ݋݊ݏ‬ℎ

(3.1)

Snow density values beyond three standard deviations from the snow course density mean
(57.1 – 613.6 kg m-3; Figure 3.13) were considered outliers and removed from the dataset. This
criterion was also applied to the snow pillow dataset.

27

Figure 3.13: Provincial snow course density data, colour-coded by 1, 2, and 3 sigma thresholds (BC,
2021a).

A preliminary analysis of the snow density data (Figure 3.14) shows the expected
seasonal increase in snow density from January to May. This analysis uses observations from
across the province. The data show no noticeable changes in the rate of snow densification
between decades, so simple empirical snow density models based on the day of the year may
be helpful.

28

Figure 3.14: Provincial snow course data spanning seven decades (BC, 2021a). Continuous lines:
decadal median snow density. Dashed lines: interquartile ranges.

To have a suitable timeframe, the snow course dataset was reduced to match the earliest
snow pillow observations available (October 10, 2003) up to a week after the last lidar survey
(Appendix C) May 20, 2020; n= 11,230). Half of the samples were randomly chosen for model
calibration (section 3.3.3 and 3.3.4; n= 5615), while the remaining values were used as a
validation dataset (section 3.3.5; n= 5615).
The provincial automated snow pillows (Automated Snow Weather stations, ASWS)
were also used for model calibration (section 3.3.3; section 3.3.4). ASWS records hourly air
temperature, snow depth, SWE, and other variables. ASWS snow density data were retrieved
from the Aquarius Time-Series database (Province of British Columbia, 2021a), filtered for
outliers, and pre-processed into a format compatible with the snow density models and snow
course data.

29

For the snow density calculations and model development, ASWS sites with SWE and
depth (Appendix D) entries were retained.
Hourly snow depth (cm) and water equivalent (mm) entries were converted into meters.
The snow pillow temporal resolution was then reduced by averaging the snow depth and water
equivalent by date so the snow pillow timestamp (year-month-day hour: minute: second; UTC)
matches the snow course observations format (year-month-day).
After the stations were merged, the data frame was transposed, maintaining the date
column for chronologic subsetting. The snow density was estimated from the average daily
snow depth and SWE following Equation 3.2.
With data from 51 stations, a daily-averaged dataset from October 10, 2003, to May
20, 2020 (n=142,451) was assembled. This timespan includes entries since the earliest snow
pillow observations were available up to a week after the last lidar survey used in this study.

3.3.2. Existing empirical snow density models
The first approach aimed to evaluate published empirical models to estimate snow
density. Three models that use time and elevation are examined here (Jonas et al., 2009;
Pistocchi, 2016; Hill et al., 2019) and tested with their calibrated parameters.
Hill et al. (2019)
The SWE model developed by Hill et al. (2019) was derived from telemetric and field
data from British Columbia and the western United States. Most of the dataset was provided

30

by the Natural Resources Conservation Service in the United States and the Automated Snow
Weather Stations operated by the Water Management Branch of British Columbia.
The model proposed by Hill et al. (2019) uses a linear equation to estimate SWE from
snow depth (h) and the number of days (DOY) since October 1st, which represents the
beginning of the snow season. This model (Equation 3.2) uses two different densification
parameters during accumulation (SWEacc) and the ablation (SWEabl) periods, which are
identified by the day of peak SWE (DOYp):
1
1
ܹܵ‫( · ܿܿܽܧܹܵ = ܧ‬1 − tan ℎ [0.01{‫ ܻܱܦ‬− ‫ )]} ܻܱܲܦ‬+ ܹܵ‫( · ݈ܾܽܧ‬1 + tan ℎ [0.01{‫ ܻܱܦ‬− ‫)]} ܻܱܲܦ‬
2
2

(3.2)

The SWE values modelled through the Hill et al. (2019) equation were converted into
snow densities to compare with other models used in this study.
Pistocchi et al. (2016)
The model proposed by Pistocchi et al. (2016) estimates the snow density based on the
depth and the number of days since the start of the snow season (defined as November 1 in this
study). The simple regression and its coefficients were derived from 16 telemetric stations
operated by the Autonomous Province of Bolzano in the Italian Alps.
The equation provided by Pistocchi et al. (2016) suggests an initial snow density (ρ0)
of 200 kg m-3 and an initial densification rate (K) of 1 kg m-3 per day. The densification
parameter multiplies the days since November 1st, which is assumed to be the start of the snow
year (DOY). This equation describes the observed densification pattern shown by snow
samples (Equation 3.3).
31

(3.3)

ߩ‫ߩ = ݓ݋݊ݏ‬0 + ‫ ܻܱܦ(ܭ‬+ 61)

Jonas et al. (2009)
The equation provided by Jonas et al. (2009) estimates snow density from snow depth
and two additional variables: altitude and time-dependent coefficients (Equation 3.4). The
model was produced from 36 snow telemetry stations operated by the WSL Institute for Snow
and Avalanche Research in Switzerland:
(3.4)

ߩ‫ݐ݌݁݀ ݓ݋݊ݏ ∗ ܽ = ݓ݋݊ݏ‬ℎ + ܾ

Coefficients a and b are fitted for three different elevation ranges and months during the
accumulation season (Table 3.3).
Table 3.3: Regression coefficients calibrated by Jonas et al. (2009; Equation 3.4)
Elevation between 0 and 1400 m.a.s.l.

a
b

Jan

Feb

Mar

April

May

JuneOct.

Nov

Dec

235
31

279
9

333
3

347
25

413
19

0
0

149
37

201
26

Nov

Dec

149
37

201
26

Nov

Dec

149
37

201
26

Jan
a
b

235
31

Jan
a
b

235
31

Elevation between 1400 and 2000 m.a.s.l.
JuneFeb
Mar
April
May
Oct.
279
9

333
3

347
25

413
19

0
0

Elevation between 2000 and 3000 m.a.s.l.
JuneFeb
Mar
April
May
Oct.
279
9

333
3

347
25

32

413
19

0
0

3.3.3. Least squared regressions
Provincial snow pillow and snow course observations of snow depth and SWE were
used to create least squares regression models between density and snow depth, a relation
explored by Pelto et al. (2019). Robust and ordinary least-square regressions were used to
describe the snow behaviour at a regional scale on any given day (Figure 3.15).

Figure 3.15: Snow depth and density subsets (March 05, 2019) from the snow pillow and snow course
observations.

Data outliers can drastically bias the estimators of a least square regression (Rousseeuw
& Leroy, 1987). Thus, once the chronologic subsets were retrieved (Figure 3.16a), the density
values beyond historically typical values (57.1 – 613.6 kg m-3; section 3.3.1; Figure 3.16b) and
statistical outliers within the subset (values beyond a 3-sigma score; Figure 3.16c) were
removed from the analysis.

33

Figure 3.16: Linear regression performance after data filters: a) unfiltered, positive data b) filtered by
density range, c) values within three sigma score. This real-time data subset (2020-03-08) needed data
filters to be trusted. Navy line: Ordinary least squares regression. Aquamarine line: Robust Huber loss
least square regression.

With each subset, both an ordinary least square regression and a robust Huber’s Mestimator least square regression (Huber, 1973) were used to retrieve a density slope from snow
depth and SWE, following Pelto et al. (2019) and Jonas et al. (2009). Linear regressions
(Equation 3.5) predict a certain amount (p) of unknown parameters (θ1, …, θp) from
observations (x1, …, xn) with correlation (Equation 3.6). A regression equation requires
coefficients (cij) and close-to-evenly distributed independent random errors (Ui) to make
estimates.

‫݌‬

‫ = ݅ݔ‬෍

݆ =1

݆ܿ݅ ∗ ߠ݅ + ܷ݅

(3.5)

Regression estimators were retrieved from the snow depth and SWE datasets through ordinary
least squares (OLS) regressions, i.e., by minimizing the sum of the squared distance between
the prediction line and the observations. OLS regressions, however, assume the absence of data
outliers, impacting the regression estimators and unevenly distributed residuals (Rasheed et
al., 2014).

34

Instead, the robust least squares (RLS) regressions require weights on their estimators
to enhance the predictions. The Huber loss function considers the mean squared error function
and the absolute value, proposes a tuning constant (k = 1.345σ) for maximum estimator
likelihood in ordinary cases, and offers protection against outliers.
The ordinary and robust linear squares regressions were calculated from the calibration
datasets through the openly available Python and shell scripts compiled by Seabold and
Perktold (2010). This snow density estimation scales density as a function of snow depth and
can be used directly with lidar-derived snow depth data.
3.3.4. Multi-parameter, non-linear equation
The third approach to estimating snow density uses a non-linear relation that expects
snow density to increase with depth and cumulative shortwave radiation.
The proposed model (Equation 3.6) uses five tuning parameters, snow depth and the
DOY, to model snow density. The parameters were fitted to provincial snow course data. The
empirical function was built assuming that the most critical factors controlling snow density
are the weight of the snow (snow depth) and energy transfers to the snowpack from sensible
heat (air temperature) and shortwave radiation.
1

ߩ‫݌ = ݓ݋݊ݏ‬1 + ‫ݐ݌݁݀ ݓ݋݊ݏ‬ℎ ∗ ‫݌‬2 + ‫݌‬3 ∗ ݁ ቀ−2‫݌‬5

2 (‫ܻܱܦ‬−‫) ݌‬2 ቁ
4

(3.6)

The relation between snow depth and density was assumed to be linear, whereas the
connection between density and energy inputs depends on the day of year.
The coefficients in Equation 3.6 were fitted using initial estimates for minimum density
(p1 = 197 kg m-3), the direct snow depth and density relation (p2 = 9 kg m-3 per meter), the
35

amplitude and spread of the function (p3 = 20, p5 = 0.015) and the timing of maximum density
due insolation (p4 = 180).
I optimized the RMSE between the modelled and the observed snow densities by
recording the function's minimum error and the p variables that minimized it. This scalar
function minimization was achieved through the open-source SciPy Python and shell scripts
compiled by Virtanen et al. (2020).
3.3.5. Model uncertainty assessment
All density models were tested with the validation dataset using the modelled densities
as response variables of observed snow depths. These modelled densities (Y’) were evaluated
against observed (n= 5,615) densities (Y). The root-mean-square error (Equation 3.7) was
calculated to estimate the errors’ average magnitude and assess the model performance.

ܴ‫ = ܧܵܯ‬ඨ

∑݊݅=1(ܻ − ܻ′)2
݊

(3.7)

The absolute mean difference between the modelled values (Y) and the observations
(MAE, Equation 3.8). MAE does not consider the errors’ direction, so it reflects the model’s
symmetric distribution.

݊

1
‫ = ܧܣܯ‬෍|ܻ − ܻ′|
݊
݅=1

36

(3.8)

MBE, alternatively, measures the average difference between the predicted and actual values,
considering the direction of the error (Equation 3.9). It identifies over- and underestimations,
systematic biases, and a skewed distribution.

݊

1
‫ = ܧܤܯ‬෍(ܻ − ܻ ′ )
݊

(3.9)

݅=1

Although the RMSE, MAE, and MBE are considered model performance metrics, the RMSE
penalizes errors more heavily, making RMSE more sensitive to outliers (Montgomery et al.,
2012; Neter et al., 1988). Thus, the RMSE was chosen as a conservative estimate of the density
model uncertainty.

3.4.

Basin-wide SWE estimations and uncertainty propagation

The best-performing snow density models were used to estimate snow density for each
survey date. With estimates of snow density for each pixel, SWE was calculated following
Equation 3.2. The average SWE for (a) non-glacierized and (b) glacierized portions of the basin
(taking each glacier individually) was then calculated, and the uncertainty in each mean was
calculated by propagating the individual errors in quadrature.
At a given point, SWE is a function of snow depth and density (Equation 3.10), with
independent random errors for snow depth (δsnow depth) and snow density (δρsnow). The error for
a single estimate of SWE was considered a product of independent and random uncertainties,
as in Taylor (1997; Equation 3.11)

37

(3.10)

ߩ‫ݓ݋݊ݏ‬
൰ ∗ ܵ‫ܦ‬
ܹܵ‫ = ܧ‬൬
ߩ‫ݎ݁ݐܽݓ‬

ߜ‫ݖ‬
ߜ‫ ݔ‬2
ߜ‫ ݕ‬2
= ඨ൬ ൰ + ൬ ൰
‫ݖ‬
‫ݔ‬
‫ݕ‬

(3.11)

Therefore, the relative uncertainty (ߜௌௐா ⁄ܹܵ‫ ) ܧ‬in the SWE estimations was estimated as
(Equations 3.12 and 3.13).
2

ߜ‫ݐ݌݁݀ ݓ݋݊ݏ‬ℎ
ߜܹܵ‫ܧ‬
ߜߩ‫ ݓ݋݊ݏ‬2
൰ +ቆ
= ඨ൬
ቇ
ܹܵ‫ܧ‬
ߩ‫ݓ݋݊ݏ‬
‫ݐ݌݁݀ ݓ݋݊ݏ‬ℎ

(3.12)

∴

2

ߜ‫ݐ݌݁݀ ݓ݋݊ݏ‬ℎ
ߜߩ‫ ݓ݋݊ݏ‬2
ඨ
൬
൰ +ቆ
ߜܹܵ‫= ݈݁݀݋݉ܧ‬
ቇ ∗ ܹܵ‫݈݁݀݋݉ܧ‬
ߩ‫ݓ݋݊ݏ‬
‫ݐ݌݁݀ ݓ݋݊ݏ‬ℎ

(3.13)

The standard deviation found in the 24 pseudo ground control points (Appendix A) served
as an estimate of the snow depth uncertainty. The RMSE for modelled snow density was used
as the snow density uncertainty.
തതതതതതത ; Equation 3.14) can be
For non-glacierized portions of the basin, the average SWE (ܹܵ‫ܧ‬
calculated by summing all SWE observations (SWE i >= 0) from non-glacierized pixels and
dividing by the number of pixels (N):

തതതതതതത
ܹܵ‫= ܧ‬

∑ ܹܵ‫݅ܧ‬
ܰ

38

(3.14)

Following the principles of error propagation, a conservative estimate of the uncertainty in
mean SWE (ߪௌௐா
തതതതതതത ; Equation 3.15) is then calculated as the square root of the sum of the squares
of the individual error terms (ߪௌௐா_௜ ) divided by the sample size (N)

∑(ߪܹܵ‫) ݅_ܧ‬2
ߪܹܵ‫ܧ‬
തതതതതതത = ඨ
ܰ

(3.15)

As described earlier, vertical changes due to glacier dynamics are zero when integrated
തതതതതതത௜ ) was averaged by its area
over a glacier. Therefore, each glacier’s overlying SWE (ܹܵ‫ܧ‬
തതതതതതത‫ ; ݁ݓ‬Equation
(‫ܽ݁ݎܣ‬௜ ; RGI Consortium, 2017). Then, the basin-wide average is weighted (ܹܵ‫ܧ‬

3.16).

തതതതതതത
ܹܵ‫ܧ‬௪௘ =

തതതതതതത௜ )
∑௡௜ୀଵ(‫ܽ݁ݎܣ‬௜ ∗ ܹܵ‫ܧ‬
ܶ‫ܽ݁ݎܽ ݈ܽݐ݋‬

(3.16)

The uncertainty of the weighted average (ߪௌௐா
തതതതതതതೢ೐ ; Equation 3.17) can be calculated by
interpreting the watershed glaciers as repeat observations (n).
To summarize, the methods used in this analysis include:
-

Snow depth retrievals from co-registered and bias-corrected lidar DEMs

-

Testing and development of empirical snow density models

-

Calculation of basin-wide SWE and its uncertainty from snow depths and modelled
snow density

Results are presented in the following chapter.

39

4.

Results
Lidar-derived, vertically corrected, and co-registered DEMs from the LaJoie basin,

Gold Bridge, British Columbia, were subtracted from the bare earth DEM to produce
georeferenced snow depth grids. These gridded snow depth values were analyzed and
compared to observed snow observations from snow courses and snow pillows.
The snow depth grids were combined with snow density models that yielded the lowest
RMSE to produce basin-wide snow water equivalent (SWE) estimations. As described in
Chapter 3, the SWE uncertainty was propagated from the independent errors in the lidarsurveying technique, identified through pseudo ground control points and RMSE from each
density model.

4.1.

Lidar data summaries

On average, the distribution of lidar-derived snow depths is positively skewed and
ranges from 0 to 5 m (Figures 4.1 and 4.2).

40

Figure 4.1: Lidar-derived snow depth over ice-free terrain, colour-coded by acquisition date.

41

Figure 4.2: Vertically corrected snow depth estimations. Colours over glaciers represent average
snow depths.

42

Mean snow depths in glacier-free terrain ranged between 1.2 and 1.7 m. (Figure 4.3;
table).

Figure 4.3: Mean snow depth by survey, colour coded by survey (Figure 4.1). Hatched bars: mean
snow depth weighted by glacier area (RGI Consortium, 2017).

4.1.1. Snow depth over glaciers
As expected, snow depth over glaciers was deeper than ice-free terrain and ranged from 2.2 to
5.4 m (Figure 4.3). One hundred and four glacier polygons in the watershed (RGI 6.0; RGI
Consortium, 2017) were computed to average snow depth weighted by glacier size (Table 4.1).
Table 4.1: Minimum, maximum, median, and weighted average depth over glaciers.

Survey date Min. depth

Max. depth Mdn. depth

Weighted
mean depth

2017-05-07

1.799

6.462

4.952

5.356

2018-02-19

1.123

4.446

3.092

3.333

2018-04-22

1.066

5.895

4.065

4.325

2019-03-18

1.149

4.412

2.413

2.498

2019-05-04

1.053

4.827

2.503

2.666

2020-03-09

0.923

4.368

2.311

2.213

2020-05-10

0.708

4.393

2.301

2.198

43

During the study, the median snow depth over glaciers ranged from 5.0 to 2.3 m., while
weighted averages ranged from 5.4 to 2.2 m. Snow depth distribution also changed from a
negative to a positive skewness (Figure 4.4). The May 2017 survey reflects 1 to 3 m—deeper
glacier snowpacks than any other sample, suggesting a high accumulation winter.

Figure 4.4: Glacier polygons by their averaged snow depth; colour-coded by survey (Figure 4.1).

44

4.1.2. Vertical snow accumulation gradients
Generally, snow depth increases with elevation for all surveys (Figure 4.5). On icefree terrain, median snow depth increases as elevation rises from the lowest elevations to
approximately 1800-1900 m. After this, the median snow depth remains consistent for a few
hundred meters and decreases beyond 2300 meters.

Figure 4.5: Vertical gradients of the observed snow depths at the LaJoie basin. 100-m box plots colourcoded by survey (Figure 4.1). The whisker plots show the median (black line), the first and third
quartiles (box limits), and the minimum and maximum values observed by elevation range.

45

The elevations above 2100 meters above sea level experienced the highest snow depth values,
peaking at 8 m (‘2015-05-07’), 7 m (‘2018-04-22’) and close to 6 m (every other survey).
4.1.3. Orientation of snow depth
Lidar-derived snow depth showed a preference for topographic orientation. The NEfacing slopes (NE, 0 to 90°) showed the deepest snowpacks, nearly half a meter deeper than
SW (180° to 270°) facing slopes, i.e., the shallowest snowpacks (Figure 4.6). In contrast, The
NW (270 to 0°) and SE (90° to 180°) facing slopes showed similar snow depth.

Figure 4.6: Azimuthal distribution of snow depths in non-glacial portions of the LaJoie Basin. Radial
axis: 0.5-meter increments. Dashed lines: mean snow depth by aspect, colour-coded by survey (Figure
4.1).

46

4.1.4. Uncertainty assessment
The pseudo-ground control points (GCPs) were chosen in snow-free pixels, allowing
any vertical bias in the lidar surveys to be identified and corrected. A total of 24 points were
collected from stable terrain (i.e., snow-free terrain between subsequent surveys).
The median elevation difference in snow-free pixels was isolated as the net elevation
bias by survey (in meters; Table 4.2) accounting for the number of measurements (N)
estimated, the standard deviation, and errors (Table 4.2). These standard errors served as a
conservative estimate of the uncertainty in the lidar-derived measurements.

Table 4.2: Bias, standard deviation, and standard error identified through 24 pseudo-GCP.

Survey date

Elevation bias (m.)
Median

Std.

Std. Error

* 2017-09-22

-1.63

0.56

0.11

2017-05-07

0.77

0.47

0.10

2018-02-19

0.12

0.18

0.04

2018-04-22

-0.10

0.35

0.07

2019-03-18

-0.02

0.16

0.03

2019-05-04

-0.05

0.18

0.04

2020-03-09

-0.03

0.20

0.04

2020-05-10

0.04

0.22

0.05

47

4.2.

Snow density modelling

Basin-wide lidar-derived snow depths require a snow density estimate to yield SWE.
Three approaches were used to model the snow density from the snow depth empirically. In
total, three snow density models from existing studies, six models calibrated with a snow
course dataset (n=5,615), and four models calibrated with snow pillow observations
(n=142,451) were tested against independent snow course data (n=5,615) to determine the best
empirical approach.
The RMSE between the modelled and the observed values served as an indicator of the
model’s uncertainty. The mean absolute (MAE) and the mean bias (MBE) errors were also
used to describe each model’s performance. Furthermore, the data and subsequent analysis
were limited to the Coast Mountains, the Lajoie Basin physiographic region.
4.2.1. Previous empirical snow density models
Empirical models found in the literature were tested with the validation dataset
(n=5,615) to evaluate the model performance (Figure 4.7). The model coefficients come from
each source; therefore, the modelled densities reflect their local calibrations. Model residuals
were examined to analyze model behaviour (Figure 4.8) further.

48

Figure 4.7: Observed snow densities compared to existing empirical snow density models. The line
represents a line of equality (1:1 line).

49

Figure 4.8: Observed snow densities compared to model residuals from existing empirical snow
density models.

The Hill et al. (2019) approach ranged snow densities between 250 and 400 kg m-3.
This model has the lowest RMSE (68.18 kg m-3) for any of the published models, but it also
produced underestimations below 300 kg m-3 and overestimations above 350 kg m-3.
Jonas et al. (2009) produced snow densities between 200 and 500 kg m-3 with an RMSE
of 84.58 kg m-3. This model showed the group's lowest MAE (52.52 kg m-3) and MBE (2.08
kg m-3). However, this model exhibited data clustering from monthly and elevation-based
parameters.
The Pistocchi et al. (2016) modelled densities ranged between 300 and 500 kg m-3, the
narrowest range in these existing snow models. This model also achieved the group’s highest
RMSE (97.77 kg m-3), mean absolute (81.29 kg m-3), and bias (-70.91 kg m-3) errors. Since the

50

model depends on the year’s day, it produces quasi-uniform predictions for a given day or
month (Figure 4.7).
4.2.2. Linear regression models from snow pillow data
The second set of models fits linear regression models to daily pillow data (Figures 4.9
and 4.10).

Figure 4.9: Observed snow densities compared to the snow pillow-derived linear regression model.
The line represents a line of equality (1:1 line).

51

Figure 4.10: Observed snow densities compared to model residuals from the snow pillow derived linear
regressions.

The pillow-derived linear regressions yielded densities that ranged between 200 and
650 kg m-3. The robust regressions outperformed the ordinary regressions, with a smaller
RMSE of 75.2 kg m-3 compared to 77.74 kg m-3. The robust regressions also achieved a lower
MAE of 56.96 kg m-3, whereas the ordinary model scored 59.49 kg m-3.
Both models displayed negatively skewed estimates with MBE values of -35.72 kg m-3 (robust)
and -38.07 kg m-3 (ordinary).
Once the pillow-derived linear regressions were restricted to the Coast Mountains, the
estimated densities ranged from 250 to 550 kg m-3. The robust model exhibited a lower RMSE
of 61.24 kg m-3 than the ordinary model at 63.25 kg m-3. However, both models demonstrated
similar MAE values of 44.90 kg m-3 (robust) and 47.21 kg m-3 (ordinary), and both models
displayed similar mean bias errors of -11.54 kg m-3 (robust) and -11.29 kg m-3 (ordinary).
52

4.2.3. Snow course-derived linear regressions
The modelled densities calibrated with snow course observations (Figure 4.11) produce
lower degrees of heteroscedasticity and RMSE values than those calibrated with snow pillow
data (Figure 4.12). The regression models using snow course observations achieved the
smallest evaluation metrics in all density models (Figure 4.15).

Figure 4.11: Observed snow densities compared to the snow course derived linear regression model.
The line represents a line of equality (1:1 line).

53

Figure 4.12: Observed snow densities compared to the model residuals from the snow course derived
linear regression model.

The snow course-derived linear regression modelled densities ranged between 180 and
600 kg m-3, the largest spread observed in this study. The robust regressions yielded marginally
higher RMSE (47.61 kg m-3) than the ordinary model (47.23 kg m-3). The robust regressions
also scored higher MBE (2.48 kg m-3) than the ordinary regressions (0.98 kg m-3), but both
showed comparable MAE (34.93 and 35.06 kg m-3; ordinary and robust, respectively).
Restricting the snow course data to the Coast Mountains returned a range of 200 to 550 kg m3

. The coastal snow course robust regressions model scored higher RMSE (36.79 and 36.31 kg

m-3; ordinary and robust, respectively) than province-wide analysis but lower MAE and MBE
(Table 4.3). The Coast Mountains’ snow course derived linear regressions scored this study's
lowest degree of heteroscedasticity.

54

4.2.4. Multi-parameter, non-linear equation
Five parameters in a non-linear equation were fitted to the validation dataset to
minimize the modelling error (Figures 4.13 and 4.14). This group achieved the second-lowest
evaluation metrics in this study and showed the second-fastest computing times and simplicity.

Figure 4.13: Observed snow densities compared to the multi-parameter, non-linear equation model.
The line represents a line of equality (1:1 line).

Figure 4.14: Observed snow densities compared to the model residuals from the multi-parameter,
non-linear equation.

55

After the error minimization, the minimum density (p1) achieved was 211.25 kg m-3,
the relation between snow depth and density (p2) was 48.48 kg m-3 per change in depth, and
the amplitude and spread of the function (p3, p5) were 325.27 and -0.02. The timing of
maximum density (p4) became the 203rd ordinary day (July 21st or 22nd).
The non-linear equation model achieved ranges between 211 and 650 kg m-3 and an RMSE of
51.93 kg m-3. The MBE (2.43 kg m-3) reflects positively skewed estimations.
Subsetting the data to the Coast Mountains replaced the minimum density (p1) with
299.41 kg m-3. The direct snow depth and density relation (p2) changed to 20.3, the amplitude
and spread of the function (p3, p5) became 255.5 and -0.02, and the timing of maximum density
(p4) became the 188th ordinary day (July 6th or 7th)
When restricted to data from the Coast Mountains, the modelled density values obtained
through the multi-parameter, non-linear equation ranged between 299 to 600 kg m-3. The model
achieved an RMSE of 52.79 kg m-3, more prominent than the province-wide analysis.

4.2.5. Model evaluation and selection
The density model performance was used to select those models to estimate SWE
(Figure 4.15; Table 4.3). The selection of the best snow density models hinged on evaluating
the model performance against snow course data through achieved RMSE, MAE, and MBE.

56

Table 4.3: Snow density model performance (RMSE, MAE and MBE are expressed in kg m-3).

Group

Literature
Equations

Snow pillow linear
regressions

Snow course linear
regressions

Multi-parameter,
non-linear equation

Model

RMSE

MAE

MBE

Hill et al. (2019).

68.18

53.05

7.7

Jonas et al. (2009)

84.58

52.52

2.08

Pistocchi et al. (2016)

97.77

81.29

-70.91

Robust - All samples

75.2

56.96

-35.72

Ordinary - All samples

77.74

59.49

-38.07

Robust - Coastal samples

61.24

44.9

-11.54

Ordinary - Coastal samples

63.25

47.21

-11.29

Robust - All samples

47.61

34.93

2.48

Ordinary - All samples

47.23

35.06

0.94

Robust - Coastal samples

36.79

27.73

1.18

Ordinary - Coastal samples

36.31

28.26

0.49

All samples

51.93

39.43

2.43

Coastal samples

52.79

40.82

-11.54

Figure 4.15: RMSE between the observed and the modelled snow densities. Samples restricted to the
Coast Mountains are denoted with an asterisk.

57

The snow course-derived linear regressions achieved the lowest performance metrics
but required at least eight field observations per regression. In contrast, existing empirical snow
density models required no field data. However, they produced overestimations (Hill et al.,
2019) or data clustering (Jonas et al., 2009), scoring the highest bias in this study.
Therefore, a) the multi-parameter, non-linear equation and b) the robust least square
regressions were chosen as the best density estimations. These models are calibrated with
different data sources, ensuring the generalization of this study.
The multi-parameter, non-linear equation achieved this study's second-to-lowest
RMSE, mean absolute and bias errors. The equation is based on stable snow densification
patterns in the province (Chapter 3.3.1), so no further field data are required.
The snow pillow-derived linear regressions require only automated data, which reduces
exposure to winter conditions and enhances temporal resolution and scalability. This model
achieved lower performance metrics than the existing empirical snow density models.
None of the calibrated density models showed significant metric improvement when
the datasets were reduced to Coast Mountains (Table 4.4). Any observed differences are likely
due to random variability rather than a real improvement in the modelled values.
Table 4.4: Comparison of RMSE values for different sample sets, expressed in kg m-3
Province-wide
analysis

Only coastal
samples

All but coastal
samples

Pillow-derived robust regressions

75.2

61.24

58.79

Pillow-derived ordinary regressions

77.74

63.25

60.29

Course-derived robust regressions

47.61

44.9

30.91

Course-derived ordinary regressions

47.23

36.79

31.11

Course-derived multi-parameter, non-linear eq.

51.93

52.79

37.66

Model

58

4.3.

Snow water equivalent estimations

The lidar-derived snow depth grids were paired with the best-performing snow density
models to produce basin-wide snow water equivalent (SWE) estimates. Each snow depth cell
served as the input value for each density model. The calculated SWE was stored in a
georeferenced grid with the same spatial extent and resolution as the snow depth grids. The
best density models were chosen based on their applicability and performance metrics against
snow course data.

4.3.1. Basin-wide summary
Basin-wide SWE estimations use the lidar-derived snow depths and density estimates
from the pillow-derived robust linear regressions (Figure 4.16) and the multi-parameter, nonlinear snow density models (Figure 4.17). These SWE estimates also include glacierized terrain
and the rest of the basin.
Unsurprisingly, SWE is highest (Figures 4.16 and 4.17) during acquisitions when snow
is most profound, and estimated snow density is high (e.g., acquisition dates ‘2018-04-22’,
‘2019-05-04’, and ‘20-05-10’).

59

Figure 4.16: Basin-wide SWE obtained with the multi-parameter, non-linear equation density model.
Colours over glaciers represent the propagated SWE produced with average snow depths.

60

Figure 4.17: Basin-wide SWE from with the pillow-derived robust linear regression density model.
Colours over glaciers represent the propagated SWE produced with average snow depths.

61

The total SWE estimates were partitioned into glacierized and non-glacierized terrain
(Tables 4.5 and 4.6; Figures 4.18 and 4.19) and integrated into a basin-wide estimate (Table
4.7).
Table 4.5: Non-glaciated terrain average snow depth (m.) and SWE (m.w.e) estimates, *Root mean squared
propagated uncertainties (± m.w.e.)
Snow water equivalence (m.w.e.)
Survey date Snow depth (m)
Non-linear equation*
Robust linear regressions*
2017-05-07

2.123

1.061 ± 0.143

0.924 ± 0.197

2018-02-19

1.553

0.555 ± 0.1

0.489 ± 0.142

2018-04-22

1.887

0.858 ± 0.125

0.773 ± 0.175

2019-03-18

1.516

0.572 ± 0.096

0.546 ± 0.137

2019-05-04

1.465

0.697 ± 0.098

0.709 ± 0.139

2020-03-09

1.579

0.587 ± 0.102

0.552 ± 0.145

2020-05-10

1.402

0.697 ± 0.099

0.509 ± 0.135

Table 4.6: Average snow depth and SWE estimates weighted by glacier size. **Standard deviation
with Bessel’s correction (± m.w.e.).
Survey date Snow depth (m)

Snow water equivalence (m.w.e.)
Non-linear equation**
Robust linear regressions**

2017-05-07

5.356

2.928 ± 0.517

1.920 ± 0.175

2018-02-19

3.333

1.266 ± 0.313

1.130 ± 0.287

2018-04-22

4.325

2.103 ± 0.488

1.813 ± 0.384

2019-03-18

2.498

0.962 ± 0.332

0.944 ± 0.376

2019-05-04

2.666

1.293 ± 0.457

1.280 ± 0.385

2020-03-09

2.213

0.818 ± 0.333

0.772 ± 0.289

2020-05-10

2.198

1.086 ± 0.462

0.794 ± 0.449

Table 4.7: Integrated basin-wide SWE averages, i.e., glacierized and non-glacierized SWE estimates
aggregated by area.
Survey date
2017-05-07
2018-02-19

Snow water equivalence (m.w.e.)
Non-linear equation
Robust linear regressions
1.416
1.113
0.690

0.611

62

2018-04-22

1.095

0.971

2019-03-18

0.646

0.622

2019-05-04

0.810

0.818

2020-03-09

0.631

0.594

2020-05-10

0.771

0.563

Figure 4.18: Mean snow water equivalent (multi-parameter, non-linear equation). Hatched bars: mean
water equivalent weighted by glacier size. Colour-coded by survey (Figure 4.1).

Figure 4.19: Mean snow water equivalent (snow pillow derived robust linear regressions). Hatched
bars: mean water equivalent over glaciers. Colour-coded by survey (Figure 4.1).

While comparable, the multi-parameter non-linear model yielded higher SWE
estimates during the 2017-05-07 and 2018-02-19 surveys than the linear regression model.
63

This can be attributed to the snow course dataset, wherein 18 of 11,230 snow course entries
represent snowpacks deeper than 5 m, with a timespan between 58 and 124 days of the year
(two months; Table 4.8). This class imbalance impacts the modelled SWE in the deepest
snowpacks.
Table 4.8: Temporal distribution of the snow course observations by depth
Snow
depth (m.)

Obs.
Count

min DOY

mdn DOY

max DOY

0
1
2
3
4
5
6

5589
3769
1344
407
102
18
1

1
1
1
1
29
58
97

63
86
86
89
91
95
97

366
366
366
364
152
124
97

In the non-glacierized terrain, the estimations’ interquartile range from 0 to 1.5 m.w.e.
and show a vertical gradient comparable to the snow depth values (Figure 4.20).

64

Figure 4.20: Vertical gradient observed in the SWE estimates (ice-free terrain) achieved with the
pillow-derived robust linear regressions (navy blue) and the optimized multi-parameter, non-linear
equation (gold bars).

The non-linear equation model achieved smaller values in the early season than the pillow
robust regressions. Inversely, the late-season pillow-derived regressions showed lower values
than the non-linear model. The April 2018 survey (‘2018-04-22’) shows the best agreement
between both models.

65

4.3.2. Point-based validation
I used two snow pillows in the LaJoie basin (Figure 4.21) to validate the SWE estimates
independently. Hourly SWE observations were excluded from the snow pillow calibration
dataset as they lacked snow depth measurements.

Figure 4.21: Snow pillows at the LaJoie Basin. The glacier polygons display the Randolph Glacier
Inventory (RGI Consortium, 2017), also used in this study.

The validation consisted of extracting snow depth from a 15-m radius at the exact
coordinates of the snow pillow sites (Figures 4.22 and 4.23) and applying the snow density
models to obtain estimates of SWE for each station. These SWE estimates were validated
against 14 pillow observations (seven dates observed in two snow pillows).

66

In half of the comparisons against the multi-parameter non-linear model, the pillow
observations are within the interquartile range, and five overlap the modelled SWE median
counterpart (Figure 4.22). Pillow measurements beyond the interquartile range indicate
disagreements smaller than 0.14 m.w.e. (Table 4.9).
Table 4.9: Snow pillow observations against SWE estimates (m.w.e.) achieved with the snow course
derived multi-parameter, non-linear snow density model.

Survey date

1C12P obs. (m.w.e.) SWE estimate (m.w.e.)
Mdn.
Std.
Min. Mdn.
Max.

2017-05-07

0.713

0.002

0.442

0.670

0.998

13.230

2018-02-19

0.690

0.000

0.369

0.509

0.681

15.616

2018-04-22

0.812

0.003

0.661

0.865

1.165

12.508

2019-03-18

0.682

0.001

0.384

0.545

0.764

14.591

2019-05-04

0.704

0.004

0.412

0.633

0.942

11.842

2020-03-09

0.537

0.000

0.400

0.555

0.719

14.972

2020-05-10

0.506

0.006

0.356

0.495

0.767

12.134

Survey date

1C38P obs. (m.w.e.)
SWE estimate (m.w.e.)
Mdn.
Std.
Min. Mdn.
Max.

2017-05-07

0.858

0.004

0.861

0.999

1.266

11.948

2018-02-19

0.604

0.001

0.182

0.572

0.700

15.421

2018-04-22

0.699

0.001

0.674

0.843

0.982

12.558

2019-03-18

0.686

0.001

0.517

0.677

0.812

14.285

2019-05-04

0.714

0.002

0.433

0.730

0.941

11.660

2020-03-09

0.598

0.000

0.535

0.681

0.824

14.642

2020-05-10

0.557

0.004

0.399

0.580

0.769

11.811

67

Relative unc.
Mdn. (%)

Relative unc.
Mdn. (%)

Figure 4.22: SWE observations compared to estimates from the multi-parameter, non-linear model
colour-coded by survey (Figure 4.1). The line represents a line of equality (1:1 line).

Alternatively, all pillow observations fell within the SWE modelled with the snow
pillow-derived robust linear regressions, and five fell inside the interquartile ranges. The rest
showed discrepancies smaller than 0.15 m.w.e from the interquartile range (Table 4.10).

Table 4.10: Snow pillow observations compared with SWE estimates (m.w.e) achieved with the snow
pillow-derived robust linear regressions model.

Survey date

SWE estimate
1C1P2 obs. (m.w.e.)
Mdn.
Std.
Min. Mdn.
Max.

2017-05-07

0.713

0.002

0.488

0.703

0.975

16.924

2018-02-19

0.690

0.000

0.419

0.554

0.704

20.505

2018-04-22

0.812

0.003

0.657

0.833

1.074

18.279

2019-03-18

0.682

0.001

0.337

0.494

0.722

22.993

2019-05-04

0.704

0.004

0.468

0.689

0.968

15.466

2020-03-09

0.537

0.000

0.390

0.535

0.684

22.184

2020-05-10

0.506

0.006

0.408

0.563

0.856

14.916

68

Relative unc.
Mdn. (%)

Survey date

1C38P obs. (m.w.e.)
SWE estimate
Mdn.
Std.
Min. Mdn.
Max.

Relative unc.
Mdn. (%)

2017-05-07

0.858

0.004

0.866

0.975

1.167

16.831

2018-02-19

0.604

0.001

0.220

0.611

0.720

20.655

2018-04-22

0.699

0.001

0.668

0.814

0.929

18.275

2019-03-18

0.686

0.001

0.467

0.630

0.773

21.997

2019-05-04

0.714

0.002

0.490

0.780

0.967

15.536

2020-03-09

0.598

0.000

0.516

0.649

0.778

21.957

2020-05-10

0.557

0.004

0.456

0.656

0.858

14.698

Figure 4.23: SWE observations compared to the snow pillow-derived robust linear regressions SWE
estimates, colour-coded by survey (Figure 4.1). The line represents a line of equality (1:1 line).

4.4.

Uncertainty propagation

The modelled SWE uncertainty was calculated for every raster, considering the independent
errors in a) snow depth (from the lidar pseudo ground control points’ standard deviation) and
b) modelled snow density (from the RMSE between the achieved modelled densities and the

69

snow course observations; Figures 4.24 and 4.26). Relative SWE estimates consider the
uncertainty magnitude against the estimated SWE values (Figures 4.25 and 4.27).

70

Figure 4.24: Snow course derived multi-parameter, non-linear model – SWE uncertainty. Colours
over glaciers represent the propagated uncertainty produced with average snow depths.

71

Figure 4.25: Snow course derived multi-parameter, non-linear model – relative SWE uncertainty.
Colours over glaciers represent the propagated uncertainty produced with average snow depths.

72

Figure 4.26: Snow pillow derived robust linear regression model – SWE uncertainty. Colours over
glaciers represent the propagated uncertainty produced with average snow depths.

73

Figure 4.27: Snow pillow derived robust linear regression model – relative SWE uncertainty. Colours
over glaciers represent the propagated uncertainty produced with average snow depths.

74

Uncertainties in gridded SWE values ranged from ±0.010 to ±0.424 m.w.e. for the nonlinear snow density model, while the SWE uncertainties from snow pillow-derived regressions
ranged from ±0.011 to ±0.599 m.w.e. respectively (Tables 4.11 and 4.12; Figures 4.28 and
4.29). SWE uncertainties display positive skewness (Figure 4.28).
Table 4.11: Modelled SWE uncertainty ranges achieved with the snow course derived multi-parameter,
non-linear equation.
Survey Date
2017-05-07
2018-02-19
2018-04-22
2019-03-18
2019-05-04
2020-03-09
2020-05-10

Min. SWE
Mdn. SWE
Max. SWE
uncertainty (m.w.e.) uncertainty (m.w.e.) uncertainty (m.w.e.)
0.054
0.011
0.027
0.008
0.020
0.013
0.012

0.151
0.103
0.126
0.101
0.096
0.104
0.085

0.599
0.599
0.599
0.596
0.599
0.598
0.514

Median relative
uncertainty (%)
12.782
16.003
13.335
14.859
12.401
15.333
13.008

Figure 4.28: Histogram of SWE uncertainty in non-glacierized terrain calculated with the snow coursederived multi-parameter non-linear equation, colour-coded by survey (Figure 4.1).

75

Table 4.12: Modelled SWE uncertainty ranges achieved with the snow pillow derived linear regression.

Survey Date

Min. SWE
uncertainty (m.w.e.)

Mdn. SWE
uncertainty (m.w.e.)

Max. SWE
uncertainty (m.w.e.)

Median relative
uncertainty (%)

2017-05-07

0.044

0.111

0.424

17.725

2018-02-19

0.012

0.073

0.420

20.776

2018-04-22

0.028

0.092

0.422

18.696

2019-03-18

0.010

0.071

0.420

23.840

2019-05-04

0.017

0.068

0.420

15.910

2020-03-09

0.013

0.074

0.420

22.431

2020-05-10

0.022

0.064

0.421

15.319

Figure 4.29: Histogram of SWE uncertainty in non-glacierized terrain calculated with the snow pillowderived linear regressions, colour-coded by survey (Figure 4.1).

76

5.

Discussion
5.1.

Lidar-derived snow depth

5.1.1. General observations
The lidar observations demonstrated a strong relation between topography and snow
depth. Deeper snowpacks were consistently found at higher elevations of 1800-1900 m, where
snow depth remains constant for a few hundred meters. Higher elevations are more prone to
erosion (Déry et al., 2010), thus above 2300 m., the snow depth decreases, possibly due to
redistribution due to wind and snow avalanches (Grünewald et al., 2014; Vionnet et al., 2021).
The measured snow depth also reflects the underlying terrain orientation. The NE-facing slopes
(NE, 0 to 90°) host the deepest snowpacks, nearly half a meter deeper than the SW (180° to
270°) counterpart, i.e., the shallowest snowpacks (Figure 4.6). While this can denote the
leeward and windward slopes, this contrast can also be attributed to incoming radiation
patterns, wherein the NE-facing sides receive less solar exposure, which helps preserve the
snowpacks.
Wind circulation patterns also influence variation in snow depth in the Coast
Mountains. The Pacific High and the Aleutian Low systems peak during winter, resulting in a
dominant southwesterly circulation (Klock et al., 2001) evidenced in low-slope areas with
wind-swept snow.
The glacierized terrain contains deeper snowpacks than the rest of the basin (Tables 4.5
and 4.6). The topographic shelter, negligible ground heat flux from glacier ice, and skin
temperatures of the ice at the melting point all contribute to deep snowpacks. Despite the
increased snow retention on glaciers, glacier shrinkage will limit the areas of deep snowpacks
77

in the Lajoie basin for the foreseeable future (Bevington & Menounos, 2022; Rounce et al.,
2023).
5.1.2. Comparative analysis of snow depth uncertainties
This study uses multiple lidar surveys for a given year (Table 3.2), unlike studies
relying on a single (Moreno et al., 2009) or yearly observations (Deems et al., 2013; Pflug &
Lundquist, 2020). These data allow for a more detailed understanding of the temporal
dynamics of snowpacks over topographically complex terrain.
Given the lack of ground observations, I estimated the vertical uncertainty in the lidarderived snow depths by calculating the variance of elevation over stable terrain, a methodology
distinct from relying on field data (López-Moreno & Nogués-Bravo, 2006; Painter et al., 2016).
A significant limitation of this approach is that the area of stable terrain required to complete
this analysis is much smaller than conventional methods used in error propagation for glacier
change (e.g., Pelto et al., 2017; Menounos et al., 2019).
The lidar-derived snow depth uncertainties ranged from ±0.03 to ±0.11 m. (Table 4.2),
with a sampling resolution of approximately one laser return per m-2. Increasing the density of
laser returns, up to 364 per m-2, significantly reduced uncertainties as low as ±0.0122 m,
particularly in open ground (Jacobs et al., 2021).
However, the largest source of uncertainty in SWE estimates lies in the applied snow
density models (Mital et al., 2022; Painter et al., 2016). Consequently, potential improvements
in lidar precision need to be weighed against the associated costs and practical implications.
Alternative methods can be used to estimate mountainous snow depth. Using
photogrammetric reconstructions, the high-resolution stereo pair imagery from the Pléiades
78

satellites can achieve 4 m ground sampling resolutions, with uncertainties from ±0.36 m (Shaw
et al., 2019) to ± 0.2 m. (Deschamps-Berger et al., 2020; Shaw et al., 2020) Meanwhile,
uncrewed aerial vehicles (UAV) have centimetre-scale uncertainties (±0.05 to ±0.11 m;
Deschamps-Berger et al., 2020). While these alternatives offer different acquisition flexibility,
the precision achieved by lidar should be duly acknowledged. Further research should weigh
precision over project requirements and terrain.

5.2.

Empirical snow density modelling

5.2.1. Major findings
The linear regressions calibrated with snow course data produced the lowest RMSE
(47.61 and 47.23 kg m-3) of this study, outperforming the regressions calibrated with snow
pillows (75.2 and 77.74 kg m-3, robust and ordinary, respectively). The multi-parameter, nonlinear density parameterization (also calibrated with snow course data) achieved the secondto-lowest RMSE of this study (51.93 kg m-3; Table 4.3),
A provincial subset of the calibration data was used to evaluate physiographic
sensitivity. Reducing the models to the Coast Mountains did not significantly improve the
statistical performance, supporting the model’s transferability. However, there were notable
differences in the parameter fitting of the multi-parameter, non-linear models. When
considering the Coast Mountains samples exclusively, the timing of maximum density shifted
from the 203rd ordinal day (July 21st to 22nd) to the 188th day (July 6th or 7th). The earlier

79

end of the snow season highlights different snowmelt patterns in the southern Coast Mountains
than elsewhere.

5.2.2. Comparative analysis of density modelling errors
The empirical density models developed here can be compared to physically-based
models. The snowmelt Snobal model used in California, US, achieved RMSE values ranging
from 52 to 79 kg m-3 (Smyth et al., 2019), and the SNOWPACK model achieved a bias between
22 to 31 kg m-3 over Antarctic snow and firn (Keenan et al., 2021). The models in this study
achieved RMSE values between 36 and 77 kg m-3, affirming their comparability to state-ofthe-art physically-based approaches.
Unlike the empirical models used in this study, a physically-based model requires the
start of the snow season, among other variables. Contrarily, Hill et al. (2019) considered the
beginning of the snow season (DOY) by October 1st, while Pistocchi (2016) chose November
1st instead. Generalizations regarding uniform start or end dates may serve practical purposes,
but the elevation, topography, and atmospheric patterns influence the initial snowfall
conditions. Analyzing snowpack timespans with an empiric perspective may further allow us
to explain densification patterns.

5.3.

Propagated basin-wide SWE and error

5.3.1. Major findings
Lidar-derived SWE estimates showed interannual variability but significant differences
between glacierized and non-glacierized terrain. Glaciers have significantly higher snow
80

retention than non-glacierized terrain (Tables 4.5 and 4.6) and need to be considered separately
in any calculation of basin-wide SWE.
Notably, the uncertainty in my modelled snow densities tends to be more pronounced
in thin, early-season snowpacks. Since the lidar samples were strategically planned to capture
near-maximum SWE in the watershed, extrapolating my findings to thinner, less dense
snowpacks may significantly induce errors. Similar studies (Helfricht et al., 2014; Painter et
al., 2016) reinforce the need to explore thinner, less dense snow cover scenarios.

5.3.2. Comparative analysis
This study’s non-linear model combined with lidar-derived snow depth resulted in a
median relative uncertainty of 12.9% of the SWE estimates (Table 4.10). Comparatively,
Guyennon et al. (2019) achieved a relative uncertainty of 15.6% using manual snow depth and
a non-linear empirical snow density model. Since the SWE variations primarily come from
depth fluctuations (Painter et al., 2016), the reduced relative uncertainty in this study is more
influenced by greater snow depth sample size and enhanced control on snow mass than the
modelled densities.
This study’s 5 m ground sampling resolution was sufficiently detailed, considering the 0.03 to
0.1 m uncertainty achieved and the 962.5 km2 sampled. At finer, 1 m resolutions, snow depth
uncertainties were achieved around 0.18 m (Pyrenees Mountains, Spain; Moreno et al., 2009)
and 0.15 m. (Tyrolean Alps, Austria; Schöber et al., 2014), However, time and operational
expenses constrained their studies to relatively small watersheds (38 and 36 km2).

81

Larger watersheds are attainable but demand an extraordinary level of effort. The
Tuolomne Watershed, USA, hosts the “most [airborne laser surveyed] region in the world”
(Pflug & Lundquist, 2020). With a surveyed area of 1625 km2 and seven years of observations,
3-meter snow depth rasters achieved uncertainties smaller than 0.8 m. These cases demonstrate
that larger areas and finer resolutions may impact the produced uncertainties, but not
significantly. However, these comparisons should be taken carefully as these studies compare
their surveys to in-situ observations.
This study's design balances resolution, area coverage, and associated uncertainties
over the basin-wide estimates. Further contributions to minimizing the uncertainty of snow
density are encouraged.

5.4.

Study limitations

Several notable limitations of this study exist. The lack of ground snow depth and
density observations prevents independent error analysis through time and space. Significant
uncertainties arise from 1) the temporal resolution achieved with lidar, 2) the complexity of
the glacierized terrain, and 3) the independent data used to calibrate the density models.
The temporal coverage in my lidar-derived snow depths captured conditions near
maximum snow depths, restricting the generalization of findings to other winter periods.
Increasing the sampling resolution would also help to characterize snowpack dynamics,
identify hiatus, and minimize errors in snow depth estimates. Meanwhile, transferring the
regression models primarily during the melt period, assuming consistent metamorphism rates
(e.g., compaction, hydraulic conductivity), is advisable.

82

Secondly, averaging glacial snow depths (or SWE) oversimplifies the complexity of
glacier accumulation patterns. Ice dynamics widely impact the winter mass balance gradients
observed in glaciers of western Canada (e.g., Mukherjee et al., 2023; World Glacier
Monitoring Service (WGMS), 2023), particularly in glaciers with high rates of mass flux and
wind erosion. Acknowledging these complexities through an ice dynamics model is necessary
to describe snow depth over glacierized terrain.
Additionally, the snow depth and accumulation rates beneath vegetation are often
subject to significant uncertainties (Zaremehrjardy et al., 2021). Recent research emphasized
the influence of canopy height and soil type on snow distribution, revealing distinct clusters or
zone groupings that impact measurement and modelling accuracy (Geissler et al., 2023).
Stratifying lidar and modelling data by representative land cover classes could enhance model
accuracy (Deems et al., 2013) and enable snow density and error maps by region.
Lastly, the data used to calibrate these empirical models are imperfect datasets. The
snow course data (Province of British Columbia, 2021a) has a coarse temporal resolution and
scarcely reflects the deepest snowpacks (i.e., more profound than 3 m). Snow pillow data,
instead, are readily available (Province of British Columbia, 2021b) but susceptible to data
outliers, as evidenced by the higher RMSE achieved with robust regressions against ordinary
regressions.
These data outliers can be identified on a pillow-scale basis. A snow column naturally
densifies over time, increasing daily. Some snow pillows, however, do not record a
progressively increasing snow density (Figure 5.2).

83

Figure 5.1: Densification pattern recorded by two snow pillows in one snow season. The snow pillow
entries denote densification inconsistencies in the bottom panel (Station 1F06P).

While noise or equipment malfunction should be considered, these densification artifacts may
be caused by pillows prone to developing ice lenses or being exposed to runoff, compromising
the column resistivity and the inferred snow depth

84

5.5.

Transferability of results

The uncertainty in lidar measurements depends upon operational factors, such as flight
conditions and sensor specifications. Nonetheless, the major findings of this study may assist
in creating representative study designs.
In large basins, the elevation and aspect dependence of snow depth can be addressed by
designing surveys that sample slope, aspect, and elevations.
A survey that also samples the most representative land cover classes (e.g., clear-cut versus
canopy), different hierarchy discharge pathways, and glacial and non-glacial terrain (if
applicable); may help us to further understand the hydraulic distribution within the catchment
area and isolate independent densification patterns.
In contrast, the snow density models exhibit should be transferrable across the province
and may be applied to other regions with representative calibration datasets. The laser altimetry
and empirical modelling approaches are interchangeable with other snow depth and density
measurement techniques.

85

6.

Conclusions
This thesis achieved three objectives: (1) acquire and process lidar snow depth data, (2)

calibrate and model empirical snow density, and (3) estimate basin-wide snow water equivalent
propagating independent uncertainties.
In the first part of this study, lidar-derived snow depths showed positively skewed
distributions, with most measurements falling within zero and five meters. Snow depths in
glacier-free terrain averaged between 1.4 and 2.1 m. In comparison, the glacierized landscape
weighted average snow depth between 2.2 and 5.4 m, emphasizing the importance of snow
retention on glaciers in basin-wide SWE calculations.
Snow depth increases with elevation up to 1900 m. Above 1900 m.a.s.l. snow depths decrease,
reflecting the importance of elevation (favourable lapse rates of precipitation and adiabatic
cooling), followed by snow redistribution by wind and gravity. Deeper snowpacks were found
on NE-facing slopes, nearly half a meter deeper than the shallowest SW-facing snowpacks.
Snow-free pixels were isolated to identify and correct vertical biases in the lidar-derived snow
depths. The median elevation difference in 24 stable terrain pixels served as the net elevation
bias, while the standard error denoted an uncertainty in the estimates from ±0.03 to ±0.11 m.
In the second part of this study, I evaluated three previously published snow density
models (Hill et al., 2019; Pistocchi, 2016; Sturm et al., 2010) and calibrated ten linear and nonlinear empirical models with provincial snow course and pillow data. The modelled density
values were tested against a snow course validation dataset (n=5,616) to evaluate the model
performance.

86

Provincial snow pillow (n=142,451) and snow course (n=5,616) observations were used to
calibrate ordinary and robust linear regressions of snow density versus depth. The robust
regressions with snow pillows achieved lower RMSE (74.12 and 76.97 kg m-3), MAE (55.76
and 58.64 kg m-3), and MBE (-34.97 and -36.97 kg m-3) compared to ordinary linear
regressions. Conversely, the snow course-derived robust and ordinal regressions demonstrated
marginal differences in RMSE (47.65 and 47.21 kg m-3), MAE (34.93 and 34.87 kg m-3) and
MBE (3.05 and 1.55 kg m-3), achieving the best performance metrics in this study.
I also calibrated a multi-parameter, non-linear snow density model to the validation dataset,
achieving RMSE (51.6 kg m-3), MAE (39.17 kg m-3), and MBE (3.29 kg m-3) values.
Observations from the Coast Mountains region were used as a subset to test the influence of
physiography on snow densification. Marginal differences between coastal and non-coastal
groups were more likely due to random variability than real estimate improvements.
The snow course-derived multi-parameter, non-linear model, and the snow pillow-derived
robust regressions provided the lowest errors in estimated snow density and were carried
forward for subsequent comparison. These models are also the most transferrable across the
province and beyond to other regions
In the third part of this study, basin-wide SWE estimates were produced from the lidarderived snow depth and the two best snow density models. Late winter surveys showed higher
SWE estimates than their respective early winter counterparts. The first survey produced the
highest average SWE (1.061 and 0.924 m.w.e.), and the last survey had the lowest (0.697 and
0.509 m.w.e.; non-linear and linear regression model respectively).

87

Independent data from two snow pillow observations were compared against lidar and
snow density model-derived SWE estimates from a 15-meter radius around each site, totalling
28 comparison points. The modelled SWE values fell within the SWE observations, from
which 12 also fell within the daily interquartile ranges.
Any differences between the pillow observations and the SWE estimates’ interquartile range
were smaller than 0.15 m.w.e. In 24 out of the 28 comparisons, the combined value of the
median SWE estimates and their uncertainties exceeded the pillow observations, fulfilling a
conservative approach and demonstrating this study’s comparability to remotely sensed data.
Finally, the SWE uncertainty was propagated from the independent errors in snow
depth (determined from repeat observations of elevation on snow-free pixels) and the RMSE
of the best snow density models. Both density models showed a higher contribution to basinwide uncertainty than the lidar-derived snow depth.
The complexities of glacier dynamics present challenges in calculating snow depth and water
equivalent on individual glacierized pixels. Without an underlying ice dynamics model,
averaging snow depth by glacier lets us approximate glacial snowpack inventories by weighing
averages (5.356 to 2.198 m.).
Nonetheless, these averages, based on bare earth models collected during a time of rapid
glacier change, may misrepresent true glacial snow depths and underscore the necessity for
considering glacier dynamics when assessing mountainous glacial snowpacks.
In non-glacierized terrain, median uncertainties in SWE from ±0.096 to ±0.143 m.w.e.
using the non-linear density model. Alternatively, the snow pillow derived regressions

88

achieved median SWE uncertainties between ±0.135 and ±0.197 m.w.e. Estimated
uncertainties showed positively skewed distributions.
Snow density models that linearly correlate station-based observations of SWE and snow
depth, and fitted non-linear empirical models that consider depth, density and day of the year
should be developed and tested in other regions. Further explorations of snow pillow datasets
may help to isolate a snowpack’s granular texture, effective permeability, residence time and
ablation rates. These factors collectively contribute to understanding the behaviour of a
snowpack over time.
In the overall backdrop of this study, continued loss of alpine glacier mass (Anderson, 1976;
Pelto, 2020) and progressively shallower snowpacks can be expected due to climate change.
This study contributes to understanding the snowpack dynamics in glacierized alpine regions,
and these efforts will improve our resource management amidst the uncertainty of evolving
water cycles and weather patterns.

89

7.

References

Anderson, E. A., & States., U. (1976). A point energy and mass balance model of a snow cover (pp.
xviii, 150 p.). U.S. Dept. of Commerce, National Oceanic and Atmospheric Administration,
National Weather Service, Office of Hydrology. file://catalog.hathitrust.org/Record/101642680
Bevington, A. R., & Menounos, B. (2022). Accelerated change in the glaciated environments of western
Canada revealed through trend analysis of optical satellite imagery. Remote Sensing of
Environment, 270, 112862. https://doi.org/10.1016/j.rse.2021.112862
Boon, S. (2009). Snow ablation energy balance in a dead forest stand. Hydrological Processes, 23(18),
2600–2610. https://doi.org/10.1002/hyp.7246
Brun, E., Vionnet, V., Boone, A., Decharme, B., Peings, Y., Valette, R., Karbou, F., & Morin, S. (2013).
Simulation of Northern Eurasian Local Snow Depth, Mass, and Density Using a Detailed
Snowpack Model and Meteorological Reanalyses. Journal of Hydrometeorology, 14(1), 203–219.
https://doi.org/10.1175/JHM-D-12-012.1
Csele, M. (2004). Fundamentals of light sources and lasers. In John Wiley & Sons, Inc.
https://doi.org/10.1097/00024382-200205000-00022
Deems, J. S., Painter, T. H., & Finnegan, D. C. (2013). Lidar measurement of snow depth: a review.
Journal of Glaciology, 59(215), 467–479. https://doi.org/10.3189/2013JoG12J154
Déry, S. J., Clifton, A., MacLeod, S., & Beedle, M. J. (2010). Blowing Snow Fluxes in the Cariboo
Mountains of British Columbia, Canada. Arctic, Antarctic, and Alpine Research, 42(2), 188–197.
https://doi.org/10.1657/1938-4246-42.2.188
Deschamps-Berger, C., Gascoin, S., Berthier, E., Deems, J., Gutmann, E., Dehecq, A., Shean, D., &
Dumont, M. (2020). Snow depth mapping from stereo satellite imagery in mountainous terrain:
evaluation using airborne laser-scanning data. The Cryosphere, 14(9), 2925–2940.
https://doi.org/10.5194/tc-14-2925-2020
Dingman, S. L. (2002). Physical Hydrology (2nd Editio). Prentice Hall.
Dozier, J., M. Melack, J., Marks, D., Elder, K., & Kattelmann, R. (1988). Snow deposition, melt, runoff,
and chemistry in a small alpine watershed, Emerald Lake Basin, Sequoia National Park. Final
report, 1 July 1984-31 March 1987.
Elder, K., Dozier, J., & Michaelsen, J. (1991). Snow accumulation and distribution in an Alpine

90

Watershed. Water Resources Research, 27(7), 1541–1552. https://doi.org/10.1029/91WR00506
Euskirchen, E. S., Goodstein, E. S., & Huntington, H. P. (2013). An estimated cost of lost climate
regulation services caused by thawing of the Arctic cryosphere. Ecological Applications, 23(8),
1869–1880. https://doi.org/10.1890/11-0858.1
Flanner, M. G., Shell, K. M., Barlage, M., Perovich, D. K., & Tschudi, M. A. (2011). Radiative forcing
and albedo feedback from the Northern Hemisphere cryosphere between 1979 and 2008. Nature
Geoscience, 4, 151. https://doi.org/10.1038/ngeo1062
Frei, A., Tedesco, M., Lee, S., Foster, J., Hall, D. K., Kelly, R., & Robinson, D. A. (2012). A review
of global satellite-derived snow products. Advances in Space Research, 50(8), 1007–1029.
https://doi.org/10.1016/j.asr.2011.12.021
Geissler, J., Rathmann, L., & Weiler, M. (2023). Combining Daily Sensor Observations and Spatial
LiDAR Data for Mapping Snow Water Equivalent in a Sub‐Alpine Forest. Water Resources
Research, 59(9). https://doi.org/10.1029/2023WR034460
Grünewald, T., Bühler, Y., & Lehning, M. (2014). Elevation dependency of mountain snow depth. The
Cryosphere, 8(6), 2381–2394. https://doi.org/10.5194/tc-8-2381-2014
Hedrick, A. R., Marks, D., Havens, S., Robertson, M., Johnson, M., Sandusky, M., Marshall, H.-P.,
Kormos, P. R., Bormann, K. J., & Painter, T. H. (2018). Direct Insertion of NASA Airborne Snow
Observatory-Derived Snow Depth Time Series Into the iSnobal Energy Balance Snow Model.
Water Resources Research, 54(10), 8045–8063. https://doi.org/10.1029/2018WR023190
Hedstrom, N. R., & Pomeroy, J. W. (1998). Measurements and modelling of snow interception in the
boreal

forest.

Hydrological

Processes,

12(10–11),

1611–1625.

https://doi.org/10.1002/(SICI)1099-1085(199808/09)12:10/11<1611::AID-HYP684>3.0.CO;2-4
Helfricht, K., Kuhn, M., Keuschnig, M., & Heilig, A. (2014). Lidar snow cover studies on glaciers in
the Ötztal Alps (Austria): comparison with snow depths calculated from GPR measurements. The
Cryosphere, 8(1), 41–57. https://doi.org/10.5194/tc-8-41-2014
Hill, D. F., Burakowski, E. A., Crumley, R. L., Keon, J., Hu, J. M., Arendt, A. A., Wikstrom Jones, K.,
& Wolken, G. J. (2019). Converting snow depth to snow water equivalent using climatological
variables. The Cryosphere, 13(7), 1767–1784. https://doi.org/10.5194/tc-13-1767-2019
Hirst, S. M. (1991). Impacts of the operation of existing hydroelectric developments on fishery
resources in British Columbia, Vol. II Inland Fisheries.

91

Huber, P. J. (1973). Robust Regression: Asymptotics, Conjectures and Monte Carlo. The Annals of
Statistics, 1(5), 799–821. http://www.jstor.org.prxy.lib.unbc.ca/stable/2958283
Huss, M., Bookhagen, B., Huggel, C., Jacobsen, D., Bradley, R. S., Clague, J. J., Vuille, M., Buytaert,
W., Cayan, D. R., Greenwood, G., Mark, B. G., Milner, A. M., Weingartner, R., & Winder, M.
(2017). Toward mountains without permanent snow and ice. Earth’s Future, 5(5), 418–435.
https://doi.org/10.1002/2016EF000514
Intergovernmental Panel on Climate Change (Ed.). (2014). Climate Change 2013 – The Physical
Science Basis. Cambridge University Press. https://doi.org/10.1017/CBO9781107415324
International Energy Agency. (2018). World Energy Balances 2018.
Jackson, M. (1993). 35th Annual Grammy Awards. CBS Broadcasting Inc.
Jacobs, J. M., Hunsaker, A. G., Sullivan, F. B., Palace, M., Burakowski, E. A., Herrick, C., & Cho, E.
(2021). Snow depth mapping with unpiloted aerial system lidar observations: a case study in
Durham,

New

Hampshire,

United

States.

The

Cryosphere,

15(3),

1485–1500.

https://doi.org/10.5194/tc-15-1485-2021
Jonas, T., Marty, C., & Magnusson, J. (2009). Estimating the snow water equivalent from snow depth
measurements

in

the

Swiss

Alps.

Journal

of

Hydrology,

378(1–2),

161–167.

https://doi.org/10.1016/j.jhydrol.2009.09.021
Keenan, E., Wever, N., Dattler, M., Lenaerts, J. T. M., Medley, B., Kuipers Munneke, P., & Reijmer,
C. (2021). Physics-based SNOWPACK model improves representation of near-surface Antarctic
snow and firn density. The Cryosphere, 15(2), 1065–1085. https://doi.org/10.5194/tc-15-10652021
Klock, R., Mullock, J., of Canada, M. S., & Canada, N. (2001). The Weather of British Columbia:
Graphic

Area

Forecast

31.

NAV

CANADA.

https://books.google.com.mx/books?id=iEvqlgEACAAJ
Libbrecht, K. G. (2005). The physics of snow crystals. Reports on Progress in Physics, 68(4), 855–
895. https://doi.org/10.1088/0034-4885/68/4/R03
López-Moreno, J. I., & Nogués-Bravo, D. (2006). Interpolating local snow depth data: an evaluation
of methods. Hydrological Processes, 20(10), 2217–2232. https://doi.org/10.1002/hyp.6199
Lutz, D. A., & Howarth, R. B. (2015). The price of snow: albedo valuation and a case study for forest

92

management. Environmental Research Letters, 10(6), 064013. https://doi.org/10.1088/17489326/10/6/064013
MacKenzie, W. H., & Meidinger, D. V. (2018). The Biogeoclimatic Ecosystem Classification
Approach: an ecological framework for vegetation classification. Phytocoenologia, 48(2), 203–
213. https://doi.org/10.1127/phyto/2017/0160
McCreight, J. L., & Small, E. E. (2014). Modelling bulk density and snow water equivalent using daily
snow depth observations. The Cryosphere, 8(2), 521–536. https://doi.org/10.5194/tc-8-521-2014
McKay, G. A., & Findlay, B. (1971). Variation of snow resources with climate and vegetation in
Canada.

In

39th

Annual

Western

Snow

Conference.

Western

Snow

Conference.

http://sites/westernsnowconference.org/PDFs/1971McKay.pdf
Meybeck, M., Green, P., & Vörösmarty, C. (2001). A New Typology for Mountains and Other Relief
Classes: An Application to Global Continental Water Resources and Population Distribution.
Mountain Research and Development, 21(1), 34–45. http://www.jstor.org/stable/3674130
Mital, U., Dwivedi, D., Özgen-Xian, I., Brown, J. B., & Steefel, C. I. (2022). Modeling Spatial
Distribution of Snow Water Equivalent by Combining Meteorological and Satellite Data with
Lidar Maps. Artificial Intelligence for the Earth Systems, 1(4). https://doi.org/10.1175/AIES-D22-0010.1
Montgomery, D. C., Peck, E. A., & Vining, G. G. (2012). Introduction to Linear Regression Analysis.
Wiley. https://books.google.com.mx/books?id=0yR4KUL4VDkC
Moreno, I., Ruiz, A., Marturia, J., Oller, P., Piña Iglesias, J., García, C., Martínez, P., & Talaya, J.
(2009). Snowpack depth modelling and water availability from LIDAR measurements in eastern
Pyrenees.
Mortezapour, M., Menounos, B., Jackson, P. L., Erler, A. R., & Pelto, B. M. (2020). The role of
meteorological forcing and snow model complexity in winter glacier mass balance estimation,
Columbia

River

basin,

Canada.

Hydrological

Processes,

34(25),

5085–5103.

https://doi.org/10.1002/hyp.13929
Mukherjee, K., Menounos, B., Shea, J., Mortezapour, M., Ednie, M., & Demuth, M. N. (2023).
Evaluation of surface mass-balance records using geodetic data and physically-based modelling,
Place and Peyto glaciers, western Canada. Journal of Glaciology, 69(276), 665–682.
https://doi.org/10.1017/jog.2022.83

93

National Energy Board. (2019). Provincial and Territorial Energy Profiles – British Columbia. Energy
Production. https://www.neb-one.gc.ca/nrg/ntgrtd/mrkt/nrgsstmprfls/bc-eng.html
National Oceanic and Atmospheric Administration (NOAA) Coastal Services. (2012). Lidar 101 : An
Introduction to Lidar Technology, Data, and Applications. NOAA Coastal Services Center,
November, 76. https://doi.org/10.5194/isprsarchives-XL-7-W3-1257-2015
Neter, J., Wasserman, W., & Whitmore, G. A. (1988). Applied Statistics. Allyn and Bacon.
https://books.google.com.cu/books?id=RtftAAAAMAAJ
Nolan, M., Larsen, C., & Sturm, M. (2015). Mapping snow depth from manned aircraft on landscape
scales at centimeter resolution using structure-from-motion photogrammetry. The Cryosphere,
9(4), 1445–1463. https://doi.org/10.5194/tc-9-1445-2015
Nuth, C., & Kääb, A. (2011). Co-registration and bias corrections of satellite elevation data sets for
quantifying glacier thickness change. The Cryosphere, 5(1), 271–290. https://doi.org/10.5194/tc5-271-2011
Painter, T. H., Berisford, D. F., Boardman, J. W., Bormann, K. J., Deems, J. S., Gehrke, F., Hedrick,
A., Joyce, M., Laidlaw, R., Marks, D., Mattmann, C., McGurk, B., Ramirez, P., Richardson, M.,
Skiles, S. M., Seidel, F. C., & Winstral, A. (2016). The Airborne Snow Observatory: Fusion of
scanning lidar, imaging spectrometer, and physically-based modelling for mapping snow water
equivalent

and

snow

albedo.

Remote

Sensing

of

Environment,

184,

139–152.

https://doi.org/10.1016/j.rse.2016.06.018
Pelto, B. M. (2020). An approach to remotely monitor glacier mass balance at seasonal to annual time
scales, Columbia and Rocky Mountains, Canada. University of Northern British Columbia.
Pelto, B. M., Menounos, B., & Marshall, S. J. (2019). Multi-year Evaluation of Airborne Geodetic
Surveys to Estimate Seasonal Mass Balance, Columbia and Rocky Mountains, Canada. The
Cryosphere Discuss., 2019, 1–30. https://doi.org/10.5194/tc-2019-30
Pflug, J. M., & Lundquist, J. D. (2020). Inferring Distributed Snow Depth by Leveraging Snow Pattern
Repeatability: Investigation Using 47 Lidar Observations in the Tuolumne Watershed, Sierra
Nevada, California. Water Resources Research, 56(9). https://doi.org/10.1029/2020WR027243
Pigeon, K. E., & Jiskoot, H. (2008). Meteorological Controls on Snowpack Formation and Dynamics
in the Southern Canadian Rocky Mountains. Arctic, Antarctic, and Alpine Research, 40(4), 716–
730. https://doi.org/10.1657/1523-0430(07-054)[PIGEON]2.0.CO;2

94

Pistocchi, A. (2016). Simple estimation of snow density in an Alpine region. Journal of Hydrology:
Regional Studies, 6, 82–89. https://doi.org/10.1016/j.ejrh.2016.03.004
Pomeroy, J. W., Parviainen, J., Hedstrom, N., & Gray, D. M. (1998). Coupled modelling of forest snow
interception

and

sublimation.

Hydrological

Processes,

12(15),

2317–2337.

https://doi.org/10.1002/(SICI)1099-1085(199812)12:15<2317::AID-HYP799>3.0.CO;2-X
Province of British Columbia. (2021a). Archive Manual Snow Survey Data - Data Catalogue. Ministry
of

Environment

and

Climate

Change

Strategy

-

Knowledge

Management.

https://catalogue.data.gov.bc.ca/dataset/705df46f-e9d6-4124-bc4a-66f54c07b228
Province

of

British

Columbia.

(2021b).

Automated

Snow

Weather

Data.

https://www2.gov.bc.ca/gov/content/environment/air-land-water/water/water-sciencedata/water-data-tools/snow-survey-data/automated-snow-weather-station-data
Rasheed, B. A., Adnan, R., Saffari, S. E., & Dano Pati, K. (2014). Robust Weighted Least Squares
Estimation of Regression Parameter in the Presence of Outliers and Heteroscedastic Errors. Jurnal
Teknologi, 71(1). https://doi.org/10.11113/jt.v71.3609
RGI Consortium. (2017). Randolph Glacier Inventory – A Dataset of Global Glacier Outlines: Version
6.0: Technical Report. https://doi.org/10.7265/N5-RGI-60
RIEGL Laser Measurement Systems GmbH. (2017). LAS Extrabytes Implementation in RIEGL
Software.
http://www.riegl.com/uploads/tx_pxpriegldownloads/Whitepaper_LASextrabytes_implementati
on_in-RIEGLSoftware_2017-12-04.pdf
Rounce, D. R., Hock, R., Maussion, F., Hugonnet, R., Kochtitzky, W., Huss, M., Berthier, E.,
Brinkerhoff, D., Compagno, L., Copland, L., Farinotti, D., Menounos, B., & McNabb, R. W.
(2023). Global glacier change in the 21st century: Every increase in temperature matters. Science,
379(6627), 78–83. https://doi.org/10.1126/science.abo1324
Rousseeuw, P. J., & Leroy, A. M. (1987). Robust Regression and Outlier Detection. John Wiley &
Sons, Inc. https://doi.org/10.1002/0471725382
Seabold, S., & Perktold, J. (2010). Statsmodels: Econometric and Statistical Modeling with Python.
Proceedings of the 9th Python in Science Conference. https://doi.org/10.25080/majora-92bf1922011
Shaw, T. E., Deschamps-Berger, C., Gascoin, S., & McPhee, J. (2020). Monitoring Spatial and

95

Temporal Differences in Andean Snow Depth Derived From Satellite Tri-Stereo
Photogrammetry. Frontiers in Earth Science, 8. https://doi.org/10.3389/feart.2020.579142
Shea, J. M. (2010). Regional-scale distributed modelling of glacier meteorology and melt, southern
Coast

Mountains,

Canada

[The

University

of

British

Columbia].

https://open.library.ubc.ca/cIRcle/collections/ubctheses/24/items/1.0069481
Smyth, E. J., Raleigh, M. S., & Small, E. E. (2019). Particle Filter Data Assimilation of Monthly Snow
Depth Observations Improves Estimation of Snow Density and SWE. Water Resources Research,
55(2), 1296–1311. https://doi.org/10.1029/2018WR023400
Steppuhn, H. (1976). Areal water equivalents for prairie snow covers by centralized sampling.
Proceedings 44th Annual Meeting, Western Snow Conference, Calgary, AB, 63–68.
Sturm, M., Taras, B., Liston, G. E., Derksen, C., Jonas, T., & Lea, J. (2010). Estimating Snow Water
Equivalent Using Snow Depth Data and Climate Classes. Journal of Hydrometeorology, 11(6),
1380–1394. https://doi.org/10.1175/2010JHM1202.1
Taylor, J. R. (1997). An Introduction to Error Analysis. Journal of the Acoustical Society of America.
Teti, P. (2008). Effects of overstory mortality on snow accumulation and ablation.
Varhola, A., Coops, N. C., Weiler, M., & Moore, R. D. (2010). Forest canopy effects on snow
accumulation and ablation: An integrative review of empirical results. Journal of Hydrology,
392(3–4), 219–233. https://doi.org/10.1016/j.jhydrol.2010.08.009
Verschuren, J. P. (1972). Handbook on the Principles of Hydrology, edited by D. M. Gray. Published
by the Canadian National Committee for the International Hydrological Decade, Ottawa, Canada.
1971.

Journal

of

the

American

Water

Resources

Association,

8(1),

220–220.

https://doi.org/10.1111/j.1752-1688.1972.tb05116.x
Vionnet, V., Marsh, C. B., Menounos, B., Gascoin, S., Wayand, N. E., Shea, J., Mukherjee, K., &
Pomeroy, J. W. (2021). Multi-scale snowdrift-permitting modelling of mountain snowpack. The
Cryosphere, 15(2), 743–769. https://doi.org/10.5194/tc-15-743-2021
Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E.,
Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J.,
Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., … Vázquez-Baeza, Y. (2020).
SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3),
261–272. https://doi.org/10.1038/s41592-019-0686-2

96

Viviroli, D., Dürr, H. H., Messerli, B., Meybeck, M., & Weingartner, R. (2007). Mountains of the
world, water towers for humanity: Typology, mapping, and global significance. Water Resources
Research, 43(7). https://doi.org/10.1029/2006WR005653
World Glacier Monitoring Service (WGMS). (2023). FLUCTUATIONS OF GLACIERS DATABASE.
World Glacier Monitoring Service (WGMS). https://doi.org/10.5904/WGMS-FOG-2023-09
Zaremehrjardy, M., Razavi, S., & Faramarzi, M. (2021). Assessment of the cascade of uncertainty in
future snow depth projections across watersheds of mountainous, foothill, and plain areas in
northern

latitudes.

Journal

of

https://doi.org/10.1016/j.jhydrol.2020.125735

97

Hydrology,

598,

125735.

8.

Appendices
Appendix A: Selected snow-free pixels (i.e., pseudo-GCPs).
Latitude
(N)

Longitude
(W)

Elevation (m.a.s.l.)

1

50.82

-123.58

1719.64

2

50.84

-123.03

807.05

3

50.78

-123.33

2328.44

4

50.84

-122.86

669.17

5

50.84

-122.85

663.60

6

50.84

-122.86

740.07

7

50.86

-123.20

2417.32

8

50.87

-123.27

2396.77

9

50.90

-123.28

2552.33

10

50.92

-123.59

2611.52

11

50.85

-123.65

2427.56

12

50.81

-123.61

1916.61

13

50.79

-123.59

2000.85

14

50.78

-122.94

2181.57

15

50.83

-122.88

859.40

16

50.83

-122.85

770.61

17

50.88

-122.98

2221.86

18

50.86

-123.10

2472.16

19

50.84

-123.03

723.66

20

50.82

-122.85

891.78

21

50.83

-122.96

728.84

22

50.83

-123.06

753.17

23

50.82

-123.12

719.20

24

50.83

-123.19

835.97

ID

Appendix B: Glacier extents used in this study, extracted from the Randolph Glacier Inventory 6.0
(RGI Consortium, 2017).
RGIId

*Latitude
(N)

*Longitude
(W)

Min.

Max.

Mdn.

RGI60-02.02636
RGI60-02.02631
RGI60-02.02037
RGI60-02.02107
RGI60-02.02041
RGI60-02.02490

50.822
50.872
50.745
50.792
50.764
50.903

-123.639
-123.621
-123.306
-123.414
-123.343
-123.588

1388
1577
1691
1945
1763
1909

2885
2841
2668
2510
2701
2648

2289
2251
2202
2248
2317
2255

Elevation (m.a.s.l.)

98

Area (km2)

Name

79.336
24.644
13.161
6.048
5.687
5.426

Bridge
Sykora
Zavisha

RGI60-02.01826
RGI60-02.02784
RGI60-02.01922
RGI60-02.02526
RGI60-02.02752
RGI60-02.01923
RGI60-02.02745
RGI60-02.02613
RGI60-02.02049
RGI60-02.02751
RGI60-02.02671
RGI60-02.02274
RGI60-02.01919
RGI60-02.02042
RGI60-02.02145
RGI60-02.01921
RGI60-02.01938
RGI60-02.01920
RGI60-02.02111
RGI60-02.02132
RGI60-02.02064
RGI60-02.01910
RGI60-02.02059
RGI60-02.01874
RGI60-02.02038
RGI60-02.01835
RGI60-02.02785
RGI60-02.02093
RGI60-02.02514
RGI60-02.02747
RGI60-02.02786
RGI60-02.02444
RGI60-02.01937
RGI60-02.02670
RGI60-02.02324
RGI60-02.02216
RGI60-02.02615
RGI60-02.02590
RGI60-02.01988
RGI60-02.02094
RGI60-02.02202
RGI60-02.02606
RGI60-02.02220
RGI60-02.01964
RGI60-02.02289
RGI60-02.02036
RGI60-02.02252
RGI60-02.02110

50.695
50.991
50.738
50.919
50.956
50.733
50.974
50.792
50.776
50.953
50.964
50.875
50.730
50.776
50.814
50.727
50.748
50.732
50.802
50.810
50.783
50.736
50.780
50.722
50.733
50.710
50.974
50.797
50.935
50.976
50.978
50.898
50.746
50.960
50.889
50.854
50.785
50.953
50.764
50.798
50.849
50.855
50.857
50.757
50.883
50.762
50.870
50.806

-123.198
-123.421
-123.158
-123.574
-123.490
-123.112
-123.446
-123.513
-123.304
-123.465
-123.408
-123.248
-123.180
-123.130
-123.335
-123.084
-123.139
-123.148
-123.358
-123.295
-123.141
-123.213
-123.323
-123.193
-123.279
-123.182
-123.421
-123.384
-123.341
-123.474
-123.403
-123.289
-123.080
-123.402
-123.611
-123.147
-123.504
-123.403
-122.970
-123.370
-123.569
-123.813
-123.536
-123.272
-123.201
-123.315
-123.235
-123.372

1866
1979
1979
1989
2093
1938
1986
1898
1934
2101
2069
2122
2052
2142
2082
2010
2086
2121
2058
2016
2120
1954
2000
2078
2048
2288
2216
2139
2241
1908
2175
2259
2043
2085
2232
2119
2016
2217
2204
2284
2094
1201
2103
2035
2277
2261
2175
2189

99

2506
2465
2498
2538
2478
2555
2446
2375
2689
2470
2455
2450
2518
2540
2515
2412
2534
2487
2488
2516
2631
2433
2639
2633
2493
2673
2394
2387
2539
2509
2389
2522
2464
2435
2467
2470
2371
2388
2457
2474
2288
2954
2216
2322
2525
2622
2445
2355

2225
2197
2209
2218
2311
2257
2169
2162
2270
2298
2235
2285
2331
2325
2254
2254
2216
2261
2265
2269
2344
2197
2350
2347
2220
2498
2296
2232
2409
2209
2292
2404
2266
2260
2374
2338
2218
2335
2359
2355
2196
2436
2166
2164
2410
2518
2346
2263

3.898
3.114
3.109
2.967
2.823
2.502
2.463
2.101
2.035
1.830
1.335
1.325
1.234
1.234
1.065
0.978
0.856
0.836
0.811
0.785
0.743
0.717
0.691
0.687
0.592
0.591
0.568
0.504
0.480
0.448
0.448
0.416
0.397
0.392
0.378
0.350
0.345
0.323
0.308
0.301
0.296
0.246
0.244
0.232
0.227
0.212
0.206
0.203

McParlon

Surfusion

RGI60-02.02011
RGI60-02.02621
RGI60-02.02068
RGI60-02.01972
RGI60-02.01980
RGI60-02.02520
RGI60-02.02211
RGI60-02.02262
RGI60-02.02231
RGI60-02.02104
RGI60-02.02627
RGI60-02.01989
RGI60-02.02044
RGI60-02.02672
RGI60-02.01900
RGI60-02.02103
RGI60-02.01982
RGI60-02.02622
RGI60-02.01932
RGI60-02.02214
RGI60-02.02746
RGI60-02.02448
RGI60-02.02633
RGI60-02.01864
RGI60-02.02750
RGI60-02.02144
RGI60-02.02624
RGI60-02.02039
RGI60-02.01825
RGI60-02.02040
RGI60-02.02106
RGI60-02.02443
RGI60-02.02749
RGI60-02.02626
RGI60-02.02491
RGI60-02.02492
RGI60-02.02105
RGI60-02.01873
RGI60-02.01824
RGI60-02.02632
RGI60-02.02792
RGI60-02.02533
RGI60-02.02617
RGI60-02.02518
RGI60-02.01914
RGI60-02.02493
RGI60-02.02614
RGI60-02.01945

50.768
50.781
50.788
50.759
50.761
50.939
50.853
50.875
50.861
50.780
50.776
50.764
50.780
50.959
50.738
50.798
50.764
50.905
50.748
50.855
50.943
50.905
50.787
50.727
50.948
50.809
50.800
50.752
50.691
50.723
50.786
50.901
50.960
50.817
50.909
50.917
50.777
50.718
50.684
50.781
51.002
50.935
50.791
50.936
50.742
50.911
50.782
50.743

-123.303
-123.507
-123.151
-123.130
-123.136
-123.335
-123.166
-123.263
-123.544
-123.414
-123.711
-122.976
-123.501
-123.419
-123.037
-123.453
-123.308
-123.713
-123.071
-123.157
-123.465
-123.282
-123.674
-123.069
-123.476
-123.329
-123.723
-123.353
-123.212
-123.301
-123.392
-123.276
-123.509
-123.770
-123.626
-123.606
-123.402
-123.173
-123.206
-123.619
-123.367
-123.566
-123.529
-123.304
-122.989
-123.649
-123.516
-123.359

2217
2167
2064
2191
2265
2170
2216
2314
2019
2283
1725
2273
2028
2282
2163
2131
2273
1685
2143
2255
2213
2309
1564
2172
2352
2336
2664
2520
2155
2160
2277
2259
2094
1373
1900
2044
2401
2026
2238
2149
2159
2007
2113
2151
2141
2136
2198
1871

100

2489
2387
2351
2530
2568
2454
2356
2494
2155
2481
2827
2442
2218
2419
2279
2398
2518
2888
2348
2407
2330
2438
2864
2244
2433
2441
2829
2670
2434
2532
2376
2538
2448
2884
2650
2607
2487
2686
2461
2353
2494
2561
2272
2552
2348
2626
2392
2636

2320
2313
2148
2388
2400
2297
2284
2399
2085
2373
2183
2377
2140
2317
2221
2233
2386
2500
2283
2310
2289
2390
2220
2214
2421
2382
2718
2566
2349
2405
2311
2361
2207
2346
2295
2376
2437
2364
2340
2317
2300
2271
2181
2387
2257
2406
2306
2331

0.201
0.200
0.182
0.180
0.163
0.152
0.151
0.145
0.143
0.136
0.106
0.105
0.101
0.085
0.080
0.080
0.072
0.065
0.064
0.062
0.062
0.054
0.054
0.053
0.053
0.052
0.052
0.051
0.046
0.043
0.041
0.040
0.040
0.034
0.028
0.021
0.018
0.011
0.008
0.007
0.007
0.006
0.006
0.005
0.002
0.002
0.002
0.001

RGI60-02.02527
RGI60-02.02616

50.932
50.765

-123.600
-123.559

2015
1704

2627
2336

2303
2036

0.001
0.001

Appendix C: Manual snow samples; n= number of days considered in this study—data from the
province of British Columbia (2021b).
ID

Name

Latitude
(N)

Longitude
(W)

Elevation (m.a.s.l.)

n

1A01

Yellowhead

52.91

-118.55

1860

7

1A02

McBride Upper

53.31

-120.32

1610

42

1A05

Longworth (Upper)

53.95

-121.44

1693

80

1A06A Hansard

54.08

-121.87

608

22

1A10

Prince George arpt.

53.88

-122.68

689

57

1A11

Pacific Lake

54.38

-121.58

755

90

1A12

Kaza Lake

56.02

-126.29

1250

77

1A14

Hedrick Lake

54.11

-121.01

1110

50

1A15

Knudsen Lake

54.30

-120.78

1580

86

1A16

Burns Lake

54.23

-54.23

800

79

1A18

Holmes River

53.28

-119.46

1900

10

1A19

Dome Mountain

53.62

-121.02

1766

42

1A21

Narrow Lake

53.58

-121.9

1650

10

1A23

Bird Creek

53.67

-124.34

1180

62

1A24

Forfar Creek (Upper)

54.98

-125.58

1410

6

1B01

Mount Wells

53.73

-126.42

1490

78

1B02

Tahtsa Lake

53.59

-127.66

1300

82

1B05

Skins Lake

53.78

-125.99

890

65

1B06

Mount Swannell

53.34

-125.27

1620

71

1B07

Nutli Lake

53.46

-126.29

1490

75

1C01

Brookmere

49.81

-120.87

980

48

1C05

McGillivray Pass

50.89

-122.62

1725

75

1C06

Pavilion

50.91

-121.82

1230

21

1C07

Lac Le Jeune Lower

50.68

-120.59

1370

18

1C08

Nazko

53.02

-123.65

1070

59

1C09A Highland Valley

50.49

-120.98

1510

35

1C13A Horsefly Mountain

52.34

-121.06

1550

51

1C14

Bralorne

50.77

-122.79

1389

71

1C17

Mount Timothy

51.89

-121.26

1660

90

1C19

Gnawed Mountain

50.43

-120.97

1580

41

1C21

Big Creek

51.71

-123.04

1140

54

1C22

Puntzi Mountain

52.13

-124.09

940

56

Penfold Creek

52.75

-120.56

1685

60

1C23

101

Lac Le Jeune Upper

50.45

-120.49

1509

80

1C28

Duffey Lake

50.38

-122.48

1200

35

1C29

Shovelnose Mountain

49.86

-120.86

1450

57

1C32

Deadman River

51.10

-120.61

1430

37

1C33

Granite Mountain

52.53

-122.27

1175

13

1C33A Granite Mountain

52.53

-122.27

1150

71

1C37

Bralorne (Upper)

50.79

-122.75

1981

64

1C38

Downton Lake Upper

50.86

-123.18

1887

65

1C39

Bridge Glacier Lower

50.85

-123.46

1390

80

1C40

Tyaughton Creek (N)

51.15

-122.79

1947

76

1C42

Caverhill Lake New

51.33

-120.38

1400

26

1D01

Coquitlam Lake

49.36

-122.79

300

1

1D06

Tenquille Lake

50.54

-122.93

1680

16

1D08

Stave Lake

49.58

-122.32

1250

75

1D09

Wahleach Lake

49.23

-121.58

1480

72

1D10

Nahatlatch River

49.83

-122.06

1550

63

1D13

Wolverine Creek

50.51

-122.98

250

16

1D16

Dickson Lake

49.32

-122.07

1160

67

1D18

Disappointment Lake

49.55

-122.76

1050

101

1E01B

Blue River

52.17

-119.24

690

1

Blue River

52.13

-119.29

670

79

1E02A Mount Cook

52.18

-119.32

1570

7

1E03A Trophy Mountain

51.81

-119.94

1860

64

1E05

Knouff Lake

50.98

-120.13

1200

40

1E06

Cook Forks

52.17

-119.32

1390

6

1E07

Adams River

51.58

-119.43

1720

100

1E08

Azure River

52.60

-119.72

1640

3

1E13

North Clemina Creek

52.56

-118.95

1860

11

1F01A

Aberdeen Lake

50.15

-119.05

1310

56

1F02

Anglemont

50.99

-119.19

1190

63

1F04

Enderby

50.65

-118.93

1900

78

2A01A Canoe River

52.78

-119.27

870

37

2A02

51.25

-117.49

1250

85

2A03A Field

51.39

-116.51

1285

56

2A07

Kicking Horse

51.44

-116.36

1650

86

2A11

Beaverfoot

51.28

-116.88

1890

80

2A14

Mount Abbot

51.24

-117.51

2010

84

2A16

Goldstream

51.70

-118.45

1920

81

2A17

Fidelity Mountain

51.24

-117.69

1870

85

2A18

Keystone Creek

51.41

-118.36

1890

80

2A19

Vermont Creek

50.96

-116.94

1520

80

2A20

Watam Lake

51.33

-118.51

1810

1

Sunbeam Lake

51.59

-117.65

2010

78

1C25

1E01B

2A22

Glacier

102

Bush River

51.72

-117.42

1920

71

2A25

Kirbyville Lake

51.57

-118.77

1750

82

2A27

Downie Slide Lower

51.51

-118.53

980

71

2A29

Downie Slide Upper

51.50

-118.55

1630

77

2B02A Farron

49.27

-118.13

1220

94

2B04

Whatshan Lower

50.13

-117.98

1250

1

2B05

Whatshan Upper

50.19

-118.04

1525

69

2B06

Barnes Creek

50.07

-118.36

1620

55

2B07

Koch Creek

49.72

-117.98

1860

64

2B08

St. Leon Creek

50.44

-117.7

1860

48

2B09

Record Mountain

49.09

-117.87

1890

112

2C01

Sinclair Pass

50.66

-115.96

1370

43

2C04

Sullivan Mine

49.73

-116.02

1550

90

2C07

Fernie East

49.51

-115.02

1250

73

2C09A Morrissey Ridge

49.45

-114.97

1800

5

2C11

Kimberley (Upper)

49.56

-116.08

2140

26

2C12

Kimberley Middle

49.56

-116.05

1680

23

2C14

Floe Lake

51.06

-116.14

2090

56

2C15

Mount Assiniboine

50.92

-115.61

2230

68

2C16

Mount Joffre

50.53

-115.13

1750

79

2C17

Thunder Creek

50.05

-115.23

2010

75

2C20

Vermilion riv. No. 3

51.20

-116.07

1570

27

2D02

Ferguson

50.68

-117.46

880

80

2D03

Sandon

49.98

-117.23

1070

29

2D04

Nelson

49.42

-117.23

930

80

2D05

Gray Creek Lower

49.61

-116.68

1650

75

2D06

Char Creek

49.10

-116.95

1310

104

2D07A Duncan Lake No. 2

50.26

-116.96

630

41

2D09

Mount Templeman

50.72

-117.2

1860

44

2D10

Gray Creek Upper

49.62

-116.66

1930

77

2E01

Monashee Pass

50.08

-118.51

1370

77

2E02

Carmi

49.50

-119.09

1250

37

2E03

Big White Mountain

49.72

-118.97

1680

90

2E06

Bluejoint Mountain

49.54

-118.53

2040

22

2F01

Trout Creek

49.74

-120.18

1430

47

2F01A

Trout Creek West

49.74

-120.19

1430

53

2F02

Summerland Res.

49.81

-120.02

1280

78

2F03

McCulloch

49.79

-119.19

1280

55

2F04

Graystoke Lake

49.99

-118.87

1840

69

2F05

Mission Creek

49.94

-118.95

1780

1

2F07

Postill Lake

49.99

-119.21

1370

66

2F08

Greyback Reservoir

49.62

-119.42

1550

74

Whiterocks Mountain

50.02

-119.75

1780

89

2A23

2F09

103

Silver Star Mountain

50.37

-119.06

1840

77

2F11

Isintok Lake

49.55

-119.97

1680

86

2F12

Mount Kobau

49.12

-119.68

1810

107

2F13

Esperon Creek upp.

50.08

-119.74

1620

56

2F14

Esperon Creek mid.

50.07

-119.7

1420

29

2F18

Brenda Mine

49.86

-119.98

1460

32

2F19

Oyama Lake

50.12

-119.28

1340

63

2F20

Vaseux Creek

49.28

-119.33

1400

49

2F21

Bouleau Lake

50.28

-119.65

1400

50

2F23

MacDonald Lake

49.88

-120.03

1740

51

2F24

Islaht Lake

49.99

-119.81

1480

73

2F25

Postill Lake Upper

49.96

-119.21

1540

33

2G04

Lost Horse Mountain

49.28

-120.13

1940

75

2G05

Missezula Mountain

49.74

-120.54

1550

76

2G06

Hamilton Hill

49.49

-120.79

1490

68

3A01

Grouse Mountain

49.38

-123.08

1100

81

3A02

Powell River (Upper)

50.27

-124.33

1040

15

3A05

Powell River Lower

50.26

-124.34

910

15

3A09

Palisade Lake

49.45

-123.03

880

44

3A10

Dog Mountain

49.37

-122.96

1080

119

3A19

Orchid Lake

49.54

-123.05

1190

112

3A20

Callaghan Creek

50.14

-123.1

1040

75

3A26

Chapman Creek

49.58

-123.59

1022

15

3A27

Edwards Lake

49.59

-123.63

1070

13

Forbidden Plateau

49.65

-125.21

1100

83

3B02A Mount Cokely

49.25

-124.59

1190

29

3B04

Elk upp. riv.

49.86

-125.82

270

36

3B10

Thelwood Lake

49.52

-125.7

990

46

3B13

Heather Mountain

48.94

-124.45

1170

5

3B18

Wolf River Middle

49.71

-125.7

990

77

3B19

Wolf River Lower

49.72

-125.7

640

65

3B22

Tennent Lake

49.58

-125.64

1134

4

3B24

Heather mt. upp.

48.94

-124.45

1190

4

3C07

Wedeene River South

54.23

-128.73

220

36

3D01C Sumallo River West

49.25

-121.25

790

55

3D02

49.05

-120.85

1220

51

3D03A Klesilkwa

49.13

-121.31

1175

65

4A01

Mcleod Lake

54.55

-122.53

980

1

4A02

Pine Pass

55.35

-122.64

1450

59

4A03

Ware (Upper)

57.39

-125.71

1575

85

4A04

Ware Lower

57.41

-125.67

970

78

4A05

Germansen Upper

55.82

-124.69

1480

84

Tutizzi Lake

56.31

-12559

1045

84

2F10

3B01

4A06

Lightning Lake

104

Lady Laurier Lake

56.70

-123.75

1440

84

4A09

Pulpit Lake

57.55

-126.75

1335

83

4A10

Fredrickson Lake

56.95

-126.52

1325

84

4A11

Trygve Lake

56.99

-127.46

1410

78

4A12

Tsaydaychi Lake

55.42

-124.76

1190

84

4A13

Philip Lake

55.11

-123.89

1035

82

4A14

Pink Mountain

57.03

-122.51

1170

6

4A16

Morfee Mountain

55.42

-123.04

1430

82

4A18

Mount Sheba

54.53

-121.8

1490

82

Mackenzie Airport

55.30

-123.14

700

3

4A19A Mackenzie Airport

55.18

-123.08

700

3

4A20

Monkman Creek

54.73

-121.23

1570

72

4A21

Mount Stearns

57.07

-123.40

1505

85

4A25

Fort St. John Airport

56.24

-120.73

690

57

4A28

Bullhead Mountain

56.03

-122.10

798

10

4B01

Kidprice Lake

53.86

-127.44

1370

79

Johanson Lake

56.60

-126.18

1420

83

4B03A Hudson Bay Mt.

54.77

-127.28

1480

116

4B04

Chapman Lake

54.88

-126.73

1460

55

4B06

Tachek Creek

54.64

-126.29

1140

60

4B07

McKendrick Creek

54.80

-126.8

1050

57

4B08

Mount Cronin

54.93

-126.81

1480

55

4B10

Ningunsaw Pass

56.78

-129.95

690

52

4B11A Bear Pass

56.10

-129.63

460

46

4B13A Terrace Airport

54.47

-128.58

180

50

4B14

Equity Mine

54.19

-126.25

1420

50

4B15

Lu Lake

54.20

-126.3

1300

50

4C01

Sikanni Lake

57.26

-124.13

1385

83

4C02

Summit Lake

58.65

-124.72

1280

46

4C03

Dease Lake

58.43

-130.01

820

60

4C05

Fort Nelson Airport

58.84

-122.57

380

59

4C09

Deadwood River

59.26

-128.33

1300

3

4C15

Jade City

59.24

-129.66

940

51

4D01

Telegraph Creek

57.95

-131.14

580

29

4D02

Iskut

57.85

-129.99

1000

40

4D10

Tumeka Creek

57.24

-129.71

1220

3

4D11

Kinaskan Lake

57.58

-130.27

1020

3

4D14

Wade Lake

58.19

-128.91

1370

3

Log Cabin

59.76

-134.97

900

57

4E02A Atlin Lake

59.59

-133.73

710

22

4E02B

59.59

-133.71

730

39

4A07

4A19

4B02

4E01

Atlin Lake

105

Appendix D: Automated snow pillows and the number of days (n) with valid SWE and depth
measurements (BC, 2021a).
Pillow Name
1A01P
1A02P
1A03P
1A05P
1A12P
1A14P
1A15P
1A17P
1A19P
1C14P
1C29P
1C41P
1D06P
1D17P
1D18P
1D19P
1E02P
1E08P
1E10P
1E14P
1F03P
1F04P
1F06P

Yellowhead Lake
McBride Upper
Barkerville
Longworth Upper
Kaza Lake
Hedrick Lake
Knudsen Lake
Revolution Creek
Dome Mountain
Bralorne
Shovelnose Mountain
Yanks Peak
Tenquille Lake
Chilliwack River
Disappointment Lake
Spuzzum
Mount Cook
Azure River
Kostal Lake
Cook Creek
Park Mountain
Enderby
Celista

2C10P
2D14P
2E07P
2F01AP
2F05P
2F08P
2F10P
2G03P
3A09P
3A20P
3A22P
3A24P
3A25P
3A28P
3B23P
3B24P
3B26P

Moyie Mountain
Redfish Creek
Grano Creek
Trout Creek West
Mission Creek
Greyback Reservoir
Silver Star Mountain
Blackwall Peak
Palisade Lake
Callaghan
Nostetuko River
Mosley Creek Upper
Squamish River Upper
Tetrahedron
Jump Creek
Heather Mountain Upper
Mount Arrowsmith

Latitude
(N)

Longitude
(W)

Elevation
(m.a.s.l.)

n

52.915
53.308
53.051
53.951
56.021
54.102
54.303
53.796
53.627
50.775
49.856
52.830
50.542
49.017
49.555
49.661
52.169
52.596
52.203
52.168
50.441
50.653
51.407
49.254
49.695
49.559
49.740
49.949
49.616
50.371
49.085
49.454
50.137
51.237
51.779
50.153
49.598
48.970
48.944
49.207

-118.542
-120.317
-121.491
-121.443
-126.298
-121.003
-120.783
-120.373
-121.017
-122.789
-120.863
-121.356
-122.932
-121.712
-122.758
-121.661
-119.322
-119.723
-120.017
-119.304
-118.627
-118.927
-118.999

1860
1611
1520
1740
1257
1100
1601
1690
1774
1382
1460
1670
1680
1600
1050
1180
1550
1652
1770
1280
1890
1950
1500
1930
2104
1860
1420
1780
1550
1839
1940
900
1017
1500
1650
1340
1420
1160
1190
1465

4472
3589
4129
1093
1305
4964
1312
3883
4036
641
497
4769
4070
4328
770
4288
4253
4007
2994
2121
2893
1331
4099

106

-115.779
-117.085
-118.678
-120.189
-118.949
-119.423
-119.062
-120.779
-123.032
-123.103
-124.458
-124.627
-123.441
-123.605
-124.276
-124.452
-124.568

4790
4965
4151
746
4502
922
1410
4962
663
230
2846
3241
4189
689
4006
1111
737

3C08P
4A02P
4A09P
4A27P
4A30P
4A36P
4B15P
4B16P
4B17P
4B18P
4D11P

Burnt Bridge Creek
Pine Pass
Pulpit Lake
Kwadacha North
Aiken Lake
Parsnip Upper
Lu Lake
Shedin Creek
Tsai Creek
Cedar-Kiteen
Kinaskan Lake

52.474
55.353
57.549
57.623
56.437
54.606
54.197
55.862
54.646
55.153
57.579

107

-126.237
-122.637
-126.748
-125.073
-125.742
-122.138
-126.304
-127.698
-127.674
-128.712
-130.263

1330
1400
1311
1554
1050
790
1300
1480
1360
885
1020

3750
1551
1762
944
1472
580
4109
4493
4825
3671
1290