ASSESSMENT OF LIDAR SNOW DEPTH MEASUREMENTS AND THE SPATIOTEMPORAL VARIABILITY OF THE SNOWPACK IN A FORESTED WATERSHED ON
VANCOUVER ISLAND

by

Alison Bishop
Bachelor of Science in Environmental Science, Dalhousie University, 2016
Advanced Diploma in GIS, Vancouver Island University, 2019

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

UNIVERSITY OF NORTHERN BRITISH COLUMBIA
May 2024
©Alison Bishop, 2024

Abstract
This study uses repeat airborne LiDAR surveys to assess snow distribution within complex forested
terrain in Russell Creek, a 33 km2 sub-basin of the Tsitika watershed, on northern Vancouver Island.
LiDAR surveys used within this study were acquired in 2022 and 2023, with 4 – 5 flights per year timed
with the intent to capture maximum snow depth, to the completion of the melt season. In addition, two
separate bare earth models were used to assess error by acquisition date (2017 vs 2023). Over 2000
manual snow depth measurements were taken over the course of this study via a) cardinal plots within
six different forest cover types and b) weather station snow courses conducted on gravel roads.
Mean difference (MD) of the averaged manual and LiDAR measurements were overall lowest for the
weather station snow courses (-29 to 13 cm) and highest in the harvested (juvenile and regenerating)
plots (-39 to 134 cm). Analysis of the LiDAR bare earth point returns showed that juvenile and
regenerating plots both had a low percentage of pixels with ground returns – average of 24 and 35%
respectively - a result of the tall (~ 2 meter) and complex ground vegetation. Error was reduced in the
juvenile forest plots (-22 to 29 cm) when LiDAR snow depth was processed with an alternative bare
earth model acquired in 2017, in which the average coverage of ground returns was much greater (96%).
Total snow storage within Russell Creek was relatively similar (<5% difference) between the two bare
earth models, and the higher elevations (>1000 m) of the watershed - comprised of old growth forest and
alpine cover types – store the majority (56 – 82%) of the snow. At high elevations (1400 – 1700 m) in
the watershed, the snow volumes were much greater (28 – 44%) when processed with the 2023 bare
earth model compared to the 2017.
An additional component to this study was a paired control-treatment approach to assess the impacts of
forest harvest on a) the bare earth model and b) snow distribution. Evaluation of pre- (2020) and postii

harvest (2023) snow free models showed a slight (0.15 m) bias within treatment sites due to changes to
the ground surface. After bias correction, the rate of increase in snow water equivalence between the
pre-and post-harvest snow years in the treatment sites (1.4 – 4.2) doubled the control (0.7 – 1.7).
Overall, these results highlighted the importance of the bare earth model to accurately measure snow,
especially within complex forested and harvested watersheds.

iii

Table of Contents
Abstract ....................................................................................................................................................... ii
Table of Contents ....................................................................................................................................... iv
List of Tables ............................................................................................................................................ vii
List of Figures ............................................................................................................................................. x
Acknowledgement .................................................................................................................................. xvii
Chapter 1 : Introduction .............................................................................................................................. 1
Chapter 2 : Background .............................................................................................................................. 4
2.1 Snow and Remote Sensing ................................................................................................................ 4
2.1.1 LiDAR systems........................................................................................................................... 6
2.1.2 Data processing........................................................................................................................... 7
2.1.3 LiDAR validation ....................................................................................................................... 8
2.1.4 LiDAR Error and Uncertainty .................................................................................................... 9
2.2 Snow distribution in forests ............................................................................................................ 11
2.3 Forest Disturbance and Hydrological Recovery ............................................................................. 14
2.4 Research Gaps and Objectives ........................................................................................................ 18
Chapter 3 : Study Area and Data ............................................................................................................. 19
3.1 Study Area ....................................................................................................................................... 19
3.2 Weather Station Network ................................................................................................................ 20
3.3 Land Cover Classification ............................................................................................................... 21
3.4 LiDAR surveys ................................................................................................................................ 23
3.4.1 Snow Free acquisition .............................................................................................................. 24
3.5 In-Situ Observations ........................................................................................................................ 27
3.5.1 Snow Sampling ......................................................................................................................... 27
3.5.2 Sampling error and mitigation .................................................................................................. 33
3.5.3 Vegetation surveys ................................................................................................................... 34
Chapter 4 : Data Processing and Analyses ............................................................................................... 35
4.1 LiDAR processing ........................................................................................................................... 36
4.2 In-Situ Observations and Supplementary Data ............................................................................... 39
4.2.1 Weather Stations ....................................................................................................................... 40
4.2.2 Game Cameras .......................................................................................................................... 40
iv

4.3 Plot Analyses ................................................................................................................................... 40
4.3.1 Test of Normality...................................................................................................................... 42
4.3.2

Snow Depth Variability ........................................................................................................ 42

4.3.3 Mean Snow Depth .................................................................................................................... 43
4.3.4 Timing error and mitigation ..................................................................................................... 44
4.4 Watershed Scale Analysis ............................................................................................................... 46
Chapter 5 : Results .................................................................................................................................... 47
5.1 Snowpack and Weather ................................................................................................................... 47
5.2 Vegetation Surveys ......................................................................................................................... 51
5.3 LiDAR Point Densities.................................................................................................................... 53
5.4 In-situ vs LiDAR Snow Depths ...................................................................................................... 55
5.4.1 Snow depths Validation ............................................................................................................ 56
5.4.2 Variability of Observed Snow Depths ...................................................................................... 70
5.5 Bare Earth DEM .............................................................................................................................. 76
5.6 Watershed Snow Storage ................................................................................................................ 79
5.6.1 Snow Volumes by Elevation .................................................................................................... 80
5.6.2 Snow Volumes by Forest Cover ............................................................................................... 82
Chapter 6 : Forest Harvest Experiment ..................................................................................................... 85
6.1 Introduction ..................................................................................................................................... 85
6.2 Methods ........................................................................................................................................... 86
6.2.1 Snow free surface assessment................................................................................................... 88
6.2.2 Snow distribution and volume change ...................................................................................... 90
6.3 Results ............................................................................................................................................. 91
6.3.1 Weather Station Overview ........................................................................................................ 91
6.3.2 Ground surface change ............................................................................................................. 93
6.3.3 Snow Distribution and water storage ........................................................................................ 96
6.4 Discussion and Recommendations .................................................................................................. 98
Chapter 7 : Discussion and Recommendations ....................................................................................... 102
7.1 LiDAR Accuracy within complex ground cover .......................................................................... 102
7.2 Impact of the Bare Earth Model .................................................................................................... 106
7.3 Snow Distribution: Plot and Watershed Scale .............................................................................. 109
7.3.1 Plot scale distribution ............................................................................................................. 109
7.3.2 Watershed snow storage ......................................................................................................... 109
v

7.4 Study Limitations and Recommendations..................................................................................... 113
Chapter 8: Summary ............................................................................................................................... 116
Data Availability Statement .................................................................................................................... 117
References ............................................................................................................................................... 118
Appendices .............................................................................................................................................. 127
Appendix A: ACO 2023 Bare Earth Model ........................................................................................ 127
lasAppendix B: LiDAR Processing..................................................................................................... 129
Appendix C: GNSS report................................................................................................................... 130
Appendix D: Snow depth variability results tables ............................................................................. 133
Appendix E: Old Growth Forest: Canopy Snow Depth Variability.................................................... 136
Appendix F: Average snow depth results tables ................................................................................. 137
Appendix G: Bias corrections for meteorological variability ............................................................. 141
Appendix H: Watershed Snow Volume .............................................................................................. 142
Appendix I: Harvest Experiment Slash Pile Analysis ......................................................................... 146

vi

List of Tables
Table 3.1 Breakdown of the land cover categories within the Russell Creek watershed, which were used
to: a) determine plot locations and b) for watershed scale analyses. ........................................................ 22
Table 3.2 LiDAR products and metadata for the Tsitika survey area from 2020 to 2023 used within this
research project ......................................................................................................................................... 26
Table 3.3 Meta data for cardinal plots used within the in-situ sampling campaign. Precise plot locations
and GNSS accuracies for the 2022 and 2023 sampling campaigns can be found in Appendix C............ 29
Table 3.4 Game camera sites and metrics for the 2022 and 2023 LiDAR campaigns. Precise position of
the game camera stakes can be found in Appendix C .............................................................................. 32
Table 5.1 Average daily air temperature across all weather stations. Highlighted values represent that the
temperature was higher between the two years, and a dash indicates there was no data for that month. 50
Table 5.2 Summary of 2022 and 2023 snow depth measurements at four weather stations, selected for
consistent snow depth data across 2022 and 2023. Metrics include: 1) peak snow depth and date, 2)
length of melt period in days (defined as sustained loss of 0.5 cm per day) and rate of snow melt, and 3)
snow free date ........................................................................................................................................... 50
Table 5.3 Vegetation composition by cover type, summarized by the average height (μ) and the standard
deviation (σ) of the shrubs and ferns, woody debris, and heather. ........................................................... 51
Table 5.4 Median ground point density per meter squared for the snow surveys, as well as the average
proportion of valid pixels per each of the cover types .............................................................................. 54
Table 5.5 Summary table of LiDAR and manual dates for 2022 and 2023 and how the snow depth
changed within that time frame. Snow depth change includes the SR50 data of all seven weather
stations in Russell Creek, provided the sensor was collecting data in that time frame ............................ 56

vii

Table 5.6 Standard Deviation (SD) of the in-situ snow depth observations at each weather station snow
course. Snow courses which were snow free are shaded blue. To assess whether the in-situ snow depth
variability changed as snow melted, the Levene’s test was conducted, where for each weather station
and year, each groups (LiDAR survey) variance were compared to each other. Snow free and no data
surveys were not included as groups. ....................................................................................................... 71
Table 5.7 Standard Deviation of the LiDAR collected snow depth at each weather station snow course.
Snow courses which were snow free are shaded blue. To assess whether the LiDAR snow depth
variability changed as snow melted, the Levene’s test was conducted, where for each weather station
and year, each groups (LiDAR survey) variance were compared to each other. Snow free and no data
surveys were not included as groups. ....................................................................................................... 71
Table 5.8 Standard deviation for in-situ snow depth measurements for the cardinal plot sampled across
2022 and 2023. Cells highlighted blue indicate snow-free. To assess whether the in-situ snow depth
variability changed as snow melted, the Levene’s test was conducted, where for each plot and year,
groups (LiDAR survey) variance were compared to each other. Snow free and no data surveys were not
included as groups..................................................................................................................................... 73
Table 5.9 Standard deviation for LiDAR cardinal plot snow depth measurements across both years. Cells
highlighted blue indicate snow-free. To assess whether the LiDAR snow depth variability changed as
snow melted, the Levene’s test was conducted, where for each plot and year, groups (LiDAR surveys)
variance were compared to each other. Snow free and no data surveys were not included as groups. .... 74
Table 5.10 Snow storage by elevation band, with percent of overall snow stored in the Russell Creek
watershed shown in brackets, between the two different bare earth models (2017 vs 2023), for the 9th of
March 2022 and 10th March 2023 acquisitions. Highlighted values indicated which model predicted the
greater snow volume per elevation band. ................................................................................................. 82

viii

Table 5.11 Snow storage by cover type, with the percent of overall snow stored in the Russell Creek
watershed shown in brackets, between the two different bare earth models (2017 vs 2023), for the 9th of
March 2022 and 10th March 2023 acquisitions. Highlighted shows which model predicted more snow
per cover type ............................................................................................................................................ 83
Table 6.1 Topographic and canopy metrics for the control, treatment and road area which was used
within this study. Metrics for slope, canopy cover, and height represent the average value. ................... 88
Table 6.2 LiDAR surveys for 2022 and 2023 paired for comparing snow distribution change pre-harvest
(2022) and post-harvest (2023). Weather station snow depth data from Lower Cain (1260 m) included
for each survey .......................................................................................................................................... 90
Table 6.3 Results from pre and post-harvest ground surface changes, in which the control, treatment, and
road plot types average and standard deviation elevation were calculated, and error metrics mean bias
error (MBE), mean absolute error (MAE) and root mean square error (RMSE) were calculated between
the 2022 (pre-harvest) and 2023 (post-harvest) bare earth models. .......................................................... 93
Table 6.4 Snow water equivalence (mm) for both the pre-harvest (2022) and post-harvest years (2023),
+/- standard error, for the three snow on LiDAR surveys. Rate of change represents the increase or
decrease in SWE from 2022 to 2023, for the control and treatment blocks ............................................. 98
Table 7.1 Error results from airborne LiDAR snow depth studie which evaluated LiDAR accuracy using
in-situ measurements and in forested terrain. ......................................................................................... 103

ix

List of Figures
Figure 3.1 Overview map of the study site: Russell Creek (32 km2) sub-basin of the Tsitika watershed
(397 km2), located Northern Vancouver Island, British Columbia .......................................................... 20
Figure 3.2 Overview of Russell Creek watershed showing A) the landcover type, B) the landcover
broken down by elevation, C) the orthoimagery of the watershed, and D) aspect ................................... 22
Figure 3.3 Images of the forest cover types where plots were sampled within this research project:
Juvenile (A), Regenerating (B), Dense Old Growth (C), Old Growth (D) and Alpine (E) ...................... 23
Figure 3.4 LiDAR survey area within the Tsitika watershed, showing the ACO 2023 bare earth
boundary and snow surveys(red), the 2017 BCTS bare earth coverage of the AOI (orange), and the ACO
2020 acquisition which was targeted to fill the missing area over Mount Cain (blue). Note that the linear
boundary to the SE was used to cover the highway and provide generally snow free stable terrain ....... 25
Figure 3.5 Location of the game camera sites (blue), weather stations (purple), and cardinal plots (red)
within the Russell Creek watershed .......................................................................................................... 27
Figure 3.6 Cardinal plot snow sampling strategy of density (open circles) and snow depth (dash)
measurements. Distances (red) are in meters ............................................................................................ 28
Figure 3.7 GNSS set-up showing the base station (A) established on a permanent ground control point,
and the rover (B) set up at the plot centre of Juvenile site ST2CC2 ......................................................... 30
Figure 3.8 Game camera snow stake set up – this imagery is of game camera site ST2GC1 as seen in
snow free and snow on conditions ............................................................................................................ 32
Figure 4.1 Hillshade of snow-on surface for the Russell Creek watershed, showing the A) interpolated
(gap-filled) raster surface and B) the non-gap filled (pixels with zero ground returns removed, and
shown in red) for March 10th 2023 LiDAR survey. For both years, the area of no ground point returns
increased through the snow melt season. .................................................................................................. 37
x

Figure 4.2 Histograms of the snow-free pixel values of the stable terrain for each snow survey
acquisition in 2023. The ground point density has been classified by greater than five points per m2
(green), 2 – 5 points per m2 (yellow) and less than two points per m2 (grey). Pixels with no ground
points were excluded from this analysis. .................................................................................................. 39
Figure 4.3 Ortho-imagery and snow depth data for dense old growth forest plot, ST1OG6, from the
April 29th 2023 survey. The in-situ (A) and LiDAR (B) snow depths are averaged separately, and mean
difference calculated. Missing pixels in the LiDAR are a result of no valid snow surface returns for the
survey. ....................................................................................................................................................... 41
Figure 5.1 Hourly air temperature (upper) and snow depth (lower) data recorded from Steph 2(1164m)
weather station, for water years 2022 and 2023 ....................................................................................... 48
Figure 5.2 Snow depth records from 2016 - 2023 for Steph 6 (upper) and Steph 2 (lower) weather
station, with the years of this study highlighted with a blue (2022) and red (2023) dashed line
respectively. Of the 7 years of data shown for Steph 6 and Steph 2 weather stations. 2022 had the latest
date on record for snow free, whereas 2023 had the second earliest ........................................................ 49
Figure 5.3 Average vegetation height by cover type broken into the categories of shrub, woody debris,
and heather. Vegetation measured 5 cm or less, as well as trees, were not included in the results. Error
bars represent the 95% confidence interval. ............................................................................................. 52
Figure 5.4 Percentage of the cardinal plots which had at least one LiDAR ground or snow surface point
return per m2, grouped by forest cover type, for each LiDAR acquisition in 2022 (upper) and 2023
(lower). Hatched bars denotes acquisitions that were snow-free for all plots within the cover type. Error
bars represent the 95% confidence interval. ............................................................................................. 53
Figure 5.5 Imagery of the juvenile forest from the (A) May 25th 2023 and (B) June 12th 2023
acquisitions showing the green up of the vegetation ................................................................................ 55

xi

Figure 5.6 Steph 2 weather station showing the change in average daily snow depth which occurred
between the manual sampling window and the LiDAR survey for the April 29th 2023 acquisition ........ 57
Figure 5.7 The mean difference for snow courses from the 2022 and 2023 snow season. A 10 day delay
between the LiDAR acquisition and manual sampling -which coincided with a period of snow melt - is
attributed to the large negative bias for June 24th 2022 ............................................................................ 58
Figure 5.8 2022 and 2023 snow-on surveys mean difference grouped by whether LiDAR was flown pre
or post manual sampling, and the snow depth change between acquisition and sampling. Snow depth
readings are compiled from the weather stations in Russell Creek, Steph 1 – 8. ..................................... 59
Figure 5.9 Alpine cardinal plots mean difference between LiDAR and in-situ snow depth measurements
for the 2022 and 2023 surveys. Meteorological corrections from the weather station snow depth data is
denoted with the darker bars, and asterisk mark if LiDAR was acquired after manual measurements.
Hatched bars indicate plot was snow free. ................................................................................................ 61
Figure 5.10 Dense old growth cardinal plots mean difference between LiDAR and in-situ snow depth
measurements for the 2022 and 2023 surveys. Meteorological correction applied using the game camera
data is denoted with the darker bars, and asterisk mark if LiDAR was acquired after manual
measurements. Hatched bars indicate plot was snow free. ....................................................................... 62
Figure 5.11 Game camera ST1GC1 showing the rapid snow melt which occurred, with a loss of
approximately 10 cm between 10 am on the 28th of April to 3 pm on the 29th of April 2023. Snow
stakes are marked by 20 cm increments. .................................................................................................. 63
Figure 5.12 Old growth cardinal plots mean difference between LiDAR and in-situ snow depth
measurements for the 2022 and 2023 surveys. Error corrected using the game camera imagery is denoted
with the darker bars, and asterisk mark if LiDAR was acquired after manual measurements. Hatched
bars indicate plot was snow free. .............................................................................................................. 64

xii

Figure 5.13 Regenerating cardinal plots mean difference between LiDAR and in-situ snow depth
measurements for the 2022 and 2023 surveys. Error corrected using the meteorological or game camera
data is denoted with the darker bars, and asterisk mark if LiDAR was acquired after manual
measurement. Hatched bars indicate plot was snow free. ........................................................................ 65
Figure 5.14 Juvenile cardinal plots mean difference between LiDAR and in-situ snow depth
measurements for the 2022 and 2023 surveys. Error corrected using SR50 snow depth sensor denoted
with the darker bars, and asterisk mark if LiDAR was acquired after manual measurement. Hatched bars
indicate plot was snow free. ...................................................................................................................... 66
Figure 5.15 Boxplots of the mean difference for each plot by cover type, with “open” referring to the
weather station snow courses. For each vegetation type, the left box represents 2022 data and the right
box 2023. For 2022, phase 4 has been excluded due to the noted offset due to delay in LiDAR
acquisition and manual sampling, and subsequent snow melt. The box shows the quartiles of the dataset,
with whiskers extending to points within 1.5 x the inter-quartile range. Observations outside of this
range are displayed as points. ................................................................................................................... 68
Figure 5.16 Scatterplot of the slope angle of the cardinal plots and the mean absolute bias error for the
2022 and 2023 LiDAR surveys. Each marker represents the mean absolute bias error of the cardinal plot,
color coded by cover type, and the vertical bars represent +/- sigma. ...................................................... 69
Figure 5.17 Snow distribution for Steph 1 weather station snow course for each LiDAR survey in 2022,
showing the manual (red) and LiDAR (blue) snow depth distribution. None of the measurements were
significantly different in snow depth variability. The offset between the LiDAR and the manual
sampling in the June 24th acquisition is due to the snow melt which occurred between manual sampling
and LiDAR survey. No in-situ measurements were taken at Steph 1 on April 11, 2022. ........................ 72

xiii

Figure 5.18 KDE plot for 2022 showing the variability of snow depth for old growth plot ST1OG7.
None of the survey dates were significantly different in snow depth variability between the in-situ and
LiDAR....................................................................................................................................................... 75
Figure 5.19 KDE plot for the regenerating forest plot ST4RF1 for 2022. The snow depth distribution
was significantly different across all phases between the in-situ and LiDAR. For the 24th of June survey,
the 33 in situ measurements not visible on the graph recorded zero snow. .............................................. 76
Figure 5.20 Boxplot comparing the mean difference of plot snow depth for the manual observations and
LiDAR across both years, with open cover representing the weather station snow courses. Blue boxes
represent results of the 2017 BCTS model, and the orange boxes are of the 2023 bare earth model. The
box shows the quartiles of the dataset, with whiskers extending to points within 1.5 x the inter-quartile
range. Observations outside of this range are displayed as points. ........................................................... 77
Figure 5.21 Bar graphs to show the comparison of the two bare earth models - 2017 BCTS (hatched)
and 2023 ACO (non-hatched) showing the total ground points per plot, broken down by point density
per meter squared (top) and the proportion of the plot which had valid ground returns (bottom). Bars
represent the average of the cardinal plots by each cover type, and error bars represent standard error . 78
Figure 5.22 Comparison of the total snow storage in the Russell Creek watershed between the 2017
BCTS model, and the 2023 ACO model for each snow survey, broken down by forest cover type. Grey
bars for the 12th June 2023 represent that the snow volumes were negative overall. ............................... 80
Figure 5.23 Snow volume by elevation band for the 2017 BCTS model (red) and the 2023 ACO model
(blue), showing the 2022 (top) and 2023 (bottom) LiDAR surveys ......................................................... 81
Figure 6.1 Site maps for the natural experiment with a) showing the AOI relative to the Tsitika
watershed, northern Vancouver Island, and b) scaled to the AOI, showing the road, control, and
treatment area used within the analysis and c) the extent of Cain creek catchment inclusive of the Lower
Cain (1260 m) and Ridge Run (1330 m) weather stations........................................................................ 87
xiv

Figure 6.2 Overview of the Lower Cain (1260 m) weather station air temperature (red) and snow depth
(blue) measurements for the 2022 and 2023 peak snow and snow melt season. ...................................... 92
Figure 6.3 Histogram of the residuals for the control, treatment, and road sites between the post-harvest
(2023) and pre-harvest (2020) bare earth. A positive skew is indicative of the post-harvest elevation
greater than the pre-harvest. ...................................................................................................................... 94
Figure 6.4 Snow depth distributions of the three snow free acquisitions, ie. no snow, comparing the
control and the treatment site processed with the post-harvest (2023) bare earth model. The 24th of June
2022 (left) shows a negative offset when using the harvested ground model to measure snow in the
treatment site which - as it is prior to harvesting – is forest covered. The 25th of May 2023 (middle) and
12th June 2023 (right) show good agreement between the post-harvest ground surface model and the
harvested treatment and forested control sites .......................................................................................... 95
Figure 6.5 Box plots of the snow depth (m) in the control and treatment sites for the pre and post-harvest
conditions. The box shows the quartiles of the dataset, with whiskers extending to points within 1.5 x the
inter-quartile range. Observations outside of this range are not shown on the graph ............................... 96
Figure 6.6 Total water storage in cubic metres for the control (green) and treatment (orange) sites for the
pre- and post-harvest years. ...................................................................................................................... 97
Figure 7.1 Images from the regenerating stand (A) and dense old growth (B), showing notable difference
in tree branch height in relation to the ground. ....................................................................................... 104
Figure 7.2 Mean difference of the cardinal plots from both 2022 and 2023, for only harvested plots (the
regenerating and juvenile), using three different bare earth models: 2017 BCTS survey acquired in
September 2017 (orange), 2023 ACO survey (blue), and the substitute bare earth model: the May 2023
snow survey which was snow-free in the juvenile and regenerating stands (green). Individual points
represent each plot validation across 2022 and 2023. ............................................................................. 108

xv

Figure 7.3 Imagery of the ground point density (GPD) within the Russell Creek watershed showing the
March 10th acquisition (A1, A2), where snowpack is covering much of the watershed and high snow
accumulation survey, and low snow accumulation May 25th 2023 (B1, B2). As snow melts in the
watershed, the number of no data pixel increases, and the ground point density decreases. ................. 113
Figure 7.4 Distribution of the cardinal (blue) and weather station snow course (orange) mean difference
for 2022 (left) and 2023 (right) processed with the 2023 ACO bare earth model. ................................. 114

xvi

Acknowledgement
This work was supported by the Tula Foundation, Ministry of Forests, and BC Timber Sales, as well as
the BC Graduate Scholarship Program. Thank you to Dr. Brian Menounos and the Airborne Coastal
Observatory (ACO) for providing the LiDAR datasets, as well as Santiago Gonzalez Arriola, Rob
White, and Steve Beffort on the Hakai geospatial team for all your assistance with the datasets. I would
like to thank my co-supervisor Dr. Joseph Shea; you were a driving force in the development of my
coding skillset, keeping me motivated, and helped me with realistic goal and scope setting. Another
special thanks to my co-supervisor Dr. Bill Floyd for taking a chance on hiring a new GIS graduate 5
years ago and for giving me the opportunities to learn and develop my skills which inspired me to
pursue my masters. A huge component to completing this project is the support and field work help from
my fellow lab mates over the last 2.5 years: Anna Kaveney, Warren Olmstead, Sergey Marchenko,
Julien Bodart, and Sara Darychuk. From the countless sled stucks, spectacular ski wipeouts, the “rarely
good” ski down from Steph 1, and the extravagant field meals, it was all the more fun with you all!
Special thanks to Rosie Bisset, both for your LiDAR expertise, and for all the help in less than
glamorous field work in the regenerating forests (“I will never eat another blueberry!”). Another special
thanks to Griffin Fisk for putting up with my chaos in the field and at home. A final thank-you to my
parents, I am extremely grateful for your endless support and enthusiasm, which has shaped who I am
today.

xvii

Chapter 1 : Introduction
Climate change directly affects coastal mountain snowpacks. Warm temperatures reduce snow water
storage, shift snow melt earlier in spring, and decrease summer soil moisture, which can increase the
magnitude of flood and drought events (Schnorbus et al. 2014; López-Moreno et al. 2021; Gillett et al.
2022; Hidalgo-Hidalgo et al. 2022). Further, numerous studies (eg. Swanson et al. 1986; Gelfan et al.
2004; Winkler et al. 2005) have shown that forest management can have an impact on snow at the plot
scale, but scaling these results to the watershed is challenging due to landscape and forest complexity.
Thus, basin wide estimates of snow water storage are critical to understand the impacts of climate
change and forest management on watershed scale hydrology.
Advances in remote sensing technology have improved the ability to map snow depths for entire
watersheds. A Light Detection and Ranging (LiDAR) sensor fixed to an airborne platform can acquire
millions of high accuracy surface elevation points over several hundred square kilometres in a time span
of just several hours (Hojatimalekshah et al. 2021). To derive a snow depth surface, a snow-free LiDAR
digital elevation model (DEM) is subtracted from a snow-on LiDAR DEM. LiDAR performs best in
open, relatively flat terrain where light pulses can easily penetrate to the ground (Hojatimalekshah et al.
2021). Uncertainty increases when using LiDAR to measure snow in forested environments, as trees and
vegetation obstruct the ground and the snow surface, which reduces the number of surface returns
(Hopkinson et al. 2004; Hopkinson et al. 2012; Zheng et al. 2016). Dense ground cover vegetation can
also increase elevation uncertainty by being mis-classified as ground in the bare earth model (Hopkinson
et al. 2005; Gould et al. 2001). Additionally, forest snow accumulation and melt patterns are much more
variable than in open areas, due to complex relationships between the canopy and energetic inputs into
1

the snowpack (Pomeroy and Gray 1995; Varhola, Coops, Bater, et al. 2010). This holds particular
significance to watershed management in regions where a significant proportion of the snowpack is
stored within forests.
As winter temperatures increase due to climate change, lower elevation coastal mountains will see an
increase in winter precipitation falling as rain (Evan and Eisenman 2021). Having accurate data for the
snow water storage in the forest is essential for sustainable water resource consumption, in particular for
managing watersheds with “at-risk” snow (Nolin and Daly 2006). Forest management and natural
disturbance lead to a mix of forest structure from clear open areas to regenerating stands of mixed
species, age and structure (Floyd 2012; Winkler et al. 2021). Furthermore, there is a need to assess how
snow storage varies with forest regeneration to contribute to our understanding of how stands and
watersheds hydrologically recover after disturbance, as well as advance models for large and complex
watersheds.
This research project will focus on the use of airborne LiDAR to measure snow depth in Russell Creek
Experimental watershed (established in 1991), a heavily forested watershed on Northern Vancouver
Island (Hudson and Anderson 2006; Floyd 2010). Beginning in the late 1970s, forest harvest began in
Russell Creek, resulting in a range of stand ages and forest stand composition. The objective of this
research is to assess LiDAR uncertainty in forested watersheds typical of the coastal mountain region of
British Columbia to improve our understanding of catchment scale snow water storage.
With repeat LiDAR acquisitions and plot scale sampling for depth, this research aims to answer the
following questions:
1. How accurate is LiDAR at capturing the spatial variability of the snowpack under different forest
stand structures?

2

2. How does the spatial variability of the snowpack vary among different forest stands, and how
does this variability change through the snow melt season?
3. How does the bare earth model impact LiDAR snow depth accuracy under different forest stand
structures?
This research project will contribute towards current literature on forest snow interactions in a region
with at-risk snow due to climate change and a history of forest management going back to the 1970s. It
will also contribute to development of hydrological recovery assessment methods in coastal British
Columbia. This part of a larger study using LiDAR snow surveys in forested watersheds in British
Columbia, which is the first of its kind in the province.
This thesis is organized as follows:
Chapter 2 provides a comprehensive review of literature related to LiDAR observations of snow depth
in forested catchments and summarizes the current state of knowledge with respect to forest disturbance
hydrology. Chapter 3 identifies the study area and data collection methods used in this research, while
Chapter 4 summarizes the methods of analysis. In Chapter 5, the main results section addresses the
research questions and objectives. Chapter 6 follows with a case study, in which harvesting activity that
occurred over the course of this research is used to quantify the impacts of harvesting on the ground
surface and subsequent DEM, and the change to snow storage. Chapter 7 discusses the key findings and
implications of this research. Chapter 8 provides the overall conclusions and recommendations for
future research.

3

Chapter 2 : Background
Airborne LiDAR has drastically up scaled the speed, size, and quantity of snow depth acquisitions, from
manually collected points or data from discrete weather stations, to millions of measurements at the
catchment scale. Despite this, using LiDAR within heavily forested terrain is challenging both due to
limitations of the LiDAR instrument, as well as the complexity of the snowpack under canopies. Section
2.1 will discuss how LiDAR is used to measure the snowpack, as well as sources of error and
uncertainty. Section 2.2 will follow with an overview of the complex snow distribution within forests.
Section 2.3 will focus on how forest disturbance changes snow distribution and can influence watershed
scale hydrology.

