Characterizing the spatial distribution of field-scale snowpack using unpiloted aerial system (UAS) lidar and structure-from-motion (SfM) photogrammetry

Unpiloted aerial system (UAS) light detection and ranging (lidar) and structure-from-motion (SfM) photogrammetry have emerged as viable methods to map high-resolution snow depths (∼1 m). These technologies enable a better understanding of snowpack spatial distribution and its evolution over time, advancing hydrological and ecological applications. This is particularly critical in mixed vegetation environments, where both forest canopy and open areas influence snow accumulation and melt patterns. In this study, a series of UAS lidar/SfM snow depth maps were collected during the 2020/2021 winter season in Durham, New Hampshire, USA, with three objectives: (1) quantifying UAS lidar/SfM snow depth retrieval performance using in situ magnaprobe measurements, (2) conducting a quantitative comparison of lidar and SfM retrievals of shallow snow depths (<35 cm) throughout the winter, and (3) understanding the spatial distribution of snow depth and its relationship with terrain features. Eight UAS surveys were conducted over approximately 0.35 km² including both open fields and a mixed forest. In the field, lidar had a slightly lower error than SfM, compared with in situ observations, with a mean absolute difference (MAD) of 3.5 cm for lidar and 4.0 cm for SfM. Snow depth maps from SfM and lidar were fairly consistent in the field, with only marginal differences on most dates. In the forest, SfM greatly overestimated in situ snow depths compared with lidar (lidar MAD = 6.3 cm, SfM MAD = 31.4 cm). There was no clear agreement between SfM and lidar snow depth values for individual 1 m² pixels in the forest (MAD = 55.7 cm). Using the concept of temporal stability, we found that the spatial distribution of snow depth captured by lidar was generally consistent throughout the period, indicating a strong influence from static land characteristics. Considering both areas (forest and field), the spatial distribution of snow depth was primarily influenced by vegetation type while also reflecting the effects of soil variables (e.g., soil organic matter). When the field and forest areas were analyzed separately, the spatial distribution was distinctly affected by slope and the shadowing effects of the forest canopy.

1 Introduction

Snowpacks are vital to hydrological, climatic, and ecological processes across multiple scales (Barnett et al., 2005; Clark et al., 2011). An understanding of snowpack distribution and its temporal evolution is important to determine snowmelt runoff, infiltration, and groundwater recharge (Carroll et al., 2019; Harpold et al., 2015; Maurer and Bowling, 2014), as well as energy partitioning processes (Lawrence and Slater, 2010; Stieglitz et al., 2001; Sturm et al., 2017). Snowpacks also exert a strong control on snow–soil interactions because the insulating capacities of snowpacks affect the underlying soil freeze–thaw state, influencing soil respiration, nutrient retention, and carbon dynamics (Anderton et al., 2002; Schlögl et al., 2018; Cho et al., 2021; Monson et al., 2006; Sorensen et al., 2018; Reinmann and Templer, 2018; Wilson et al., 2020; Yi et al., 2015).

The spatial variability of a snowpack is a function of static (e.g., slope, aspect, vegetation type, soil properties) and dynamic (e.g., solar radiation, wind direction and speed, temperature) variables and fluxes over a range of spatial scales (Clark et al., 2011; Grayson et al., 2002; Mott et al., 2011; Trujillo et al., 2007). Over time, spatial patterns may evolve and change, but many hydrological patterns persist until they are modified by weather conditions. Spatial snowpack patterns and their consistency, or repeatability, play a crucial role in various applications, including operational snowmelt predictions, the downscaling of remotely sensed or model outputs, the integration of in situ observations through upscaling, the assimilation of data to enhance model simulations, and the utilization of snowpack characteristics as proxies or analogs for similar hydrological units (Pflug and Lundquist, 2020; Cho et al., 2023). They can also provide insight into underlying landscape features, biogeochemical processes, and wildlife habitats (Boelman et al., 2019; Pflug et al., 2023).

Previous studies have proposed diverse approaches to characterize snow distribution patterns and their temporal evolution across a range of climatic and topographic settings (Sturm and Wagner, 2010; Vögeli et al., 2016; Revuelto et al., 2020; Pflug et al., 2021). For example, Sturm and Wagner (2010) found that snow depth patterns in an Arctic region remain stable across years due to persistent topographic and vegetation influences, highlighting the value of empirical snow distribution patterns for improving snow model accuracy. Vögeli et al. (2016) used high-resolution airborne digital sensors to refine precipitation scaling in a snow distribution model (Alpine3D), demonstrating the potential of remote sensing data to better simulate complex snow dynamics in alpine regions. Pflug et al. (2021) examined the interannual consistency of snow patterns and proposed a downscaling approach based on historical snow patterns in the California Tuolumne River Watershed; this approach is particularly useful for predicting snow distribution in years with limited observations. Revuelto et al. (2020) introduced a method combining in situ snow depth measurements with terrestrial laser scanner and time-lapse photography to produce temporal snow depth distribution patterns in a subalpine mountain environment, offering a transferable approach for deriving spatial snow data from limited ground observations.

Traditionally, field (approximately 100 m) or local-scale (approximately 1 m) snow features are captured through in situ observations and field campaigns (Clark et al., 2011; Trujillo et al., 2007), whereas regional or continental-scale patterns are typically observed using airborne and satellite remote sensing techniques (Lievens et al., 2022; Painter et al., 2016; Derksen et al., 2005). Airborne and satellite remote sensing methods have provided the ability to collect snowpack data over a large spatial extent, thus expanding the understanding of snow distribution (Cho et al., 2020; Lievens et al., 2022; Painter et al., 2016; Tsang et al., 2021). However, the ability to capture small-scale snow patterns, discerned through field campaigns or less frequent routine operational collections, is often hindered by challenges, such as weather conditions, tree canopies, and site accessibility, which can lead to infrequent sampling during the winter season.