2.1 Snow and Remote Sensing
In-situ methods to measure snow depth and density are generally limited to single point measurements,
either through manual depths with snow probes or tubes, imagery from snow stakes or remote cameras,
or automated instruments, like snow depth sensors and snow pillows associated with weather stations
(Hopkinson et al 2001; Kinar and Pomeroy 2015). To quantify catchment wide snow water storage,
these snow depth and density point measurements can be used to validate physically based hydrological
models (e.g. Pomeroy et al. 2007), however the sparse sampling leads to high uncertainty in estimates.
Snow depth measurements predominantly occur in relatively flat, un-forested and sheltered areas which
a) are recommended for representative meteorological station locations that aim to measure incoming
precipitation and snow storage that in unaffected by wind (World Meterological Organization 2008) and,
b) avoid exposure to avalanches and other hazards endemic to complex terrain. Despite efforts to
establish weather stations in more remote and under-represented regions within British Columbia (Foord
et al. 2014), as well as the extensive Snow Telemetry Network (SNOTEL) network in the U.S. (U.S.
4

Department of Agriculture), point measurements are limited in their ability to describe catchment-wide
snow depth variability, particularly within forested watersheds.
Major advantages to remote sensing methods are the ability to capture spatially continuous snow depth
data, and coverage of areas too remote or complex to access with manual measurements. Snow covered
area has been observed through satellite remote sensing since the 1960s, where the high reflectance in
visible wavelengths, and low reflectance in the SWIR band clearly distinguishes snow from other
surfaces (Hall et al. 2005), however these data provide no information about depth and are limited to
cloud free days. Currently, satellites like Sentinel-1 synthetic aperture radar (SAR) data and ICESat-2
can collectively determine snowpack properties such as the depth, wetness, and water equivalence
(Karbou et al. 2021), though coverage can be limited and both forests and wet snow surfaces complicate
interpretation of data. While recent studies have shown that ultra-high-resolution satellite imagery can
derive high resolution (1-4 m) and decimeter accuracy snow depth mapping (Marti et al. 2016;
Deschamps-Berger et al. 2023), the majority focus on large scale snow measurements in alpine terrain
due to limitations of spatial resolution and complications associated with forest cover (Lievens et al.
2019; Karbou et al. 2021). In contrast, Structure from Motion (SfM), in which a camera is attached to an
unpiloted aircraft (UAV), is capable of centimeter resolution, but generally focused on relatively small
study areas (<5 km2) due to aircraft battery limitations and acquisition time (Harder et al. 2016). Studies
have shown that UAV SfM can accurately map snow depths (RMSE: 8.0 – 13.0 cm) in open
environments, but data are poor in areas of exposed vegetation and forests due to the obscuration of the
snow surface (Harder et al. 2020; Dickinson 2022; Landry et al, 2021).
Airborne LiDAR has been used to measure snow since 2001 and is now a well-established technology to
map snow depth across expansive alpine terrain (Hopkinson et al. 2001; Painter et al. 2017; Currier et al.
2019). Applications of LiDAR snow depth data cover a diverse range, inclusive of evaluation of
hydrological models (Dickerson-Lange, Lutz, Gersonde, et al. 2015; Vionnet et al. 2020), improved
5

knowledge of snow distribution patterns (Deems et al. 2006; Currier and Lundquist 2018; Pelto et al
2019), and quantification of derived snow water equivalence (Tedesco et al. 2014; Floyd et al. 2020;
Geissler et al. 2023). Despite widespread application in alpine environments, there are limited examples
in forested watersheds. This is especially important in the low elevation coastal mountains, where
typically snowpack accumulation is high but alpine terrain limited, thus a considerable proportion of
snow is stored in the forest.

2.1.1 LiDAR systems
LiDAR is a remoting sensing technology which emits light pulses and measures the time of flight for a
given photon to calculate distance from the sensor. A LiDAR system is composed of a laser scanner, a
global navigation satellite system (GNSS), and dedicated inertial measurement unit (IMU). The laser
scans perpendicular to the line of flight and measures the return time of the laser pulse to the unit to
determine a three-dimensional elevation point. To calculate snow depth, LiDAR data are acquired in
both snow-on and snow-free conditions, processed into a separate DEMs, co-registered and then
differenced. The LiDAR scanner can be either terrestrial laser system (TLS) or airborne laser system
(ALS); each hold advantages and disadvantages. TLS provides high density returns- in the order of
thousands of returns per metre squared - therefore suited for producing accurate, high resolution datasets
(< 10 cm mean bias error; Prokop et al. 2008; Hartzell et al. 2015; Maté-González et al. 2022). A TLS is
well suited to describe fine scale forest-canopy snow relationships, as the ground-based system provides
the ideal scan angle to capture vegetation structure, as well not be obscured by the canopy (Disney 2019;
Hojatimalekshah et al. 2021). For example, Revuelto et al (2016) used TLS to quantify the forest
management practice of pruning lower branches reduced canopy area by 36% and increased total snow
depth by 14%. Handheld terrestrial LiDAR sensors have become more widely available, a recent study
by King et al (2022) conducted with an iPhone 12 Pro LiDAR module produced very high accuracy
6

(RMSE ~ 6mm) snow depth measurements. Studies have found snow depth accuracy decreases with
distance and high (>60 degrees) incidence angles, therefore TLS are best suited for short-to-medium
range (<500 m2) areas and recommended several scans per site (Prokop et al. 2008; Revuelto et al. 2014;
Hartzell et al. 2015).
Airborne LiDAR configurations are commonly mounted to aircraft such as helicopters or fixed wing
airplanes, or Uncrewed Aerial Vehicles (UAVs). Plane mounted systems are well suited for larger area
(order of hundred kilometers) campaigns at the watershed scale (e.g. Floyd et al. 2020), and under ideal
circumstances, accurate to the decimetre level at a 1 metre resolution (Deems et al. 2013). The drawback
for a more extensive survey area is a generally a lower density of ground point returns (~ < 10pts/m2)
than helicopter or UAV mount LiDAR, which can be a source of uncertainty in the measurements
(Menounos et al. 2020). UAV-mounted LiDAR systems provide a compromise between the terrestrial
and plane mounted LiDAR systems in terms of point density (~100 pts/m2) and area coverage, and are
thus well-suited to stand level studies (Harder et al. 2020). LiDAR sensor capabilities continue to
expand as technologies rapidly improve. For example, geiger mode LiDAR – a photon counting system
– has become more prominent, and can capture high resolution datasets (on par with UAV LiDAR) with
the coverage of plane mounted systems (Deems et al. 2013; Lin et al. 2022).

2.1.2 Data processing
Imagery - either ortho-imagery collected coincidentally with a LiDAR survey, or satellite - can be an
important part of the processing workflow for identifying snow free and stable terrain in the survey area
to help with co-registration and confirm snow cover when snowpacks are shallow (Deems et al. 2013;
Floyd et al. 2020; Vionnet et al. 2020). Raw LiDAR data are represented as a point cloud, comprised of
single x,y,(horizontal) and z (vertical) points, of which z may represent either a discrete return or the full
waveform (Deems et al. 2013). Discrete returns are the more common approach, as full waveform
7

requires a complex processing methodology, higher computational power, and a greater cost (Anderson
et al. 2016). Full waveform LiDAR is popular for forestry applications as it allows for much more
detailed description of vegetation, including the canopy and ground vegetation characteristics, as
discrete points only separate sufficiently spaced objects into separate returns (Lindberg et al. 2012).
Anderson et al (2016) found bare earth was captured with greater accuracy in vegetated terrain by full
waveform, and conversely in non-vegetated terrain (roads) by discrete returns. Thus, full waveform data
shows great potential in the snow depth mapping application to better distinguish vegetation from
ground, and thereby improve accuracy of the ground surface.
From the raw LiDAR point cloud, algorithms are used to a) filter and remove noise and outlier data
points, and b) for classification to identify the ground surface and other attributes like vegetation and
structures. A digital elevation model (DEM) is rendered from the ground point class to represent either
the snow or snow-free surface. The highest resolution a DEM can support is usually a function of the
ground point density, for example an average ground point density between 1 - 2 points per m2 could
support a 1 m2 gridded surface (Guo et al. 2010).

2.1.3 LiDAR validation
Typically, LiDAR snow depth accuracy is assessed against in-situ measurements, which are
georeferenced with a GNSS system to ensure high accuracy and precision. A common sampling strategy
is by a linear transect with GNSS measurements at the start and end with manual measurements
between, or a radius approach with transects at cardinal directions from a centre GNSS position
(Hopkinson et al. 2012; Floyd et al. 2020; Proulx et al. 2022). Generally, studies compare a one-to-one
manual-to- LiDAR snow depth measurement, which can be a source of geo-positional error, especially
under the canopy where GNSS accuracy is reduced (Murgaš et al. 2018). To address this, Broxton et al
(2019) used canopy height maps in the field combined with GNSS to ensure trees were correctly
8

identified for under canopy snow depth measurements. In-situ measurements are generally taken with
graduated avalanche probe or an automated version known as a magnaprobe – or a Federal snow
sampler, which also measures density. A review of the different manual sampling methods by Proulx et
al (2022) showed a slight positive snow depth bias (average 1.7 cm) between the magnaprobe and the
snow tube, attributed to penetration of soil by the probe, with a greater difference in the forest than open
sites. As snow tube measurements are considerably slower and more strenuous than a snow probe, a
hybrid approach (snow tube and probe) is recommended for in-situ measurements to obtain both depth
and density at a selection of sampling points (Proulx et al. 2022).

2.1.4 LiDAR Error and Uncertainty
Sources of error in airborne LiDAR system can be classified as either systematic, errors which can be
reproduced and potentially can be removed, or random, which reflects the measurement uncertainty.
Sensor based geo-positional errors are a result of poor accuracy within the GPS, INS, IMU systems,
with generally horizontal error greater than vertical (Hodgson and Bresnahan 2004a; Deems et al. 2013).
As well, an increase in laser scan angle, flight height, and laser beam incidence angle tends to
correspond to an increase in LiDAR error (Goulden et al. 2016)
Terrain features that can affect the precision of LiDAR surveys include forest canopy cover, ground
vegetation, and slope (Hopkinson et al. 2012; Zheng et al. 2016). Slope introduces error as a small
horizontal geo-positional error will lead to a large error in the vertical position on steep terrain (Nuth
and Kääb 2011). Several studies have found that the vertical error on slope angles greater than 30
degrees nearly doubles (Hodgson and Bresnahan 2004b; Tinkham et al. 2012).
Total LiDAR ground returns are reduced under the canopy as the trees interact with emitted photons,
which results in data gaps in the ground point cloud. Areas with very high forest density are likely to
have larger areas with no LiDAR ground returns, in which interpolation of pixel values is required.
9

Zheng et al (2016) found that under-sampling under mixed deciduous forest resulted in an
overestimation of snow depth of an average of 10 cm under the canopy due to the interpolation of data
with the adjacent open pixels, which had a higher snow depth than under canopy.
Dense ground vegetation has been shown to cause a systematic underestimation of snow depth as a
result of both the inability of the laser to penetrate through to ground as well as misclassification of
vegetation as ground from the bare earth acquisition (Hopkinson et al. 2005; Gould et al. 2013; Currier
et al. 2019). Gould et al (2013) found that airborne LiDAR ground models in areas with ceanothus, a
dense ground shrub species, had elevation errors (RMSE: 17 to 37 cm) that were attributed to the
inability of LiDAR pulses to penetrate the shrub, and misclassified as ground. An accuracy assessment
of LiDAR by vegetation types conducted by Hodgson and Bresnahan (2004a) showed the highest error
in brush and low trees (RMSE 23.3 cm) and the lowest error in pavement and evergreen forests (RMSE
of 14.9 cm and 12.9 cm). The low error within the evergreen forest is attributed to the minimal ground
vegetation and generally uniform height of forest canopy, resulting in two distinct levels of which
ground can be easily distinguished. The opposite effect is true for the brush and low trees: the multistory
environment, and variable height of the low vegetation, increases the likelihood of misclassifying
ground. Reutebuch et al (2003) found minor difference in mean ground surface error (16 – 31 cm)
between four forest categories (cutblock, heavily thinned, lightly thinned, and uncut), attributed to the
relatively lower flying height (200 m), which results in stronger return signals due to the shorter distance
the laser travels. Thus, the LiDAR errors in snow depth mapping attributed to vegetation may be largely
dependent on the accuracy of the ground surface representation.
Recent studies with UAV LiDAR systems generally show lower snow depth error in forest (RMSE: 9 –
16 cm) than compared to plane mounted systems (RMSE: 5 - 35 cm), likely due to their ability to fly
lower to the ground surface, though generally errors were worse in forest compared to the open
(Mazzotti et al. 2019; Harder et al. 2020; Jacobs et al. 2021; Geissler et al. 2023). Menounos et al (2020)
10

suggested the use of lower flight altitudes (to increase the density of the point cloud) and in-situ
measurements (to correct detected bias) to improve LiDAR snow depth accuracy in forested terrain.
Other studies have highlighted the importance of in-situ snow depth measurements which encompass
forested terrain in addition to the open to better understand snow depth uncertainties (Painter et al.
2016).

2.2 Snow distribution in forests
Snow distribution is governed by energy inputs into the snowpack, as well the interactions of vegetation,
topography and wind. The snow energy balance of the snowpack can be represented as
[2.1]

∆ܳ = ܴ௡௘௧ + ‫ܪ‬ௌ + ‫ܪ‬௅ + +‫ ீܪ‬+ ‫ܪ‬ோ

where ∆ܳ is change in snowpack energy, and the energetic fluxes are: net radiation (ܴ௡௘௧ ), sensible heat
(‫ܪ‬ௌ ), latent heat (‫ܪ‬௅ ), ground heat (‫) ீܪ‬, and advective heat transfer from precipitation (‫ܪ‬ோ ). The
contributions of net radiation to the snow energy balance are five to ten times greater than all other
components, though the relative inputs of the short and longwave radiation vary considerably (Marks
and Dozier 1992). Net radiation is the sum of the incoming and outgoing radiant fluxes, and can be
expressed as
[2.2]

ܴ௡௘௧ = (1 − ߙ)ܴ௦ + ܴ௟௜ + ܴ௟௢

Where α is albedo, Rs is direct and diffuse incoming radiation, Rli is the incoming longwave, and Rlo is
the outgoing longwave radiation. Longwave radiation is emitted by the ground surface and objects, and
by the atmosphere, and is generally higher: 1) under cloudy skies as opposed to clear, 2) in higher air
temperatures and humidity, and 3) within a forest vs in a clearing (Dingman 2015). Shortwave radiation
is emitted directly from the sun and is a function of the declination of the sun, latitude, and time of year
(Dingman 2015). Albedo measures the proportion of short-wave radiation which is reflected; fresh snow
11

with a small grain size has a higher albedo value (0.9) than older and wetter snow (0.5) (Jones et al.
2001). As ground and vegetation are exposed, albedo will decrease and longwave and conductive energy
from the ground and vegetation will generally increase in snow melt rate (Pomeroy et al. 2003).
Snow distribution under forests differs from in the open due to the interactions with the canopy altering
snow redistribution and components of the snow energy balance. Generally, snow will accumulate
slower and reach shallower depths within forests, due to first, snow interception, and then a combination
of sublimation, meltwater drip, evaporation and mass unloading from the canopy. Mechanics for snow
loss to the canopy vary across climatic regions: in cold and dry climates, up to 30 to 45% of annual
snowfall may be lost to canopy sublimation (Pomeroy and Gray 1995), versus in maritime climates
where sublimation generally only account for approximately 10% of loss, with meltwater drip and mass
unloading from the canopy are the dominant processes (Storck et al. 2002). Methods to assess how
intercepted snow is removed from the canopy have included:
-

Individual branch experiments: warmer branches will bend easier, therefore snow will unload
more rapidly at warmer temperatures, or in maritime climates (Schmidt and Pomeroy 1990)

-

Cut-Tree experiments: Storck et al (2002) found that the variability of snow interception was
minimal between the species compared (douglas fir, white fir, and ponderosa pine) in a maritime
environment, and interception variability was mostly a function of the snow event itself. A novel
approach using TLS by Russell et al. (2020) shows potential, with initial results showing good
agreement between the TLS and cut-tree measured snow mass.

-

Time lapse photography: Floyd and Weiler (2008) captured mass release events from the canopy,
and observed that snow release occurred over a broad area under the tree, and snow depth
remained the same in areas which snow unloaded, but snow density increased.

Trees can either enhance or reduce ablation rates of the snowpack, but generally melt rates are
slower under forests than in the open due to the sheltering effects of the trees from short-wave radiation
12

and turbulent energy fluxes associated with wind (Hedstrom and Pomeroy 1998; Varhola, Coops,
Bater, et al. 2010). Recent studies conducted in maritime climates, which typically have deeper
snowpacks, warmer temperatures, and larger trees, have shown that in many cases the pattern of snow
melt under forest and in the open is reversed. In maritime climates, an increase in long-wave radiation
from the trees can dominate the energy balance, which in turn accelerates the melt rate within the
forest, and results in longer snowpack duration in adjacent open sites (Lundquist et al. 2013;
Dickerson-Lange, Lutz, Martin, et al. 2015; Roth and Nolin 2016). This pattern is more pronounced
below tree-line, in warm dense forest stands sheltered from the wind (Dickerson-Lange, Lutz, Martin,
et al. 2015). In colder and more arid climates, trees alter the energy balance by reducing incoming
shortwave radiation which reduces melt beneath forest canopies.
Wind also impacts snow distributions through a) transportation of snow from exposed to sheltered areas,
b) turbulent fluxes of energy into the snowpack causing melt or c) sublimation of the snow into the
atmosphere, though this tends to be a small component in coastal watersheds (Currier and Lundquist
2018). Blowing snow is more prevalent in cold climates due to a lower wind speed threshold required to
transport cold and low density snow, than the warm, and higher density snow typical of the coastal
climate (Li and Pomeroy 1997). A study by Currier and Lundquist (2018) showed significant differences
to snow depth between the windward and leeward side of a forest in a high wind cold climate site in
Colorado, and no significant different in a maritime location with equal or greater wind speeds. High
wind speed sites (> 8 – 17 m/s) in a study in the Pacific Northwest showed accumulation rates were
equal between forest and open due to the depositional differences from the wind transported snow
(Dickerson-Lange et al. 2017).
Large rain-on-snow (ROS) events occur with greater frequency in temperate climates and can rapidly
change the snowpack and its distribution (Marks et al. 1998; López-Moreno et al. 2021; Gillett et al.
2022). When the snowpack is near isothermal, there is a low energy requirement to begin snowmelt
13

during a ROS event, which can rapidly result in flood and landslide events (Harr 1981). Melt rates
during ROS events have been found to be significantly higher in openings and clear cuts than in forests
(Berris and Harr 1987; Marks et al. 1998). The major difference in the energy balance between forested
and un-forested sites is the turbulent fluxes of sensible and latent heat, as the trees shelter the surface of
the snow from wind. Slope and aspect also affect melt rates, as areas positioned towards the incoming
wind will be subject to greater turbulent transfer and increased melt rate (Harr 1981). In areas where
snowpacks are generally shallow, such as the transient snow zones, the snowpack can melt entirely
during a single large rain-on snow event (Berris and Harr 1987). Musselman et al (2018) investigated
the impacts of climate change to ROS events, and in a maritime environment frequency at lower
elevations will reduce due to declining snowpack and increase at mid elevations, as well as the major
ROS events shifting earlier in the season. While studies have shown these differences at the stand level
during rain-on-snow, linking these stand level effects to watershed scale stream flows has been
challenging due a sparsity of data.

2.3 Forest Disturbance and Hydrological Recovery
Forest disturbance can cause major shifts to the snow energy balance, by both a) changes to
accumulation resulting from loss of canopy interception and b) changes to ablation from loss of canopy
cover altering the long and shortwave radiation inputs and turbulent fluxes (Boon 2009; Marks et al.
1998; Varhola, Coops, Weiler, et al. 2010; Winkler et al. 2015). Numerous studies have been conducted
to assess how forest and open site melt rates differ as it has clear implications for forest management
practices (eg. Hubbart et al. 2006). Gelfan et al (2004) found melt rate in the open to be 3.2mm day-1
faster than in the forest. Increased snow accumulation mostly occurs in the openings, but even with a
higher antecedent snowpack, averaged snow melt duration were longer in the forest (30 days) than in the
open (22 days) (Gelfan et al 2004).
14

The size and aspect of the harvested block affects the snow retention, which is an important
consideration for water resource management. Swanson et al. (1986) found that cut blocks that were 1-2
tree height in width, and 60 – 80 meters in length, with the long side perpendicular to the prevailing
wind, were ideal for snow retention. Clear cuts which are south facing and situated on level terrain have
shown an increase in melt rate than north aspect (Murray and Buttle 2003; Pomeroy et al. 2012). Forest
thinning also impact snowmelt timing; Ellis et al (2013) found that on southern aspects melt rate
increased due to greater shortwave radiation, and on northern aspects decreased, due to loss of longwave
radiation from the trees.
Stand level hydrological recovery describes to what extent a regenerating stand is behaving as an old
growth or mature forest, and may be evaluated on the following metrics: rainfall interception, snow
melt, snow accumulation and interception, ROS, and evapotranspiration (ET), although in coastal
watersheds ET recovery is not deemed to be a major component (Hudson and Horel 2007). Generally in
British Columbia, hydrological recovery focuses on snow accumulation and melt, and is estimated via
canopy characteristics, such as height, density and species ( Hudson and Horel 2007; Winkler and Boon
2017). Floyd (2012) investigated hydrological recovery through ROS events within the coastal Russell
Creek watershed (the location of this study) and found recovery to be variable; low elevation
(approximately 500 m) and sheltered forest plots with a canopy height of 13 m were 100% recovered in
relation to energy inputs during ROS, but recovery of forests of a similar height at higher elevation was
lower likely due to higher exposure to turbulent fluxes.
Since the 1990’s, equivalent clear cut area (ECA) has been the primary metric within British Columbia
to quantify the extent of disturbance, and is defined as the area that functions hydrologically as a clear
cut in a disturbed area (Winkler and Boon 2017). To measure ECA, recovery curves as described above
are used to assess each harvested area and its regenerating forests. For example, an area which is
considered permanently disturbed, such as roads or powerlines, has an ECA of 100%, while a forest
15

15m in height may be considered 50% recovered. The next step is to calculate the area disturbed at the
watershed level based on the ECA and use that resulting metric to assess the potential for current and
additional harvest to negatively impact stream flow. A major shortcoming of the use of ECA is it is a
highly simplified approach based on empirical relationships from a handful of studies, which fails to
consider the impact of important watershed characteristics such as aspect, size, and elevation gradient
(Winkler and Boon 2017). Despite these limitations, ECA along with a suite of other indicators is widely
used to assess overall risk for harvesting to impacts peak flows and downstream values.
There is a long history of using paired watershed studies to assess the impacts of forest harvest on
stream flows (Richardson et al. 2023). Here, a control and treatment watershed which are of similar size
and characteristics, are monitored pre and post forest harvesting to determine changes to stream flow,
including water yield, low flows, and peak flows. The long-standing nature of many of the experimental
paired watersheds (mean duration: 27 years), as well as their spatial distribution across North America
(34 operational as of 2023) provides both a ) a method to assess relative change due to disturbance
which is not obscured by the effects of climate change, and b) comparison of watershed responses
between different hydrological regimes and c) insights into the effects of forest regeneration on stream
flow (Richardson et al. 2023). For example, Perry and Jones (2017) evaluated the impact of forest
plantation on summer streamflow across 8 paired watersheds in the Pacific Northwest over six decades.
They found that relatively young stands (25 – 45 years) of Douglas fir within the treatment watersheds
resulted in summer stream flow deficits of 50% compared to the control basin, which was comprised of
150 – 500 year old forests. Coble et al (2020) found these effects can reverse at larger scales, when a
broad range of regenerating forest are represented in the catchment. Despite widespread usage,
criticisms and limitations of the paired catchment approach include:1) challenges to find representative
paired watersheds, as well as basins which have pre-disturbance stream flow records, 2) the uncertainty
of scaling for small basin to large, 3) problems with data quality and accessibility between watersheds,
16

and 4) the cost of maintaining long-term observations (Pike et al. 2010; Neary 2016; Richardson et al.
2023).
In theory, a greater area of disturbance results in greater potential for hydrologic change in terms of
snow melt rate, peak SWE timing, and overall water yield. However, there is disagreement within the
forest hydrology community on the effects of forest harvesting on the magnitude of the hydrological
response, and whether disturbance impacts only average flows, or the frequency of peak flows (Alila et
al. 2009; Bathurst et al. 2020; Pham and Alila 2024). Disagreement largely stems from a) differences in
disturbance distributions within basins and hydro-climatic regimes, and b) the analysis methods which
have been employed to assess change, specifically flood frequency vs chronological pairing (Alila et al.
2009). Some studies - both paired catchments and through other methods - have shown watersheds that
are more than 30% logged have increased snow melt rates and peak flows, as well as peak flows that
occur earlier in the season and low flows that occur later (Bosch and Hewlett 1982; Moore and Scott
2005; Winkler et al. 2015; Winkler et al. 2021). However, a recent study by Johnson and Alia (2023),
which employed a flood frequency approach, found a lower threshold, in which 21% harvested
watershed resulted in flood events of all magnitudes increasing in frequency, with 7, 20, 50, and 100
year flood events becoming 2, 4, 6 and 10 times more frequent respectively.
Much of the current understanding of forest disturbance has been gained through plot scale data, paired
watersheds, or hydrological models which have utilized the paired watershed and plot data (Alila et al.
2009). Compounded with forest harvesting are the effects of climate change, which can result in snow
melt occurring earlier, and can increase risk of early melt events, and summer droughts (Winkler et al.
2015). This highlights the need for robust snow water measurements at the catchment scale to improve
our understanding of watershed processes and forest management practices.

17

2.4 Research Gaps and Objectives
Based on the background review presented above, there remain some key gaps in our understanding of
the ability of LiDAR to evaluate snow distributions in forested catchments where logging has occurred.
These gaps, and the research questions identified in Chapter 1, will be addressed in this thesis through
the following objectives:
1) Establish a manual sampling network to assess accuracy of snow depth derived from LiDAR
samples and calculate error metrics at the plot scale (Chapter 4)
2) Characterize forests and vegetation attributes at manual sampling network, and ground return
point density in different forest cover types, to determine if these contribute to LiDAR error
(Chapter 5)
3) Compare the plot results using an alternative bare earth DEM with different ground sampling
point densities and acquisition times (Chapter 5)
4) Assess the spatial distribution of snow at the watershed scale, and the differences between the
two bare earth models (Chapter 5)
5) Evaluate the impact of harvesting activities on bare earth models and snow depth distributions in
a control-treatment forest harvest experiment (Chapter 6)

18

Chapter 3 : Study Area and Data
3.1 Study Area
This study is situated within the Russell Creek Experimental Watershed, a sub-basin of the Tsitika River
watershed (397 km2), located on Northern Vancouver Island (Figure 3.1), British Columbia, Canada.
Access to the watershed is via a network of resource roads. The watershed encompasses an area of
approximately 32 km2, with predominantly north to northwest aspects, and spans an elevation range
from 220 to 1670 m. The area is 85% forested and with limited alpine terrain. Three biogeo-climatic
zones are within the watershed: coastal mountain-heather alpine, coastal western hemlock, and mountain
hemlock (Meidinger and Pojar 1991). The coastal mountain-heather alpine zone is characterized by most
vegetation being low growing white and pink mountain-heathers, and treeline tree species are
predominantly mountain hemlock (Tsuga mertensiana). It encompasses only 12% of the watershed and
is generally above 1500 m. The mountain hemlock zone occurs at the approximately 900 to 1400 metre
range, with extensive parkland from 1100 to 1400 depending on aspect, with Mountain Hemlock,
amabilis fir (Abies amalabis) and yellow cedar (Chamaecyparis nootkatensis) being the dominant tree
species. The largest zone is the coastal western hemlock, which is approximately 65% of the basin,
situated in the low to mid elevation range and bordered by the mountain hemlock zone. Dominant tree
species are western hemlock (Tsuga heterophylla), western red cedar (Thuja plicata), and large areas of
regenerating stands of douglas fir (Pseudotsuga menziesii). The watershed is characterized by high
precipitation (approx. 2300 mm per year), rain-on-snow events, and deep winter snowpacks from 3 to 5+
m in the alpine (Floyd 2010). Average temperatures are mild, with annual average temperatures of 7.1,
5.2 and 2.9 degrees Celsius at 400, 840, and 1500 m respectively (Floyd 2010). It is typical for seasonal
snowpacks to develop above 700 m, with transient snowpacks below that elevation, but this can vary
19

from year to year. Rain-on-snow events are typical in this watershed and are associated with the highest
peak flows on record.

Figure 3.1 Overview map of the study site: Russell Creek (32 km2) sub-basin of the Tsitika watershed (397 km2), located Northern
Vancouver Island, British Columbia

Active logging has occurred in the Russell Creek Watershed since 1978, which presents the opportunity
to study a spectrum of regenerating and old growth forest stands within the same hydrological zone. The
regenerating forest at lower elevations is a mix of planted western hemlock, western red cedar and
douglas fir (Pseudotsuga menziesi), while at higher elevation regenerating forests western hemlock and
douglas fir are more common. As of 2023, approximately 40% of the watershed has been harvested, all
of which occurs below 1300 metres.

3.2 Weather Station Network
Russell Creek Experimental Watershed was established in 1991, originally to address concerns that
timber harvesting was causing ecological harm to a nearby rubbing beach for Orcas (Floyd 2010). Since
20

2009, ten weather stations have been operated by the Ministry of Forests Coast Area Research, though
the construction and data acquisition of several of the stations precede this. The weather stations are
located at an elevation gradient from 280 to 1516 meters and located to capture the elevation gradient of
the watershed and in areas with different exposure to energy fluxes. Each station measures at a
minimum rain via a tipping bucket, temperature and relative humidity, while most measure wind speed
and direction, solar radiation and snow depth, with three stations located at 493, 820 and 1164 m
measuring total precipitation in a standpipe. Data can be found at http://graph.viu-hydromet-wx.ca/.
Meteorological data from these stations, specifically snow depth and air temperature, were used to
describe the seasonal snowpack, as well as assess change when manual measurements and the LiDAR
acquisition occurred.

3.3 Land Cover Classification
To assess whether forest cover had an impact on LiDAR error, broad land cover types based on age and
height were created in Russell Creek watershed (Figure 3.2). Roads within the watershed, both the
highway and resource roads, were digitized from the bare earth model, and alpine terrain was classified
by a combination of the BC bio geo-climatic zones, and ortho-imagery. The remaining forested portion
of the watershed was first categorized by whether the forest was previously harvested, or not (old
growth). Of the harvested forest, age was a relatively good proxy to distinguish between the different
stand characteristics, and was broken into four categories: mature regenerating, regenerating, juvenile,
and new harvest (Table 3.1; Figure 3.2). Figure 3.3 provides images of the typical cover types in the
watershed.

21

Table 3.1 Breakdown of the land cover categories within the Russell Creek watershed, which were used to: a) determine plot locations and
b) for watershed scale analyses.

Cover Type
Alpine
Old Growth
Mature Regenerating
Regenerating
Juvenile
New Harvest
Roads

Canopy Height
Canopy Density
(Median, m)
[%]
1.1
35%
21.1
84%
14.4
95%
7.9
74%
2.7
25%
0.4
0%
2.8
33%

Harvest year

1986 or before
1987 - 2005
2006 - 2017
2018 - 2023

Area
[% of watershed]
12
45
16
21.8
5
1.1
0.7

Figure 3.2 Overview of Russell Creek watershed showing A) the landcover type, B) the landcover broken down by elevation, C) the
orthoimagery of the watershed, and D) aspect

22

Figure 3.3 Images of the forest cover types where plots were sampled within this research project: Juvenile (A), Regenerating (B), Dense
Old Growth (C), Old Growth (D) and Alpine (E)

3.4 LiDAR surveys
Airborne LiDAR data used in this study was collected by the Airborne Coastal Observatory (ACO), a
partnership between the Hakai Institute, the University of Northern British Columbia, and Kisik Aerial
Survey. As part of a greater project to quantify seasonal snow water storage, LiDAR surveys were flown
over portions of four key watersheds in southern coastal British Columbia: Tsitika, Cruickshank,
Englishman and Seymour. Five flights were targeted between March and June to capture peak snow
accumulation and the subsequent melt for each year between 2020 and 2023. The exact timing of the
flights is subject to weather and equipment availability, as well as the snow conditions. The focus of this
23

study is the 2022 and 2023 snow years within the Tsitika watershed; an approximately 80 km2 LiDAR
survey area (Figure 3.4). The Hakai geospatial team processed the raw LiDAR data and produced a
cleaned point cloud, ortho-imagery, and derived raster products (Table 3.3).

3.4.1 Snow Free acquisition
The LiDAR bare earth DEM for the Tsitika watershed was acquired in 2023. The motivation for an
updated DEM was to have one acquisition for the full watershed, acquired with the same LiDAR system
used to create the snow surface DEMs. The prior bare earth DEM was a mosaic of two separate
acquisitions: the majority of the area was a BC Timber Sales LiDAR acquisition from 2017, and a small
area was flown by the ACO in 2020 to cover the southern Mount Cain region (Figure 3.4). This
provided an opportunity to compare two different bare earth DEMs created from separate acquisitions
and subsequent impacts on resulting calculated snow depths and volumes. The updated LiDAR bare
earth survey was captured by the ACO in August 2023 (see Appendix A: ACO 2023 Bare Earth Model),
and is the primary DEM used within this research project.

24

Figure 3.4 LiDAR survey area within the Tsitika watershed, showing the ACO 2023 bare earth boundary and snow surveys(red), the 2017
BCTS bare earth coverage of the AOI (orange), and the ACO 2020 acquisition which was targeted to fill the missing area over Mount Cain
(blue). Note that the linear boundary to the SE was used to cover the highway and provide generally snow free stable terrain

25

Table 3.2 LiDAR products and metadata for the Tsitika survey area from 2020 to 2023 used within this research project

LiDAR product

Format

Acquisition date

Source

Raster
Resolution
(m)

Average ground/snow
surface point density
(pt/m2)

Bare Earth:
2017 BCTS

raster (.tiff);
point cloud (.laz)

September 3 2017

BC Timber
Sales

1

n/a

Bare Earth:
2020 ACO (Mount
Cain)

Raster(.tiff);
Point cloud (.laz)

September 29 2020

ACO

1

n/a

Bare Earth:
2023 ACO

Raster(.tiff);
Point cloud (.laz)

August 19 2023

ACO

1

1.19

Snow surveys: 2022

raster (.tiff);
Point cloud (.laz)

Phase 1: March 9
Phase 2: April 11
Phase 3: May 21
Phase 4: June 24
Phase 5: cancelled

ACO

1

1.24
1.12
1.09
0.86

Snow Surveys: 2023

raster (.tiff);
Point cloud (.laz)

Phase 1: March 10
Phase 2: April 4
Phase 3: April 29
Phase 4: May 25
Phase 5: June 12

ACO

1

2.88
n/a
1.08
1.01
0.97

Ortho Imagery

Georeferenced
imagery(.tiff)

Produced for each
snow survey
acquisition

ACO

10 – 20 cm

-

Canopy Density

raster (.tiff)

Produced for each
bare earth

ACO

1

-

Canopy Height

raster (.tiff)

Produced for each
bare earth

ACO

1

-

Canopy Cover

raster (.tiff)

Produced for each
bare earth

ACO

1

-

Ground returns

raster (.tiff)

Produced for each
bare earth

ACO

1

-

All
returns

raster (.tiff)

Produced for each
snow survey
acquisition

ACO

1

-

26

3.5 In-Situ Observations
Field validation measurements occurred at cardinal plots and weather station snow courses within an
approximately 4 km2 area in the southern portion of the Russell Creek watershed (Figure 3.5)

Figure 3.5 Location of the game camera sites (blue), weather stations (purple), and cardinal plots (red) within the Russell Creek watershed

3.5.1 Snow Sampling
Thirteen cardinal plots were located at elevations ranging from 809 m to 1505 m. Sites were selected to
incorporate the different forest cover types captured in the watershed (Figure 3.5), while also allowing
for summer and winter access. Each plot had a GNSS-sampled plot centre, with cardinal and intercardinal arms extending 10 m (Figure 3.6). Every 2.5 m, a snow depth measurement was taken, with a

27

total of 33 snow depth measurements per plot. A federal snow sampler was used to measure depth, snow
water equivalent and density at the plots enter and at the 10 m of each cardinal arm (n=5).

Figure 3.6 Cardinal plot snow sampling strategy of density (open circles) and snow depth (dash) measurements. Distances (red) are in
meters

Manual snow depth measurements were recorded using an iPad via the Device Magic field form
application (Device Magic, 2022). For each probed snow depth, canopy cover was observed (open, edge
or under canopy, tree well) as well as characteristics of the ground hit (wood, rock, or other notes). Plot
naming conventions started with the prefix of the next highest weather station: for example, ST1 is
located at an elevation below Steph 1 and above Steph 2 weather station. The suffix denotes the cover
type of the plot (A: alpine, OG: old growth, CC: clearcut, RF: regenerating forest), and numbered from
highest to lowest elevation. For the plot-scale study, old growth with denser canopy closure (>70% of
plot) was further sub-categorized as dense old growth. Slope and aspect metrics were derived from the
bare earth DEM.

28

Table 3.3 Meta data for cardinal plots used within the in-situ sampling campaign. Precise plot locations and GNSS accuracies for the
2022 and 2023 sampling campaigns can be found in Appendix C

Cover
(harvest year)

Alpine

Old
Growth

Dense
Old
Growth

Plot ID

Average Canopy
height (m)

Canopy
Cover
(%)

Average
Slope
(degrees)

Elevation
(m, ellipsoid)

Aspect

ST1A1

2

6

18

1505

E

ST1A2

2

27

10

1507

SW

ST1OG3

8

30

33

1379

W

ST1OG7

15

40

20

1254

N

ST1OG8

18

41

25

1228

NE

ST1OG4

18

84

30

1331

W

ST1OG5

15

81

26

1311

W

ST1OG6

15

72

18

1296

W

ST4RF1

6

76

30

832

NE

ST4RF2

5

53

20

830

NE

ST6GC1

7

88

18

809

W

ST2CC1

2

9

25

1007

N

ST2CC2

2

2

12

1076

NE

Regenerating
(1988)

Juvenile
(2008)

Accurate geo-referencing of the snow validation site locations is essential to compare in-situ snow depth
measurements to LiDAR depths. Two GNSS units were used to establish the plot position: one base
head set up on a previously established and known point, and the second rover point at the plot centre
(Figure 3.7). Precise position and accuracy of the GNSS plot centres is provided in Appendix C. The
horizontal position of the plot location was of higher importance than the vertical, as purpose was to
29

geo-position the probed snow depths, rather than the elevation of the snow surface. Generally, plots
were surveyed only the first validation trip of the season, as the plot centres remained in place to mark
location for the future sampling campaigns.
Additional GNSS data points were collected in 2023 snow season as a means to improve the process of
LiDAR survey alignment. This included elevation transects of the snow surface, as well as stable terrain
points such as a cabin roof. High accuracy data points in the x,y, and z planes were required, therefore
only FIXED solution points (corrections from the base received; 1 – 3 cm accuracy) were collected.
Technicians also ensured best practices were followed, for example careful measurement of snow
surface to rover, as well as rover was level and stable when data point collected.

Figure 3.7 GNSS set-up showing the base station (A) established on a permanent ground control point, and the rover (B) set up at the plot
centre of Juvenile site ST2CC2

30

As part of long-term snow data collection in the watershed, manual snow survey measurements using a
federal snow sampler were taken at seven weather stations (Steph 1 through Steph 8). For each of the
stations, 10 measurements are taken at intervals of approximately 1 m, with GNSS positions for the start
and end of the transect. With the exception of Steph 1, all weather station snow courses are conducted
on relatively flat (less than 10 degree slope) gravel roads. Steph 1 weather station is in the alpine, where
the ground is relatively flat and is either rock, moss or low-lying heather. The weather station sites
provide a dataset free from steep slopes, canopy cover and ground vegetation, which are known to
increase LiDAR uncertainty (see sec 2.1.3).
In addition to the manual snow depth plots, nine camera sites captured automated continuous data
throughout the snow season. Each site was equipped with a wildlife camera, and five PVC snow stakes
with associated GNSS position, either 2 or 3 m in length with increments of 20 cm (Figure 3.8).
Cameras were programmed to take two images per day, which provided a continuous dataset for the
snow accumulation and melt period, as well as the snow disappearance date. These data were used to
measure change between surveys.

31

Table 3.4 Game camera sites and metrics for the 2022 and 2023 LiDAR campaigns. Precise position of the game camera stakes can be
found in Appendix C

Cover
(harvest year)

Old
Growth

Regenerating
(1988)
Mature
Regenerating
(1986)

Plot ID

Average
Canopy height
(m)

Canopy
Cover
(%)

Average
Slope
(degrees)

Elevation
(ellipsoid)

Aspect

ST1_GC1

18

95

17

1290

W

ST1_GC2

18

93

15

1268

W

ST1_GC3

22

93

16

1211

N

ST2_GC1

25

97

21

1094

NW

ST6_GC2

41

99

18

787

W

ST6_GC3

36

98

24

771

NW

ST6_GC1

7

88

18

809

W

ST7_GC1

22

100

19

606

NW

ST7_GC2

15

100

23

584

NW

Figure 3.8 Game camera snow stake set up – this imagery is of game camera site ST2GC1 as seen in snow free and snow on conditions

32

3.5.2 Sampling error and mitigation
A known error in manual snow depth measurements can occur by over-probing into ground, with
increased likelihood in unfrozen soils, complex terrain or areas with dense leaf litter. Though
challenging to assess, other research has found probing error to be relatively small and infrequent, with
an average error sub-10 cm, occurring in less than 10% of measurements (Hiemstra et al. 2020). Depth
measurements when using a snow sampler will have less “over-probe” error as a soil plug would be in
the tube and accounted for. For the probed measurements, technicians minimized this error by retaking a
snow depth measurement if uncertain of ground strike (dirt on the probe end would be indicative of
this).
Minor positional errors are likely to occur within the sampling plot, as a result of:
1. Plot cord: Variation in the cord tension between each transects, bends in the cord when
navigating through trees within the plot, as well as the slope of the plot can also contribute to
positional errors.
2. Compass: Though care is taken by technicians, compass error can result in slight offsets in
direction. Even with correct use, situational choices, for example which side of tree to measure,
are a factor in error.
3. Movement of centre position marker: Snow creep can move or break the centre position post.
Nearest trees were flagged with a bearing and distance to centre point to ensure plot center is
returned to the same position if moved.
Since we compared the plot averages to the LiDAR averages over the same area i.e. not point to point,
these errors were expected to be minor.

33

3.5.3 Vegetation surveys
Low lying vegetation can be a source of error when classifying bare earth by reducing the number of
points that reach the ground, requiring interpolation over larger distances, and when especially dense,
can be misclassified as ground which can result in a systematic bias. To account for this error the
cardinal plots were surveyed in snow free conditions to evaluate the height and cover of ground
vegetation in the different forest types. These data were used to evaluate whether a systematic bias in
snow depth measurements was a result of ground vegetation.
A point sampling transect method was used to evaluate ground cover to mirror the cardinal snow plot
methodology. At every 2.5 m interval on the cardinal and inter-cardinal directions, a 50 cm radius area
was sampled. The following metrics were gathered: total number of species present, height of each
species, and classification of species. Broad categories of species identification were used as the focus
was to characterize the height of ground cover, not conduct a full species inventory. An avalanche snow
probe was used to measure vegetation up to 3.2 m, exceeding that was noted as greater than 3.2 metres.
The vegetation categories were broken down into major groupings: Ground cover, heather, shrub, and
woody debris. Ground cover was composed of mosses and lichens, which were less than 5 cm in height,
and woody debris encompassed logs, stumps, and slash piles. The cardinal plots were grouped by cover
type, and the average height and variability each vegetation category were calculated.

34

Chapter 4 : Data Processing and Analyses
The following chapter describes the processing methods of each dataset, including LiDAR, manual
observations, and weather station data. As this study is part of a larger study, some stages of the
processing workflows were completed by other members of the team.
To assess snow depth variability and LiDAR accuracy within the different forest cover types, the
following steps were completed and described in detail below:
i)

Analysis of weather and snowpack conditions for 2022 and 2023 snow years

ii)

Summary of vegetation height and broad species overview for each cardinal plot cover type

iii)

Calculation of LiDAR snow surface point density for each forest cover type and changes
through the snow melt season

iv)

Comparison of the average and variability of snow depths between the manual and LiDAR
observations

v)

Bare earth error analysis: Comparison of bias error when LiDAR snow depth is processed
with alternate bare earth model (ACO 2023 vs BCTS 2017)

The final portion of analysis focused on watershed scale snow distribution using the LiDAR derived
snow depth layer. Total snow storage in Russell Creek was calculated with both the 2017 and 2023 bare
earth models, for each LiDAR acquisition, to assess how the bias and error analysis affected the final
volumes.

35

4.1 LiDAR processing
The initial processing of the LiDAR data was completed by the Hakai Geospatial team. Processing steps
closely follow the methodology described by Deems et al. (2013):
I.
II.
III.
IV.
V.

Flight strip alignment
Ground classification and filtering of noise
Integration of GNSS ground control points
A manual co-registration using GCP’s and pseudo-points within the LiDAR
Generation of snow surface raster, and additional raster products

Ground classification of the point cloud is an important step in the processing, as the choice of algorithm
can be a significant factor in the elevation accuracy (Tinkham et al. 2011). All ACO LiDAR snow-on
acquisitions and the 2023 bare earth model were classified with LASground (Isenberg 2014), which is
based on Axellson’s (2000) progressive Triangulated Irregular Network (TIN) densification algorithm.
First, a subset of low points are selected and used to form a coarse TIN surface. The algorithm then
iterates through the points, and either adds to the TIN ground surface or filters as non-ground based on
whether they fall within assigned threshold parameters. Specific LAStools processing parameters used
are found for ground classification, noise filtering, and generation of the rasters are found in Appendix
B: LiDAR Processing. The 2017 BCTS - which was merged with the 2020 ACO acquisition over Mount
Cain - was re-classified by the Hakai geo-spatial team with the proprietary Terrasolid software product
TerraScan, which also uses a TIN densification algorithm. All raster products were then aligned to the
bare earth model using GDAL warp, to ensure consistency between the processing extents. LiDAR
datasets were primarily processed in python, with some supplementary workflows completed with the
GIS software ArcGIS Pro and QGIS.

36

4.1.1 Non-gap Filled Rasters
The creation of LiDAR bare earth and snow surface rasters requires data interpolation between pixels
when there are gaps in the ground returns. Where there are no ground or snow surface returns in a pixel,
a value can be derived through a nearest neighbour algorithm (Deems et al. 2013).
To eliminate the smoothing effect from interpolation between pixels only the “true” LiDAR derived
pixels in the plot-scale analyses were used in the comparisons between in-situ observations and LiDARderived snow depths. Non-gap filled (NGF; zero point return pixels removed) rasters were generated for
each snow survey (Figure 4.1). The interpolated, or gap filled, snow surfaces were utilized in the
watershed snow volume analyses, where the priority shifted to equal area (ie. number of pixels) between
surveys, and a viable operational workflow.

Figure 4.1 Hillshade of snow-on surface for the Russell Creek watershed, showing the A) interpolated (gap-filled) raster surface and B) the
non-gap filled (pixels with zero ground returns removed, and shown in red) for March 10th 2023 LiDAR survey. For both years, the area
of no ground point returns increased through the snow melt season.

37

4.1.2 Snow Depth Rasters
Co-registered snow depth rasters were produced by subtracting the bare earth DEM from the snow
surface DEM, then converted to cm of snow (Eq. 1),
(1)

ܵ݊‫ݐ݌݁ܦ ݓ݋‬ℎ [ܿ݉] = (Z snow ) − (Z bare ) ‫ ݔ‬100

where Z snow is the surface elevations with snow cover, and Z bare is the surface elevation without
snow cover. This was repeated for the gap filled and non-gap filled snow surfaces, for both the primary
2023 ACO bare earth acquisition, and the 2017 BCTS.
4.1.3 Bias Assessment
Snow-free and stable terrain within the LiDAR boundaries were isolated to assess each snow survey for
systematic error. A section of the highway within the LiDAR survey boundary was ideal for this
analysis, as it generally remained snow-free, had a high point return count, and was free of vegetation.
Calculated snow depths on snow-free terrain should be zero, and can be used to assess bias. Orthoimagery collected contemporaneously with the LiDAR data was used to confirm the highway was snowfree, and a hillshade applied to the bare earth model guided the digitization of the highway. The nongap-filled snow depth and point density rasters were then clipped to the extent of the highway.
Snow depth values were plotted in histograms, and categorized as high (greater than 5), medium (2 – 5)
and low (less than two) ground point density (points/m2; Figure 4.2). Unbiased datasets will have peak
in the snow depth values at or close to zero. Histograms that skew above or below zero indicate of a
systematic error in the acquisition, and prompted a survey re-alignment to remove bias.

38

Figure 4.2 Histograms of the snow-free pixel values of the stable terrain for each snow survey acquisition in 2023. The ground point
density has been classified by greater than five points per m 2 (green), 2 – 5 points per m2 (yellow) and less than two points per m2 (grey).
Pixels with no ground points were excluded from this analysis.

4.2 In-Situ Observations and Supplementary Data
Over 2000 manually probed snow depth data points were collected over the course of the LiDAR
campaigns. A python script was developed to clean, summarize, and generate scatter plots from the data
output by Device magic. The plot summaries and scatterplots were reviewed to identify and correct any
errors within the dataset. In-situ snow depth data was then geo-positioned to the GNSS cardinal plot
centre.

39

4.2.1 Weather Stations
With the exception of Steph 3 which transmits via an Iridium Satellite Modem, weather station data
were manually downloaded through the season. The meteorological data then underwent a QA/QC
process to remove noise and outlier data points. With the cleaned datasets, analyses completed were:
i) Seasonal review: Overview of the snow seasons in 2022 and 2023 in terms of air temperature
and snow depth
ii) Short-term change: An assessment of air temperature and snow depth data in the period between
the manual sampling and LiDAR acquisition to determine if depth had changed
iii) Historical averages: comparison of how 2022 and 2023 snow melt seasons and average monthly
air temperature compare to previous years weather station data, going back to 2016

4.2.2 Game Cameras
To assess the relative change in snow depth between the LiDAR acquisition and manual observations,
game cameras captured daily images of five snow stakes within the forest at nine locations. Imagery for
2022 and 2023 was processed using the image manipulation software GIMP to create a workflow to
measure pixels of snow, using the snow stakes increment of 20 cm to calibrate a pixel to snow ratio. A
snow depth measurement per snow stake was acquired on the date of the LiDAR survey, as well as 10
days before and after each survey, to calculate the relative change in snow depth. Game camera
observations were removed from this analysis if the snow stake had shifted or was not straight, and
changes in snow depth were averaged between the five stakes.

4.3 Plot Analyses
To compare the manual and LiDAR snow depths at the plot scale, the mean snow depth per plot, rather
than a point-to-point comparison, was employed. The use of mean snow depth reduces errors that may
40

be introduced due to minor (less than 1 meter) offsets in position, the accuracy (<1m) of GNSS plot
centres, as well as field sampling practices (see 3.5.2 Sampling error and mitigation).
LiDAR snow depth data was extracted for each plot using a 10 meter radius buffer from the GNSS plot
centre for the cardinals, and a 5 m width line from the start and end GNSS points for the weather station
snow courses. The total number of pixels captured in the LiDAR were variable due to a) the number of
full and partial pixels contained within the plot buffer and b) the number of valid data points within the
plot buffer for the non-gap filled rasters (Figure 4.3). The maximum number of LiDAR data points were
312 for the cardinal plots and 33 for the weather station snow courses. Manual snow depth
measurements were generally 33 per cardinal plot, and 10 per a weather station snow course. Cardinal
plot results were analysed individually, as well as aggregated by their respective cover type: juvenile,
regenerating, old growth, or alpine (see sec 3.5.1)

Figure 4.3 Ortho-imagery and snow depth data for dense old growth forest plot, ST1OG6, from the April 29th 2023 survey. The in-situ (A)
and LiDAR (B) snow depths are averaged separately, and mean difference calculated. Missing pixels in the LiDAR are a result of no valid
snow surface returns for the survey.

Additional raster datasets used in the plot analysis include LiDAR-derived aspect, slope, and ground
point density, and the canopy metrics of height and density. Canopy coverage and slope are both known
topographical components which impact both LiDAR error and snow distribution (see sec 2.1 and 2.2).
The ground point density is indicative of LiDAR error, as poor ground point coverage results in more
41

interpolated data points. The ground point density rasters were used to determine the total ground points,
the proportion of pixels with no data (data gaps), and average point density per each plot. The point
density analysis was completed for all snow-on surveys, as well as both (2017 BCTS and 2023 ACO)
bare earth models.

4.3.1 Test of Normality
Snow depth distributions for each plot were assessed for normality, as this dictated which statistical tests
were valid for testing equality of variability and central tendency. Three metrics which were used to
determine whether the plot snow depths followed a normal distribution were 1) a q-q plot, 2) a
histogram, and 3) the sample Kolmogorov-Smirnov (KS) Goodness-of-Fit Test. The KS test is used to
test the null hypothesis that a sample is distributed according to a normal distribution. Using a
confidence level of 95%, the null hypothesis is rejected if the p-value is less than 0.05, and was
calculated using the built in KS function within the SciPy stats package.

4.3.2 Snow Depth Variability
Snow depth variability is known to be influenced by factors such as canopy cover and wind exposure. A
highly variable snowpack poses more of a challenge to LiDAR, as the conventional gridded format will
naturally smooth a dataset. A variability analysis was conducted for each of the different forest covers to
assess a) how the variability of the snowpack changes through the snow melt season and b) whether
LiDAR captures the same variability as the manual measurements.
For each plot, the standard deviation of snow depth was calculated independently for the manual
observations and LiDAR. The coefficient of variation was not used due to misleading results as snow
depths approach zero. Standard deviation was first compared separately for manual and LiDAR
separately for each cover type, to establish whether the snow depth variability changes as snow melts.
42

Then, to assess whether LiDAR captures the same variability as the manual observations, the standard
deviations were compared between the manual and LiDAR measurements.
The Levene’s test was selected to assess for homogeneity of variance, as it does not require normality of
data or equal sample sizes. The Levene’s test was conducted for two different null hypotheses for each
plot:
i) Snow depths (LiDAR and manual observations separately) are equal in variance through the
snow melt period, with number of groups (n) ranging between 2 and 5, dependant on number of
surveys and/or manual observations per snow season
ii) LiDAR and manual measurements captured equal variance (n=2)
Assumptions for the Levene’s tests were that observations are independent, and the test variable is
quantitative (not nominal or ordinal). Levene’s test results were calculated with the SciPy python
package, where a p-value of less than 0.05 at 95% confidence level rejects the null hypothesis for
homogeneity of variance.

4.3.3 Mean Snow Depth
To assess the accuracy of LiDAR within the different forest cover types, the average snow depth from
manual sampling and LiDAR were compared for each plot. Mean difference (MD) was then calculated
using the in-situ observations as the true measure of snow depth (Eq. 2). A positive value for MD
indicates LiDAR over-estimation of snow depths and vice versa.
(2)

‫ = ܦܯ‬µௌ஽ ௅௜஽஺ோ − µௌ஽ ௠௔௡௨௔௟

43

The Mann Whitney U test, a non-parametric statistical test, was selected to compare whether the LiDAR
and manual measurements are equal in median. Again, this was selected for the less stringent
assumptions, including:
i)

Two independent, non-paired samples

ii)

A minimum of 10 samples per group

iii)

Independence of observations

As with the previous statistical test, the manual and LiDAR snow depth measurements were not required
to be normally distributed, or of equal sample size. The Mann Whitney test uses the rank sum of the
snow depths, and the test statistic is calculated as the smaller of ܷ௠௔௡௨௔௟ (Eq. 3) and ܷ௅௜஽஺ோ (Eq. 4)
௡ (௡೘ ାଵ)

ܷ௠௔௡௨௔௟ = ݊௠ ݊௅ + ೘

ଶ

௡ (௡ ାଵ)

ܷ௅௜஽஺ோ = ݊௅ ݊௠ + ಽ ଶಽ

(3)

− ܴ௠

(4)

− ܴ௅

where n is the number of samples and R is the sum of the ranks, respective for the manual (m) and
LiDAR (L) snow depth measurements. The test statistic was then compared to the critical value at the
5% significance level, and null hypothesis rejected if the test statistic exceeds the critical value.
A consistent MD within a particular cover type, ie. if LiDAR over- or under-estimates snow depth, may
be the result of a systematic bias within the bare earth model. To establish the role of the bare earth
model on LiDAR snow depth accuracy, the workflow described above was repeated with the 2017
BCTS bare earth model.

4.3.4 Timing error and mitigation
A delay between LiDAR acquisition and manual sampling can be a source of error, as changes to the
snowpack may occur between the two observations. A longer period of time between LiDAR and
44

manual sampling will generally increase the likelihood of snowpack change between the two
observations, especially when in the snow ablation period. However, Geissler et al. (2023) noted a
systematic bias due to snowpack settling in as short as a four-hour delay between manual sampling and
LiDAR collection. To address this source of error, the weather station snow depth sensors and snow
depth derived from the game camera stations were used to detect snow depth changes between the
LiDAR and manual measurements. Plots which were forested (regenerating and old growth) used the
nearest game camera imagery, and the more open plots (juvenile, alpine, and weather station snow
courses) used the nearest weather station snow depth data. Relative change to the snowpack was
calculated by subtracting snow depth observed on the date of the manual sampling, from the snow depth
on the date of the LiDAR acquisition. The snow depth change value is referred to in the results as a
meteorological bias correction.
A second timing-based error may incur when plots are manually sampled prior to the LiDAR
acquisition. Here, there is potential for snow compaction from the survey team activities, which may
result in an under prediction of the snowpack by LiDAR. Snow compaction is difficult to account for, as
compaction may be non-existent when the snowpack is firm, or substantial in softer conditions. Ideally,
approximate snow compaction was documented in field notes, however this was not a consistent
protocol for each plot. Plots which were sampled prior to the LiDAR acquisition are noted in the results,
but no error correction applied.

45

4.4 Watershed Scale Analysis
The final portion of the analysis describes total snow storage within the Russell Creek watershed. Here,
the objectives were to assess:
i) Snow distribution in the watershed, in terms of the forest cover types and elevation bands
ii) How the snow distribution changes through the snow melt season
iii) Whether biases observed at the plot scale were present when scaled to the watershed
iv) The impact of the bare earth model to the proportion of snow stored by cover type and elevation
band, as well as the total snow storage in the watershed
The interpolated (gap-filled) snow depth grids were used for the watershed scale analysis. The Russell
Creek watershed was delineated from a 30 m SRTM DEM using ArcGIS Pro hydro tools, and the gapfilled snow depth grids were masked to the watershed extent. As the pixel size of the raster was 1 meter,
snow volume for each pixel was calculated as:
ܵ݊‫݉[ ݁݉ݑ݈݋ܸ ݓ݋‬ଷ ] = (ܵ݊‫ݐ݌݁ܦ ݓ݋‬ℎ [݉]) × (1݉ଶ )

(5)

This process was completed using the snow depth grids generated with both the 2023 ACO and the 2017
BCTS bare earth model. For each snow survey, snow volume was aggregated by 200 meter elevation
bands, as well as the cover type classifications.

46

Chapter 5 : Results
This chapter summarizes meteorological and snowpack observations collected in 2022 and 2023,
vegetation characteristics of the snow survey plots, and the plot scale mean and variability in snow
depths from (a) ground surveys, (b) LiDAR measurements made with the 2023 ACO bare earth model,
and (c) LiDAR measurements made with the 2017 BCTS bare earth model. The final portion of this
chapter examines watershed-scale snow storage, the impacts of bare earth model selection, and the
variability of snow storage across stand types.

5.1 Snowpack and Weather
Observations of air temperature from a selection of weather stations located in the Russell Creek
watershed (Figure 3.5) were used to compare the 2022 and 2023 accumulation and melt seasons, and
identify changes in snow depth (snow accumulation or melt) that occurred in the days between LiDAR
acquisitions and in-situ snow depth sampling. Automated observations of snow depth were available at
four weather stations.
Visual inspection of 2022 and 2023 snow depth and temperature records (Figure 5.1) show distinct
differences in the snow accumulation and melt periods between seasons. Snow accumulation was
slightly greater in 2022, and the timing of maximum snow depth was approximately 1 week earlier in
the alpine in 2022 than in 2023. The snowpack melted rapidly in 2023, with a short period (26 – 42
days) between melt onset and snow-free conditions. In 2022, the snowpack persisted approximately one
month later.

47

Figure 5.1 Hourly air temperature (upper) and snow depth (lower) data recorded from Steph 2(1164m) weather station, for water years
2022 and 2023

The 2022 snow season was the longest since 2017 at both the Steph 6 (827 m) and Steph 2 (1164 m)
weather stations (Figure 5.2). The 2022 snowpack development started slow, but by the time of the first
ACO flight was completed, snow depth was near average for the rest of the season and the second latest
seasonal peak snow, before a rapid melt at Steph 2 and 6.

48

Figure 5.2 Snow depth records from 2016 - 2023 for Steph 6 (upper) and Steph 2 (lower) weather station, with the years of this study
highlighted with a blue (2022) and red (2023) dashed line respectively. Of the 7 years of data shown for Steph 6 and Steph 2 weather
stations. 2022 had the latest date on record for snow free, whereas 2023 had the second earliest

Air temperatures in March 2022 were warmer than March 2023, with average temperatures ranging
from 0.4 to 2 °C across all stations (Table 5.1; Table 5.2). For April, May, and June the temperatures
were warmer in 2023 than 2022. Notably, May 2023 had the greatest divergence in average monthly air
49

temperatures. When compared to the previous 7 years of weather station data, average monthly air
temperatures in the alpine were -3.3 °C cooler in 2022, versus 3.8 °C warmer in 2023.
Table 5.1 Average daily air temperature across all weather stations. Highlighted values represent that the temperature was higher between
the two years, and a dash indicates there was no data for that month.

station
Steph 8
Steph 7
Steph 6
Steph 4
Steph 3
Steph 2
Steph 1

elevation
576
707
813
935
1004
1158
1508

March
2022
2023
2.6
1.4
1.4
1.0
1.1
0.1
1.4
1.0
0.9
-0.4
-0.3
-1.8
-2.3
-4.3

Average Air Temperature
.
April
May
2022
2023
2022
2023
2.3
3.5
6.0
11.2
2.9
6.2
10.6
0.6
2.3
3.8
5.6
2.9
6.2
10.6
0.1
0.9
3.3
8.6
-1.3
0.2
2.0
8.9
-3.8
-2.0
0.5
7.6

June
2022
2023
11.0
12.1
10.6
11.4
9.4 10.6
11.4
9.2
9.7
7.6
10.3
5.4
8.4

Peak snow depths in 2022 were greater than 2023 at all weather stations, ranging from 8-26 % deeper
(Table 5.2). With the exception of Steph 3, peak snow depth occurred 3 – 16 days later in 2022. The
cooler temperatures in May 2022 resulted in the snowpack persisting 41 days longer at the alpine
stations in 2022 than in 2023. Generally, snow melt periods were longer in 2022 than 2023 weather due
to the greater antecedent snowpack.
Table 5.2 Summary of 2022 and 2023 snow depth measurements at four weather stations, selected for consistent snow depth data across
2022 and 2023. Metrics include: 1) peak snow depth and date, 2) length of melt period in days (defined as sustained loss of 0.5 cm per
day) and rate of snow melt, and 3) snow free date

station
Steph 6
Steph 3
Steph 2
Steph 1

peak snow depth
(cm, date)
elevation
2022
2023
813
125.8 (5-Apr)
91.9 (21-Apr)
1004
106.1 (5-Apr)
97 (3-Mar)
1158
212.2 (5-Apr)
170.3 (21-Apr)
1508
334.5 (13-Apr) 246.8 (16-Apr)

50

snow loss rate
[cm/day, (days of melt)]

2022
4.3 (42)
6.1 (8)
5.1 (35)
6.1 (47)

2023
4.3 (36)
5.7 (13)
6.1 (26)
4.7 (42)