Unpiloted aerial systems (UASs) have been used to provide spatially continuous opportunistic snow-covered area and snow depth observations at scales between in situ and airborne and satellite remote sensing (Bühler et al., 2016; De Michele et al., 2016; Harder et al., 2016, 2020; Meyer et al., 2022; Revuelto et al., 2021a; Geissler et al., 2023). UAS-based remote sensing enables the acquisition of data at fine spatial resolutions, reaching scales as precise as centimeters for a designated area. UAS platforms also offer a cost-effective alternative to aerial surveys, facilitating routine monitoring of snow conditions (Gaffey and Bhardwaj, 2020). Hence, the capabilities of UAS platforms equipped with diverse sensors can observe snowpack properties and support analyses of field-scale physical interactions between snowpacks and land/soil characteristics (Cho et al., 2021).

UAS light detection and ranging (lidar) and structure-from-motion (SfM) photogrammetry have emerged as viable methods for mapping high-resolution snow depths (∼1 m), enabling a better understanding of snowpack spatial distribution and its evolution over time at the field scale (Feng et al., 2023; Harder et al., 2016; Jacobs et al., 2021; Koutantou, 2022; Geissler et al., 2023). There is a growing need for new multi-temporal snow datasets, especially during transition periods characterized by shallow and patchy snow conditions, which pose measurement challenges but are critical for understanding snowpack dynamics and their hydrological, ecological, and energy implications (Harrison et al., 2021; Harpold et al., 2017; Grogan et al., 2020; López-Moreno et al., 2024). Different UAS-based sensors offer complementary strengths and weaknesses that warrant further investigation across diverse environmental settings.

This study aims to achieve three main objectives using a series of UAS lidar/SfM snow depth maps over a mixed-use temperate forest landscape: (1) to quantify UAS snow depth retrieval performance by comparing it with in situ measurements; (2) to conduct a quantitative comparison of lidar and SfM snow depths throughout the snow period for a range of depths that reflect the specific conditions observed in our dataset (i.e., 0 to 35 cm); and (3) to gain a better understanding of the spatial distribution of snow depth, its stability over time, and its relationship with physical terrain features. This paper is organized as follows. Section 2 provides an overview of the study area, including its land characteristics. Section 3 describes the datasets utilized in the study, including UAS lidar, SfM photogrammetry, and field observations, and the methods employed, such as the relative difference concept. Section 4 presents the results, with subsections detailing comparisons between UAS snow depth and in situ measurements (Sect. 4.1), as well as comparisons between lidar and SfM snow depth (Sect. 4.2). Sections 4.3 and 4.4 further examine the spatial patterns and temporal dynamics of snow depth, along with the physical variables influencing these patterns. Section 5 discusses new insights derived from the comparison results and spatial patterns of snow depth, along with the limitations of this study and future perspectives. Finally, conclusions are drawn in Sect. 6.

2 Study area

This study was conducted at the University of New Hampshire Thompson Farm Research Station in southeastern New Hampshire, the United States (43.10892° N, 70.94853° W, 35 m above sea level), which was chosen for its mixed hardwood forest and open field land covers (Cho et al., 2021; Burakowski et al., 2015; Jacobs et al., 2021), which are characteristic of the region (Fig. 1). Thompson Farm has a rich history of forest ecology research and data collection. Thompson Farm has an area of 0.83 km² and little topographic relief (18 to 36 m a.s.l.). The agricultural fields are actively managed for pasture grass. The deciduous, mixed, and coniferous forest is composed primarily of white pine (Pinus strobus), northern red oak (Quercus rubra), red maple (Acer rubrum), shagbark hickory (Carya ovata), and white oak (Quercus alba). The forest soils are classified as Hollis/Charlton very stony fine sandy loam and well-drained; field soils are characterized as Scantic silt loam and poorly drained. There are two logging access roads running north–south through the pasture and western forest section. The winter climate at Thompson Farm has a mean winter air temperature of −3.0 °C and an annual snowfall of 114 cm, with 3 weeks to over 3 months of days with snow cover (Burakowski and Hamilton, 2020; Johnston et al., 2024). Snow depth can range from a trace up to 94 cm and typical snow density ranges from 100 to 400 kg m⁻³ (Burakowski and Hamilton, 2020). The snowpack at Thompson Farm is short-lived and warm, and snow climatologies from Sturm and Liston (2021) and Johnston et al. (2024) both classify the area as ephemeral. A review of existing research on the snow classes defined by Sturm and Liston (2021) and Johnston et al. (2024) determined that, despite covering large areas in the northern hemisphere, the ephemeral snow class is largely understudied, making new research on ephemeral snowpacks valuable. López-Moreno et al. (2024) also emphasized the importance of studying these transient snow conditions, highlighting their sensitivity to climate variability and their implications for hydrological and ecological processes.

Figure 1 Thompson Farm survey area located in Durham, NH, USA. (a) In situ sample sites and field and forest boundaries are overlain on the snow-off imagery. The pond, section of dense shrubs, outbuildings, and USCRN station in the western field were not representative of the field and were removed. (b–g) Maps of (b) vegetation type, (c) slope, (d) shadow hours, (e) soil hydraulic conductivity (K_sat), (f) organic matter, and (g) aspect are shown for the field and forest areas. The derivation of each of these variables is explained in Sect. 3.

3 Datasets and methods

Series of UAS lidar surveys, UAS SfM photogrammetry surveys, and in situ sample campaigns were conducted at Thompson Farm during the winter 2020–2021. Eight snow-on campaigns were conducted between 10 February and 11 March 2021, during which period UAS lidar, UAS SfM photogrammetry, and in situ data were collected (Table S1). The UAS snow-on surveys were conducted prior to in situ sampling on each of the campaign dates. Because compaction of underlying vegetation during snow-on periods can result in negative snow depths when compared with the snow-off baseline (Masný et al., 2021), the snow-off baseline survey was conducted on 2 April 2021 following snowmelt.

3.1 UAS lidar and SfM photogrammetry

The lidar sensor payload consisted of a Velodyne VLP-16 laser scanner and an Applanix APX-15 inertial navigation system (INS: global navigation satellite system (GNSS) + inertial measurement unit (IMU)). The VLP-16 is a lightweight (≈830 g) low-power (≈8 W) sensor; this makes it ideal for UAS deployment. The sensor incorporates 16 rotating infrared (IR) lasers that are arranged and oriented on the payload to provide a 30° along-track field of view with a cross-track field of view, limited only by the range of the sensor (approximately 100 m). At an altitude of 65 m, the sensor range produces an effective cross-track field of view of approximately 98°. Each laser operates at a wavelength of 903 nm.