Snow-free date
2022
2023
28-May
6-May
26-May
3-May
23-Jun
18-May
7-Jul
27-May

5.2 Vegetation Surveys
To evaluate the role of ground vegetation in LiDAR ground return point density, characteristics of the
vegetation at each cardinal plot were collected in October, 2022. The height of ground vegetation is used
here as a proxy for the vigour and density of ground vegetation that might affect LiDAR point densities.
Between 70 to 232 individual plant height measurements were collected at each cardinal plot (Figure
5.3; Table 5.3). Though a thorough species inventory was not collected, the most numerous species
within the forested sites belonged to the genus Vaccinium (blueberries and huckleberries) and Asplenium
(ferns), and Phyllodoce (mountain heathers) in the alpine. Vegetation measuring less than 5 cm in height
and trees were not included in the results. Species were broadly categorized by whether they were a)
woody debris, for example logs and stumps b) heather, or c) shrubs and ferns.
Table 5.3 Vegetation composition by cover type, summarized by the average height (μ) and the standard deviation (σ) of the
shrubs and ferns, woody debris, and heather.

Vegetation

Shrubs and Ferns

Woody Debris

Heather

Cover Type

μ
[cm]

σ
[cm]

%
of total

μ
[cm]

σ
[cm]

%
of total

μ
[cm]

σ
[cm]

%
of total

Juvenile

116.3

44.9

84.7

71.2

66.6

15.3

n/a

n/a

n/a

Regenerating

78.8

87.8

88.1

72.8

55.1

11.9

n/a

n/a

n/a

Dense Old Growth

36.7

29.8

70.3

28.0

38.1

12.9

10.0

8.0

16.8

Old Growth

52.4

37.0

71.9

35.7

15.8

1.3

32.3

15.2

26.8

Alpine

26.5

26.3

57.1

n/a

n/a

n/a

19.5

10.0

42.9

The average shrub height was greatest in the cardinal plots located in the juvenile and regeneratingforests, at 116 and 78 cm, respectively. The regenerating and juvenile forests had a high proportion of
woody debris – stumps, logs, and slash piles – that are bi-products of the forest harvest, and they also
51

had the greatest variability in shrub height of the cover types. Large canopy gaps and relatively small
trees (less than 5 m) in juvenile forest stands allow for the ground vegetation to grow to greater heights.
The tree and shrub layer were the least distinct within the juvenile forest, and in some cases the shrubs
exceeded the tree heights.

Figure 5.3 Average vegetation height by cover type broken into the categories of shrub, woody debris, and heather. Vegetation measured 5
cm or less, as well as trees, were not included in the results. Error bars represent the 95% confidence interval.

The old growth forest was predominantly comprised of a single shrub layer, which was highest in the
stands with large canopy gaps. As the canopy density increased in the cardinal plots, the average height
of the shrub layer decreased. Alpine vegetation was predominantly low-lying mountain heather, with
clusters of stunted trees and sparsely spaced shrubs. The alpine shrub layer was the lowest average
height and the lowest variability of all the cover types.
Though ground vegetation layers were present at all cardinal plots, the harvested (juvenile and
regenerating) cardinal plots were considerably more complex in terms of the height, density and
52

variability of the vegetation. As a result, the LiDAR ground earth model is likely to be less accurate in
the harvested forests due to the a) potential misclassification of vegetation as ground, and b) reduced
ground points due to the lasers inability to penetrate through the vegetation to the ground.

5.3 LiDAR Point Densities
The LiDAR snow surface point density returns varied by both the forest cover type, and the date of
acquisition (Figure 5.4). Overall, as the snowmelt season progressed, the average point density
decreased and the number of data gaps (pixels with zero ground returns) also increased (Figure 5.4). No
changes were made to the planned snow survey density between the acquisitions.

Figure 5.4 Percentage of the cardinal plots which had at least one LiDAR ground or snow surface point return per m2, grouped by forest
cover type, for each LiDAR acquisition in 2022 (upper) and 2023 (lower). Hatched bars denotes acquisitions that were snow-free for all
plots within the cover type. Error bars represent the 95% confidence interval.

53

With the exception of the regenerating forest, all plots in 2022 were snow-covered for each acquisition
due to the late snow melt. As a result, 2022 had higher average point density per m2, as well as
proportion of snow surface returns per plot than 2023. In contrast, the earlier snowpack meltout in 2023
led to reduced ground point densities in the snow-free acquisitions. This is interpreted as being the result
of ground vegetation leaf-out, particularly in the juvenile and regenerating forests and between the May
25th and June 12th acquisition (Figure 5.5). An anomaly to this pattern occurs in April 2023 acquisition
in the dense old growth forest, as the average proportion of ground point returns drops from 92% in the
March 10th acquisition to 18% in the April 4th, despite snow conditions being stable. A less pronounced
decrease in point density is also observed within the old growth in April 2023, but the overall decrease
in returns through the melt season is consistent between the two acquisitions. Changes to the snow
surface, i.e. reduced reflectivity of snow with vegetation emerging or forest debris, may also contribute
to variability in the snow surface point density.
Table 5.4 Median ground point density per meter squared for the snow surveys, as well as the average proportion of valid pixels per each
of the cover types

LiDAR Acquisition
Date

2023

2022

Cover
Mar-9
Apr-11
May-21
Jun-24
Mar-10
Apr-4
Apr-29
May-25
Jun-12

Median GPD/m2 (Mean proportion of ground returns per plot [%])

juvenile
5 (99.0%)
5 (99.4%)
5 (97.4%)
3 (84.0%)
5 (98.5%)
3 (95.4%)
3 (96.6%)
1 (53.6%)
1 (14.2%)

regenerating
2 (66.6%)
3 (71.5%)
2 (58.6%)
1 (37.5%)
2 (62.9%)
1 (50.5%)
2 (49.6%)
1 (32.9%)
1 (12.4%)

dense old
growth
2 (83.3%)
2 (85.8%)
1 (69.1%)
1 (69.2%)
2 (92.2%)
1 (18.7%)
2 (75.8%)
1 (57.4%)
1 (57.4%)

old growth
2 (91.6%)
3 (93.6%)
2 (89.4%)
2 (85.0%)
3 (94.0%)
1 (57.4%)
2 (87.1%)
1 (71.5%)
1 (65.1%)

alpine
3 (97.6%)
3 (99.6%)
3 (98.2%)
2 (85.8%)
4 (97.8%)
3 (98.6%)
3 (98.3%)
2 (64.2%)
2 (66.8%)

Despite highest canopy cover occurring in the dense old growth forest, the overall lowest proportion of
snow surface and ground returns per plot was found in the regenerating forest, ranging from an average
54

high of 71.5% in high snow accumulation conditions, to a low of 12.4% when snow free. The alpine had
the fewest gaps in ground returns, and median point returns of 3 - 4 points m-2 in high snow
accumulation conditions, and 2 point m-2 in intermittent to snow-free conditions.
The greatest variability in the point density occurred in the juvenile forest. In snow-on conditions,
minimal canopy cover resulted in high point density and less than 10% of data points missing per plot.
As the snowpack thinned and melted, and vegetation exposed, the median point density dropped to 1
point per m2, and only 14.2% of ground points per plot.

Figure 5.5 Imagery of the juvenile forest from the (A) May 25 th 2023 and (B) June 12th 2023 acquisitions showing the green up of the
vegetation

5.4 In-situ vs LiDAR Snow Depths
LiDAR-derived snow depths, were compared with in-situ observations made at the weather stations and
cardinal plots. Weather station data was used to establish to what extent the snowpack conditions
changed between the time of the LiDAR survey and in-situ measurements. The mean difference (MD)
and standard deviation (SD) was compared between the different forest cover types.

55

5.4.1 Snow depths Validation
In-situ snow depth sampling was targeted to overlap with the LiDAR acquisitions, however this was not
always possible due to weather conditions and other logistics. In general, the manual observations took
place within 3 days pre- or post-LiDAR (Table 5.5). An exception to this was the June 24th 2022
acquisition, which occurred 6-10 days after the manual measurements. As this survey occurred during
the snowpack ablation period, automated snow depth observations indicate that a large loss of snow (35
– 55 cm) occurred between manual sampling and the LiDAR acquisition (Table 5.5).
Table 5.5 Summary table of LiDAR and manual dates for 2022 and 2023 and how the snow depth changed within that time frame. Snow
depth change includes the SR50 data of all seven weather stations in Russell Creek, provided the sensor was collecting data in that time
frame

Year

2022

2023

LiDAR
acquisition
date
March 9th

Manual sampling

SR50 snow depth change

March 8th - March 10th

-1 cm to 2 cm change from manual
sampling to LiDAR acquisition

April 11th

April 5th - April 10th

18 – 25 cm increase from the start of the
manual sampling to LiDAR acquisition

May 21st

May 21st - May 25th

June 24th

June 14th - June 18th

16 – 24 cm loss from LiDAR acquisition
to end of manual sampling
35 – 55 cm loss from start of manual
sampling to LiDAR acquisition

March 10th

March 9th - March 13th

April 4th

April 5th - April 7th

April 29th

April 26th - April 27th

May 25th

May 23rd - May 25th

June 12th

n/a

-2 to 2 cm change from manual sampling
period to LiDAR acquisition
-1 to 1 cm change from LiDAR
acquisition to manual sampling
13 – 22 cm of loss from the manual
sampling period to LiDAR acquisition
7 cm loss at Steph 1 recorded from start of
manual sampling period to LiDAR
acquisition
Snow free at plots and weather stations

LiDAR surveys in 2023 were generally aligned closer to the manual sampling windows than in 2022.
Manual sampling occurred in the days immediately prior to the April 29th LiDAR survey, yet a rapid rise

56

in temperature resulted in an average of 17 cm snow loss between start of sampling and the LiDAR
survey (Figure 5.6).

Figure 5.6 Steph 2 weather station showing the change in average daily snow depth which occurred between the manual sampling window
and the LiDAR survey for the April 29 th 2023 acquisition

5.4.1.1 Weather Station Snow Courses
The weather station snow courses (with the exception of Steph 1) are conducted on gravel resource
roads, and are therefore free from bias introduced by ground vegetation. Across the 2022 snow year, the
3rd of March survey had the lowest bias error between the LiDAR and manual sampling, with mean
difference (MD) ranging from – 6 to 8 cm (Figure 5.7). This snow course represents the ‘ideal’
comparison between manual and LiDAR-derived snow depths, as it was conducted shortly after the
March 9th LiDAR acquisition, with no disturbance to the snow surface from manual sampling, and the
weather was stable with no accumulation or melt (Table 5.5). The greatest MD occurred in the June 24th
acquisition, in which snow depth was underestimated by 37 to 56 cm due to melt between the snow
course and LiDAR acquisition. A comparison of the snow-on and snow-free plots from the June survey
57

show a much reduced, though positive bias for the snow free plots, confirming the large negative bias is
a result of the changes to the snowpack (Figure 5.7).

Figure 5.7 The mean difference for snow courses from the 2022 and 2023 snow season. A 10 day delay between the LiDAR acquisition and
manual sampling -which coincided with a period of snow melt - is attributed to the large negative bias for June 24th 2022

In 2023, the greatest deviation between the manual and LiDAR snow depths occurred in April, where
the April 29th LiDAR acquisition underestimated snow depth by 14 to 24 cm. This again coincided with
a melt period between the manual sampling and LiDAR acquisition (Figure 5.6). The snow courses
which were snow-free, Steph 7 and Steph 8, had a MD of -2 and 6 cm respectively, which indicates that
snowpack changes between manual and LiDAR surveys were responsible for the observed bias.
In general, the systematic error was greater in plots in which LiDAR was acquired after the manual
sampling occurred as opposed to before, as well as if it had snowed or melted in between the manual
surveys and LiDAR acquisition (Figure 5.8). This suggests that both meteorological variability as well
as potential snowpack interference from manual sampling can explain some of the systematic error for
specific surveys.

58

Figure 5.8 2022 and 2023 snow-on surveys mean difference grouped by whether LiDAR was flown pre or post manual sampling, and the
snow depth change between acquisition and sampling. Snow depth readings are compiled from the weather stations in Russell Creek,
Steph 1 – 8.

5.4.1.2 Cardinal Plots
The cardinal plots were situated in more complex terrain than the weather station snow courses, with
varying degrees of canopy cover and ground vegetation (see Table 3.3 and Figure 3.3). Over each of the
314 m2 plot, the average snow depth of 33 in-situ measurements was compared to a maximum of 312
LiDAR measurements, dependant of pixels with at least one ground return point within the plot. Overall,
systematic differences between LiDAR and manual snow depth observations at the cardinal plots were
found to be related to the forest cover type.
As identified previously, changes to the snowpack that occur between LiDAR acquisition and manual
sampling will affect the mean difference. To account for sources of accumulation or melt error, weather
station snow depth sensors at nearby meteorological stations were used to estimate the snowpack change
in the alpine and juvenile sites between manual and LiDAR surveys, while game camera images of snow
stakes were used to evaluate snow depth changes in the forested sites. When the LiDAR was flown after
59

the manual measurements, compaction of the snow surface during manual snow sampling may introduce
additional errors, depending on snow conditions. Plots in which the LiDAR was captured after manual
sampling occurred are marked with a star, but no offset applied.
Despite a significant melt period between manual sampling and the June 24 th 2022 LiDAR acquisition
(Table 5.5), these data will still be used within the validation of snow depths, as the meteorological bias
correction applied should counteract the melt. However, data from this survey will be omitted from the
aggregated results.
5.4.1.2.1 Alpine
Differences in manual and LiDAR derived snow depths were inconsistent between acquisitions (Figure
5.9), however not all plots were sampled every time due to elevated avalanche hazard. In 2022, alpine
plots were sampled for three of the four LiDAR surveys. The lowest mean difference (-7 cm) was
observed during the first acquisition. While the May 21st LiDAR survey showed average snow depths 40
and 50 cm greater than the manual survey, this was reduced by a correction of 10 cm to account for
snow melt that occurred between the LiDAR and manual measurements. In contrast, the negative bias
for the June 24th survey is consistent with the snow melt which occurred between LiDAR flight and
manual sampling (see Table 5.5).
For 2023, all five of the acquisitions were sampled, but the final acquisition were all verified snow free.
Overall, MD was lower in 2023 for the alpine plots, ranging from -20 to +11 cm across all surveys with
snowpack conditions that were stable between LiDAR and manual surveys. In snow-free conditions,
MD was reduced to -4 and -2 cm, showing good agreement between the bare earth surfaces within the
plots.

60

Figure 5.9 Alpine cardinal plots mean difference between LiDAR and in-situ snow depth measurements for the 2022 and 2023 surveys.
Meteorological corrections from the weather station snow depth data is denoted with the darker bars, and asterisk mark if LiDAR was
acquired after manual measurements. Hatched bars indicate plot was snow free.

5.4.1.2.2 Old Growth
Bias was similar across acquisitions and plots for the dense old growth forest (bias across all
acquisitions and plots. For both years, LiDAR snow depths were 3 to 33 cm lower than manual
observations (Figure 5.10). After the snow melt correction was applied to the June 24th 2022 acquisition,
the LiDAR snow depths were slightly higher than manual observations.

61

Figure 5.10 Dense old growth cardinal plots mean difference between LiDAR and in-situ snow depth measurements for the 2022 and 2023
surveys. Meteorological correction applied using the game camera data is denoted with the darker bars, and asterisk mark if LiDAR was
acquired after manual measurements. Hatched bars indicate plot was snow free.

The highest MD (29 to 33 cm) for 2023 occurred on the April 29th survey, which may be attributed to
the manual measurements occurring prior to LiDAR acquisition. Field notes documented approximately
10 cm compaction - as well as rapid warming and melt (see Table 5.5). In the one day between manual
sampling and the LiDAR survey on April 29th, the forest game camera situated in nearby forest stand
(ST1GC1) show an average 10 cm of melt (Figure 5.11). When conditions were snow-free, the bias
remained negative, suggesting a consistent bias to the ground surface model in the dense canopy forest.

62

Figure 5.11 Game camera ST1GC1 showing the rapid snow melt which occurred, with a loss of approximately 10 cm between 10 am on
the 28th of April to 3 pm on the 29th of April 2023. Snow stakes are marked by 20 cm increments.

Manual and LiDAR-derived snow depths in the old growth plots were characterized by higher snow
depths and variability (Figure 5.12). The median of the mean difference in snow depth was -17 cm. No
significant differences in snow depths were observed for the March 9th 2022 LiDAR survey. For the
May 21st 2022 survey, the LiDAR overestimated the snow-pack by an average of 82 cm for plot
ST1OG3, which was reduced to 74 cm after a correction for melt was applied. Notably, a) ST1OG3 is
situated in a steeper terrain, and is the highest of the old growth sites, and b) the May 21st survey also
saw an over-estimation of snow depth in the two alpine plots. The greatest bias for old growth plots
OG7 and OG8 – which are spatially closest together – occurred on the April 4th 2023, where LiDAR
underestimated the snowpack by average of 57 and 103 cm respectively. The snow conditions were
stable, sampling occurred post-LiDAR surveys, and manual measurements have been checked for error.
However, this is consistent with the reduced proportion of snow surface returns shown in Figure 5.4.

63

Figure 5.12 Old growth cardinal plots mean difference between LiDAR and in-situ snow depth measurements for the 2022 and 2023
surveys. Error corrected using the game camera imagery is denoted with the darker bars, and asterisk mark if LiDAR was acquired after
manual measurements. Hatched bars indicate plot was snow free.

5.4.1.2.3 Regenerating
In the regenerating forest plots LiDAR consistently over-estimated snow depths by 25 to 106 cm (Figure
5.13). Biases in snow depth during three snow-free acquisitions (i.e. snow depth should be zero) ranged
from +160 to -30 cm, which indicates poor agreement between the snow-free elevations and the ground
surface model. Plot ST6GC1 generally had the lowest MD (29 to 56 cm) in snow-on conditions, and is a
denser forest canopy than both ST4RF1 and ST4RF2.

64

Figure 5.13 Regenerating cardinal plots mean difference between LiDAR and in-situ snow depth measurements for the 2022 and 2023
surveys. Error corrected using the meteorological or game camera data is denoted with the darker bars, and asterisk mark if LiDAR was
acquired after manual measurement. Hatched bars indicate plot was snow free.

5.4.1.2.4 Juvenile
LiDAR over predicted snow depth in the juvenile forest cover in both snow years, though at much
different magnitudes between the two plots (Figure 5.14). Of all the cover types, the juvenile plots
showed the greatest error variability in snow-on conditions (32 to 127 cm), despite being similar in
topography. Plot ST2CC2 had a MD which was consistently in a range three times greater than that of
ST2CC1. Notably, as the vegetation greened up, the snow free error changed from positive on the May
25th survey (31 and 127 cm) to negative on the June 12th (-9 and -27 cm).

65

Figure 5.14 Juvenile cardinal plots mean difference between LiDAR and in-situ snow depth measurements for the 2022 and 2023 surveys.
Error corrected using SR50 snow depth sensor denoted with the darker bars, and asterisk mark if LiDAR was acquired after manual
measurement. Hatched bars indicate plot was snow free.

66

5.4.1.3 Summary
Overall, differences between the LIDAR and manual snow depth measurements were greater in the
cardinal plots (characterized by complex terrain and a range of ground and canopy covers) than the
weather station snow courses (characterized by graded surfaces with no ground vegetation or canopy
cover; Figure 5.15). Though the meteorological corrections proved beneficial to some plot biases, the
results were mixed, and in some cases bias increased. Therefore, MD ranges will be reported for the
non-meteorological corrected version, with the June 24th 2022 survey excluded due to the significant
melt which occurred (Table 5.5).
The Median MD for the weather station snow courses were -4 and -6 cm for 2022 and 2023
respectively. The greatest MD occurred in the regenerating (median MD: 55 cm) and juvenile (median
MD: 76 cm) stands, with LiDAR generally over-predicting snow depth in the harvested stands. Ground
vegetation in the regenerating and juvenile stands was notably much higher and more complex than the
other cover types. Snow depths in regenerating plots were consistently overestimated by LiDAR across
all three of the plots, with MD between +12 to +106 cm (Figure 5.14). The two juvenile plots were
similar in topography, canopy cover, and ground vegetation, however they had two distinctly different
ranges in MD (2 to 41 cm and 121 to 150 cm). Old growth forest plots had a general negative bias, with,
a smaller range of MD observed in the denser canopy sites (-25 to -6 cm) than the more open sites (-58
to 4 cm; excluding ST1OG3, complex, sub-alpine plot). Alpine plots were inconsistent between snow
years, with a mostly positive, and wider MD range in 2022 (-6 to 50 cm), and a negative range in 2023 (18 to -4 cm). However, only three pairs of LiDAR and manual measurements are included for 2022. As
observed with the weather station snow courses, the snow depth bias observed at the cardinal plots was
generally reduced when sampling took place shortly after the LiDAR acquisition, when minimal
changes to the snowpack had occurred.

67

Figure 5.15 Boxplots of the mean difference for each plot by cover type, with “open” referring to the weather station snow courses. For
each vegetation type, the left box represents 2022 data and the right box 2023. For 2022, phase 4 has been excluded due to the noted offset
due to delay in LiDAR acquisition and manual sampling, and subsequent snow melt. The box shows the quartiles of the dataset, with
whiskers extending to points within 1.5 x the inter-quartile range. Observations outside of this range are displayed as points.

In addition to the forest cover and ground vegetation, slope is a topographic parameter which will
impact LiDAR accuracy (see sec 2.1.3). The sampling strategy of this thesis focused on forest
composition, and did not target slope angle to assess LiDAR accuracy due to both the well-established
relationship between slope angle and LiDAR error (eg. Painter et al. 2017), as well as challenges of
manual sampling in steep terrain. Therefore the average slope angle of the plots did not exceed 35
degrees, and, as a result, there was no clear relationship between the gradient of the plot, and the mean
absolute bias (Figure 5.16).

68

Figure 5.16 Scatterplot of the slope angle of the cardinal plots and the mean absolute bias error for the 2022 and 2023 LiDAR surveys.
Each marker represents the mean absolute bias error of the cardinal plot, color coded by cover type, and the vertical bars represent +/sigma.

69

5.4.2 Variability of Observed Snow Depths
The spatial variability of snow depth is an important parameter (Lopez-Moreno et al., 2011) that may
vary through the accumulation and melt seasons (Grunewald et al., 2010). However, there are few
comparisons between LiDAR-derived and manual snow depth that specifically investigate the snow
depth variability at the plot scale. This section summarizes the variability of snow depths measured (a)
manually and (b) with airborne LiDAR scans at the snow courses and cardinal plots throughout the snow
melt season.
5.4.1 Weather Station Snow Courses
Snow depth variability from both manual (Table 5.6) and LiDAR-derived (Table 5.7) snow depth
observations were greatest at Steph 1, where standard deviation ranged from 16 to 31 cm across all three
acquisitions. Steph 1 is an alpine site with reduced canopy and high wind exposure which results in
greater wind transport of snow. As well, the ground surface of the Steph 1 snow course is much more
variable than the rest of the weather station snow courses.
The lowest snow depth variabilities for both manual and LiDAR observations occurred in the low
elevation and sheltered snow courses (Steph 8 and Steph 7), which do not receive much direct sun or
wind. Manual observations yielded standard deviations between 3.0 and 5.2 cm at these sites, while
LiDAR observations gave standard deviations between 7.4 and 12.1 cm. Overall, the snow depth
variability did not change through the snow melt period for the manual nor LiDAR, though Steph 2
showed statistically different variance between sampling periods in 2022.
For the weather station snow courses, the manual sampling and LiDAR observations captured similar
snow depth variability across both melt seasons (Figure 5.17). Of the 46 instances of manual and LiDAR
comparison measurements, only 6 were statistically different in variability (full results in Appendix D).
70

These six occurred at lower elevations and were spread over all acquisitions. Both manual and LiDAR
snow depth measurements identified the same sites as having the greatest snow depth variability (Steph
1 and Steph 3).
Table 5.6 Standard Deviation (SD) of the in-situ snow depth observations at each weather station snow course. Snow courses which were
snow free are shaded blue. To assess whether the in-situ snow depth variability changed as snow melted, the Levene’s test was conducted,
where for each weather station and year, each groups (LiDAR survey) variance were compared to each other. Snow free and no data
surveys were not included as groups.

2023

2022

In-situ Observations - Snow Depth Standard Deviation (cm)
Weather Station
Steph 1
Steph 2 Steph 3 Steph 4 Steph 6 Steph 7
9-Mar
33.0
4.5*
5.3
7.5
4.7
5.2
11-Apr
4.5*
14.9
8.4
3.3
4.8
21-May
29.9
10.7*
8.6
10.5
4.4
0.0
24-Jun
30.3
4.6*
0.0
0.0
0.0
0.0
10-Mar
16.9
3.9
5.4
5.0
3.1
2.3
4-Apr
29.5
4.2
6.3
5.2
5.5
1.8
29-Apr
7.3
7.7
5.9
4.4
0.0
25-May
22.9
0.0
0.0
0.0
0.0
0.0
12-Jun
0.0
0.0
0.0
0.0
0.0
0.0
*
Statistically significant at a p-value of 0.05.

Steph 8
4.1
3.0
0.0
0.0
1.7
1.6
0.0
0.0
0.0

Table 5.7 Standard Deviation of the LiDAR collected snow depth at each weather station snow course. Snow courses which were snow free
are shaded blue. To assess whether the LiDAR snow depth variability changed as snow melted, the Levene’s test was conducted, where for
each weather station and year, each groups (LiDAR survey) variance were compared to each other. Snow free and no data surveys were
not included as groups.

LiDAR - Snow Depth Standard Deviation (cm)

2023

2022

Weather Station
9-Mar
11-Apr
21-May
24-Jun
10-Mar
4-Apr
29-Apr
25-May
12-Jun
*

Steph 1
25.0
27.5
33.5
35.7
16.8
25.5
22.3
21.1
2.7

Steph 2 Steph 3
7.1*
9.1*
5.9*
15.5*
6*
13.7*
8.2*
5.1
6.5
10.8*
6.3
10.5*
6
9.6*
5.1
1.6
3.5
2.8

Steph 4
4.1
6.7
7.5
4.5
3.8
4.1
7.8
4.8
3.3

Statistically significant at a p-value of 0.05

71

Steph 6
6.9*
5.8*
9.1*
6.5
5.4
7.0
8.6
5.4
9.3

Steph 7
9.3
12.1
8.3
8.3
7.5
8.3
6.3
7.5
8.6

Steph 8
7.4*
8.7*
6.0
8.5
4.9
4.1
4.2
4.4
5.1

Figure 5.17 Snow distribution for Steph 1 weather station snow course for each LiDAR survey in 2022, showing the manual (red) and
LiDAR (blue) snow depth distribution. None of the measurements were significantly different in snow depth variability. The offset between
the LiDAR and the manual sampling in the June 24th acquisition is due to the snow melt which occurred between manual sampling and
LiDAR survey. No in-situ measurements were taken at Steph 1 on April 11, 2022.

5.4.2 Cardinal plots
Snow depth variability was greater at the cardinal plots than at the weather stations for both the manual
and LiDAR observation (Table 5.8; Table 5.9). While alpine sites were inconsistently sampled due to
varying avalanche conditions, snow depth variability calculated using both manual and LiDAR
observations was relatively stable through the 2022 and 2023 snow melt seasons.

72

Table 5.8 Standard deviation for in-situ snow depth measurements for the cardinal plot sampled across 2022 and 2023. Cells highlighted
blue indicate snow-free. To assess whether the in-situ snow depth variability changed as snow melted, the Levene’s test was conducted,
where for each plot and year, groups (LiDAR survey) variance were compared to each other. Snow free and no data surveys were not
included as groups.

12-Jun

25-May

29-Apr

4-Apr

10-Mar

34.6

27.8

30.0

37.3

-

28.2

0

ST2A2

-

-

44.5

46.9

30.8

38

-

34.9

0

O Growth

24-Jun

11-Apr
-

ST1OG3

-

-

57.8

44.1

39.2

42.3

-

26.2

0

ST1OG7

34.7

46.9

41.7

36.8

26.0* 37.1* 40.2* 21.1*

0

ST1OG8

44.8

46.4

44.3

39.9

33.2* 51.4* 43.4* 16.3*

0

D O Growth

9-Mar

Alpine

25.5

ST1OG4

-

30*

17.3*

13.6* 21.3* 15.8* 21.3* 0

0

ST1OG5

24.6

29.1

28.0

19.8

21.3* 29.4* 25.9* 5.1*

0

ST1OG6

21.3

28.7

25.5

20.1

22.4* 19.1* 18.2* 5.8*

0

ST4RF1

54.2*

40.7* 30.1*

0

24*

30.9* 33.2* 0

0

ST4RF2

37.6

40.3

46.7

0

27.5

27.2

0

0

ST6GC1

33.2*

55.3* 29*

0

29.7* 35.5* 27.9* 0

0

ST2CC1

-

51.2* 46.4*

5.4*

25.4

34.3

40.8

0

0

17.6* 23*

26*

34.6* 0

0

Regen

ST1A1

Juvenile

*

2023

21-May

2022
In-Situ
measurements
SD (cm)

ST2CC2
30.2* 46.5*
Statistically significant at a p-value of 0.05

34.3

Open canopy old growth forest plots had the greatest snow depth variability among the manual
observations (SD = 16.3 to 57.8) but not in the LiDAR (SD = 21.1 to 61.5 cm) observations, and the
highest snow depth variability coincided with the timing of the greatest average snow depth in both
years (April acquisitions). With both snow depth observation types, plot ST1OG3 stands out with higher
variability in snow depths due to its location in the sub-alpine and greater exposure to wind. At all old
growth sites and in both datasets and seasons, snow depth variability decreases from the maximum snow
depth through to the end of the melt season.
Dense old growth forest overall had less variable snow depths than open canopy sites in both manual
(SD = 5.1 to 30.0 cm) and LiDAR (SD = 7.1 to 26.1 cm) observations, and was the least variable
73

snowpack across all the forest cover types. Similar to the old growth, snow depth variability was
consistent across all acquisitions and observation types for 2022, and more variable in 2023 as the plots
were largely snow free by the May 25th acquisitions.
Regenerating and juvenile plots showed no change in snow depth variability as the melt seasons
progressed, and had the highest variability observed with the LiDAR observations (SD = 46 to 93 cm).
This variability is likely due to ground features endemic to harvested terrain, such as the stumps, slash
piles and compressed shrubs. Leaf-out in the spring would also affect the variability observed through
the LiDAR observations (but not the manual probing).

25-May

12-Jun

ST2A2

35.2

34.2

28.5

32.5

26.9

22.9

ST1OG3

51.8* 61.5* 52.3*

41* 34.7* 27.2* 33.5* 25.5*

10.8

ST1OG7

33.5* 44.3* 39.7* 36.7* 26.4* 28.4* 29.3* 23.3*

10.9

ST1OG8

34.7* 41.3* 36.4* 32.6* 26.3*

19.5

ST1OG4

20.5* 23.7* 23.7* 18.5* 13.5* 7.1*

ST1OG5

23.6* 26.1*

28.8

30.5

29-Apr

10.8

4-Apr

21-May

34.4 32.1* 32.9* 33.5* 17.2*

10-Mar

11-Apr

36.0

9-Mar

36.2

Alpine

32.1

Old Growth

2023

ST1A1

32.8

31* 25.6* 21.1*

9.8

26* 16.8* 17.9* 12.3* 20.1* 10.2*

10.2

ST1OG6

18.2* 22.7* 20.7* 16.4* 14.0* 12.8* 16.4* 11.4*

9.1

ST4RF1

80.0

85.2

78.4

93.1

83.3

76.1

72.2

72.2

32.5

ST4RF2

66.0

62.6

67.7

88.1

61.1

62.3

63.6

65.0

59.7

ST6GC1

75.6

75.8

70.0

81.7

71.4

69.9

74.2

67.7

41.2

ST2CC1

58.2

62.3

61.0

60.4

58.1

57.7

59.4

55.0

80.5

ST2CC2
53.5 50.9 51.1
Statistically significant at a p-value of 0.05

63.6

57.4

54.4

52.6

58.0

46.0

Regen

10.8

Juvenile

*

2022

Dense Old
Growth

Standard Deviation
(cm)
LiDAR

24-Jun

Table 5.9 Standard deviation for LiDAR cardinal plot snow depth measurements across both years. Cells highlighted blue indicate snowfree. To assess whether the LiDAR snow depth variability changed as snow melted, the Levene’s test was conducted, where for each plot
and year, groups (LiDAR surveys) variance were compared to each other. Snow free and no data surveys were not included as groups.

74

13.9*

Only two of the cardinal plots were equal in variance between the manual and LiDAR measurements for
every acquisition: the alpine plot ST1A1, and old growth plot ST1OG7 (Figure 5.18). Of the cover
types, the old growth forest had the most cases of equal variance between the manual and LiDAR, with
only 13% significantly different in variance through the snow melt period (full results table in appendix
D).

Figure 5.18 KDE plot for 2022 showing the variability of snow depth for old growth plot ST1OG7. None of the survey dates were
significantly different in snow depth variability between the in-situ and LiDAR

The majority of the regenerating and juvenile sites were significantly different in variability between the
manual and LiDAR measurements (89% and 84% of all plots respectively). In each case, the LiDAR
snow depth variability was greater than the manual observations (Figure 5.19).

75

Figure 5.19 KDE plot for the regenerating forest plot ST4RF1 for 2022. The snow depth distribution was significantly different across all
phases between the in-situ and LiDAR. For the 24th of June survey, the 33 in-situ measurements not visible on the graph recorded zero
snow.

5.5 Bare Earth DEM
To assess bias which may be inherent to the bare earth acquisition, the snow depth analysis was re-run
with an alternative bare earth model for Russell Creek, a BCTS acquisition flown in 2017. Here, the
manual to LiDAR snow depth validation results are compared between both models to assess whether
the errors differ by cover type. Additionally, the LiDAR ground point coverage and density per each plot
is compared between the bare earth models.
Validation results for the weather station snow courses, conducted on stable ground with no canopy
cover (open), are relatively uniform using either the 2017 or 2023 bare earth DEMs (Figure 5.20).
However, distinctly different errors were produced within the juvenile, regenerating and dense old
growth forest when processed with the 2017 bare earth model.

76

Figure 5.20 Boxplot comparing the mean difference of plot snow depth for the manual observations and LiDAR across both years, with
open cover representing the weather station snow courses. Blue boxes represent results of the 2017 BCTS model, and the orange boxes are
of the 2023 bare earth model. The box shows the quartiles of the dataset, with whiskers extending to points within 1.5 x the inter-quartile
range. Observations outside of this range are displayed as points.

The juvenile plots had the greatest MD range between plots with the 2023 BCTS model (-27 to 134 cm),
whereas with the 2017 BCTS model, it had the lowest MD ranges across the cardinal cover types (-13 to
29 cm). The juvenile forest low MBE was consistent for both of the plots for the 2017 model, versus the
inconsistent results with the 2023 model. Within the regenerating forest plots, bias shifted from
generally positive with the 2023 bare earth model (median MD = 45 cm), to predominantly negative
(median MD = -30 cm) with the 2017 model, though both had a wide range in error (>80 cm).
Error increased in the dense old growth forest with the 2017 bare earth model, with LiDAR overestimating the snow depth in the range of 25 – 86 cm across all acquisitions. However, the old growth
with larger canopy gaps did not over-estimate snow, and the MD was the most consistent between
models of all the cardinal cover types.

77

Overall, the percentage of the plots with valid ground return points is greater in the 2017 model than the
2023, but to varying degrees across the different cover types (Figure 5.21). The greatest discrepancy
occurs in the juvenile forest plots, with ground point returns is on average 95% coverage with the 2017
model – akin to the alpine - and drops to 35% with the 2023 model. The standard error for the juvenile
is also much greater, due to the inconsistency between the two plots. The regenerating forest plots are
the lowest proportion (24 – 25%) and most consistent result between both models.

Figure 5.21 Bar graphs to show the comparison of the two bare earth models - 2017 BCTS (hatched) and 2023 ACO (non-hatched)
showing the total ground points per plot, broken down by point density per meter squared (top) and the proportion of the plot which had
valid ground returns (bottom). Bars represent the average of the cardinal plots by each cover type, and error bars represent standard
error

With the exception of the juvenile forest, the 2023 bare earth model had more total ground points per
plot than the 2017 model. Despite a greater proportion of data gaps, the point density was higher in the
2023 model, resulting in overall a higher number of ground returns per plot. Notably, the dense old
growth forest, which over-predicts snow depth in the 2017 model, is almost entirely comprised of single
78

point returns in the 2017 model. The juvenile forest is the exception; total ground point returns for the
juvenile forest plots are approximately double within the 2017 model versus the 2023.

5.6 Watershed Snow Storage
Total snow volumes were calculated for the Russell Creek watershed with both the 2017 BCTS and
2023 ACO bare earth models using interpolated (gap-filled) rasters (Figure 5.22). Using the 2023 bare
earth model resulted in greater snow volumes (3-6%) across all LiDAR acquisitions.
Consistent with the meteorological observations (Sec. 5.1), peak snow storage in 2022 was captured in
the April 4th acquisition, and substantial snow volumes persisted in the watershed until the last survey on
June 24th. In 2023, peak snow storage was captured on the March 3rd survey, with significant melt (76 –
82% of snow volume) occurring between the April 4th and May 25th acquisitions. The watershed was
largely snow free, with intermittent snow in the alpine, for the final acquisition on June 12th 2023.
Overall, peak snow storage captured in 2022 was 28% greater than in 2023. For the 2023 snow surveys,

79

Figure 5.22 Comparison of the total snow storage in the Russell Creek watershed between the 2017 BCTS model, and the 2023 ACO model
for each snow survey, broken down by forest cover type. Grey bars for the 12th June 2023 represent that the snow volumes were negative
overall.

5.6.1 Snow Volumes by Elevation
Across all acquisitions, peak snow volumes were observed between 1200 and 1400 meters (Figure 5.23).
In 2022, more snow was stored below 800 meters, (20 – 33%) and stayed through until the June
acquisition. Contrasted with 2023, from the April 4th survey onwards, less than 9% of total snow is
stored in the lower 800 meters. The choice of bare earth model had negligible impacts on the distribution
of snow volume by elevation, with the exception of the 1400 – 1700 m range.
80

Figure 5.23 Snow volume by elevation band for the 2017 BCTS model (red) and the 2023 ACO model (blue), showing the 2022 (top) and
2023 (bottom) LiDAR surveys

The pattern was consistent between all acquisitions as to which bare earth model predicted more snow
per elevation band. With the exception of two elevation bands, the 2017 bare earth model predicted
more snow than the 2023 model across all acquisitions, but not by a substantial volume overall.
However, as the 2023 bare earth model predicted a much greater (28 – 44%) volume of snow in a subalpine and alpine elevation band (1400 – 1700 m), overall the 2023 total snow storage was consistently
higher.

81

Table 5.10 Snow storage by elevation band, with percent of overall snow stored in the Russell Creek watershed shown in brackets,
between the two different bare earth models (2017 vs 2023), for the 9th of March 2022 and 10th March 2023 acquisitions. Highlighted
values indicated which model predicted the greater snow volume per elevation band.

Date Flown
Bare Earth

3/9/2022
2017, BCTS
2023, ACO

200 - 400
400 - 600
600 - 800
800 - 1000
1000 - 1200
1200 - 1400
1400 - 1700

Total

3/10/2023
2017, BCTS
2023, ACO

m³ (percent of total snow volume)

m³ (percent of total snow volume)

432,489(1.3%)

438,489(1.3%)

103,995(0.3%)

120,090(0.3%)

3,327,715(10.2%)

2,879,070(8.4%)

2,417,108(7.2%)

1,829,453(5.3%)

3,425,212(10.4%)

3,314,036(9.7%)

3,254,608(9.8%)

2,925,779(8.5%)

4,706,684(14.3%)

5,126,971(15.0%)

5,490,746(16.4%)

5,883,739(17.1%)

5,641,676(17.2%)

5,231,501(15.3%)

7,120,985(21.3%)

6,715,517(19.5%)

9,030,259(27.5%)

8,791,304(25.7%)

10,044,987(30.1%)

9,755,978(28.3%)

6,237,263(19.0%)

8,404,252(24.6%)

4,952,145(14.8%)

7,211,324(20.9%)

32,801,298

34,185,623

33,384,574

34,441,880

5.6.2 Snow Volumes by Forest Cover
The majority of snow in Russell Creek is stored in the alpine and old growth forest (Figure 5.22). The
alpine makes up only 6.8% of the watershed by area, but contains 32 – 58% of the total snow storage,
regardless of the bare earth model used (Figure 5.22). Though total snow volume decreases, the
proportion of snow stored in the alpine increases steadily by acquisition date, as the lower elevation
portions of the watershed loses snow. Old growth forest is the dominant forest type at higher elevation
where a deeper snowpack develops (see Figure 3.3), so it contains a higher relative proportion of the
total snow volume.

82

Table 5.11 Snow storage by cover type, with the percent of overall snow stored in the Russell Creek watershed shown in brackets, between
the two different bare earth models (2017 vs 2023), for the 9th of March 2022 and 10 th March 2023 acquisitions. Highlighted shows which
model predicted more snow per cover type

Date Flown
Bare Earth

3/9/2022
2017, BCTS
2023, ACO

3/10/2023
2017, BCTS
2023, ACO

m³ (percent of total snow volume)

Alpine
Old Growth
Mature Regen.
Regenerating
Juvenile
New Harvest
FSR
Total

m³ (percent of total snow volume)

11,814,549(36.0%) 13,712,941(40.1%) 11,053,398(33.1%)
13,007,977(39.7%) 11,230,885(32.9%) 14,129,254(42.3%)
1,299,083(4.0%)
1,070,746(3.1%)
1,131,221(3.4%)
4,823,783(14.7%) 5,947,794(17.4%)
4,816,056(14.4%)
1,572,433(4.8%)
2,059,789(6.0%)
1,883,514(5.6%)
195,834(0.6%)
73,387(0.2%)
269,670(0.8%)
87,639(0.3%)
90,081(0.3%)
101,461(0.3%)
32,801,298
34,185,623
33,384,574

12,968,040(37.6%)
12,200,555(35.4%)
763,779(2.2%)
5,934,958(17.2%)
2,358,069(6.8%)
116,560(0.3%)
99,919(0.3%)
34,441,880

Differences in snow cover volumes and percentages of total snow volume are consistent between all
acquisitions in terms of which bare earth model is used (Table 5.11). The biases observed at the plot
scale (Sec 5.4.1) were present when scaled to the watershed. In juvenile forest plots, the use of the ACO
2023 bare earth model results in overestimated snow depths, whereas with the BCTS 2017 model there
is good agreement between manual and LiDAR-derived observations of snow depth. When scaled to the
watershed, juvenile forests appear to contain 33% more snow (e.g. 1.5 x 106 m3 vs 2.0 x 106 m3 in
March 9 2022 survey) using the 2023 ACO bare earth model, though the overall impact on total snow
volume is ~1%.
In contrast, snow depths in old growth forest are likely overestimated when using the 2017 BCTS bare
earth model versus the 2023 ACO bare earth model, which had a much lower MBE (Sec 5.6). At the
watershed scale, the use of the less accurate 2017 bare earth DEM in old growth forests results in
substantially greater relative contributions to overall snow volumes (39.8 vs 33.1% in 2022). The choice
of bare earth DEM, and their respective accuracies in different forest types, thus plays an important role
in determining the relative importance of different forest cover types for total snow storage.

83

In other regions of the Tsitika watershed outside of Russell Creek, forest harvesting has recently
occurred at higher elevations, where snow accumulates to greater volumes. The following chapter will
use this recent forest harvesting as an experiment to assess the changes which occurred to both the
ground surface model and snow distribution pre- and post-harvest.

84

Chapter 6 : Forest Harvest Experiment
6.1 Introduction
Forest harvesting can change snow accumulation and melt due to removal of the canopy which alters the
energy balance, generally increases wind and eliminates interception of snow (Boon 2009; Varhola,
Coops, Weiler, et al. 2010; Winkler et al. 2015). Current forest practices also change the ground
characteristics, with slash piles, woody debris and roads potentially changing the elevation. Further,
after harvest low lying shrubs can become a dominant ground cover as the stand regenerates. Though
numerous studies have investigated the impact of forest harvesting on snow storage (e.g. Swanson et al.
1986, Winkler et al. 2015), most focus on the post-harvest condition, rather than the change which
occurs before and after disturbance. As well, studies which have looked and pre- and post-disturbance
changes to the snowpack and hydrological response pre-dominantly are based in cold climate
watersheds which develop shallower snowpacks than coastal British Columbia (Richardson et al. 2023).
Forest harvesting in coastal British Columbia is occurring at higher elevations, of which the impacts to
snow accumulation and melt are poorly understood.
Active forest harvesting within the Tsitika watershed which occurred during the span of the LiDAR
campaign presents a unique opportunity to assess how snow distribution changed pre-and postharvesting, and through the accumulation season. As described in section 2.3, paired catchment studies
which compare water yield pre and post disturbance through control and treatment watersheds, have
been the longstanding method to assess impacts of forest disturbance (Richardson et al. 2023). Here, a
similar control-experiment before-after approach was taken, in which cut blocks harvested in 2022 were
designated “treatment” sites, and adjacent forest stands, not harvested, were assigned “control”. Then,
the change in snow storage and water equivalence between 2022 and 2023 for both the control and
treatment areas were compared to evaluate the impact of forest harvesting.
85

As demonstrated within the previous chapter, an accurate snow-free surface model is critical to quantify
snow depth variability and snow storage (see sec 5.5 Bare Earth DEM). An advantage to recently
harvested areas is vegetation is yet to regrow and trees have been removed and therefore provide
relatively good conditions for LiDAR. However, the debris which remains after forest harvesting may
pose challenges. This study area has ACO acquired bare earth models pre- and post-harvest, which
provides the opportunity to assess whether the bare earth model changes from canopy covered to
harvested conditions. The first objective of this analysis was to assess whether the LiDAR snow-free
ground surface model changed after harvesting. The results of the ground surface analysis were then
integrated into the second objective: quantification of snow distribution change that occurred after forest
harvesting. This section focused on snow volume and water storage changes at the plot scale, and used
LiDAR derived snow depths to quantify the change. This also provided an opportunity to assess how
these changes impact overall snow storage at the watershed scale.

6.2 Methods
In June to August 2022, forest harvesting (16 hectares) occurred in the southern portion of the Tsitika
watershed adjacent to the Mount Cain ski hill, within the Cain Creek catchment (634 hectares). LiDAR
snow surveys acquired in 2022 therefore reflect the pre-harvest conditions, and surveys acquired in 2023
reflect post-harvest conditions. Additionally, the ACO acquired bare earth models for this region in 2020
(pre-harvest), and in 2023 (post-harvest). Meteorological data were collected at the Lower Cain (50.225
° N, -126.356 ° W, 1260 m) and Cain Ridge Run (50.227 ° N, -126.352 ° W, 1330 m) weather stations,
situated approximately 1.2 km south and 100 - 200 meters higher in elevation (Figure 6.1).

86

(b)

Figure 6.1 Site maps for the natural experiment with a) showing the AOI relative to the Tsitika watershed, northern Vancouver Island, and
b) scaled to the AOI, showing the road, control, and treatment area used within the analysis and c) the extent of Cain creek catchment
inclusive of the Lower Cain (1260 m) and Ridge Run (1330 m) weather stations

87

A relatively homogenous 30ha area of interest (AOI) where harvesting occurred between the 2022 and
2023 snow survey was delineated (Figure 6.1b). Harvested areas within the AOI were identified through
ortho-imagery collected by the ACO, and were designated as the treatment sites. Only areas which were
harvested after June 24th 2022 - the last LiDAR acquisition of the season - were included. Unharvested
areas were the control sites. A minimum 10 meter buffer was applied to the boundaries of the control
and treatment to each other, as well as the road (Figure 6.1c). The buffer distance was based on current
literature on forest edge effects as well as canopy height on snow distribution (Currier and Lundquist
2018). Treatment and control areas were summarised by the metrics of slope, elevation, aspect, forest
canopy density and stand area (Table 6.1). One section of road was built prior to the pre-harvest bare
earth collection in 2020, and was used to assess any change to stable terrain between the two bare earth
surfaces.
Table 6.1 Topographic and canopy metrics for the control, treatment and road area which was used within this study. Metrics for slope,
canopy cover, and height represent the average value.

Plot

Control
Treatment
road

Elevation
range

Slope

Canopy Cover

Canopy Height

Area

[m]

[degrees]

[%]

[m]

[ha]

1039 - 1361

27

88

18

14

1080 - 1353

28

91

25

16

1077 - 1124

4

0

0

0.54

Aspect
South
South
South

6.2.1 Snow free surface assessment
A primary goal of the snow free surface assessment was to identify whether the LiDAR ground models
changed after harvesting. Post-harvest surfaces are a complex mix of slash piles, stumps, and debris
which may influence the ground surface model if the LiDAR pulse cannot penetrate to the ground. This
88

will have direct implications as to a) whether a bare earth LiDAR model captured prior to harvesting is a
viable product to measure snow once harvesting has occurred or vice versa and b) whether the bare earth
model post-harvest is capturing the ground, or surface debris.
The road constructed prior to the 2020 bare earth acquisition was used to first establish whether the two
bare earth models had good agreement on the stable terrain. If there was no significant change between
the elevation values between the two bare earth models, this would be indicative of no systematic error
between the models due to timing of acquisition and potential instrument bias. Average and standard
deviation of the road surface elevation values were calculated for both the 2020 (pre-harvest) and 2023
(post-harvest) models. As this portion of the analysis employs a point to point approach to compare
LiDAR pixels, different error metrics were able to be calculated which were not possible with the
unequal sample sizes of the plot-scale analyses. Error metrics were mean bias error (Eq. 6.1), mean
absolute error (6.2), and root mean squared error (6.3), where z is the elevation raster value.
ଵ

∑௡௜ୀଵ(‫ݖ‬ଶ଴ଶଷ − ‫ݖ‬ଶ଴ଶଶ )

(6.1)

‫ = ܧܣܯ‬௡ ∑௡௜ୀଵ|‫ݖ‬ଶ଴ଶଷ − ‫ݖ‬ଶ଴ଶଶ |

(6.2)

‫= ܧܤܯ‬

௡
ଵ

ଶ

ଵ

ܴ‫ = ܧܵܯ‬ට ∑௡௜ୀଵ(‫ݖ‬ଶ଴ଶଷ − ‫ݖ‬ଶ଴ଶଶ )
௡

(6.3)

Then, the same process was applied to the control and treatment sites. As the ground surface
distributions were not normally distributed, the Levene’s and Mann-Whitney U test were selected for
testing whether the variability and mean of the surfaces were significantly different.
Additionally, the snow depth rasters in which the AOI was snow-free were used to confirm the existence
or absence of changes in ground surface elevations pre- and post-harvest within the treatment and
89

control areas. Ortho-imagery captured coincidently with the LiDAR was used to verify whether there
was snow in the AOI. Then the “snow depth” i.e. no snow pixel values were extracted for the control
and treatment sites. If the treatment site snow depth distributions were centered at zero in the pre- and
post-harvest snow free surveys, this would indicate the LiDAR bare earth DEM is unchanged by forest
harvesting. If the treatment site elevations differ after the treatment, indicates a change to the ground
surface due to harvesting.

6.2.2 Snow distribution and volume change
Pre-harvest snow depth distributions were compared between control and treatment plots to confirm
snow depths were similar before harvesting. Snow depths were then converted to snow water equivalent
(SWE mm) by multiplying the snow depth pixels by a modelled snow density pixel value. Snow
densities were modelled using a random forest regression with the key factors being snow depth,
elevation, and meteorological variability, which is a proxy using hourly weather station data (Bisset et
al, 2024 (in prep)). A relative change assessment was conducted by comparing the snow-on acquisitions
from 2022 and 2023. This analysis is based on the control area remaining stable, despite the date of the
LiDAR acquisitions and the snowpack changing between the two years. Three surveys in both 2022 and
2023 were acquired when the average snow depth exceeded 1 m at the weather stations and snow visible
in the ortho-images. These three surveys were defined as phases, identified to compare change between
the two years (Table 6.2). A ratio was calculated for the control and treatment sites, to reflect the
relative increase or decrease in SWE from 2022 to 2023.
Table 6.2 LiDAR surveys for 2022 and 2023 paired for comparing snow distribution change pre-harvest (2022) and post-harvest (2023).
Weather station snow depth data from Lower Cain (1260 m) included for each survey

Phase
1
2
3

LiDAR Acquisition Date
2022
2023
9th March
10th March
11th April
4th April
st
21 May
29th April

Lower Cain: snow depth [m]
2022
2023
1.69
2.21
2.27

90

1.74
1.94
1.76

6.3 Results
6.3.1 Weather Station Overview
The two weather stations located on Mount Cain - Lower Cain (1260 m) and Ridge Run (1330 m) provide important context to the snow distribution results and the differences in the 2022 and 2023 snow
seasons. In 2022, peak snow depth occurred on the 12th of April, 266 cm at Lower Cain and 303 cm at
Ridge Run (Figure 6.2). The start of snow melt was defined as the date at which an average loss of .5 cm
of snow depth occurred for a minimum of 14 day period. Sustained snow melt started on the 19th of May
for both stations lasting 42 -45 days, with an average 6 cm snow loss per day. Both weather stations
snow depth sensors were recording no snow as of July 3rd. Relative to the LiDAR acquisitions, peak
snow occurred 2 days after the 10th of April acquisition. Based on the orthoimagery from the June 24th
2022 LiDAR acquisition, the snowpack within the area of interest was intermittent and mostly snow
free, with snow patches still present in natural openings within forested areas.

91

Figure 6.2 Overview of the Lower Cain (1260 m) weather station air temperature (red) and snow depth (blue) measurements for the 2022
and 2023 peak snow and snow melt season.

The 2023 snow season was characterized by an earlier melt, with both stations being snow free by May
27th, over a month earlier than 2022. The weather stations recorded a shorter sustained snowmelt period
and lower peak snow depths, 34 and 47 cm less at Lower Cain and Ridge Run stations respectively. By
the May 25th LiDAR survey, the natural experiment AOI was snow-free. Peak snow occurred
approximately halfway between two LiDAR acquisitions, on April 18th 2023.
The first three LiDAR acquisitions of both the pre-harvest 2022 acquisitions, and the post-harvest 2023
acquisitions had a measurable snowpack, so will be referred to in this analysis as the snow-on surveys.
There were three surveys near or completely snow-free within the experiment AOI: the 24th of June
2022 and the last two surveys in 2023 (25th May and 12th June). All of the LiDAR surveys will be
included in the snow distribution assessment, and only the snow-on surveys will be incorporated into the
water storage analysis.

92

6.3.2 Ground surface change
The road within the harvested block showed good agreement in elevation between the 2022 and 2023
bare earth models, with an RMSE value of 0.02 m (Table 6.3). A histogram of the residuals showed the
peak close to zero, indicating minimal change between the two bare earth surfaces (Figure 6.3). Thus,
the stable terrain within the AOI was consistent between both surface models.
Table 6.3 Results from pre and post-harvest ground surface changes, in which the control, treatment, and road plot types average and
standard deviation elevation were calculated, and error metrics mean bias error (MBE), mean absolute error (MAE) and root mean square
error (RMSE) were calculated between the 2022 (pre-harvest) and 2023 (post-harvest) bare earth models.

Plot type
2022
control
treatment

Average [m]
2023

Standard Deviation [m]
2022
2023

1171.93
1172.08
1179.03
1179.33
road
1107.31
1107.32
*
Statistically significant at a p-value of 0.05

97.38
69.38
13.04

97.43
69.25
12.95

MBE
[m]

MAE
[m]
0.15
0.3
0.01

0.37
0.4
0.12

RMSE
[m]
0.48
0.49
0.2

The variability of the ground surface showed good agreement for both the control and treatment sites.
The average elevation for the control and treatment area showed an increase from 2022 to 2023, with a
mean bias error of 15 and 30 cm, respectively. The RMSE was marginally higher for the treatment (0.49
m) versus the control (0.48m), though both were quite high overall (Table 6.3). The histogram of the
residuals shows a slight positive skew for the treatment. Neither the Levene’s test for equal variance nor
Mann-Whitney U test produced significant results for the control, treatment or road groups (Figure 6.3).

93

Figure 6.3 Histogram of the residuals for the control, treatment, and road sites between the post-harvest (2023) and pre-harvest (2020)
bare earth. A positive skew is indicative of the post-harvest elevation greater than the pre-harvest.

The second component to the ground surface analysis compared the calculated snow depth distributions
for the snow free acquisitions, which, if unbiased, should be centred around zero (Figure 6.4). The
distributions of “no snow" snow depth for the snow free acquisitions exhibit a negative bias for the preharvest treatment areas, suggesting that the process of forest clearing introduced a significant change to
the ground surface. For the 2023 post-harvest surveys, the control and treatment average snow depth
were within .05 m of zero.

94

Surface elevation differencing for two snow-free and post-harvest acquisitions (25th May 2023 and 12th
June 2023) with the average depth ranging from -0.05 to 0.03 m. The standard deviation for the forested
control sites was greater for both surveys than the harvested treatment.

Figure 6.4 Snow depth distributions of the three snow free acquisitions, ie. no snow, comparing the control and the treatment site
processed with the post-harvest (2023) bare earth model. The 24 th of June 2022 (left) shows a negative offset when using the harvested
ground model to measure snow in the treatment site which - as it is prior to harvesting – is forest covered. The 25th of May 2023 (middle)
and 12th June 2023 (right) show good agreement between the post-harvest ground surface model and the harvested treatment and forested
control sites

As the treatment sites were already harvested in the bare earth model, the 2022 treatment distribution
was a harvested ground model used to measure snow within forest. This resulted in a negative offset (0.21 m) in comparison to the control (-0.05 m).
When using the post-harvest 2023 bare earth model, snow depth was underestimated in the pre-harvest
treatment sites. This was consistent with the results of the snow-free bias assessment, which showed a
positive bias to the 2023 ground surface model. To accurately quantify snow storage change from
harvesting, a bias correction must be applied to the pre-harvest treatment snow depths. The difference
between the mean bias error for the control (0.30m) and treatment (0.15m) was differenced, and a 0.15
m increase was applied to the treatment sites snow depth measurements for 2022.

95

6.3.3 Snow Distribution and water storage
Prior to harvesting, the average snow depth for the control and treatment were closely aligned, with 0.04
m or less difference across all the acquisitions (Figure 6.5). The average snow depths ranged from 0.31
to 0.92 cm for both control and treatment. The snow depths increased for the first three acquisitions,
with the average snow depth peaking on the 21st May 2022 (0.91 – 0.92 cm).

Figure 6.5 Box plots of the snow depth (m) in the control and treatment sites for the pre and post-harvest conditions. The boxes (green
control and orange treatment) shows the quartiles of the dataset, with whiskers extending to points within 1.5 x the inter-quartile range.
Observations outside of this range are not shown on the graph

After harvest, the average snow depth was consistently higher for the treatment versus control sites. The
largest discrepancy in average snow depth occurred in the March 10th, 2024 survey, where the average
snow depth was 0.75 m greater in the treatment than the control. Between the April 29th and May 25th
2023 surveys, all the snow in both groups melted. Sentinel data were examined for dates between April
9th and May 25th, and while the approximate date for snow free status identified in the treatment area, it
was too difficult to see the snow cover in the forested treatment areas, thus it was not possible to
compare the bulk ablation rates. However, as the average snow depth for the 29th April survey was 76%
greater in the treatment than control sites, a greater quantity of snow was available to melt than prior to
harvest. Despite the observed similarities in distribution of the control and treatment sites pre-harvest,
96

and differences post-harvest, statistical testing produced significant differences in both the before and
after harvest scenarios.
After forest harvest, there was a notable shift in both the distribution and the total water storage between
the control and treatment groups (Figure 6.6). Isolating the control group, peak water storage across both
snow years surveys occurred on 21st May 2022. For the treatment sites, all of the post-harvest surveys
exceeded the water storage of the pre-harvest surveys. Therefore, despite the 2022 exceeding water
storage within the forest, the overall water storage was greater for the April 29th survey due to the
harvested areas.

Figure 6.6 Total water storage in cubic metres for the control (green) and treatment (orange) sites for the pre- and post-harvest years.

When comparing 2022 to 2023 snow years, the change in SWE was positive and greater for the
treatment group than control across all phases (Table 6.4). Phase 1 (March 9, 2022 and March 10, 2023)
had the greatest increase in SWE for both groups, though the treatment increased by over double that of
the control. The third surveys of each year resulted in the largest difference between the control and
treatment, as well as a direction change: whereas the control group decreased in SWE, the treatment
sites still increased compared to 2022.
97

Table 6.4 Snow water equivalence (mm) for both the pre-harvest (2022) and post-harvest years (2023), +/- standard error, for the three
snow on LiDAR surveys. Rate of change represents the increase or decrease in SWE from 2022 to 2023, for the control and treatment
blocks

SWE [mm]
Snow year
2022
2023
Rate of change
(percent)

Phase 1
(Mar 9 2022 : Mar 10 2023)
control
treatment

Phase 2
(Apr 11 2022 : Apr 4 2023)
control
treatment

Phase 3
(May 21 2022: April 29 2023)
control
treatment

125 +/- 0.44
216 +/- 0.31

95 +/- 0.32
398 +/- 0.34

192 +/- 0.45
224 +/- 0.36

187 +/- 0.38
375 +/- 0.36

340 +/- 0.52
244 +/- 0.43

317 +/- 0.54
428 +/- 0.52

1.7 (72%)

4.2 (318%)

1.2 (15%)

2.0 (100%)

0.7 (-28%)

1.4 (35%)

6.4 Discussion and Recommendations
Through a comparative experiment approach, this case study used LiDAR to quantify snow depth and
snow water volume change that occurred after forest harvesting. An important finding from this study is
the shift to the bare earth model after harvesting. The increase to the surface elevation in the cut block
suggests the debris from the harvesting caused the ground surface to overall register higher. A potential
source of this bias is the slash piles or logs are being classified as ground, causing an overall increase to
the ground elevation distribution. However, an analysis of the logs and slash compared to ground within
the treatment site showed the bias was consistent, and not caused specifically by log piles (appendix I),
suggesting that forest harvest caused a general increase to the ground surface, which is being reflected
with the bare earth LiDAR model. It is critical to address this change, as the underestimation of the preharvest treatment snow depth which occurred due to the use of the harvested ground surface model
would cause the relative change to snow depth after harvesting to be amplified. For example, without
the bias correction applied, snow depth change within the treatment sites would have been overestimated by 12 – 33%. This highlights the importance of using a ground surface model which reflects
the current state of the land surface to accurately measure snow. Any additional changes to the ground
surface which may occur (removal of slash for example) or on a greater temporal scale, the regeneration
of forest, may introduce error to measuring snow with LiDAR. Furthermore, this supports and adds to
results from the previous chapter, in which the regrowth of vegetation post-harvest reduced the accuracy
98