For these acquisition missions, the VLP-16 was hard-mounted to a DJI Matrice 600 to maintain constant lever arm offsets between the inertial navigation system (INS) GNSS antenna, the lidar sensor, and the INS board. As opposed to a gimbal-mounted system, this hard-mounted configuration achieves a more tightly coupled system, resulting in improved point cloud geolocation accuracy. The lidar sensor was set to dual-return mode to improve ground detection in the forested areas of the field site. The system was flown at an altitude of 65 m with a flight speed of 3 m s⁻¹ and ≈40 m spacing between flight lines. Flights produced ≈70–140 million returns per mission, depending on site ground conditions.

Lidar observations were georeferenced using position and attitude measurements acquired with the Applanix APX-15 INS. The INS produced 2–5 cm positional, 0.025° roll and pitch, and 0.08° true heading uncertainties following post-processing. Post-processing of INS data was performed using POSPac UAV (v. 8.2.1, Applanix Corporation 2018), correcting differentially against a permanent continuously operating reference station (CORS) at the University of New Hampshire in Durham, NH (NHUN). Position and attitude data were output as a smoothed best estimate of trajectory (SBET), then time synchronized with lidar returns to produce a georeferenced point cloud using LidarTools (v. 3.1.4, Headwall Photonics, Inc.).

Three-dimensional point clouds were processed using a progressive morphological filter (PMF) within the R programming language package “lidR” to identify ground returns. For ground classification, point clouds were chunked into 100 m square tiles with a 15 m buffer on all sides using catalog options in lidR to ensure that returns near tile edges were classified. The PMF was parameterized using a set of window sizes of 1, 3, 5, and 9 m, and elevation thresholds of 0.2, 1.5, 3, and 7 m, which were determined by varying value sets and assessing digital terrain models (DTMs) to determine the parameter sets that produced a visually smooth surface over a dense grid (Muir et al., 2017). Following ground classification for each tile, returns within the 15 m tile buffers were removed, and all resulting 100 m square ground classified tiles were merged. The result of using the PMF is that non-ground returns (i.e., trees, shrubs, and noise) were filtered out of the point cloud data sets, so that only returns from ground surfaces remained. The two data sets, non-ground returns and ground returns from the original point clouds, were coded according to LAS specifications and merged. Lidar snow depths were calculated as the difference between the ground classified snow-on and snow-off elevations within each pixel. For comparison with in situ observations, ground returns were extracted for the 1×1 m square sampling sites, corresponding to the alignment and orientation of the transect. The lidar snow depth was calculated as the difference between the mean snow-on and mean snow-off elevations within each sampling grid.

Photogrammetry bare-earth and snow-on elevation models were constructed from UAS-borne optical imagery. RGB images were collected using the DJI Phantom 4 real time kinematic (RTK) UAS platform equipped with a 20 megapixel complementary metal oxide semiconductor (CMOS) sensor. The RTK system integrates a static base station that relays GNSS corrections to the UAS, enabling approximately 3 cm accuracy of image geotags. To ensure that photogrammetry snow depth products aligned correctly with the UAS lidar products, the RTK base station was placed over a monument with known coordinates, which were entered into the DJI flight app. Flights were conducted at an altitude of 65 m AGL and a flight speed of 8 m s⁻¹. The shutter triggering interval was set to achieve a forward overlap of 80 % between image pairs and the flight lines were spaced to achieve 80 % side overlap. Three ground control points (GCPs) were placed within the area of interest to verify the accuracy of the photogrammetry products. The GCPs were surveyed using a Trimble© Geo7X GNSS positioning unit and Zephyr™ antenna with sub-centimeter accuracy.

The acquired image datasets were processed through the basic photogrammetry workflow using Agisoft Metashape (v. 1.8.4). Sparse clouds were constructed using the default key point and tie point limits of 40 000 and 4000, respectively. Points with high errors within the sparse clouds were then removed using the gradual selection tool. This included points exceeding the following thresholds: reprojection error >0.5, reconstruction uncertainty >50, and projection accuracy >5. The camera intrinsic/extrinsic parameters were optimized following the removal of the poorly localized points. Dense clouds were produced with the quality setting set to high and depth filtering set to moderate. Ground returns were classified using the ground classification tool within Metashape. A first pass at establishing a ground surface was done by triangulating the lowest point elevation within 50 m grid cells. The default thresholds for maximum distance and angle (1 m and 15°, respectively) of all points relative to the triangulated surface were used to determine points forming part of the ground surface. Finally, digital elevation models (DEMs) were derived from the ground classified points within the dense clouds. Snow depth products were derived following the same procedure as for lidar by calculating the difference between the ground classified snow-on and snow-off elevations within each pixel. Additional filtering based on the point confidence metric was completed for the 20 and 24 February snow depth maps to remove points with high uncertainty. GCPs surveyed using the base/rover equipment were used to co-register the UAS data. Linear, horizontal, and vertical shifts were applied to align all SfM and lidar DEMs to the GCPs.

3.2 Field observations

In situ snow depth sampling was conducted in the field and forest using a Snow-Hydro LLC magnaprobe (Sturm and Holmgren, 2018) and three Moultrie Wingscapes Birdcam Pro field cameras. The magnaprobe sampling followed a single long transect (18 points) and two short transects (3 points each). The long transect was approximately 145 m long and laid out from east to west (Fig. 1). From east to west, the transect started in the open field area, then transitioned to the coniferous, then mixed, and finally, deciduous forested areas. The two short transects were located in the open field, one in the northwestern portion and the other in the southeast. At each point, nine evenly spaced measurements were taken within 1 m × 1 m grid cells. It is worth noting that some dates were missing sample points due to disturbance of the sample area, either by collection on previous days or recreational use at the site or due to personnel and equipment limitations (Table S1). All sampling locations were geolocated using a Trimble Geo7X GNSS positioning unit and Zephyr antenna with an estimated horizontal uncertainty of 2.51 cm (standard deviation, 0.95 cm) in the field and 4.17 cm (standard deviation, 4.60 cm) in the forest after differential correction. Field camera snow depths were acquired following the method used in NASA's 2020 SnowEx field campaign in Grand Mesa, CO (personal communication, 16 November 2020). The three cameras were placed in different land cover types; one in the open field, one in the coniferous forest, and one in the deciduous forest. Each camera was mounted approximately 0.85 m above the ground and placed approximately 5.5 m from its respective 1.5 m marked PVC pole. Each PVC pole was spray-painted red and marked in 1 and 10 cm increments. The cameras captured images of the poles every 15 min for the duration of the study period. Snow depth was derived by manual inspection of the photos and recorded to the nearest centimeter. Precipitation equivalent and mean temperature were measured by a NOAA Office of Oceanic and Atmospheric Research U.S. Climate Reference Network (USCRN) station (NH Durham 2 SSW) located in the western portion of the field. Hourly air temperatures at the USCRN station were averaged from 2 s readings from three independent thermometers. Hourly precipitation is computed from 5 min readings of depth change, recorded by a weighing precipitation gauge.