of the bare earth model, and subsequently the snow depth measurements. In this case, the removal of
vegetation reduced the accuracy of retrospective snow depth measurements by the introduction of a
positive (.15m) bias, with the post-harvest open ground surface being higher than the pre-harvest
forested surface. With the change to the ground surface addressed, the overall confidence in the LiDAR
snow depth measurements was high. First, both bare earth models were flown by the ACO within 4
years of each other, using the same LiDAR sensor and systems, with similar ground classification
algorithms employed (see sec 4.1). As discussed in the previous chapters, LiDAR snow depth
measurements in the open and old growth plots were the two most accurate cover types, with median
MBE of -0.9 m and -0.16 m respectively, as well as captured the variability. Despite no in-situ snow
depth data, the snow-free validation provided confidence that the snow on and bare earth surfaces were
in good agreement pre- and post-harvest.
Results of the snow distribution analysis follows the most prevalent finding from numerous studies of
this nature: snow will accumulate to greater snowpacks in the clear cuts than in the adjacent forest
(Berris and Harr 1987; Hudson 2000). Exceptions to this can occur include larger clear cuts under high
winds, due to evaporative losses in the opening and snow may be blown in to and stored within the
forest (Troendle and King 1985; Swanson et al. 1986). Studies have found maximum snow
accumulation occurs when the width of the cutblock is 2 to 5 times the height of the trees (Swanson and
Stevenson 1971; Troendle and King 1985), though these results were found in cold climate conditions.
In this case study, cut blocks fell within this threshold, and were on average between 4 times in width
(77 m) than surrounding tree height (18m).
Hydrological effects of forest harvesting are a combination of the harvested area and recovery, the
spatial distribution of the cut blocks, as well as the specific topographic characteristics of the watershed
(Ellis et al. 2013; Winkler et al. 2015; Johnson and Alila 2023). The 2022 harvested sites encompassed
2.5% of the Cain Creek catchment, though with prior harvesting which occurred, mostly within the last
99

10 years, 9% of the catchment was harvested overall. Snow water storage between the plots increased at
a higher rate in the treatment sites than control between 2022 and 2023, with increased rates respective
for control and treatment were 72% vs 318% for phase 1, 15% vs 100% for phase 2, and -28% to 35%
for phase 3. To infer the impact of the harvested sites on snow storage at the watershed scale, the
relative difference which occurred between the plots was scaled to the Cain Creek watershed, for
example, for the phase 1 (Mar 9 2022 – Mar 10 2023), the treatment SWE increased 40% more than the
control. This only resulted in a 1.6 , 1.02, and 1.36 % overall increase in water storage for the march
10th, April 4th, and April 29th LiDAR surveys within the Cain Creek watershed as a result of the forest
removal. This suggests that snow at this elevation makes up a proportionately smaller amount of snow at
the watershed scale. The relatively high elevation of the study area was beneficial to quantifying snow
depth change. At a lower elevation, we could expect less snow to accumulate both within the open and
the dense canopy, as warmer temperatures increase the amount of precipitation falling as rain rather than
snow.
One limitation of this case study is the inability to describe and compare the snow melt rates with
LiDAR, as for both 2022 and 2023, the snowpack within the AOI melted out in entirety between the
third and fourth LiDAR surveys. Based on the available sentinel imagery, weather station data, and
understanding of the snow energy balance, we can infer some of the likely impacts on melt rate. As the
entire AOI is southerly facing, it would receive greater net shortwave from the solar incident angle and
increased longwave inputs from the heating of tree trunks, therefore likely a higher melt rate than a
comparable site on a northerly aspect (Seyednasrollah and Kumar 2019). As captured with sentinel
imagery, snow melt rate in 2023 appeared to be quite rapid in the treatment sites, with full snow
coverage May 10th, intermittent snow coverage May 13th, close to snow free as of May 18th, and no
snow as of the LiDAR acquisition on May 25th. Notably, a 2021 cut block (3 hectares) on the northern
aspect at a similar (1054 – 1117 m) elevation, was still fully snow covered as of May 18th. Higher
100

antecedent snowpack was captured on the northern cut block vs southern aspect treatment, with a 1.45 m
and 1.2 m average snow depth respectively as of the April 29th 2023 LiDAR survey. Based on recent
studies in a maritime environment, the melt rate in the forest may have exceeded that in the open, as this
site was relatively warm, sheltered, and below treeline. Notably though, the average snow depth was
44% less under the canopy, therefore less snow was available to melt. The accelerated melt on south
aspect forest harvesting may increase the synchronicity of snowmelt within the watershed, ie. the higher
and lower elevation areas melting at the same time, which Pomeroy et al (2012) found contributed to
higher peak flows.
One recommendation is to increase the temporal scope of this analysis to include additional years of
LiDAR survey data, where harvesting occurred over a broader area within the watershed. As LiDAR
surveys in the Tsitika span back to 2020, an inclusion of the 2021 harvested areas may strengthen these
findings, as well as provide greater understanding of how the size, aspect, and spatial distribution of
forest harvesting impacts snow storage within a coastal watershed. An additional recommendation is to
introduce more frequent LiDAR surveys to capture melt rate, and compare these rates between the
forested and open sites. This case study demonstrated the utility of airborne LiDAR to provide robust
figures for the change in snow storage which occurs after forest harvesting.

101

Chapter 7 : Discussion and Recommendations
This study focused on snow distributions and potential sources of LiDAR error at the plot and watershed
scale, through analysis of vegetation, ground point returns, and the bare earth model. This discussion is
centred on several key themes:
1) The vegetation and surface characteristics of the juvenile and regenerating plots, which lead to
greater observed errors in LiDAR-derived snow depths in the complex surface and canopy of
harvested forest;
2) The choice of bare earth model is important: two independent bare earth models are similar in
stable terrain but show differences between forest cover types;
3) Confidence in overall snow volumes is high, despite the differences in the bare earth models, but
there are differences and compensating errors.

7.1 LiDAR Accuracy within complex ground cover
This study assessed LiDAR accuracy by comparing plot averages of LiDAR and manual snow depth
measurements - as opposed to point-to-point which is prevalent across most studies - in order to not
introduce known geo-positional error associated with a) imprecise position of manual measurement and
b) accuracy of the GNSS in the forest (Murgaš et al. 2018; Broxton et al. 2019). Mazzotti et al (2020)
precisely referenced their snow probes in the forest with specific trees which were identified on orthoimagery, but this approach would have been impractical within this study due to the density of canopy.
A drawback of the average snow depth approach used in this study is the inability to calculate and
compare error metrics reported across other literature, commonly root mean square error (RMSE), mean
absolute error (MAE), and mean absolute difference (MAD)(Table 7.1). Median of the mean differences
for this study were lowest for the open sites, the weather station snow courses and alpine plots (-9 cm for
102

both), followed by the old growth (-16 cm) and dense old growth (-15 cm) plots, and highest for the
juvenile (76 cm) and regenerating (55 cm) plots. Generally, error values from this study were similar or
slightly greater in magnitude in the open and old growth forest sites to other LiDAR studies, and
substantially greater in the harvested. However, it is important to note that no published studies have
endeavored to measure snow with airborne LiDAR in the complexity of forest stands which were done
so in this study.
Table 7.1 Error results from airborne LiDAR snow depth studies which evaluated LiDAR accuracy using in-situ measurements and in

Study

Airborne LiDAR type Error Evaluation

Currier et al, 2019

Airplane

Geissler et al., 2023

UAV

Jacobs et al, 2021

UAV

Hopkinson et al., 2012

Airplane

Mazzotti et al, 2019

Airplane

Forest error metrics:
- Mean difference: 6 cm
- RMSD: 8 cm
- MAD: 7 cm
Equal error for open and forested terrain:
- RMSE: 2 cm, and MAE:1 cm
Greater error values in forest than open.
Field (open): MAD: 0.96 cm, RMSE: 1.22 cm
Forest: MAD: 9.6 cm, RMSE:10.5 cm
Higher mean difference in forest than open.
Mean differences were:
- Open area, high elevation: 7 cm
- Open area, alpine: -4 cm
- Forest: 13 cm
Forest metrics:
RMSE: 6 cm mean, 3cm standard deviation

forested terrain.

When compared with ground observations, LiDAR estimates of snow depth were most accurate in stable
terrain with no ground vegetation or canopy cover, demonstrated by the consistent low MD attributed to
the weather station snow courses. However, in terms of the snow storage, roads are an inconsequential
area of Russell Creek (<1%) and contributed less than half a percent of snow volume. The weather
station snow course validations confirmed that bias found at the cardinal plots was likely due to
vegetation misclassified as ground in the bare earth model.
103

Uncertainty increased in complex ground cover and vegetation, with notably the harvested plots
(regenerating and juvenile) responsible for the greatest mean differences. Ground point coverage within
the regenerating plots was low within both of the bare earth models. Even in snow-on conditions, ground
point density within the regenerating forest plots was lower than the other cardinal cover types, with 29
– 50 % of plot with no ground point returns, and less coverage later in the season. This was not
necessarily a function of the canopy coverage, as the dense old growth forest plots had the highest
canopy coverage of all the cardinal plots (90%), yet the snow-on point coverage of the plots were
considerably higher (69 - 97%). A notable difference between the two cover types were the tree
characteristics: tall (~ 20 m) trees, with a small ratio of branches to tree height, and no branches close to
ground in the dense old growth, versus the smaller (~ 8 m) trees with a comparatively dense crown and
low to ground branches in the regenerating plots (Figure 7.1). Similar results were found in a study by
Mazzoti et al (2019), in which dense spruce tree stands had more data gaps in the LiDAR snow surface
and ground returns in comparison to less dense stands. Hopkinson et al. (2004) suggested a speciesspecific bias correction, but within the Russell Creek watershed, which encompasses a range of
harvested and old growth forest, a stand age component may render more accurate results.

Figure 7.1 Images from the regenerating stand (A) and dense old growth (B), showing notable difference in tree branch height in relation
to the ground.

104

In addition, complex ground vegetation increased LiDAR error at the plot scale by obscuration of the
ground surface, which resulted in both less and low quality LiDAR ground point returns. In the
harvested plots, the ground vegetation was tall (avg: 116 cm) and variable (std: 45 cm), and as the
average tree height was ~ 3 meters, there was not a clear distinction between the two. Again, this
resulted in fewer ground point densities than in the over cover types. The ground vegetation present in
the old growth forest appeared to pose less of a problem to LiDAR penetration to ground, based on the
relatively high ground point density and coverage, despite the vegetation being not inconsequential in
the larger canopy gaps (approx. 50 cm). This may be a factor of both a) the vegetation being relatively
uniform in height, and b) the large difference in height between the canopy and the ground vegetation in
the old growth forest plots.
The ambiguity of ground in the harvested forests introduced uncertainty in both the LiDAR and in situ
measurements. Large air pockets from vegetation, as well as slash piles and stumps were sources of
error within the in situ measurements. Air pockets could be identified when using the snow tube, as the
core of snow may be substantially less than the depth measured, and measurements were retaken. The
LiDAR interpretation of ground is also important to consider, as it may not capture how snow
accumulates in the landscape. For example, if the LiDAR model classifies a slash pile as vegetation, this
will result in a substantial over-estimation of snow. Geissler et al (2020) observed fallen logs classified
as vegetation in the snow-off point cloud resulted in a large systematic error in snow depth data. Though
in this study the over-estimation of snow in the juvenile forests was not limited to areas of stumps and
slash piles, this may be a contributing factor to the error observed.
Overall, this study found that a good proxy for LiDAR snow depth accuracy was the ground point
density of the snow-free surface. Plots with both a) good ground point coverage and b) higher density
(>1) of returns generally had a lower mean difference. An exception to this is the alpine plots in 2022,
which despite being free from complex vegetation and canopy coverage, had greater range in mean
105

difference than anticipated (Figure 5.9). However, this may be an artefact of the less frequent validation
of the plots and the more variable snowpack when compared to the other cardinal plots. A similar result
was found in Tinkham et al (2014), with the highest error attributed where the snowpack was highly
variable, which in their case was the edge of the canopy.

7.2 Impact of the Bare Earth Model
The bare earth DEM and how it classifies ground and vegetation appears to be a significant source of
error. This is demonstrated when comparing the 2023 ACO bare earth to the 2017 BCTS acquisition.
Both bare earth models were accurate on stable terrain, demonstrated by the weather station snow
courses low mean difference, but biases in the cardinal plots differ (Figure 5.20). The juvenile forest
plots had the greatest difference: the highest MD when processed with the ACO 2023 model, and the
lowest MD with BCTS 2017 model.
As discussed in the previous section, the juvenile plots were characterized by tall and complex ground
vegetation, and small trees. This resulted in low ground point density and, most importantly, large data
gaps in the 2023 bare earth model, though not in the 2017 BCTS model. In snow-on conditions, juvenile
plots were consistent between both bare earth models as very high snow surface point density, due to
snow coverage of all vegetation and debris.
Several possible theories as to why the 2023 model was substantially lower in bare earth ground point
returns in the juvenile plots than the 2017 model are: 1) more vegetation has grown in the 6 years
between the bare earth acquisitions, and 2) vegetation was greened up at the time of the 2023
acquisition, thus more challenging for the laser pulse to penetrate thorough to ground. However, this is
not a product of the top of vegetation being misclassified as ground as observed in other literature
(Kucharczyk et al. 2018), as in this scenario, the overestimation of the snow-free surface elevation
would result in LiDAR under-estimation of snow depth in the juvenile forest plots. It is not clear why
106

the juvenile plots under-estimated surface elevation in the 2023 model, and this result has not been
produced in other studies.
Support for the vegetation green up theory is provided with an additional snow-free surface experiment:
the May 23rd, 2023 LiDAR acquisition. At the time of the May 2023 survey, the juvenile and
regenerating plots were snow free, and the vegetation was free of foliage (Figure 5.5). To test whether
vegetation green-up is a factor in error, the May 2023 survey was used in the place of a bare earth model
to validate only the harvested plots. With the experimental May bare earth, MD was much reduced for
both the juvenile and regenerating plots, and was in fact the lowest of the bare earth models tested
(Figure 7.2). For the Juvenile, range of MBE with the experimental May 2023 bare earth model (-16cm
to +18 cm) was consistent with the lower bias error value from the 2017 BCTS model. The regenerating
plots also were the lowest bias error recorded with the May 2023 snow-free surface. This suggests that
for areas of high and dense ground vegetation, such as regenerating and juvenile forests, timing the bare
earth acquisition with leaf off conditions will significantly reduce error. It also points to the challenge of
using a static bare earth model when vegetation not only changes seasonally, but also annually. One
possible solution to this problem would be creating a mosaic bare earth using multiple surveys and
selecting areas with a minimum number of returns. This would not reduce the problems associated with
misclassification of vegetation as ground, but would limit errors associated with interpolation over large
distances in complex terrain.

107

Figure 7.2 Mean difference of the cardinal plots from both 2022 and 2023, for only harvested plots (the regenerating and juvenile), using
three different bare earth models: 2017 BCTS survey acquired in September 2017 (orange), 2023 ACO survey (blue), and the substitute
bare earth model: the May 2023 snow survey which was snow-free in the juvenile and regenerating stands (green). Individual points
represent each plot validation across 2022 and 2023.

The algorithm used to classify LiDAR ground returns can introduce error due to misclassification,
especially within more complex terrain and vegetation (Tinkham et al. 2011). Both bare earth models
used an iterative TIN densification approach to classify ground, but the specific algorithms used did
differ between the two (see sec 3.4). A study by Moudry et al. (2018) found the LASground algorithm to
perform well overall, but in complex areas, was more likely to omit a ground point rather than falsely
classify vegetation, which resulted in low ground point density. Within this study, the omission of
points through ground classification may be the cause of the low ground point returns in the 2023 ACO
model harvested forests, although it is not clear whether the LiDAR pulse was able to penetrate the
vegetation. Tinkham et al. (2011) found between a comparison of two open source algorithms (MCC
and BCAL) the greatest deviation (40 cm) in the ceanothus, (with the BCAL being the more accurate),
108

and suggested a multiple algorithm approach by cover types may produce the most accurate ground
surface.

7.3 Snow Distribution: Plot and Watershed Scale
7.3.1 Plot scale distribution
LiDAR and ground-based observations of snow depth distributions are similar for all forest cover types
with the exception of harvested forests. LiDAR snow depths are more variable in open and wind
exposed sites (e.g. alpine) than the sheltered sites (e.g. dense old growth). In the old growth forests with
larger canopy gaps, both manual and LiDAR snow depth observations show deeper snow in the
openings, and shallower snow under and at the edge of canopy (see Appendix E: Old Growth Forest:
Canopy Snow Depth Variability).
Snow depth variability did not appear to change through the melt season, for manual or LiDAR
measurements. Notably, neither the 2022 nor 2023 LiDAR campaigns were well suited to answer
questions pertinent to snow melt rates: in 2022 the snowpack sustained until the last survey with
minimal melt, and 2023 had fairly consistent snowpack through the first three surveys, and rapid melt
between the third and fourth acquisition. Despite being an important metric, there are limited studies
which evaluate LiDAR ability to capture snow depth variability.

7.3.2 Watershed snow storage
In Russell Creek, the majority of the snow was stored above 1000 m (58 – 82%), with 14 - 28% of the
total volume of snow found in the 1400 – 1700 m elevation range. On Vancouver Island, the low to mid
elevation zones typically develop a shallow, intermittent snowpack due to the mild winter temperatures
and higher frequency of rain-on-snow events. Therefore, it is critical to understand not only LiDAR

109

snow depth accuracy at higher elevations, which in the case of Russell Creek, were the old growth forest
and the alpine, but the mid-range elevations too.
The comparison of watershed snow storage between both bare earth models produced relatively similar
results overall, with less than 5% difference in the total snow volume estimates. This is likely an effect
of the two bare earth models countering biases (over and under predictions of snow depth) by different
cover types which, to a degree, cancelled at the watershed scale. The pattern of which bare earth model
predicted more snow within each cover type was a) consistent with each survey acquisition and b)
generally reflective of the comparative biases seen at the plot scale. This gives some confidence to
which bare earth model may be more accurate within the different cover types, but still difficult to assert
which model represents the “true” snow water storage.
In terms of the cover types which store the most snow - the old growth and alpine - there is evidence
within both the plot error comparison and snow free watershed storage to suggest the ACO 2023 model
may be the more accurate representation within the old growth, but likely the “true” snow volume is
between the 2017 and 2023 models for the alpine. For the dense old growth, the reverse pattern was
observed in the snow free acquisition in 2023: median MBE of +62 cm for BCTS 2017 bare earth
model, vs -13 cm for the ACO 2023 bare earth model. The inconsistencies in biases between the bare
earth models may be attributed to differences in: 1) the acquisitions and how well the LiDAR pulses
were able to penetrate to ground, 2) the alignment of the model, and 3) the LiDAR processing and
classification. As detailed in the following section (re: study limitations), caution should be applied for
upscaling the errors observed at the plots to the watershed. However, due to the consistent relationship
between the two bare earth models - which also were consistent with the plot scale bias errors - it is
probable that the BCTS model under-predicted snow in the alpine, and over-predicted snow in the old
growth forest.
110

Within the alpine plots, the BCTS model had a MBE of – 50 and -5 cm error in the snow free acquisition
which may be a result of a positive elevation bias in the bare earth model ground surface. In support of
this, the watershed snow storage was entirely negative for the alpine on the 12th of June 2023 survey for
the 2017 BCTS model (see Appendix H: Watershed Snow Volume), despite visible snow in the orthoimagery. The 2023 ACO model generally showed good agreement on the snow free ground surface
within the plots (MBE -3 and -4 cm), but there are several substantial data gaps (450 to 13,900 m2)
isolated to steep (average slope: 51 degrees) rocky peaks, which resulted in extremely high (average: 20
m) and clearly not representative snow depth values due to the interpolated data. Though these areas
make up a relatively small proportion (4%) of the alpine, the snow volume accounted between 22 – 32%
of the total snow in the alpine. In the same areas of missing data in the ACO 2023 version, the 2017
BCTS model under-estimated snow depth, notably in 2023(average: -8 m), and snow volumes were
negative, which may be indicative of an alignment issue. This demonstrates both a) LiDAR inaccuracies
on steep slopes, which are well documented within other studies (Hodgson and Bresnahan 2004a;
Tinkham et al. 2012), and b) a limitation of the plots to assess accuracy at the watershed scale.
The over estimation in the harvested forests observed at the plot-scale with the ACO 2023 model did not
have a significant impact (<5%) on total snow volume within the watershed, due to the proportion and
elevation of the regenerating and juvenile stands within Russell Creek. The juvenile and regenerating
forest stands had the greatest uncertainties at the plot scale, therefore if represented a larger portion of
the watershed, or higher elevation areas which accumulate a deeper snowpack, could have more of an
impact on estimates of basin wide snow storage.
For the 12th of June survey 2023, in which all plots surveyed were snow-free, both bare earth models
measured negative snow depths. Though clearly erroneous, this result is consistent with the mean
difference observed at the plot scale within that survey, in which for all cover types the bias was
negative, ranging from on average a slight negative on the road (-4 +/- 3 cm) and alpine (-3 +/- 1 cm),
111

to more negative in the dense old growth (-12 +/- 3 cm), old growth (-17 +/- 4 cm), juvenile (-18 +/- 12
cm) and regenerating stands (-31 +/-7 cm). As Russell Creek appeared to only hold snow in the higher
alpine portions of the watershed outside of the sampling extent, the accumulation of the negative snow
depths would result in a negative snow volume overall. For the old growth, the negative bias was
consistent with the biases observed with the snow on, versus for the regenerating and juvenile, this was
the opposite, as all other surveys –including other snow free surveys – the bias was positive, with
LiDAR over-predicting snow. As vegetation was greened up in the regenerating and juvenile stands by
June 12th survey, the poor ground point coverage of the snow-on surface combined with the poor ground
point coverage of the bare earth model increased the uncertainty overall.
The watershed scale component of the analysis utilized the interpolated rasters, in order to be a viable
operational workflow with equal pixel count between surveys. At the watershed scale, the proportion of
missing data was substantial (42 – 68%), particularly in the lower elevation portions of the watershed
and regenerating stands (Figure 7.3). As the snowpack thinned and melted, the number of no data pixels
increased across all cover types. The most significant data gaps occur in the mature regenerating stands,
where the range of missing data was consistently high (76 to 83%) across Russell Creek regardless of
whether snow on or snow free conditions. As observed within the alpine, small, isolated areas of
extreme snow depth errors (>15 m snow depth at 600 m elevation) coincide with areas in which there
are no ground returns. Compounded with the proportion of data missing in the bare earth model (76%),
and absence of manual observations, additional analysis outside the scope of this study is required to
assess overall confidence in the mature regenerating forest. For example, the game camera sites
ST7GC1 and ST7GC2 could be used to determine approximate snow depth within the regenerating
forest for each LiDAR survey.

112

Figure 7.3 Imagery of the ground point density (GPD) within the Russell Creek watershed showing the March 10th acquisition (A1, A2),
where snowpack is covering much of the watershed and high snow accumulation survey, and low snow accumulation May 25 th 2023 (B1,
B2). As snow melts in the watershed, the number of no data pixel increases, and the ground point density decreases.

7.4 Study Limitations and Recommendations
This study highlights the importance of in-situ measurements which are spatially distributed over the
different cover types represented in a watershed to validate LiDAR snow depth. Through the sampling
strategy employed in this study, biases inherent to specific stand types at the plot scale were detected. If
the manual measurement campaign were more constrained - for example only the weather station snow
courses - error would been misconstrued as a much narrower range than when inclusive of the cardinal
plots in complex terrain (Figure 7.4).

113

Figure 7.4 Distribution of the cardinal (blue) and weather station snow course (orange) mean difference for 2022 (left) and 2023 (right)
processed with the 2023 ACO bare earth model.

Although the manual snow depth measurements encompassed a range of cover types, the location was
limited to a fairly small (~4 km2) portion of the Russell Creek watershed (33 km2). Additionally, plots of
the same cover type generally were close together. The watershed cover types were fairly broad, as well
not all cover types classified were able to be sampled (new harvest and the mature regenerating, see sec
3.4). In addition to the spatial distribution of the plots, the timing of in-situ measurements and LiDAR
acquisitions represents another limitation. In general, when in-situ measurements were timed shortly
after LiDAR acquisition, systematic error was reduced. The meteorological correction (based on the
change to snow depth observed at either the game camera snow stake plots or weather station snow
depth sensors) did seem effective at correcting error in some cases (e.g. dense old growth: April 29th
2023 survey), but not in others (e.g. dense old growth: June 24th 2022 survey; Figure 5.10). This may be
due to a thin and intermittent snowpack which results in over-compensation in the bias correction.
An accurate snow free model is crucial to measure snow in complex forested terrain. LiDAR bare earth
collection timed as soon after snow melt as possible to capture leaf-off and potentially still compressed
vegetation is one approach to reduce error. However, this may be impractical for a plane based LiDAR
system, due to the variability of snow free conditions in a watershed, and the high cost of repeat
acquisitions. For example, in 2023 there was a 46 day difference between the first (ST7GC1 at 540 m)
114

and last (ST1GC3 at 1210 m) snow-free game camera plots within Russell Creek. UAV LiDAR has
potential for improving the ground surface model, due to the flexibility of repeat acquisitions, ability to
fly low to the ground, and higher ground point density. Thus, a mosaic approach to a bare earth model,
in which particular elevations or forest stands are flown as per their optimal condition, ie. before leaf-on,
may improve LiDAR error.
The processing methods for the raw LiDAR data were outside the scope of this study. Suggested steps to
consider within the processing of LiDAR data which may improve the bare earth model include:
i)

A comparison of classification algorithms used to distinguish ground and vegetation, which
may “save” more points in the harvested forests

ii)

Investigation of full waveform data, particularly in areas of complex vegetation

iii)

Comparison of interpolation algorithms by different cover types

This study presented several key findings and novel results for measuring snow in predominantly forestcovered and heavily harvested watersheds. The bare earth model proved to be a crucial component in the
ability to measure snow depth accurately, in particular for forests with complex canopy and vegetation
structures such as the harvested stands. In the case of Russell Creek, forest harvesting has occurred
below 1300 meters, which coincides with the approximate elevation threshold where snow accumulates
to a lesser degree. Therefore, the overall impact of the LiDAR accuracy in harvested terrain is dictated
by the spatial distribution within the watershed.

115

Chapter 8: Summary
This study evaluated airborne LiDAR as a method to measure snow in a densely forested watershed,
with great variability and complexity to the stands as a result of over 40 years of forest harvesting.
Though the methodology to evaluate accuracy with LiDAR with manual snow depth measurements is
well established, there were several components to this study which were novel: 1) the complexity of
forest stands in which manual observations were collected, 2) the comparison of LiDAR accuracy in two
bare earth models, and 3) LiDAR snow depth observations for pre- and post- forest harvest. Through
both the manual plots and forest harvest experiment components of this study, it is shown that airborne
LiDAR can quantify snow storage, even in dense forest, provided a high accuracy ground surface model
is used. Important considerations for the bare earth model which may impact snow depth accuracy are:
the characteristics of the trees, the height and complexity of the ground vegetation, the seasonality, and
changes to forest cover. The need for timely and accurate snow depth data at the watershed scale has
never been greater, as forest harvesting occurs at increasingly higher elevations as well as the
compounded effects of climate in coastal watersheds is poorly understood (Bathurst et al. 2020; Gillett
et al. 2022; Johnson and Alila 2023). Results from this study support the validity of LiDAR as a tool to
increase our understanding of hydrological recovery in coastal watersheds, however there remains high
uncertainty in some regenerating stands making it difficult to quantify differences.

116

Data Availability Statement
Data and code used in this thesis are not currently available in public repositories. Please contact Dr. Bill
Floyd (Bill.Floyd@viu.ca) for data access.

117

References
Alila Y, Kuraś PK, Schnorbus M, Hudson R. 2009. Forests and floods: A new paradigm sheds light on
age-old controversies. Water Resour Res. 45(8). doi:10.1029/2008wr007207.
http://doi.wiley.com/10.1029/2008WR007207.
Anderson K, Hancock S, Disney M, Gaston KJ. 2016. Is waveform worth it? A comparison of Li DAR
approaches for vegetation and landscape characterization. Remote Sens Ecol Conserv. 2(1):5–15.
Axelsson P. 2000. DEM generation from laser scanner data using adaptive tin models. Int Arch
Photogramm Remote Sens. XXXIII(Part B4).
Bathurst JC, Fahey B, Iroumé A, Jones J. 2020. Forests and floods: Using field evidence to reconcile
analysis methods. Hydrol Process. 34(15):3295–3310.
Berris SN, Harr RD. 1987. Comparative snow accumulation and melt during rainfall in forested and
clear‐cut plots in the western Cascades of Oregon. Water Resour Res.
Boon S. 2009. Snow ablation energy balance in a dead forest stand. Hydrol Process. 23(18):2600–2610.
Bosch JM, Hewlett JD. 1982. A review of catchment experiments to determine the effect of vegetation
changes on water yield and evapotranspiration. J Hydrol. 55(1):3–23.
Broxton PD, J. D. van Leeuwen W, Biederman J. 2019. Improving Snow Water Equivalent Maps With
Machine. Water Resources Res. 55(5):3739–3757. https://doi.org/10.1029/2018WR024146
Currier WR, Lundquist JD. 2018. Snow depth variability at the Forest edge in multiple climates in the
western United States. Water Resour Res. 54(11):8756–8773.
Currier WR, Pflug J, Mazzotti G, Jonas T, Deems JS, Bormann KJ, Painter TH, Hiemstra CA, Gelvin A,
Uhlmann Z, et al. 2019. Comparing aerial lidar observations with terrestrial lidar and snow‐probe
transects from NASA’s 2017 SnowEx campaign. Water Resour Res. 55(7):6285–6294.
Deems JS, Fassnacht SR, Elder KJ. 2006. Fractal Distribution of Snow Depth from Lidar Data. J
Hydrometeorol. 7(2):285–297.
Deems JS, Painter TH, Finnegan DC. 2013. Lidar measurement of snow depth: a review. J Glaciol.
59(215):467–479.
Deschamps-Berger C, Gascoin S, Shean D, Besso H, Guiot A, López-Moreno JI. 2023. Evaluation of
snow depth retrievals from ICESat-2 using airborne laser-scanning data. The Cryosphere. 17(7):2779–
2792.
Dickerson-Lange SE, Gersonde RF, Hubbart JA, Link TE, Nolin AW, Perry GH, Roth TR, Wayand NE,
Lundquist JD. 2017. Snow disappearance timing is dominated by forest effects on snow accumulation in
warm winter climates of the Pacific Northwest, United States. Hydrol Process. 31(10):1846–1862.

118

Dickerson-Lange SE, Lutz JA, Gersonde R, Martin KA, Forsyth JE, Lundquist JD. 2015. Observations
of distributed snow depth and snow duration within diverse forest structures in a maritime mountain
watershed. Water Resour Res. 51(11):9353–9366.
Dickerson-Lange SE, Lutz JA, Martin KA, Raleigh MS, Gersonde R, Lundquist JD. 2015. Evaluating
observational methods to quantify snow duration under diverse forest canopies. Water Resour Res.
51(2):1203–1224.
Dickinson T. 2022. RPAS-SfM snow depth and snow density mapping in disturbed vegetated
mountainous environments of Coastal British Columbia.
Dingman LS. 2015. Physical Hydrology: Third Edition. Waveland Press.
Disney M. 2019. Terrestrial LiDAR: a three-dimensional revolution in how we look at trees. New
Phytol. 222(4):1736–1741.
Ellis CR, Pomeroy JW, Link TE. 2013. Modeling increases in snowmelt yield and desynchronization
resulting from forest gap-thinning treatments in a northern mountain headwater basin. Water Resour
Res. 49(2):936–949.
Evan A, Eisenman I. 2021. A mechanism for regional variations in snowpack melt under rising
temperature. Nat Clim Chang. 11(4):326–330.
Floyd B. 2010. Russell Creek Watershed. Streamline. 14(1).
Floyd B, Bishop A, Menounos B, Heathfield D, Beffort S, Gonzalez Arriola S. 2020. Fixed Wing
LiDAR Snow Mapping in the Upper Englishman and Little Qualicum River Watershed: 2020 Progress
Report. Vancouver Island University Report No.: 25.
Floyd W, Weiler M. 2008. Measuring snow accumulation and ablation dynamics during rain‐on‐snow
events: innovative measurement techniques. Hydrol Process. 22(24):4805–4812.
Floyd WC. 2012. Snowmelt Energy Flux Recovery During Rain-on-Snow in Regenerating Forests.
University of British Columbia.
Foord V, Floyd B, Weick T, Others. 2014. Enhancing provincial climate monitoring: the Forest
Ecosystem Research Network and the climate related monitoring program. Extension Note-British
Columbia Ministry of Forests and Range.(111).
https://www.cabdirect.org/cabdirect/abstract/20143249021.
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 Resour Res. 59(9). doi:10.1029
Gelfan AN, Pomeroy JW, Kuchment LS. 2004. Modeling Forest Cover Influences on Snow
Accumulation, Sublimation, and Melt. J Hydrometeorol. 5(5):785–803.
Gillett NP, Cannon AJ, Malinina E, Schnorbus M, Anslow F, Sun Q, Kirchmeier-Young M, Zwiers F,
Seiler C, Zhang X, et al. 2022. Human influence on the 2021 British Columbia floods. Weather and
Climate Extremes. 36:100441.
119