3.3 Physical land characteristics

Land and soil characteristic variables were investigated as physical drivers of the field-scale spatial distribution of snow depth. The variables used in this study were plant functional type, slope, aspect, shadow hours, saturated hydraulic conductivity (K_sat), and soil organic matter (SOM) (Fig. 1). Mapped at a 1 m scale, all variables were derived from UAS snow-off observations except the two soil variables. The two soil variables, K_sat and SOM, were obtained at soil depths of 0–5 cm from probabilistic remapping of SSURGO (POLARIS) maps at 30 m spatial resolution (Chaney et al., 2016, 2019). The soil maps were disaggregated to 1 m spatial resolution without employing interpolation methods to mitigate additional uncertainties. Vegetation cover type (field/forest) was manually delineated in geographic information system (GIS) software based on the image orthomosaics created during SfM processing. The forested area was further classified as coniferous or deciduous for the study region by applying a green leaf index (GLI) (Louhaichi et al., 2001),

$$\begin{array}{}\text{(1)}& \mathrm{GLI}={\displaystyle \frac{\left(\text{green}-\text{red}\right)+\left(\text{green}-\text{blue}\right)}{\left(\mathrm{2}\cdot \mathrm{green}\right)+\mathrm{red}+\mathrm{blue}}},\end{array}$$

to the optical three-band (red, green, and blue) orthomosaics derived from the snow-off DJI Phantom 4 RTK survey.

The GLI algorithm delineated the dense vegetation (conifer trees) from the less dense vegetation (leaf-off deciduous trees). The direct application of the GLI algorithm to the three-band orthomosaics was further filtered and refined as follows. The output was clustered using the k-means algorithm with the number of k classes equal to 2: 1 class for coniferous trees (high GLI) and 1 class for deciduous trees (low GLI). Noise within the clustered GLI map was removed by convolution with a median filter. To establish continuous delineations of coniferous regions, morphological closing was applied to the map to fill in any interior holes within the delineated regions. Forest classifications for each of the magnaprobe sample locations were estimated for a 10 m × 10 m area centered at each sampling grid, based on the percentage of coniferous pixels (< 40% = deciduous, 40 %–60 % = mixed, > 60% = coniferous). Results from the binary forest classification and the coarsened 10 m × 10 m classification are shown in Fig. S1. The slope and aspect were derived from the UAS lidar 1 m snow-off DEM using Horn's method (Horn, 1981). The shadow hours were calculated using the unfiltered UAS lidar digital terrain model and a static sun incidence angle based on the average of 4 February and 7 March. Given the minor variation in solar angles between these dates, any change in shadow hours was considered negligible for this study.

3.4 Relative difference concept

The relative difference concept, first introduced by Vachaud et al. (1985), has been widely used in the soil moisture remote sensing community to quantify the spatiotemporal variability (or stability) of soil moisture at field or regional scales (Cho and Choi, 2014; Cosh et al., 2004; Jacobs et al., 2004; Mohanty and Skaggs, 2001; Starks et al., 2006). In this study, we apply this concept to the UAS lidar snow depth measurements. The relative difference in the snow depth measurements can be expressed as

$$\begin{array}{}\text{(2)}& {\mathrm{RD}}_{i,t}={\displaystyle \frac{{\mathrm{SND}}_{i,t}-\mathrm{spatial}\mathrm{mean}\left({\mathrm{SND}}_t\right)}{\mathrm{spatial}\mathrm{mean}\left({\mathrm{SND}}_t\right)}},\end{array}$$

where SND_i,t is the individual snow depth measurement at grid i and date t, and spatial mean (SND_t) is the spatial mean value of snow depth at date t. For each grid i, the mean relative difference (MRD_i) is the average relative difference from each of the N flights and can be calculated as

$$\begin{array}{}\text{(3)}& {\mathrm{MRD}}_i={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{t=\mathrm{1}}^N{\mathrm{RD}}_{i,t}.\end{array}$$

4 Results

4.1 In situ vs. UAS-measured snow depths

Daily temperature, daily precipitation, and measured snow depths in the field and forest for the study period are shown in Fig. 2. Between 15 December 2020 and 12 March 2021, the average daily temperature was −2 °C, the maximum daily temperature was 19 °C on 31 January, and the minimum daily temperature was −19 °C on 11 March. Average wind speed was 1.4 m s⁻¹ for the study period. The cumulative precipitation for the same period was 20.4 cm. The largest precipitation event was 11 cm and occurred on 16 January when temperatures were above freezing (2–9 °C). December and early January had ephemeral snowpacks of less than 10 cm, which melted within a week. A snowpack was continuously present from late January through the middle of March. The maximum snow depth measured by the cameras occurred on 10 February, with 21 cm in the field, 21 cm in the deciduous forest, and 18 cm in the coniferous forest. A second peak snow depth occurred on 20 February, with 20 cm in the field, 12 cm in the deciduous forest, and 19 cm in the coniferous forest. A sustained period of warm temperatures in late February and early March corresponded to a decrease in snow depth due to the warming temperatures and two rain-on-snow events. The snowpack was depleted from the entire study area by 10 March.

Figure 2 Time series of conditions at the Thompson Farm, Durham, NH, study area during the 2021 winter season including (a) hourly precipitation equivalent (mm) and temperature (°C) measured by a USCRN station and (b–c) daily camera snow depths and median UAS-measured snow depths in (b) the field and (c) forest. Dates corresponding to the in situ and UAS sampling campaigns are marked by the vertical dotted lines. Periods where the temperature was colder than 0 °C are indicated by the blue plot background in (a) and periods warmer than 0 °C are indicated in pink. Median SfM-measured snow depths in the forest on 4 and 28 February 2021 exceeded 35 cm (172 and 86 cm, respectively) and are not shown.

Figure 2b and c show that the eight UAS-based SfM and lidar flights captured both the snow depth peaks and the ablation period following the last peak on 20 February, the latter referring to the phase in the seasonal snow patterns when the snowpack begins to melt and decrease in depth. During the ablation period, UAS measurements typically showed a decreasing snow depth that matched the progression from the field camera measurements. The final UAS surveys on 7 March captured the transition from a snow-cover-dominated field area on 3 March to predominantly bare ground. In the field, the UAS-based measurement techniques yielded similar snow depths throughout the winter and were able to capture snow depth changes of the order of 5 cm or less. However, the forest performance often differed between the two UAS methods. In the forest, the lidar snow depths tracked the camera observations within 5 cm. However, the SfM method had inconsistent performance. Notably, on 14 and 20 February, the SfM snow depths exceeded the lidar depths by more than 50 cm.

UAS-based snow depth retrievals were compared with in situ snow depths measured by the magnaprobe (Fig. 3; Table S2). Field camera observations were not used for validation due to their limited spatial extent compared with the magnaprobe sample locations (Fig. 3). In the field, both UAS lidar and UAS SfM snow depths were approximately 1.5 cm deeper on average than the magnaprobe measurements. While both UAS methods tended to follow the 1:1 line, SfM had several outliers, in which the SfM snow depth overestimated the in situ measurements. Overall, the UAS lidar performance was modestly better in the open field than the UAS SfM, as compared with the magnaprobe measurements based on the mean absolute difference (MAD) (SfM = 4.0 cm, lidar = 3.5 cm) and the fitted linear regression line r² values (SfM = 0.51, lidar = 0.73). Samples from the field had better agreement between UAS and magnaprobe measurements than those from the forest. In the forest, the MAD values increased modestly for the lidar, but sharply for the SfM snow depths when compared with the magnaprobe measurements (SfM = 31.4 cm, lidar = 6.3 cm). In addition, lidar measurements had a much higher r² value than SfM in the forest (SfM = 0.02, lidar = 0.70). The UAS-measured snow depths in the forest are also shifted to the right of the 1:1 line, indicating that the UAS tended to measure snow depths as being shallower than the magnaprobe. This does not necessarily indicate an error in the UAS measurements because a previous study by Proulx et al. (2023) at this study site demonstrated that the magnaprobe tends to overprobe in the forest.

Figure 3 UAS-based snow depth measurements compared with in situ snow depth measurements from the magnaprobe. UAS-measured snow depths are shown for pixels overlapping the magnaprobe sample locations. Measurements collected on all dates are shown for (a) the field and (b) forest. n values indicate the number of samples shown within the plot axes. UAS-measured snow depths that exceeded 50 cm are not shown (1 SfM field, 6 SfM forest). Magnaprobe sample locations corresponding to UAS pixels with missing snow depth values are not shown (1 SfM field, 6 SfM forest, 7 lidar forest). UAV denotes unpiloted aerial vehicle.

4.2 Comparison between lidar and SfM snow depth

Figure 4 shows a direct comparison between the snow depths at individual 1 m × 1 m pixels measured by lidar and SfM over the entire study area. The field was divided into three areas based on an early study (Cho et al., 2021) showing distinct topographic and soil characteristics in each section. While there is considerable scatter for the individual pixels in all field areas, the SfM and lidar snow depths tend to agree fairly well based on MAD (east = 4.4 cm, west = 2.9 cm, northwest = 4.3 cm). Compared with the northwestern and eastern fields, the western field had the most similar snow depth values for SfM and lidar. In that field, both lidar and SfM techniques captured relatively deeper snow depths, ranging from 50 to 100 cm. In the northwestern and eastern fields, SfM snow depths are frequently much deeper than the lidar snow depths. In contrast, there was no clear agreement between SfM and lidar snow depths at the 1 m resolution in the forest (MAD = 55.7 cm), largely due to extensive regions in which SfM snow depths were anomalously high.

Figure 4 Comparison of SfM and lidar measurements for all sample dates by location. Scatterplots (left) compare snow depth for (a, c, e) the three field areas and (g) the forest. Probability density plots (right) show the distribution of snow depth values by UAS technique for (b, d, f) the three field areas and (h) the forest. Measured snow depths exceeding 50 cm are not shown.

Time series of lidar and SfM snow depth maps over the entire study area are shown in Fig. 5. While the previous section found that individual locations may have differences, these maps show that both techniques capture the differences in snow depth between flights in both the field and forest areas. The difference between the snow depth maps by sample date shows where the UAS snow depths tend to agree and disagree. The difference between SfM and lidar snow depths was fairly consistent in the field and close to 0 cm on most dates. However, on 20 and 24 February the southeastern field showed a negative difference, indicating that the SfM snow depths were considerably deeper than the lidar snow depths in parts of the field. On other dates, including 4 February, the SfM snow depth map was missing data from extended areas in the field. In the forest, missing or patchy SfM data occurred on many days (e.g., 4, 20, and 28 February). Lidar and SfM ground return point count statistics by land cover type are summarized in Tables S3 to S6. Despite the limited overlap in the forest, there is limited agreement between the two methods. Most maps show that the SfM snow depths were much deeper than the lidar snow depth through most of the forest, except at the forest and field edge.

Figure 5 Time series of snow depths for UAS lidar and SfM in the field and forest. Difference is calculated as lidar snow depth minus SfM snow depth. All values are shown in centimeters.

4.3 Spatial distribution of snowpack and its temporal changes

Understanding and quantifying the spatiotemporal variability – or stability – of snowpack is essential for identifying the physical drivers that influence snow accumulation and ablation across heterogeneous landscapes. To explore these dynamics in detail, MRD values were mapped to reveal spatial patterns in snow depth across survey dates (Fig. 6). The average spatial distribution of snow depth across the study domain based on UAS lidar-derived maps from eight survey dates shows spatially distinct patterns. The field had relatively deeper snow, by up to 70 % greater than the spatial mean. Within the field, the snow in northern areas was generally deeper than that in the southern areas, except for near the northern edges, where it was shallower. In forested areas, the snow was shallower by up to −70 % relative to the spatial mean. Snow in northeastern forest areas was generally shallower than that in other forest areas. Distinct differences in the transition zones between field and forest show edge effects. Immediately south of the tree line at the northeastern extent of the field, the snowpack is noticeably shallower. Moving away from the forest edge, there is a transition zone where snow becomes progressively deeper. In contrast, the southern portion of the northwestern field exhibits deep snow.

Figure 6 The snow depth mean relative difference (MRD) map generated by averaging the eight relative difference maps from the UAS lidar-based snow depth maps from 4 February to 7 March.

The relative difference snow depth maps for each date show that the spatial patterns of relative differences were fairly consistent throughout the study period (Fig. 7). Generally, there was deeper snow in the northern part of the field (about 50 % larger than the spatial mean) and shallower snow was found in forested areas as well as the central part of the field. Also, there was shallower snow along the northeastern boundaries of the field. These patterns were very clear during the accumulation period before the peak snow depth around 22 February. During the ablation period, consistent spatial patterns of the relative difference were still observed, particularly in the northern/eastern field, where snow remained relatively deep, in contrast to the forest areas, which continued to show shallow snow or exposed ground. These patterns persisted despite the increased presence of patchy snow cover in some regions. These gaps emerged primarily due to differential melting, leading to sections with no remaining snow.

Figure 7 Relative difference maps generated from the UAS lidar-based snow depth maps from 4 February to 7 March. The white areas in the figures indicate either masked areas (e.g., ponds and facilities) or areas with no snow.

4.4 Physical variables characterizing spatial distribution of snow depth

To evaluate the effect of physical land characteristics on the spatial distribution of snow depth, the MRD values were analyzed with respect to five land and soil characteristic values (namely, vegetation type, slope, shadow hours, K_sat, and SOM) over the study domain. Boxplots of MRD by physical feature are shown for the combined forest and field areas, field only, and forest only (Fig. 8). Statistical significance results among groups, based on Kruskal–Wallis and Tukey tests, are summarized in Tables S7 and S8. In the combined areas (i.e., forest + field), relative snow depth significantly differs by vegetation type. Coniferous forests have low MRDs (mean: −0.36), which indicate that snow in those areas is shallower relative to the spatial mean of snow depth by around 36 %. For the deciduous forest, the mean MRD is −0.2, with a wide interquartile range from −0.23 to 0.19. MRD values in the field are higher compared with the two forest types, ranging from −0.11 to 0.22 (mean: 0.08). For the combined areas as well as field and forest only, slope contributes to snowpack spatial patterns, even though the study area has a gentle slope (less than 20 %). High MRDs are found in flat areas (0 %–5 % slope) and gradually decrease with increasing slope. The effect of slope for the forest-only area is relatively modest. The shadow hours show a clear but contradictory contribution to snow depth patterns in the field- and forest-only areas as compared with the combined area. When the field and forest are separated out, low MRDs are found in areas where shadow hours are short (e.g., less than 2 h), and the MRDs gradually increase with increasing shadow hours. For the combined area, the highest shadow hours had the lowest snow depth, but this is likely to be the result of a mixed effect due to the dense shading in the coniferous forest. K_sat shows little evidence of contributing to the spatial distribution of snow depth in the field, but there are distinct differences in MRDs among lower K_sat groups in the forest. In the combined areas, MRDs tend to decrease with increasing K_sat values, except for the highest K_sat group. Compared with K_sat, SOM exhibits a clearer pattern, with MRD decreasing as SOM increases in both the combined areas and field analysis. The forest area does not display consistent MRD patterns with changes in SOM.

Figure 8 Boxplots of the snow depth MRD for each physical feature (vegetation type, slope, shadow hours, K_sat, and soil organic matter) for (a) the combined areas (forest and field), (b) field, and (c) forest only. The 1–5 for each boxplot except for vegetation type represents the relative range of each physical variable in each area (for example, for slope in the combined areas, 1: 0 %–5 %, 2: 5 %–10 %, 3: 10 %–15 %, 4: 15 %–20 %, and 5: 20 %–25 %).

5 Discussion

5.1 Comparison with previous findings: UAS SfM and lidar snow depths

The value of lidar datasets for capturing the horizontal and vertical structure of forests and snow cover depth is well established (Harder et al., 2020; Jacobs et al., 2021; Donager et al., 2021; Dharmadasa et al., 2022). However, the technology remains expensive and data processing is complex. The lower cost of SfM techniques compared with lidar make them a valuable tool for conducting surveys of snowpack change over time (Fernandes et al., 2018). Post-processing of RGB imagery is often less complex than lidar data processing and a variety of SfM software is now available, including some open-source options. However, our study concurs with early findings that SfM accuracy for measuring forest snow depths still cannot match that of lidar (Donager et al., 2021).

In SfM processing, an insufficient number of valid tie points, used to stitch together overlapping images, may degrade the accuracy of SfM snow depth data (Harder et al., 2016). Due to the reliance on RGB optical imagery, overcast skies and poor lighting over relatively homogeneous snowpacks (e.g., fresh snow) make it difficult for SfM post-processing software to identify a sufficient number of valid tie points (Bühler et al., 2016, 2017; Harder et al., 2020; Revuelto et al., 2021b; Miller et al., 2022). A lower number of valid tie points and higher point uncertainty result in large data gaps and poor estimation of surface elevations. In this study, when there was relatively fresh snow and few features, there were gaps in the SfM snow maps in the eastern field on 4, 20, and 24 February. Areas in the field with a sufficient number of unique tie points showed better agreement between SfM and lidar-measured ground surface elevations than those with fewer tie points (Fig. S2).

SfM performance is also lower in areas with dense vegetation where the ground surface is blocked by tree branches and canopies (Harder et al., 2020). Poor penetration of the forest canopy results in fewer overall ground returns and sparser point clouds. In Thompson Farm's mixed forest with dense underbrush, SfM post-processing also exhibited numerous tie point errors due to the lack of unique tie points. This is likely to be due to a plethora of thin tree branches and brush that appear as repetitive forest features. The features in these forests differ from forests with distinct trees and limited vegetation, such as the Harder et al. (2020) mixed lodgepole pine and subalpine fir forests. On all dates in the forest, the resampled 1 m × 1 m lidar data had fewer pixels with missing data than SfM (lidar = 0.4 %–1.6 %, SfM = 3.7 %–70.5 %). This error is apparent in the profile views of the ground surface returns for SfM and lidar from 10 February (Fig. S2) where SfM-measured ground elevations had greater variability compared with lidar in the forest. Methods for reducing the errors associated with SfM in the forest are limited. Adjustments to survey techniques, such as changing the camera angle and flying at lower altitudes and speeds, and post-processing workflows, such as making selection criteria less restrictive, may improve the number of points somewhat (Lendzioch et al., 2019). However, active remote sensing techniques, including lidar, have better penetration of the forest canopy than those that rely on passive sensing (Harder et al., 2020; Bühler et al., 2016).

SfM-derived error values from our study were 4.0 cm MAD and 6.8 cm root mean square deviation (RMSD) for the field, and 31 cm MAD and 71 cm RMSD for forested areas, highlighting a clear vegetation-dependent variation in accuracy. These findings are consistent with previous studies comparing UAS SfM and snow probe measurements, which report RMSD values typically below 31 cm in sparsely vegetated or alpine environments, increasing up to 37 cm in areas with denser vegetation, such as bushes, tall grass, or forests (De Michele et al., 2016; Bühler et al., 2016; Avanzi et al., 2018; Belmonte et al., 2021). Studies using UAS SfM alongside rulers and snow stakes had a root mean square error (RMSE) of less than 14 cm in both forested and prairie land cover types (Fernandes et al., 2018; Harder et al., 2016). Findings were also similar for studies comparing UAS lidar with rulers or snow stakes, which measured a RMSE less than 17 cm, with even lower RMSE values in shallow snow and in sunny areas (Harder et al., 2020; Feng et al., 2023; Koutantou et al., 2022). Lidar RMSE also tends to increase in vegetated areas, regardless of the vegetation class or type (Harder et al., 2020). Much like Harder et al.'s (2016) observations of erroneously high SfM snow depth measurements several meters above the snow surface, we observed SfM-measured snow depths greater than 150 cm in some forest locations. Conversely, our lidar-measured snow depths never exceeded 25 cm, indicating a more consistent performance in forested areas. We also found that the SfM snow depths did not consistently agree with the lidar snow depths over the entire field and on all dates. On most dates, the difference in UAS-measured snow depths was close to 0 cm in the field. The best agreement occurred in the western portion of the field, while the southeastern and northwestern portions had a larger amount of variability in measured values. The shadow hours and land cover type in the eastern and western fields are similar; however, the eastern field has a more gentle and less variable slope and fewer unique features (e.g., access road, USCRN station, pond, dirt piles, footprints) than the western field. The relatively homogeneous features in the eastern field indicate that the difference between techniques is likely to be due to a lack of sufficient valid tie points for SfM. It is not clear what caused the differences between SfM and lidar in the NW field that were not present in the other field areas. Unique features in the NW field are prevalent drainage patterns and shadowing that could be further investigated in the future.

While it is apparent that the accuracy of SfM-derived snow depth estimates cannot match that of lidar, the results of this study indicate that both techniques provide sufficient accuracy for monitoring the median change in shallow snow depths over time in flat unforested land covers when there are a sufficient number of unique characteristics for SfM post-processing. It is clear from the results of this study and previous ones that, compared with in situ data, UAS lidar techniques produce lower errors and fewer data gaps than SfM, especially in forested land cover and over homogeneous snowpacks (Bühler et al., 2016, 2017; Harder et al., 2020; Revuelto et al., 2021b; Miller et al., 2022). While UAS lidar may be the preferred technique in most cases, UAS SfM can still provide valuable information on changes in median snowpack depth across unforested areas at a relatively low cost and with less complex post-processing compared with UAS lidar. Regardless of the sampling technique used, the unique capability of UASs for measuring snowpack properties at the field scale and at a high temporal resolution makes them useful for observing snowpack evolution over time (Fernandes et al., 2018; Harder et al., 2020). Collection of in situ snow depth time series data is often time- and cost-prohibitive and may be especially challenging in complex or avalanche-prone terrain (Bühler et al., 2016; Harder et al., 2020). Monitoring snow depth changes at the field scale provided insights into accumulation and ablation patterns across the entire study area, as well as between different land cover types (e.g., forest and field).

5.2 Physical variables at the field scale

With limited wind redistribution in the study area, the time stability analysis indicated that relative differences in the snowpack were generally consistent throughout both the accumulation and ablation periods. In addition to previous findings that snowpack patterns are relatively consistent from year to year (Pflug and Lundquist, 2020; Revuelto et al., 2014), this study demonstrates that fixed physical variables, such as vegetation, topography, and soil characteristics, can sufficiently control the spatial variations of snowpack throughout a winter period. By comparing maps of mean relative difference (MRD) with maps of physical variables at the site, specific factors influencing snowpack dynamics over the winter season were identified. Our findings highlighted that vegetation type is a dominant factor shaping snow depth patterns. In both combined and field-only areas, SOM showed a statistically significant relationship, with snow depth decreasing as SOM increased. Furthermore, shadow hours and slope were found to contribute to the spatial variability of snowpack, even though the study area features relatively gentle slopes. The findings regarding the influence of vegetation and topographical factors on the snowpack's spatial variability align with previously conducted studies (Currier and Lundquist, 2018; Deems et al., 2006; Trujillo et al., 2007).

As compared with vegetation and terrain characteristics, few studies have examined the influence of soil characteristics on the snowpack. Our results indicate that snowpack depth decreases statistically significantly with increasing SOM (significance level <0.01). This finding aligns with our previous study, which utilized maximum entropy modeling to analyze spatial variations of shallow snowpack over the same domain but during different periods (Cho et al., 2021). Even though a clear relationship between K_sat and snowpack was not found in this study (Fig. 8b), it is acknowledged that soil thermal properties, such as the thermal conductivity of the soil underneath the snowpack, generally influence the rate of heat transfer between the snow and soil layers (Kane et al., 2001; Zhang, 2005). Also, the moisture content of the soil can affect the distribution of soil frost (Bay et al., 1952) and snowpack because the energy transfer at the snow–soil interface is controlled by the wetness of the soil (Bay et al., 1952; Fu et al., 2018). Although the spatial distribution of soil moisture is typically considered to be constant (frozen) during winter, intermediate rainfall events and freeze–thaw cycles can dramatically change the spatial patterns of soil moisture and freeze–thaw states in regions having ephemeral snowpacks. This can be important because the thermal conductivity in the frozen state is more sensitive to soil type than the non-frozen condition, because the thermal conductivity of ice is four times larger than that of the liquid phase (Penner, 1970). However, few studies have investigated how soil moisture patterns may control the spatial distribution of snowpack. This is likely to be because of the difficulty of measuring spatial distributions of soil moisture and freeze–thaw states beneath the snowpack. Even though this study did not focus on it, future investigations could measure spatial patterns of soil moisture and freeze–thaw states beneath the snowpack to quantify their interactions with the snowpack. A better understanding of the soil characteristics and their impact on the snowpack in various environments would help the snow community to accurately predict and model snow distribution and snowmelt processes.

Although there are numerous studies that characterize temporal changes in the spatial distribution of snowpack across topographically uniform landscapes, encompassing both open field and forested environments (Hannula et al., 2016; Clark et al., 2011; Currier and Lundquist, 2018; Mazzotti et al., 2023), the concept of “time stability” (or “temporal stability”) using relative difference values, as implemented in this study, has not been previously applied in snow hydrology. However, in the soil moisture community, numerous investigations that have examined temporal variability have been instrumental in developing robust validation sites and sampling strategies for satellite-based soil moisture assessments (Grayson and Western, 1998; Cosh et al., 2008; Brocca et al., 2009; Mohanty and Skaggs, 2001; Jacobs et al., 2004). Similar to its utility in soil moisture studies, the integration of the time stability concept into snowpack analysis at the field scale could facilitate the identification of representative sampling locations and inform the design of sampling protocols for optimal spatial extrapolation. Extending this approach to diverse snow environments will contribute to quantifying spatiotemporal variability in snowpack, thereby enhancing the establishment of core validation sites for potential snow missions such as the Canadian Terrestrial Snow Mass Mission (TSMM; Derksen et al., 2021).

5.3 Limitations and future perspectives

Given that our investigation was conducted in a relatively uniform landscape characterized by a shallow snowpack, it is imperative to extend the analysis to encompass diverse plant functional types, climatic zones, or snow classes (Johnston et al., 2024; Sturm and Liston, 2021) to ascertain the generalizability of the findings. This is essential as snow depth distributions are influenced by terrain attributes and snow regimes (Clark et al., 2011; Currier and Lundquist, 2018). In contrast to the present study area, where spatial heterogeneity in snowpack is predominantly influenced by static terrain characteristics and vegetation cover, alpine and prairie regions experience variability due to wind-driven processes (Elder et al., 1991). Further investigation incorporating additional analyses of energy fluxes and meteorological parameters, including solar radiation, soil temperature, and wind speed/direction, would enhance the comprehensiveness of the findings concerning the primary determinants of snowpack spatial variability across both static and dynamic variables.

Even though we analyzed the spatiotemporal variability of the snowpack using well-validated UAS-based snow depth observations, this may not guarantee that the current findings also capture the snow water equivalent (SWE) variations necessary for hydrological applications. Snow density, needed to calculate SWE from snow depth, is affected by snow metamorphosis differently than snow depth. Snow density may change during snowmelt as water percolates into the snowpack and refreezes. Also, vegetation and soil characteristics strongly control turbulent and ground heat fluxes and impact snow properties including snow density (Pomeroy and Brun, 2001). Studies have used physics-based snow models, such as SnowModel (Liston and Elder, 2006), Crocus (Vionnet et al., 2012), and Flexible Snow Model (FSM2; Essery et al., 2024), to understand those physical processes that cause snowpack spatial variations. However, observational approaches focusing on spatial distribution of SWE are quite limited because only now are sensing techniques emerging that directly observe the spatial distribution of snow density with a UAS (McGrath et al., 2022). A potential future direction is to develop reliable high-resolution SWE maps by integrating emerging techniques such as lidar and gamma-ray spectrometry (Harder et al., 2024), enabling the quantification of SWE spatial distribution across diverse snow environments. Another direction could involve employing an integrative approach using physical models to maximize in situ and UAS snow observations through data assimilation or novel interpolation methods, utilizing machine (or deep) learning approaches.

6 Conclusions

In this study, UAS lidar and SfM snow depth measurements were assessed using ground-based magnaprobe data and then used to confirm that spatial patterns of snowpack depth are temporally stable. Lidar demonstrated superior performance compared with SfM when evaluated against in situ observations, exhibiting lower errors. Both UAS techniques exhibited lower errors in field settings (lidar MAD = 3.5 cm, SfM MAD = 4.0 cm) than in forested environments (lidar MAD = 6.3 cm, SfM MAD = 31.4 cm). As expected, differences between lidar and SfM snow depths were more pronounced in forested regions (MAD = 55.7 cm), with SfM often registering anomalously deep snow depth values. The spatial distribution of snow depth captured by lidar remained consistent throughout the study period. For the entire study area, deeper snow was found in the field, in locations with shallow slopes and lower soil organic matter. Within the field, snow deepened with increasing shadow hours. When examining combined landscapes including forests and fields, we observed that the spatial distribution of snow depth was predominantly shaped by the type of vegetation present. Within the field, the spatial distribution of snow depth tracked relatively modest local slope variations and shadowing effects at the forest–field edge. As ephemeral snow conditions expand in a warming climate, these results are valuable for effectively comparing UAS and in situ sampling techniques for ephemeral, shallow, seasonal snowpacks. It is also expected that this study contributes to the enhancement of land surface and snow models by offering insights into parameterizing subgrid-scale snow depths, downscaling coarse-scale remotely sensed snow observations, and comprehending snowpack evolution at the field scale, particularly in ephemeral snow environments.