Gould SB, Glenn NF, Sankey TT, McNamara JP. 2013. Influence of a Dense, Low-height Shrub Species
on the Accuracy of a Lidar-derived DEM. Photogrammetric Engineering & Remote Sensing. 79(5):421–
431.
Goulden T, Hopkinson C, Jamieson R, Sterling S. 2016. Sensitivity of DEM, slope, aspect and
watershed attributes to LiDAR measurement uncertainty. Remote Sens Environ. 179:23–35.
Guo Q, Li W, Yu H, Alvarez O. 2010. Effects of topographic variability and lidar sampling density on
several DEM interpolation methods. Photogramm Eng Remote Sens. 76:701–712.
Hall DK, Kelly REJ, Foster JL, Chang ATC. 2005. 55: Estimation of Snow Extent and Snow Properties.
In: Encyclopedia of Hydrological Sciences.
Harder P, Pomeroy JW, Helgason WD. 2020. Improving sub-canopy snow depth mapping with
unmanned aerial vehicles: lidar versus structure-from-motion techniques. cryosphere. 14(6):1919–1935.
Harder P, Schirmer M, Pomeroy J, Helgason W. 2016. Accuracy of snow depth estimation in mountain
and prairie environments by an unmanned aerial vehicle. The Cryosphere Discussions.:1–22.
Harr RD. 1981. Some characteristics and consequences of snowmelt during rainfall in western Oregon.
Journal of Hydrology. 53(3–4):277–304. doi:10.1016/0022-1694(81)90006-8.
Hartzell PJ, Gadomski PJ, Glennie CL, Finnegan DC, Deems JS. 2015. Rigorous error propagation for
terrestrial laser scanning with application to snow volume uncertainty. J Glaciol. 61(230):1147–1158.
Hedstrom NR, Pomeroy JW. 1998. Measurements and modelling of snow interception in the boreal
forest. Hydrol Process. 12(10–11):1611–1625.
Hidalgo-Hidalgo J-D, Collados-Lara A-J, Pulido-Velazquez D, Rueda FJ, Pardo-Igúzquiza E. 2022.
Analysis of the Potential Impact of Climate Change on Climatic Droughts, Snow Dynamics, and the
Correlation between Them. Water. 14(7):1081.
Hodgson ME, Bresnahan P. 2004a. Accuracy of Airborne LIDAR-Derived Elevation: Empirical
Assessment and Error Budget. Photogrammetric Engineering & Remote Sensing. 70(3):331–339.
Hodgson ME, Bresnahan P. 2004b. Accuracy of airborne lidar-derived elevation. Photogramm Eng
Remote Sensing. 70(3):331–339.
Hojatimalekshah A, Uhlmann Z, Glenn NF, Hiemstra CA, Tennant CJ, Graham JD, Spaete L, Gelvin A,
Marshall H-P, McNamara JP, et al. 2021. Tree canopy and snow depth relationships at fine scales with
terrestrial laser scanning. The Cryosphere. 15(5):2187–2209.
Hopkinson C, Chasmer LE, Sass G, Creed IF, Sitar M, Kalbfleisch W, Treitz P. 2005. Vegetation class
dependent errors in lidar ground elevation and canopy height estimates in a boreal wetland environment.
Can J Remote Sens. 31(2):191–206.
Hopkinson Chris, Collins T, Anderson A, Pomeroy J, Spooner I. 2012. Spatial Snow Depth Assessment
Using LiDAR Transect Samples and Public GIS Data Layers in the Elbow River Watershed, Alberta.
Canadian Water Resources Journal / Revue canadienne des ressources hydriques. 37(2):69–87.
120

Hopkinson C., Pomeroy J, Debeer C, Ellis C, Anderson A. 2012. Relationships between snowpack depth
and primary LiDAR point cloud derivatives in a mountainous environment. 352:354–358.
Hopkinson C, Sitar M, Chasmer L, Gynan C, Agro D, Enter R, Foster J, Heels N, Hoffman C, Nillson J,
et al. 2001. Mapping the spatial distribution of snowpack depth beneath a variable forest canopy using
airborne laser altimetry. In: Proceedings of the 58th Eastern Snow Conference, Ottawa, Ontario,
Canada: USA, Eastern Snow Conference. p. 253–264.
Hopkinson C, Sitar M, Chasmer L, Treitz P. 2004. Mapping Snowpack Depth beneath Forest Canopies
Using Airborne Lidar. Photogrammetric Engineering & Remote Sensing. 70(3):323–330.
Hudson R. 2000. Snowpack recovery in regenerating coastal British Columbia clearcuts. Can J For Res.
30(4):548–556.
Hudson R, Anderson A. 2006. Russell Creek: Summary of Research and Implications for Professional
Practice. BCMOF Report No.: EN-022.
Hudson R, Horel G. 2007. An operational method of assessing hydrologic recovery for Vancouver
Island and south coastal BC. Forest Service British Columbia Report No.: 032.
Jacobs JM, Hunsaker AG, Sullivan FB, Palace M, Burakowski EA, Herrick C, Cho E. 2021. Snow depth
mapping with unpiloted aerial system lidar observations: a case study in Durham, New Hampshire,
United States. cryosphere. 15(3):1485–1500.
Johnson RSH, Alila Y. 2023. Nonstationary stochastic paired watershed approach: Investigating forest
harvesting effects on floods in two large, nested, and snow-dominated watersheds in British Columbia,
Canada. J Hydrol. 625:129970.
Jones HG, Walker DA, Pomeroy JW, Hoham R. 2001. Snow Ecology. Cambridge University Press.
Karbou F, Veyssière G, Coleou C, Dufour A, Gouttevin I, Durand P, Gascoin S, Grizonnet M. 2021.
Monitoring Wet Snow Over an Alpine Region Using Sentinel-1 Observations. Remote Sensing.
13(3):381.
Kinar NJ, Pomeroy JW. 2015. Measurement of the physical properties of the snowpack. Rev Geophys.
53(2):481–544.
King, F., Kelly, R. & Fletcher, C. G. 2022. Evaluation of LiDAR-Derived Snow Depth Estimates From
the iPhone 12 Pro. IEEE Geoscience and Remote Sensing Letters. 19: 1–5
Kucharczyk M, Hugenholtz CH, Zou X. 2018. UAV–LiDAR accuracy in vegetated terrain.
doi:10.1139/juvs-2017-0030.
Landry A, Floyd W, Mc Innes W, Holmes K, Arriola SG, Cebulski A, Dickinson T, Butler S, Heathfield
D, Menounos B. 2021. Measuring snow depth using RPAS photogrammetry in a subalpine coastal
region of British Columbia. In: Western Snow Conference.
Li L, Pomeroy J. 1997. Probability of Blowing snow occurrence by wind. J Geophys Res. 102(D18).

121

Lievens H, Demuzere M, Marshall H-P, Reichle RH, Brucker L, Brangers I, de Rosnay P, Dumont M,
Girotto M, Immerzeel WW, et al. 2019. Snow depth variability in the Northern Hemisphere mountains
observed from space. Nat Commun. 10(1):4629.
Lin Y-C, Shao J, Shin S-Y, Saka Z, Joseph M, Manish R, Fei S, Habib A. 2022. Comparative Analysis
of Multi-Platform, Multi-Resolution, Multi-Temporal LiDAR Data for Forest Inventory. Remote Sens.
14(3):649.
Lindberg E, Olofsson K, Holmgren J, Olsson H. 2012. Estimation of 3D vegetation structure from
waveform and discrete return airborne laser scanning data. Remote Sens Environ. 118:151–161.
López-Moreno JI, Pomeroy JW, Morán-Tejeda E, Revuelto J, Alonso-González E. 2021. Changes in the
frequency of global high mountain rain-on-snow events due to climate warming. Environ Res Lett.
16(9). doi:10.1088/1748-9326/ac0dde.
Lundquist JD, Dickerson-Lange SE, Lutz JA, Cristea NC. 2013. Lower forest density enhances snow
retention in regions with warmer winters: A global framework developed from plot-scale observations
and modeling. Water Resour Res. 49(10):6356–6370.
Marks D, Dozier J. 1992. Climate and energy exchange at the snow surface in the Alpine Region of the
Sierra Nevada: 2. Snow cover energy balance. Water Resour Res. 28(11):3043–3054.
Marks D, Kimball J, Tingey D, Link T. 1998. The sensitivity of snowmelt processes to climate
conditions and forest cover during rain-on-snow: a case study of the 1996 Pacific Northwest flood.
Hydrol Process. doi:10.1002/(SICI)1099-1085(199808/09)12:10/11.
Marti R, Gascoin S, Berthier E, de Pinel M, Houet T, Laffly D. 2016. Mapping snow depth in open
alpine terrain from stereo satellite imagery. Cryosphere. 10(4):1361–1380.
Maté-González MÁ, Di Pietra V, Piras M. 2022. Evaluation of Different LiDAR Technologies for the
Documentation of Forgotten Cultural Heritage under Forest Environments. Sensors . 22(16).
doi:10.3390/s22166314.
Mazzotti G, Currier WR, Deems JS, Pflug JM, Lundquist JD, Jonas T. 2019. Revisiting snow cover
variability and canopy structure within forest stands: Insights from airborne lidar data. Water Resour
Res. 55(7):6198–6216.
Meidinger DV, Pojar J. 1991. Ecosystems of British Columbia. B.C Ministry of Forests.
Menounos B, Pelto BM, Fleming SW, Dan Moore R, Weber F, Hutchinson D, Brahney J. 2020. Glaciers
in the Canadian Columbia Basin, Technical Report. Canadian Columbia Basin Glacier and Snow
Research Network Watershed Sciences Faculty Publications.
Moore RD, Scott DF. 2005. Camp Creek Revisited: Streamflow Changes Following Salvage Harvesting
in a Medium-Sized, Snowmelt-Dominated Catchment. Canadian Water Resources Journal / Revue
canadienne des ressources hydriques. 30(4):331–344.

122

Murgaš V, Sačkov I, Sedliak M, Tunák D, Chudý F. 2018. Assessing horizontal accuracy of inventory
plots in forests with different mix of tree species composition and development stage. Journal of Forest
Science. doi:10.17221/92/2018-JFS.
Murray CD, Buttle JM. 2003. Impacts of clearcut harvesting on snow accumulation and melt in a
northern hardwood forest. J Hydrol. 271(1):197–212.
Neary D. 2016. Long-term forest paired catchment studies: What do they tell us that landscape-level
monitoring does not? For Trees Livelihoods. 7(12):164.
Nolin AW, Daly C. 2006. Mapping “At Risk” Snow in the Pacific Northwest. J Hydrometeorol.
7(5):1164–1171.
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.
Painter TH, Berisford DF, Boardman JW, Bormann KJ, Deems JS, Gehrke F, Hedrick A, Joyce M,
Laidlaw R, Marks D, et al. 2016. The Airborne Snow Observatory: Fusion of scanning lidar, imaging
spectrometer, and physically-based modeling for mapping snow water equivalent and snow albedo.
Remote Sens Environ. 184:139–152.
Painter TH, Bormann KJ, Deems JS, Berisford DF. 2017. The airborne snow observatory during NASA
snow experiment (SnowEx) year 1: Mapping of snow water equivalent and snow albedo and
constraining understanding of the physical environment. In: 2017 IEEE International Geoscience and
Remote Sensing Symposium (IGARSS). IEEE. p. 1399–1402.
Pelto BM, Menounos B, Marshall SJ. 2019. Multi-year evaluation of airborne geodetic surveys to
estimate seasonal mass balance, Columbia and Rocky Mountains, Canada. The Cryosphere. 13(6):1709–
1727.
Pham HC, Alila Y. 2024. Science of forests and floods: The quantum leap forward needed, literally and
metaphorically. Sci Total Environ. 912:169646.
Pike R, Redding TE, Wilford D, Toews DAA. 2010. Chapter 16 — Detecting and Predicting Changes in
Watersheds. In: Compendium of Forest Hydrology and Geomorphology in British Columbia. unknown.
Pomeroy J, Fang X, Ellis C. 2012. Sensitivity of snowmelt hydrology in Marmot Creek, Alberta, to
forest cover disturbance. Hydrol Process. 26(12):1891–1904.
Pomeroy JW, Gray DM. 1995. Snowcover Accumulation, Relocation, and Management. Saskatoon,
Saskatchewan : National Hydrology Research Institute (NHRI science report ).
Pomeroy JW, Gray DM, Brown T, Hedstrom NR, Quinton WL, Granger RJ, Carey SK. 2007. The cold
regions hydrological model: a platform for basing process representation and model structure on
physical evidence. Hydrol Process. 21(19):2650–2667.
Pomeroy JW, Toth B, Granger RJ, Hedstrom NR, Essery RLH. 2003. Variation in Surface Energetics
during Snowmelt in a Subarctic Mountain Catchment. J Hydrometeorol. 4(4):702–719.

123

Prokop A, Schirmer M, Rub M, Lehning M, Stocker M. 2008. A comparison of measurement methods:
terrestrial laser scanning, tachymetry and snow probing for the determination of the spatial snow-depth
distribution on slopes. Ann Glaciol. 49:210–216.
Proulx H, Jacobs JM, Burakowski EA, Cho E, G. A, Sullivan FB, Palace M, Wagner C. 2022.
Comparison of in-situ snow depth measurements and impacts 2 on validation of unpiloted aerial system
lidar over a mixed-use 3 temperate forest landscape. doi:10.5194/tc-2022-7.
Reutebuch SE, McGaughey RJ, Andersen H-E, Carson WW. 2003. Accuracy of a high-resolution lidar
terrain model under a conifer forest canopy. Can J Remote Sens. 29(5):527–535.
Revuelto J, López-Moreno J-I, Azorin-Molina C, Alonso-González E, Sanmiguel-Vallelado A. 2016.
Small-Scale Effect of Pine Stand Pruning on Snowpack Distribution in the Pyrenees Observed with a
Terrestrial Laser Scanner. For Trees Livelihoods. 7(8):166.
Revuelto J, López-Moreno JI, Azorin-Molina C, Zabalza J, Arguedas G, Vicente-Serrano SM. 2014.
Mapping the annual evolution of snow depth in a small catchment in the Pyrenees using the long-range
terrestrial laser scanning. J Maps. 10(3):379–393.
Richardson PW, Cafferata PH, Dymond SF, Keppeler ET, Wagenbrenner JW, Whiting JA. 2023. Past
and future roles of paired watersheds: a North American inventory and anecdotes from the Caspar Creek
Experimental Watersheds. Frontiers in Forests and Global Change. 6. doi:10.3389/ffgc.2023.1275392.
Roth TR, Nolin AW. 2016 Nov 1. Forest impacts on snow accumulation and ablation across an
elevation gradient in a temperate montane environment. Hydrol Earth Syst Sci Discuss.:1–24.
Russell M, U. H. Eitel J, J. Maguire A, E. Link T. 2020. Toward a Novel Laser-Based Approach for
Estimating Snow Interception. Remote Sensing. 12(7):1146.
Schmidt RA, Pomeroy J. 1990. Bending of a conifer branch at subfreezing temperatures: implications
for snow interception. Can J For Res. 20: 1250-1253.
Schnorbus M, Werner A, Bennett K. 2014. Impacts of climate change in three hydrologic regimes in
British Columbia, Canada. Hydrol Process. 28(3):1170–1189.
Seyednasrollah B, Kumar M. 2019. How Surface Radiation on Forested Snowpack Changes across a
Latitudinal Gradient. Hydrology. 6(3):62.
Storck P, Lettenmaier DP, Bolton SM. 2002. Measurement of snow interception and canopy effects on
snow accumulation and melt in a mountainous maritime climate, Oregon, United States. Water Resour.
Res. 38(11):5–1. doi:10.1029/2002wr001281.
Swanson R, Stevenson DR. 1971. Managing snow accumulation and melt under leafless aspen to
enhance watershed value. :63.
Swanson RH, Golding! DL, Rothwe RL, Bernier PY. 1986. Hydrologic effects of clear-cutting at
marmot creek and Streeter watersheds, Alberta. Can. For. Serv., Edmonton, Alberta. Inf. Rep. NOR-X278.

124

Tedesco M, Derksen C, Deems JS, Foster JL. 2014. Remote sensing of snow depth and snow water
equivalent. In: Remote Sensing of the Cryosphere. Chichester, UK: John Wiley & Sons, Ltd. p. 73–98.
Tinkham WT, Huang H, Smith AMS, Shrestha R, Falkowski MJ, Hudak AT, Link TE, Glenn NF,
Marks DG. 2011. A Comparison of Two Open Source LiDAR Surface Classification Algorithms.
Remote Sensing. 3(3):638–649.
Tinkham WT, Smith AMS, Hoffman C, Hudak AT, Gessler PE. 2012. Investigating the influence of
LiDAR ground surface errors on the utility of derived forest inventories. Can J For Res. 42(3):413–422.
Tinkham WT, Smith AMS, Marshall H-P, Link TE, Falkowski MJ, Winstral AH. 2014. Quantifying
spatial distribution of snow depth errors from LiDAR using Random Forest. Remote Sens Environ.
141:105–115.
Troendle CA, King RM. 1985. The Effect of Timber Harvest on the Fool Creek Watershed, 30 years
later. Water Resour. Res. 21: 1915-1022. doi:10.1029/WR021i012p01915
U.S. Department of Agriculture. Snow Telemetry (SNOTEL) Network. Automated Snow Monitoring.
Monitoring Programs. [accessed 2024 May 24].
https://www.nrcs.usda.gov/wps/portal/wcc/home/aboutUs/monitoringPrograms/automatedSnowMonitor
ing/.
Varhola A, Coops NC, Bater CW, Teti P, Boon S, Weiler M. 2010. The influence of ground- and lidarderived forest structure metrics on snow accumulation and ablation in disturbed forests. Can J For Res..
40(4):812–821. doi:10.1139/x10-008.
Varhola A, Coops NC, Weiler M, Moore RD. 2010. Forest canopy effects on snow accumulation and
ablation: An integrative review of empirical results. J Hydrol. 392(3):219–233.
Vionnet V, Marsh CB, Menounos B, Gascoin S, Wayand NE, Shea J, Mukherjee K, Pomeroy JW. 2020.
Multi-scale snowdrift-permitting modelling of mountain snowpack. The Cryosphere. doi:10.5194/tc2020-187.
Winkler R, Boon S. 2017. Equivalent clearcut area as an indicator of hydrologic change in snowdominated watersheds of Southern British Columbia. Prov BC, Victoria, BC Exten Note 118.
https://www.for.gov.bc.ca/hfd/pubs/Docs/En/En118.htm.
Winkler R, Spittlehouse D, Boon S, Zimonick B. 2015. Forest disturbance effects on snow and water
yield in interior British Columbia. Hydrol Res. 46(4):521–532.
Winkler RD, Allen DM, Giles TR, Heise BA, Moore RD, Redding TE, Spittlehouse DL, Wei X. 2021.
Approaching four decades of forest watershed research at Upper Penticton Creek, British Columbia: A
synthesis. Hydrol Process. 35:e14123.
Winkler RD, Spittlehouse DL, Golding DL. 2005. Measured differences in snow accumulation and melt
among clearcut, juvenile, and mature forests in southern British Columbia. Hydrol Process. 19(1):51–62.
World Meterological Organization. 2008. Guide to Meterological Instruments and Methods of
Observation. World Meterological Organization.
125

Zheng Z, Kirchner PB, Bales RC. 2016. Topographic and vegetation effects on snow accumulation in
the southern Sierra Nevada: a statistical summary from lidar data. The Cryosphere. 10:257–269.

126

Appendices
Appendix A: ACO 2023 Bare Earth Model
In order to have the full AOI an ACO acquisition (previously BC Timber Sales), a new bare earth was
initially flown on July 14th 2023. Due to high errors concentrated within the harvested areas of the
watershed, a reacquisition occurred on August 19th 2023. Although overall an improved acquisition, the
August flight has a critical data gap (~860 m2) in the old growth forest, partially within the Russell
Creek watershed AOI. Though it makes up a small portion of the watershed, six of the cardinal plots ST1OG7, ST1OG6, ST1OG5, ST1OG4, ST1OG3 and half of ST1A2 – are situated within this data gap.

Figure A 1 Hillshade of the August LiDAR acquisition, showing the data gap in red, and the Russell creek watershed in blue

127

An analysis of similar and adjacent old growth forest stands showed good agreement in the ground
surface between the July and August acquisitions. A histogram of the difference in pixel values between
the two bare earth models showed no bias within the old growth forest adjacent to the data (Figure A- 1).
Based on this result, I filled the data gap has with the pixels from the July 14th acquisition bare earth
model.

Figure A- 1 A histogram (A) of the difference in pixel values between July and August bare earth models from the Old Growth AOI
boundary (B), which is similar forest stand to the data gap in the August acquisition

128

Appendix B: LiDAR Processing
Table B 1 LAS tools processing commands used to take the raw LAS files to a gridded raster product

LAStools

Command

1 - lastile

lastile -cpu64 -i "(input).laz" -o "tile.laz" -tile_size 225 -buffer 25 -flag_as_withheld faf -cores 12 -odir “(output)\1 Tiles" –olaz

2 - lassort

-i tiles_raw\*.las ^ -odir tiles_sorted -olaz ^-cores 7

3 - lasnoise

lasnoise -cpu64 -i " tile.laz" -step_xy 1 -step_z 0.1 -isolated 100 -remove_noise -odir
"3 Noise" -o "noise5.laz"

4 – lasground

lasground -i LAS\tiles_isolated\*.laz ^ -city -ultra_fine -ignore_class 0 7 ^

(temporary)

-odir LAS\tiles_temp_ground -olaz ^ -cores 7

5 – lasheight

lasheight -cpu64 –lof *.txt -classify_below -2 7 -classify_between 0.3 20 3 -olaz

6 – las2las

-cpu64 -i "noise5.laz" -keep_classification 2 -odir "las2las de-noise" -olaz

7 - lasthin

lasthin -cpu64 -i "*.laz" -step 0.25 -adaptive 0.1 -odir " \7 Thin" -olaz

8 – lasground

lasground_new -cpu64 -i " \noise5.laz" -city -hyper_fine

(final)

-odir " \8 Final Ground" -olaz

129

Year

9-Mar

11-Mar

11-Mar

11-Mar

11-Mar

11-Mar

9-Mar

9-Mar

ST1OG7

ST1OG6

ST1OG5

ST1OG4

ST1OG3

ST1A1

ST1A2

9-Mar

ST2CC1

11-Mar

9-Mar

ST4RF2

ST1OG8

9-Mar

ST2CC2

23-May

ST4RF1

6-Mar

ST1OG5

ST1A2

6-Mar

ST1OG6

23-May

6-Mar

ST1OG7

ST1A1

6-Mar

ST1OG8

6-Mar

6-Apr

ST2CC2

23-May

6-Apr

ST2CC1

ST1OG3

9-Mar

ST1OG4

9-Mar

ST4RF2

Date

ST4RF1

Plot

688197.309

688249.46

687934.342

687840.936

687801.632

687765.507

687698.028

687709.303

687829.778

687952.609

688271.451

688443.783

688200.887

688249.878

687932.934

687843.532

687795.912

687764.006

687700.971

687704.431

5577335.367

5577217.991

688444.08

Easting
(m)
688270.471

5575975.985

5575967.128

5576153.66

5576193.032

5576230.787

5576273.717

5576397.792

5576464.042

5577335.367

5577217.991

5577506.572

5577397.82

5575971.527

5575967.973

5576154.788

5576192.942

5576230.94

5576273.928

5576400.245

5576461.385

687829.778

687952.609

5577398.062

Northing
(m)
5577507.661

1509.788

1506.961

1382.505

1336.213

1311.976

1301.808

1258.292

1232.814

1009.134

1077.448

832.733

830.827

1509.596

1507.3

1380.392141

1351.823

1313.128

1301.229

1265.579

1230.75

1009.134

1077.448

831.905

Ellipsoid Height
(m)
832.584

0.004

0.004

0.234

0.234

0.234

0.234

0.024

0.009

0.004

0.004

0.018

0.019

0.007

0.004

2.222

n/a

0.046

0.046

0.046

0.041

0.006

0.004

1.982

HRMS
(m)
0.012

FIXED

FIXED

FIXED

FLOAT

FLOAT

FLOAT

FLOAT

FLOAT

FIXED

FIXED

FIXED

FIXED

FIXED

FIXED

AUTONOMOUS

FLOAT

FLOAT

FLOAT

FLOAT

FLOAT

FIXED

FIXED

AUTONOMOUS

FIXED

STATUS

130

0.007

0.007

0.464

0.464

0.464

0.464

0.035

0.011

0.005

0.006

0.043

0.043

0.009

0.006

3.034

n/a

0.114

0.114

0.114

0.111

0.008

0.005

2.935

VRMS
(m)
0.016

n/a

12

12

15

8

9

8

9

10

13

15

12

15

13

13

12

9

9

5

6

13

13

10

11

SATS

1.7

1.7

3.3

3.4

3.4

3.5

2.9

2.6

1.6

1.8

2.1

2

1.752

1.766

1.88

n/a

3.119

3.035

3.744

2.847

1.377

2.045

1.771

PDOP
(m)
2.324

0.9

0.9

0.9

1.2

1.2

1.5

1.4

1.8

0.9

0.9

1

0.8

0.992

0.982

1.111

n/a

1.4

1.2

2.5

1.386

0.81

1.084

0.991

HDOP
(m)
1.261

1.5

1.5

3.2

3.2

3.2

3.2

2.5

1.9

1.3

1.5

1.8

1.8

1.443

1.468

1.517

n/a

2.788

2.788

2.788

2.487

1.113

1.734

1.468

VDOP
(m)
1.952

1.1

1.1

2.9

2.8

2.8

2.7

2

1

0.9

1.1

1.1

1.2

1.084

1.103

1.19

n/a

1.941

2.07

0.724

1.823

0.74

1.284

1.026

TDOP
(m)
1.483

2

2.1

4.4

4.4

4.4

4.4

3.5

2.8

1.8

2.1

2.3

2.3

2.06

2.083

2.225

n/a

3.674

3.674

3.674

3.381

1.563

2.414

2.047

GDOP
(m)
2.757

C- 1 GNSS data results for plot centres for 2022 and 2023,with metrics for horizontal root mean square (HRMS), vertical root mean square (VRMS), number of satellietes (SATS),DOP
stands for Dilution of Precision, and PDOP is Position of DOP, HDOP is Horizontal of DOP, VDOP is vertical of DOP, TDP is time DOP, and GDOP is geometric DOP,all measured im
meters

Appendix C: GNSS report

2022

2023

ST6GC3

ST6GC2

ST6GC1

ST2GC1

ST1GC3

ST1GC2

ST1GC1

Site

Date

10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-14-2020
10-13-2020
10-13-2020
10-13-2020
10-13-2020
10-13-2020
10-13-2020
10-13-2020
10-13-2020
10-13-2020
10-13-2020
10-13-2020
10-13-2020
10-13-2020
10-13-2020
10-13-2020
10-13-2020

stake

1
2
3
4
5
1
2
3
4
5
1
2
3
4
5
1
2
3
4
5
center
1
2
3
4
5
center
1
2
3
4
5
center
1
2
3

687756.159
687756.647
687758.94
687757.723
687757.085
687675.179
687673.941
687674.384
687676.216
687677.142
687712.453
687711.199
687714.662
687715.551
687715.452
687882.409
687881.895
687883.255
687883.362
687886.047
688602.465
688604.53
688604.858
688602.853
688600.971
688602.096
688528.976
688527.825
688528.784
688530.239
688531.413
688531.369
688468.36
688466.55
688468.953
688470.688

Easting
(m)
5576264.934
5576266.466
5576272.695
5576267.703
5576268.909
5576346.562
5576346.804
5576345.63
5576344.805
5576344.076
5576503.843
5576503.309
5576503.29
5576502.808
5576501.767
5577010.497
5577006.75
5577005.755
5577012.553
5577009.871
5577488.953
5577490.94
5577494.589
5577490.035
5577491.99
5577488.531
5577927.947
5577925.15
5577926.406
5577924.145
5577919.442
5577917.953
5577915.376
5577914.554
5577916.45
5577915.316

Northing
(m)

Ellipsoid
Height
(m)
1280.234
1282.879
1298.229
1295.317
1294.923
1268.239
1268.175
1267.554
1270.011
1267.936
1210.491
1207.21
1209.902
1204.947
1221.352
1090.714
1097.898
1096.33
1095.81
1090.854
809.301
809.456
811.383
805.754
814.149
807.072
786.941
784.622
781.162
787.053
787.526
785.828
771.457
773.199
771.341
774.779
4.850
1.608
3.958
2.396
1.374
0.198
0.198
0.198
0.198
0.198
0.099
0.403
0.063
0.335
0.148
1.170
1.185
1.318
0.010
1.319
0.478
1.836
1.901
0.107
0.328
0.144
2.439
0.226
0.885
0.267
0.345
0.445
1.419
1.675
5.193
3.404

VRMS
6.595
2.192
4.971
3.633
1.998
0.381
0.381
0.381
0.381
0.381
0.699
1.132
0.290
0.815
0.556
1.674
1.816
1.932
0.020
2.168
0.649
0.963
1.230
0.206
0.639
0.154
3.466
0.140
0.259
0.243
0.413
0.612
0.903
3.185
8.603
5.428

131

HRMS
DGPS
DGPS
DGPS
DGPS
DGPS
FLOAT
FLOAT
FLOAT
FLOAT
DGPS
FLOAT
FLOAT
FLOAT
FLOAT
FLOAT
DGPS
DGPS
DGPS
FIXED
DGPS
FLOAT
DGPS
DGPS
FLOAT
FLOAT
FLOAT
FLOAT
FLOAT
FLOAT
FLOAT
FLOAT
FLOAT
FLOAT
FLOAT
FLOAT
FLOAT

STATUS
7
9
8
8
10
7
12
12
9
8
7
12
5
11
9
7
11
6
6
7
7
9
9
10
6
8
9
5
11
7
10
9
7
10
10
9

SATS
5.102
3.151
3.255
2.622
2.085
3.640
3.178
3.102
3.138
3.581
3.476
3.389
2.448
2.369
2.398
2.312
2.136
2.167
3.480
2.590
2.751
2.396
3.022
1.973
1.985
3.612
1.843
2.217
1.978
1.936
1.834
2.248
1.940
2.351
2.148
2.596

PDOP
3.041
2.044
2.058
1.466
1.363
2.200
1.300
1.100
1.200
2.100
1.440
1.112
0.931
0.904
0.955
1.397
1.184
1.233
1.787
1.484
1.632
1.457
1.827
1.228
1.263
2.114
1.061
1.401
1.206
1.178
1.148
1.303
0.988
1.094
1.110
1.379

HDOP
4.097
2.398
2.522
2.174
1.578
2.900
2.900
2.900
2.900
2.900
3.163
3.201
2.264
2.190
2.200
1.842
1.778
1.783
2.986
2.123
2.214
1.902
2.408
1.544
1.531
2.929
1.507
1.719
1.567
1.537
1.430
1.831
1.670
2.081
1.839
2.199

VDOP
3.724
2.093
2.153
1.809
1.259
2.132
2.774
2.859
2.819
2.231
2.393
2.365
1.573
1.528
1.577
1.447
1.341
1.371
2.294
1.553
1.653
1.532
1.983
1.223
1.226
2.831
1.119
1.432
1.195
1.181
1.088
1.490
1.248
1.582
1.372
1.723

TDOP

6.316
3.783
3.902
3.186
2.435
4.219
4.219
4.219
4.219
4.219
4.220
4.133
2.910
2.820
2.870
2.728
2.522
2.565
4.168
3.020
3.209
2.844
3.615
2.321
2.332
4.590
2.157
2.639
2.311
2.268
2.132
2.697
2.307
2.834
2.549
3.116

GDOP

C-2 1 C- 2 GNSS data results for plot centres for 2022 and 2023,with metrics for horizontal root mean square (HRMS), vertical root mean square (VRMS), number of satellites (SATS),DOP
stands for Dilution of Precision, and PDOP is Position of DOP, HDOP is Horizontal of DOP, VDOP is vertical of DOP, TDP is time DOP, and GDOP is geometric DOP,all in meter

ST7GC2

ST7GC1

4
5
1
2
3
4
5
1
2
3
4
5

10-13-2020
10-13-2020
10-15-2020
10-15-2020
10-15-2020
10-15-2020
10-15-2020
10-15-2020
10-15-2020
10-15-2020
10-15-2020
10-15-2020

688467.947
688471.034
687504.747
687503.797
687503.883
687504.467
687503.659
687489.534
687491.202
687488.521
687486.949
687488.39

5577912.534
5577911.828
5578850.857
5578851.073
5578849.846
5578850.379
5578849.166
5578932.5
5578931.256
5578930.363
5578931.798
5578931.243

775.4
772.787
618.892
607.447
599.308
602.312
600.519
586.577
577.589
584.985
590.954
582.918

2.034
4.615
6.332
6.332
6.332
6.332
6.332
6.141
6.196
6.235
6.278
6.328

132

3.339
8.683
13.398
13.398
13.398
13.398
13.398
12.679
12.895
13.041
13.195
13.363

FLOAT
FLOAT
FLOAT
FLOAT
FLOAT
FLOAT
FLOAT
DGPS
DGPS
DGPS
DGPS
DGPS

8
9
11
5
12
11
10
12
11
11
11
10

2.632
2.495
2.629
2.669
2.593
2.629
2.629
2.098
2.696
2.250
2.507
2.667

1.369
1.171
1.000
1.100
0.900
1.000
1.000
1.257
1.277
1.247
1.230
1.257

2.248
2.204
2.432
2.432
2.432
2.432
2.432
1.680
2.375
1.873
2.185
2.353

1.824
1.700
2.020
1.967
2.066
2.020
2.020
1.322
1.838
1.341
1.622
1.807

3.202
3.020
3.315
3.315
3.315
3.315
3.315
2.480
3.263
2.620
2.986
3.222

Appendix D: Snow depth variability results tables
Table D 1 snow depth variability in cm from the weather station snow courses for 2022

STD. DEV.
(cm)

11-Apr

21-May

24-Jun

9-Mar

11-Apr

21-May

24-Jun

33.0
4.5
5.3
7.5*
4.7
5.22
4.1

12.6
33.3
4.5
14.9
8.4
3.3
4.8*
3.0

29.9
15.2
11.5
10.7*
8.6
10.5
4.4*

30.3
20.5
10.0
4.6*

25.0
18.6
13
7.1
9.1
4.1*
7.0
9.3
7.4

27.5
17.3
13.7
5.9
15.5
6.7
5.8
12.1*
8.7

33.5
15.9
8.8
6.0*
13.7
7.5
9.1*
8.4
6.0

35.7
17.8
16.5
8.2*
5.1
4.5
6.5
8.3
8.5

2022
Steph 1
Ridge Run
Lower Cain
Steph 2
Steph 3
Steph 4
Steph 6
Steph 7
Steph 8

LiDAR

9-Mar

FIELD

Table D 2 Snow depth variability in cm from the weather station snow courses for 2023

STD DEV
(cm)

10-Mar

4-Apr

29-Apr

25-May

12-Jun

22.9
12.4
6.9

16.8
18.9
7.7
6.5
10.8
3.8
5.4
7.5*
4.9*

25.5
15.8
7.7
6.3
10.5
4.1
7
8.3*
4.1

22.3
16.5
8.5
6
9.6
7.8
8.6
6.3
4.2

21.1
19.8
11.4
5.1
1.6
4.8
5.4
7.5
4.4

2.7
5.7
7.6
3.5
2.8
3.3
9.3
8.6
5.1

12-Jun

29.5 - - 4.2
7.3
6.3
7.7
5.2
5.9
5.5
4.4
1.8*
1.6

4-Apr

25-May

16.9
3.9
5.4
5
3.1
2.3*
1.7*

LiDAR

29-Apr

Steph1
ridge run
lower cain
steph2
Steph3
steph4
steph6
steph7
steph8

10-Mar

2023

FIELD

133

Table D 3 Snow depth variability in cm from cardinal plots for 2022

21-May

24-Jun

9-Mar

11-Apr

21-May

24-Jun

ST1A1
ST2A2
Alpine RRA1
ST1OG3
ST1OG7
ST1OG8
Old
RROG3
Growth RROG4
ST1OG4
D. Old ST1OG5
Growth ST1OG6
ST4RF1
ST4RF2
Regen. ST6GC1
ST2CC1
Juvenile ST2CC2

LiDAR

11-Apr

2022

Field

9-Mar

STD. DEV. (cm)

25.5
34.7
44.8*
-

46.9
46.4
30
29.1
28.7
40.7*
40.3*
55.3
51.2
30.2*

34.6
44.5*
57.8
41.7
44.3
17.3*
28
25.5
30.1*
46.7*
29*
46.4*
46.5

27.8
46.9*
68
44.1
36.8
39.9
56.4
31.7
13.6
19.8*
20.1

32.1
35.2
43.3
51.8
33.5
34.7*
42
30.7
20.5
23.6
18.2
80*
66*
75.6*
58.2
53.5

36.2
34.2
44.2
61.5
44.3
41.3
44.2
31.9
23.7
26.1
22.7
85.2*
62.6*
75.8
62.3
50.9*

36
28.5*
58.1
52.3
39.7
36.4
40.6
34.7
23.7*
26
20.7
78.4*
67.7*
70*
61*
51.1

34.4
32.5*
66.7
41
36.7
32.6
49
34.8
18.5
16.8*
16.4
93.1
88.1
81.7
60.4*
63.6*

24.6
21.3
54.2*
37.6*
33.2*
-

*5.4
17.6*

Table D 4

Table D 5 Snow depth variability in cm from cardinal plots for 2023

134

12-Jun

40.2
43.4*

25-May

28.2
34.9*
64.9
26.2
21.1
16.3

29-Apr

-

4-Apr

37.3
38
42.3*
37.1
51.4*

10-Mar

30
30.8
39.2
26
33.2

12-Jun

25-May

Old
Growth

29-Apr

Alpine

ST1A1
ST2A2
RRA1
ST1OG3
ST1OG7
ST1OG8

LiDAR

4-Apr

2023

Field

10-Mar

STD. DEV. (cm)

32.1
28.8
45
34.7
26.4
26.3

32.9
30.5
43.8
27.2*
28.4
31*

33.5
32.8
56.9
33.5
29.3
25.6*

17.2
26.9*
54.7
25.5
23.3
21.1

10.8
22.9
38.3
10.8
10.9
19.5

ST1OG4
D. Old ST1OG5
Growth ST1OG6
ST4RF1
ST4RF2
Regen. ST6GC1
ST2CC1
Juvenile ST2CC2

21.3*
21.3
22.4*
24*
27.5*
29.7*
25.4*
23*

15.8
29.4*
19.1*
30.9*
27.2*
35.5
34.3*
26*

21.3*
25.9
18.2
33.2*
34.3*
27.9*
40.8*
34.6*

5.1*
5.8*

135

13.5*
17.9
14*
83.3*
61.1*
71.4*
58.1*
57.4*

7.1
12.3*
12.8*
76.1*
62.3*
69.9
57.7*
54.4*

13.9*
20.1
16.4
72.2*
63.6*
74.2*
59.4*
52.6*

10.8
10.2*
11.4*
72.2
65
67.7
55
58

9.8
10.2
9.1
32.5
59.7
41.2
80.5
46

Appendix E: Old Growth Forest: Canopy Snow Depth Variability
Generally, canopy gaps in forests will accumulate more snow than under the trees, due to the tree
intercepting snow and subsequent sublimation. For plots such as old growth forests with large canopy
gaps, this was captured in both the field and LiDAR measurements, with higher depths in the open, and
lower depths under the canopy.

Figure E 1 Snow depth distribution of old growth forest plot ST1OG7 showing the field (orange) and LiDAR (purple) snow depths
categorized by open, edge, or under canopy for 2022 snow surveys

136

Appendix F: Average snow depth results tables

Table F 1 2022 mean snow depth from snow courses for field and LiDAR measurements. Blue highlight indicates the plot was snow free

Mean Snow
Depth (cm)

21-May

24-Jun

9-Mar

11-Apr

21-May

24-Jun

LiDAR

11-Apr

Steph 1
Ridge Run
Lower Cain
Steph 2
Steph 3
Steph 4
Steph 6
Steph 7
Steph 8

212.4 290.7
186
175.8
245.2
78.4
128.8
89.6
130.2
91.7
122.3
51.8
41.4
45.4
18.6

308.6
281
206.6
212.8
65.4
69.6
52.4

181
146.1
80.7
79.9

212.6
223
168.5
175.7
77.7
83.6
86.4
60.1
46

294.6
284
216.8
224.8
114.8
118.3
101.5
35.2
7.6

293.5
274.5
210.8
201.8
58.7
66.3
42.8
-2.6
-11.8

124.5
108.4
38.8
27.5
13.9
0.9
5.7
24.8
4.6

9-Mar

2022

FIELD

137

4-Apr

29-Apr

25-May

12-Jun

81
65.9
22.9

12-Jun

194.4 158.4
169.6
69.5
66.6
81.2
74.9
76.7
73.2
21.6
27.8

LiDAR

10-Mar

204.4
163.9
80.5
91.4
95
52.6
52.6

25-May

Steph1
Ridge Run
Lower Cain
Steph 2
Steph 3
Steph 4
Steph 6
Steph 7
Steph 8

29-Apr

2023

FIELD

4-Apr

Mean Snow
Depth (cm)

10-Mar

Table F 22023 mean snow depth from snow courses for field and LiDAR measurements. Blue highlight indicates the plot was snow free

189.5
209.7
171.5
160.6
76.2
75.1
78
43.7
42.1

207.4
217.2
173.1
156
70.8
66.9
64.5
21.4
17.3

190.4
226.9
172
154.7
48.9
45.5
49.1
6.5
-2.4

56.4
45.4
21.5
1.3
5.4
0.1
-0.5
5.3
-8.1

-4
-4.1
6.1
-8.3
-3
-4.5
-4.9
1.9
-10.1

Table F 3 Mean difference in centimeters for snow courses snow depth between field and LiDARs. Grey highlight indicates bias +/- 10 cm,
red is LiDAR under-estimated by more than 10 cm, and green highlight indicates LiDAR overestimated by more than 10 cm. Blue highlight
indicates snow free. Significantly different average snow depths between the field and LiDAR (p<0.05) are bolded with an *.

Mean
difference

11-Apr

21-May

24-Jun

10-Mar

4-Apr

29-Apr

25-May

12-Jun

2023

9-Mar

2022

Steph 1

0.2

-

-15.1

-56.5*

-14.9*

13

-

-24.6*

-4

Ridge Run

-

-6.7

-6.5

-38*

-

-

-

-20.5*

-4.1

Lower Cain

-

30.8*

4.2

-41.9*

-

-

-

-1.4

6.1

Steph 2

0.1

-20.4*

-11*

-52.4*

-3.3

-2.4

-14.9*

1.3

-8.3

Steph 3

0.7

-14*

-6.7

13.9

-4.3

1.3

-17.7*

5.4

-3

Steph 4

-6*

-11.9*

-3.3

0.9

-16.3*

-14.3*

-29.4*

0.1

-4.5

Steph 6
Steph 7

5.3
8.3*

-20.8*
-6.2

-9.6*
-2.6

5.7
24.8

-17*
-8.9*

-12.2*
-0.2

-24.1*
6.5

-0.5
5.3

-4.9
1.9

Steph 8

0.6

-11*

-11.8

4.6

-10.5

-10.5

-2.4

-8.1

-10.1

Table F 4 2022 mean snow depth from the cardinal plots for field and LiDAR measurements. Blue highlight indicates the plot was snow
free

21-May

24-Jun

Juvenile

11-Apr

Regen.

9-Mar

D. Old
Growth

ST1OG8
146.8
ST1OG4
ST1OG5
79.3
ST1OG6
93.1
ST4RF1
103.6
ST4RF2
84.2
ST6GC1
64.4
ST2CC1 ST2CC2 -

24-Jun

Old
Growth

LiDAR

21-May

Alpine

ST1A1
167.4
ST2A2
RRA1
ST1OG3 ST1OG7
163.6

172.8
191.5

229.2

157.6
201.2

69.3
71.1
206.9
88.1
94

160.4
171.3
218.4
170.9
158.8

213.7
222.9
302
223.1
205.9

222.8
232
406.1
240.1
193.2

34.8
61.1
157.7
61.6
59.8

223.9
112.6
115.3
125.6
91.4
150.5
104.8
167.3
194.7

210.4
120.1
137.6
140.4
28.1
48
14.8
147.8
117.5

94.4
16.5
23
25.5
0
0
0
27.6
13.9

151.6
69
72
79.2
143.4
191.1
120.4
169.8
284

198.3
96.1
105
107.7
147.4
216.6
117.3
196.8
316.2

185
116.1
116.9
112.6
65.4
147.6
44.1
150.5
251.5

47.2
11.8
15.5
19
162.9
118.3
76.8
68.7
164.8

11-Apr

2022

Field

9-Mar

Mean Snow Depth
(cm)

-

-

138

Table F 5 2023 mean snow depth (cm) from cardinal plots for field and LiDAR measurements. Blue highlight indicates the plot was snow
free

139

12-Jun

Juvenile

198.2
171.2
158.8
156.1
46.4
75.8
29.8
155
137.6

25-May

Regen.

178
155.2
148.2
136.5
62.9
72.8
36.2
151.5
141.6

29-Apr

22.1
0
1.6
2.2
0
0
0
0
0

175.9
143.9
143.5
130.3
90.8
106
59.6
156.4
152.4

4-Apr

D. Old
Growth

ST1OG8
ST1OG4
ST1OG5
ST1OG6
ST4RF1
ST4RF2
ST6GC1
ST2CC1
ST2CC2

10-Mar

25-May

Old
Growth

165
208.9

LiDAR

12-Jun

29-Apr

15.7
35.6
198.6
24.4
25.8

4-Apr

161.2 216.7 228.2
218.6 169.7
172.5
177.1

10-Mar
Alpine

ST1A1
ST2A2
RRA1
ST1OG3
ST1OG7

Field

0
0
0
0

148.4
183.3
291
221.9
156.5

152.5
197.6
302.8
207.7
114.9

140
172.5
335.6
224.1
164.4

11.8
33.3
134.9
3.3
11.2

-4.9
0.5
30.9
-16.4
-14

0
0
0
0
0
0
0
0
0

151.5
138.3
128.2
118.9
132.8
175.1
98.4
178.7
275.6

74.9
138.6
129.3
116.4
115.9
158.3
80.2
182
276.4

164.1
141.2
131.8
122.2
87.7
131.3
59
172.2
261.8

-7.2
-11.9
-11.8
-10.7
45.3
75.2
25.4
31.5
127.5

-22.9
-9.1
-12.2
-15.6
-27.6
-40
-25
-9.8
-27.8

Table F 6 Mean difference in centimeters for cardinal plots for snow depth between field and LiDARs. Grey highlight indicates bias +/- 10
cm, red is LiDAR under-estimated by more than 10 cm, and green highlight indicates LiDAR overestimated by more than 10 cm. Blue
highlight indicates snow free. Significantly different average snow depths between the field and LiDAR (p<0.05) are bolded with an *.

11-Apr

21-May

24-Jun

10-Mar

4-Apr

29-Apr

25-May

12-Jun

2023

9-Mar

2022

Mean difference
(cm)
ST1A1
Alpine ST2A2
RRA1
ST1OG3
Old
ST1OG7
Growth ST1OG8
ST1OG4
D. Old
Growth ST1OG5
ST1OG6
ST4RF1
Regen. ST4RF2
ST6GC1
ST2CC1
Juvenile
ST2CC2

-6.97
-4.8 -23.4*
4.8 -25.7*
- -16.6*
-7.3 -10.3*
-13.9* -17.9*
39.8* 56.1*
106.9* 66.1*
56.0
12.5
29.6
- 121.5*

50*
40.5*
82.5*
-8.1
-25.3*
-3.97
-20.7*
-27.9*
37.2
99.7
29.2
2.7
134.0*

-34.5*
-9.9*
-49.26*
-26.6*
-34.3*
-47.2*
-4.63
-7.5
-6.53
162.9
118.3
76.8
41.2
151.0*

-16.7*
-25.5*
-6.3
-13.2*
-24.4*
-5.6*
-15.3*
-11.4*
42*
69.1*
38.8*
22.3
123.2*

-8.7
-19.06*
-10.9*
-57.6*
-103.1*
-16.6*
-18.9*
-20.1*
53*
85.5*
44*
30.5*
134.8*

-12.7*
-34.1*
-30*
-27*
-33.9*
41.3*
55.5*
29.2
17.2
124.2*

-3.8*
-2.35
-63.7*
-21.1*
-14.6*
-29.3*
-11.9
-13.4*
-12.9*
45.3
75.2*
25.4
31.5
127.5*

-4.86
0.49

140

-16.4
-14
-22.9
-9.1
-12.2
-15.6
-27.6
-40
-25
-9.8
-27.8

Appendix G: Bias corrections for meteorological variability
Table G 1 Meterological bias corrections applied in the bias assessment of the cardinal plots, as shown in results section 5.4.1.2

Year

Acquisition Cover
Date

Mar-9 Regenerating
Apr-11
May-21

Juvenile
Alpine

2022

Alpine

Dense old
growth
Jun-24

ST4RF1

Snow
Weather Station (Wx) /
depth
Game Camera (GC)
change station
[cm]
GC: ST6GC1
+5

ST4RF2
ST2CC1

+11
+16

ST2CC2
ST1A1

+16
-14

ST1A2
ST1A1
ST1A2
ST1OG4
ST1OG5

-14
+43

GC: ST6GC1
Wx: Steph 2
Wx: Steph 1

dense old
growth

ST1OG3
ST1OG7
ST1OG8
ST1CC1
ST1CC2
ST1OG4
ST1OG5

old growth

ST1OG6
ST1OG3
ST1OG7

Juvenile

+10

GC: ST1GC1

+31

Average of Wx: Steph 1 &
GC: ST1GC1

+19

GC: ST1GC3

+30

Wx: Steph 2

+12.5

GC: ST1GC1

GC: ST1GC3
+13

ST1OG8
ST4RF1

Apr-29
regenerating

ST4RF2
ST6GC1

GC: ST1GC1

ST2CC2

Snowfall in plots
sampled after
LiDAR
LiDAR flown after,
melt occurred at
plots
LiDAR flown, then
warming prior to
sampling

Delayed LiDAR
flight, melt occurred;
Snow entirely
melted at ST1GC1,
approximation of
melt prior
Delayed LiDAR
flight, melt occurred

Delayed LiDAR
flight, melt occurred
Rapid snow melt
occurred between
sampling and
LiDAR
Rapid snow melt
occurred between
sampling and
LiDAR
no enough snow at
closest station
ST6GC1 for relative
change

+12.5

ST2CC1
juvenile

Notes

Wx: Steph 1

ST1OG6
Old growth

2023

Plot ID

Wx: Steph 2
+16

141

Rapid snow melt
occurred between
sampling and
LiDAR

Total

FSR

New Harvest

Juvenile

Regenerating

Mature Regen.

Old Growth

Alpine

Bare Earth

Date Flown
2023, ACO

32,801,298

11,814,549
(36.0%)
13,007,977
(39.7%)
1,299,083
(4.0%)
4,823,783
(14.7%)
1,572,433
(4.8%)
195,834
(0.6%)
87,639
(0.3%)
34,185,623

13,712,941
(40.1%)
11,230,885
(32.9%)
1,070,746
(3.1%)
5,947,794
(17.4%)
2,059,789
(6.0%)
73,387
(0.2%)
90,081
(0.3%)

m³ (percent of total snow
volume)

2017,
BCTS

3/9/2022
2023, ACO

42,734,885

14,587,480
(34.1%)
16,280,662
(38.1%)
4,669,245
(10.9%)
5,330,758
(12.5%)
1,600,730
(3.8%)
183,021
(0.4%)
82,989
(0.2%)

142

44,305,496

16,484,398
(37.2%)
14,552,546
(32.9%)
4,516,964
(10.2%)
6,512,690
(14.7%)
2,072,641
(4.7%)
75,633
(0.2%)
90,624
(0.2%)

m³ (percent of total snow
volume)

2017, BCTS

4/11/2022
2023, ACO

36,696,751

15,258,007
(41.6%)
13,875,470
(37.8%)
3,077,733
(8.4%)
3,402,653
(9.3%)
910,329
(2.5%)
126,205
(0.3%)
46,354
(0.1%)
38,622,767

17,163,290
(44.4%)
12,308,954
(31.9%)
2,973,875
(7.7%)
4,651,516
(12.0%)
1,410,716
(3.6%)
58,915
(0.1%)
55,501
(0.1%)

m³ (percent of total snow
volume)

2017, BCTS

5/21/2022

2023, ACO

23,730,500

8,963,691
(37.8%)
9,027,503
(38.0%)
1,392,165
(5.9%)
3,301,467
(13.9%)
850,146
(3.6%)
169,395
(0.7%)
26,133
(0.1%)

24,970,084

10,849,364
(43.5%)
7,172,312
(28.7%)
1,161,106
(4.7%)
4,384,568
(17.6%)
1,318,570
(5.3%)
50,404
(0.2%)
33,760
(0.1%)

m³ (percent of total snow
volume)

2017, BCTS

6/24/2022

Table H- 1 Snow Storage in the Russell Creek watershed by cover type for the 2022 LiDAR acquisitions. Snow is shown in cubic meters, with the proportion of the watershed snow shown in
brackets

Appendix H: Watershed Snow Volume

Total

1400 - 1700

1200 - 1400

1000 - 1200

800 - 1000

600 - 800

400 - 600

200 - 400

Bare Earth
Elevation
Band [m]

Date Flown
2023, ACO

432,489
(1.3%)
3,327,715
(10.2%)
3,425,212
(10.4%)
4,706,684
(14.3%)
5,641,676
(17.2%)
9,030,259
(27.5%)
6,237,263
(19.0%)
32,801,298

438,489
(1.3%)
2,879,070
(8.4%)
3,314,036
(9.7%)
5,126,971
(15.0%)
5,231,501
(15.3%)
8,791,304
(25.7%)
8,404,252
(24.6%)
34,185,623

m³ (percent of total snow volume)

2017, BCTS

3/9/2022
2023, ACO

388,107
(0.9%)
6,057,999
(14.2%)
4,394,563
(10.3%)
5,577,343
(13.1%)
7,237,320
(16.9%)
11,333,253
(26.5%)
7,746,300
(18.1%)
42,734,885

143

402,362
(0.9%)
5,864,185
(13.2%)
4,270,943
(9.6%)
5,922,151
(13.4%)
6,833,379
(15.4%)
11,102,671
(25.1%)
9,909,805
(22.4%)
44,305,496

m³ (percent of total snow volume)

2017, BCTS

4/11/2022
2023, ACO

430,814
(1.2%)
3,513,020
(9.6%)
2,911,511
(7.9%)
3,502,837
(9.6%)
6,424,254
(17.5%)
11,773,919
(32.1%)
8,140,396
(22.2%)
36,696,751

440,917
(1.1%)
3,464,926
(9.0%)
2,983,642
(7.7%)
3,874,085
(10.0%)
6,013,327
(15.6%)
11,538,073
(29.9%)
10,307,797
(26.7%)
38,622,767

m³ (percent of total snow volume)

2017, BCTS

5/21/2022

2023, ACO

526,784
(2.2%)
3,478,685
(14.7%)
3,452,784
(14.6%)
2,817,824
(11.9%)
2,921,229
(12.3%)
5,806,577
(24.5%)
4,726,617
(19.9%)
23,730,500

534,366
(2.1%)
3,047,727
(12.2%)
3,329,063
(13.3%)
3,147,159
(12.6%)
2,477,293
(9.9%)
5,542,991
(22.2%)
6,891,485
(27.6%)
24,970,084

m³ (percent of total snow volume)

2017, BCTS

6/24/2022

Table H- 2 Snow Storage in the Russell Creek watershed by elevation band for the 2022 LiDAR acquisitions. Snow is shown in cubic meters, with the proportion of the watershed snow
shown in brackets

763,779
(2.2%)
5,934,958
(17.2%)
2,358,069
(6.8%)
116,560
(0.3%)
99,919
(0.3%)

1,131,221
(3.4%)

4,816,056
(14.4%)

1,883,514
(5.6%)

269,670
(0.8%)

101,461
(0.3%)

33,384,574

M. Regen.

Regenerati
ng

Juvenile

New
Harvest

FSR

Total
34,441,880

29,581,989

61,290
(0.2%)

139,569
(0.5%)

1,449,571
(4.9%)

3,208,766
(10.8%)

305,586
(1.0%)

12,656,902
(42.8%)

12,200,555
(35.4%)

14,129,254
(42.3%)

31,280,394

70,567
(0.2%)

51,545
(0.2%)

1,934,539
(6.2%)

4,424,920
(14.2%)

175,790
(0.6%)

10,953,036
(35.0%)

13,669,997
(43.7%)

Old
Growth

12,968,040
(37.6%)

11,760,305
(39.8%)

2023,
ACO

11,053,398
(33.1%)

2017,
BCTS

Alpine

2023,
ACO

4/4/2023

m³ (percent of total snow
volume)

2017,
BCTS

3/10/2023

m³ (percent of total snow
volume)

Bare Earth

Date
Flown

2023,
ACO

144

28,232,610

47,641
(0.2%)

129,859
(0.5%)

1,120,349
(4.0%)

2,523,775
(8.9%)

165,973
(0.6%)

12,217,025
(43.3%)

12,027,988
(42.6%)

30,116,482

58,079
(0.2%)

39,308
(0.1%)

1,611,113
(5.3%)

3,782,503
(12.6%)

56,104
(0.2%)

10,635,542
(35.3%)

13,933,833
(46.3%)

m³ (percent of total snow
volume)

2017,
BCTS

4/29/2023
2023,
ACO

4,813,295

6,659
(0.1%)

63,155
(1.3%)

244,117
(5.1%)

527,137
(10.9%)

178,223
(3.7%)

1,784,975
(37.1%)

2,009,029
(41.7%)

7,018,933

18,846
(0.3%)

-4,459
(-0.1%)

789,014
(11.2%)

1,893,150
(27.0%)

50,965
(0.7%)

340,755
(4.8%)

3,930,662
(56.0%)

m³ (percent of total snow
volume)

2017,
BCTS

5/25/2023

2023,
ACO

-5,991,728

-47,433
(0.8%)

56,506
(-0.9%)

-1,342,496
(22.4%)

-3,191,502
(53.3%)

-382,670
(6.4%)

-998,707
(16.7%)

-85,426
(1.4%)

-3,589,788

-38,031
(1.1%)

-21,484
(0.6%)

-826,744
(23.0%)

-1,822,858
(50.8%)

-450,824
(12.6%)

-2,345,842
(65.3%)

1,915,995
(-53.4%)

m³ (percent of total snow
volume)

2017,
BCTS

6/12/2023

Table H- 3 Snow Storage in the Russell Creek watershed by cover type for the 2022 LiDAR acquisitions. Snow is shown in cubic kilometers, with the proportion of the watershed snow
shown in brackets

2023, ACO

1,829,453
(5.3%)

2,925,779
(8.5%)

5,883,739
(17.1%)

6,715,517
(19.5%)

9,755,978
(28.3%)

7,211,324
(20.9%)

2,417,108
(7.2%)

3,254,608
(9.8%)

5,490,746
(16.4%)

7,120,985
(21.3%)

10,044,987
(30.1%)

4,952,145
(14.8%)

33,384,574

Total
34,441,880

120,090
(0.3%)

103,995
(0.3%)

m³ (percent of total snow
volume)

2017,
BCTS

3/10/2023

200 400
400 600
600 800
800 1000
1000 1200
1200 1400
1400 1700

Bare Earth

Date
Flown

2023,
ACO

29,581,989

5,237,139
(18.6%)

11,223,485
(39.8%)

6,592,089
(23.4%)

3,397,311
(12.0%)

1,235,261
(4.4%)

541,998
(1.9%)

5,327
(0.0%)

31,280,394

7,549,572
(25.1%)

10,858,248
(36.0%)

6,186,740
(20.5%)

3,736,937
(12.4%)

1,242,975
(4.1%)

500,897
(1.7%)

41,113
(0.1%)

m³ (percent of total snow
volume)

2017,
BCTS

4/4/2023
2023,
ACO

145

28,232,610

5,587,598
(18.9%)

10,389,886
(35.1%)

6,747,699
(22.8%)

4,262,719
(14.4%)

1,891,989
(6.4%)

774,689
(2.6%)

-72,591
(-0.2%)

30,116,482

7,857,542
(25.1%)

10,061,308
(32.2%)

6,360,540
(20.3%)

4,590,145
(14.7%)

1,775,079
(5.7%)

681,567
(2.2%)

-45,787
(-0.1%)

m³ (percent of total snow
volume)

2017,
BCTS

4/29/2023
2023,
ACO

4,813,295

-449,085
(-9.3%)

2,417,601
(50.2%)

1,144,933
(23.8%)

453,504
(9.4%)

754,800
(15.7%)

517,407
(10.8%)

-25,865
(-0.5%)

7,018,933

1,950,138
(27.8%)

1,973,042
(28.1%)

818,808
(11.7%)

1,014,587
(14.5%)

787,269
(11.2%)

471,920
(6.7%)

3,169
(0.1%)

m³ (percent of total snow
volume)

2017,
BCTS

5/25/2023

2023, ACO

-5,991,728

-1,056,905
(17.6%)

762,614
(-12.7%)

-371,459
(6.2%)

-1,976,815
(33.0%)

-1,479,970
(24.7%)

-1,582,683
(26.4%)

-286,510
(4.8%)

-3,589,788

1,385,770
(-38.6%)

396,250
(-11.0%)

-668,139
(18.6%)

-1,416,145
(39.5%)

-1,385,459
(38.6%)

-1,630,707
(45.4%)

-271,358
(7.6%)

m³ (percent of total snow
volume)

2017,
BCTS

6/12/2023

Table H- 4 Snow Storage in the Russell Creek watershed by cover by elevation band for the 2023 LiDAR acquisitions. Snow is shown in cubic meters, with the proportion of the watershed
snow shown in brackets

Appendix I: Harvest Experiment Slash Pile Analysis
Further analysis of the cut block surface debris was prompted by the result of a positive increase in
ground surface post-harvesting. A potential source of error would be a misclassification of logs and
slash piles as a ground. If this occurred in numerous parts of the cut block, this would result in a false
positive bias to the ground surface.
Classification of the cut block into slash and ground using a basic supervised classification of the
orthoimagery. Shadow caused by the trees were classified separately and not included in the analysis.
The classified raster was resampled to the LiDAR 1 meter resolution, and surface elevations were
binned according to ground or slash.

Figure I 1Result of classification of cut block as ground (yellow) or slash (brown)

The summary statistics and plotted histogram do not show evidence of a bias on logs compared to the
ground. Both the slash and ground groups show a similar distribution and positive bias when comparing
the 2023 ground surface to the 2022.
146

Figure I 2 Histogram of the ground surface difference between the post-harvest (2023) and pre-harvest (2022) pixels, showing the peak for
both the slash and ground pixels are aligned, and slightly above zero

147

148