Structure from Motion of Multi-Angle RPAS Imagery Complements Larger-Scale Airborne Lidar Data for Cost-Effective Snow Monitoring in Mountain Forests

Snowmelt from mountain forests is critically important for water resources and hydropower generation. More than 75% of surface water supply originates as snowmelt in mountainous regions, such as the western U.S. Remote sensing has the potential to measure snowpack in these areas accurately. In this research, we combine light detection and ranging (lidar) from crewed aircraft (currently, the most reliable way of measuring snow depth in mountain forests) and structure from motion (SfM) remotely piloted aircraft systems (RPAS) for cost-effective multi-temporal monitoring of snowpack in mountain forests. In sparsely forested areas, both technologies give similar snow depth maps, with a comparable agreement with ground-based snow depth observations (RMSE ~10 cm). In densely forested areas, airborne lidar is better able to represent snow depth than RPAS-SfM (RMSE ~10 cm vs ~10–20 cm). In addition, we find the relationship between RPAS-SfM and previous lidar snow depth data can be used to estimate snow depth conditions outside of relatively small RPAS-SfM monitoring plots, with RMSE’s between these observed and estimated snow depths on the order of 10–15 cm for the larger lidar coverages. This suggests that when a single airborne lidar snow survey exists, RPAS-SfM may provide useful multi-temporal snow monitoring that can estimate basin-scale snowpack, at a much lower cost than multiple airborne lidar surveys. Doing so requires a pre-existing mid-winter or peak-snowpack airborne lidar snow survey, and subsequent well-designed paired SfM and field snow surveys that accurately capture substantial snow depth variability.


Introduction
Snowpack in mountain forests is a major source of water for reservoirs that provide water and hydropower for many urban and agricultural communities [1][2][3]. Mountain snowpacks are affected by many climatic, topographic and ecological variables, and are sensitive to forest disturbance such as thinning, prescribed fires, wildfire, and tree die-off [4][5][6][7][8][9][10][11][12][13]. It is important to monitor how snowpacks in these areas respond to changing environmental conditions in order to understand and forecast available water resources for both natural and human consumption.
Characterization of snowpack in these mountain forests is challenging because of the large amount of small-scale spatial variability of the snowpack, due to topographic heterogeneity and variable forest structure [14][15][16][17][18]. This extreme heterogeneity makes it difficult to monitor snowpack with point snow measurements [19]. Remote sensing of snowpack is better able to characterize this heterogeneity, though it suffers in forested environments because trees can interfere with the remote sensing of snowpack below the canopy [20]. Furthermore, remote sensing of snowpack is hindered by clouds representing the diverse forest and snowpack conditions in the region. All surveys have multi-angle aerial photography acquisitions from a RPAS paired with precisely geolocated ground-based snow surveys, while seven of these surveys are coincident with acquisitions of high point density airborne lidar. The combination of the RPAS surveys with the airborne lidar and ground surveys provides an unprecedented opportunity not only to evaluate the performance of the SfM remote sensing under a variety of forest and snowpack conditions, but also to explore opportunities to combine all of these measurements for cost-effective, high quality multi-temporal snowpack monitoring in mountain forests.

Field Data
Data in this study come from four intensively studied snow-research plots near the Mogollon Rim, a mid-to high-elevation (~2000-3000 m) forested region in central Arizona (AZ; Figure 1). Two "montane" plots, located~50 km east of Show Low, AZ, are in close proximity to each other at a high-elevation (~2800 m) field site. These plots have low to moderate relief slopes and contain a mixture of burned (from a 2011 wildfire [51]) and live mixed conifer forests and montane grasslands. One plot (the "montane meadow" plot) contains large snow variability, due to heterogeneous forest conditions (i.e., tree stands separated by meadows). In contrast, the other (the "montane valley" plot) contains large snow variability due mostly to heterogeneous topography (i.e., north-vs. south-facing slopes). A "dense forest" plot, located~20 km north of Payson, AZ, has gentle slopes and is densely forested, primarily with Ponderosa Pine (pinus ponderosa) trees. Finally,~10 km to the north of the dense forest plot, is another plot (called the "thinning comparison" plot) that is located along with the transition between recently thinned and unthinned ponderosa pine forest patches near Clints Well, AZ. The dense forest and thinning comparison plots have shallower snowpack than the montane plots. These four plots span a range of forest densities, tree heights, and topographic conditions (Table 1), allowing us to evaluate the performance of snow measurement with SfM and airborne lidar under a variety of forest conditions.  Figure 1. Locations of the four study plots (a -overview, b -montane meadow and montane valley plots, c -dense forest plot, d -thinning comparison plot). Canopy cover data are derived from 2019 airborne lidar data [52].At these field sites, we collected rich snowpack datasets, including a variety of airborne and ground-measured snowpack data that are ideal for creating and evaluating highquality SfM models. Snow surveys for these plots consisted of multi-angle RPAS imagery, airborne lidar, and ground measurements of snow depth and SWE. In total, there were 11 snow surveys (Table  1), with the most surveys occurring at the montane plots (as this area is important for snowpack monitoring for local water supply managers), with fewer surveys occurring at the dense forest and thinning comparison plots. All surveys included RPAS flights to acquire aerial imagery (using a DJI Phantom 3 Pro drone with a 12 megapixel camera-for 2017 surveys-and a DJI Phantom 4 Pro drone with a 20 megapixel camera-for 2018-2020 surveys) to generate snow-on 3-D models of the study plots using SfM (see Section 3). All RPAS flights occurred on clear, sunny days between 10:00 AM and 2:00 PM to minimize the length of shadows. For the flights in 2017, the RPAS was controlled manually with a combination of nadir and oblique (~30% nadir) imagery along back and forth flight lines spanning each domain. For the flights in 2018-2020, flight planning software called Altizure ® was used to collect the imagery in a regular flight pattern (using back and forth flight lines) with a set of nadir imagery, as well as two sets of oblique images (with 75% overlap for all imagery). All flights occurred at ~75 m above ground level with horizontal velocities of ~5-7 m/s resulting in ~4 cm and Figure 1. Locations of the four study plots (a-overview, b-montane meadow and montane valley plots, c-dense forest plot, d-thinning comparison plot). Canopy cover data are derived from 2019 airborne lidar data [52]. At these field sites, we collected rich snowpack datasets, including a variety of airborne and ground-measured snowpack data that are ideal for creating and evaluating high-quality SfM models. Snow surveys for these plots consisted of multi-angle RPAS imagery, airborne lidar, and ground measurements of snow depth and SWE. In total, there were 11 snow surveys (Table 1), with the most surveys occurring at the montane plots (as this area is important for snowpack monitoring for local water supply managers), with fewer surveys occurring at the dense forest and thinning comparison plots. All surveys included RPAS flights to acquire aerial imagery (using a DJI Phantom 3 Pro drone with a 12 megapixel camera-for 2017 surveys-and a DJI Phantom 4 Pro drone with a 20 megapixel camera-for 2018-2020 surveys) to generate snow-on 3-D models of the study plots using SfM (see Section 3). All RPAS flights occurred on clear, sunny days between 10:00 AM and 2:00 PM to minimize the length of shadows. For the flights in 2017, the RPAS was controlled manually with a combination of nadir and oblique (~30% nadir) imagery along back and forth flight lines spanning each domain. For the flights in 2018-2020, flight planning software called Altizure ® was used to collect the imagery in a regular flight pattern (using back and forth flight lines) with a set of nadir imagery, as well as two sets of oblique images (with 75% overlap for all imagery). All flights occurred at~75 m above ground level with horizontal velocities of~5-7 m/s resulting in~4 cm and~2.1 cm spatial resolution on the ground for the Phantom 3 and Phantom 4 systems, respectively. The cameras used automatic exposure settings. These surveys also included ground measurements of snow depth and snow density. A majority of these measurements were collected by more than 25 employees and volunteers from five public and private partner entities (see Acknowledgements) along carefully pre-marked snow survey transects with precisely geolocated (using a differential GPS) endpoints to enable the precise geolocation of snow measurements in order to relate them to the airborne data. Following the sampling methodology described in Reference [53], five snow depth samples were taken every 5 m along the transects in a star pattern to ensure spatial representativeness of the snow depth samples, and snow density was sampled every 10-20 m along the transects. There were also some additional measurements along transects (also positioned with a differential GPS) where snow depth is monitored using remote cameras for temporal snowpack monitoring (however, this study only uses the snow depth data, as well as manual measurements of snow density along these transects for the survey dates in Table 1). Not all dates had identical ground surveys for each study plot, depending on logistical considerations (e.g., available workers). Overall, individual surveys included ground measurements from 22-108 individual locations ( Table 1). Note that this translates to >3000 measurements as there were five snow depth measurements at each manual sample location. Further details about the snow survey methodology can be found in Reference [54].

Generation of Lidar Snow Depth, Density, and SWE Maps
For seven of the surveys, we generated state-of-the-art airborne lidar-based maps of snow depth, snow density, and SWE, as they included airborne lidar data. These maps were created for the twõ 100 km 2 lidar footprint areas (which each have three separate lidar acquisitions), that encompass the study plots. The two montane plots are within a "high-elevation" lidar footprint, and the dense forest and thinning comparison plots are within a "mid-elevation" lidar footprint. These maps were created using high point density (~10-15 points/m 2 ) snow-on airborne lidar data (which are differenced from similar snow-off lidar data) and are bias-corrected with the field-measured snow depths. All lidar data was flown by Quantum Spatial Inc. [52,[55][56][57][58].
The snow depth maps were then combined with snow density maps generated using artificial neural network (ANN) machine learning of the field-measured snow density measurements, using various lidar-derived physiographic attributes as predictors. The methodology for creating these maps (and their validation) is detailed in Reference [54], but is briefly summarized here. We use an ensemble of ANNs with relatively simple network structure (1 hidden layer with ten neurons) with Levenberg-Marquardt (L-M) backpropagation (Marquardt, 1963) to model the field snow density measurements using lidar-derived metrics as predictor variables. The predictor variables include physiographic variables (elevation, slope, northness-or sin(slope) × sin(aspect), canopy height, and canopy closure), other GIS-derived quantities using the lidar data (skyview factor, below canopy solar forcing index-or the ratio of incident solar radiation reaching the ground surface over a given period to that hitting a flat surface with no obstructions), and lidar snow depth (from 1 February 2017). These snow density maps are multiplied by the lidar snow depth maps to generate maps of SWE.

Generation of SfM Snow Depth Maps
Snow-on point clouds were generated from the RPAS imagery using Agisoft Metashape software. After the point clouds were created, the ground points were isolated by using a combination of Metashape's built-in ground filtering algorithm (for a first cut separation of the point cloud) and a Cloth Simulation Filter algorithm [59] implemented in R (to remove remaining debris from the ground point cloud).
Next, we used an automated post-processing workflow to accurately georeference the point clouds using available snow-off point cloud data. In this study, we use snow-off lidar point cloud data (from summer 2014) for this purpose (though we also tried, with similar results, applying the workflow using snow-off SfM data that we collected-see Section 3.1 below). This workflow consists of two general steps. First, the snow-on SfM point clouds were co-registered with the snow-off point cloud canopy layer by removing any vertical offsets between the snow-on and snow-off point clouds and then using an Iterative Closest Point algorithm implemented in CloudCompare software to match the canopy layers. This successfully adjusted the snow-on point clouds to match the canopy elements in the snow-off point clouds, even when there was sometimes substantial vegetation change between the snow-on and snow-off data ( Figure 2). Next, warping and tilting of the point clouds were removed using low-order polynomial filters (a 1st order polynomial correction to correct tilt followed by a 2nd order polynomial correction to correct warping) by comparing them to the lidar ground surface plus a first guess of snow depth, based on the ground survey data. This first guess was generated using the same ANN methodology as the creation of snow density maps from field snow density measurements in Reference [54] and described in Section 2.2, except that field measured snow depths, rather than field measured snow densities, are used as the predicted variable.
Remote Sens. 2020, 12, 2311 6 of 24 workflow using snow-off SfM data that we collected-see Section 3.1 below). This workflow consists of two general steps. First, the snow-on SfM point clouds were co-registered with the snow-off point cloud canopy layer by removing any vertical offsets between the snow-on and snow-off point clouds and then using an Iterative Closest Point algorithm implemented in CloudCompare software to match the canopy layers. This successfully adjusted the snow-on point clouds to match the canopy elements in the snow-off point clouds, even when there was sometimes substantial vegetation change between the snow-on and snow-off data ( Figure 2). Next, warping and tilting of the point clouds were removed using low-order polynomial filters (a 1 st order polynomial correction to correct tilt followed by a 2 nd order polynomial correction to correct warping) by comparing them to the lidar ground surface plus a first guess of snow depth, based on the ground survey data. This first guess was generated using the same ANN methodology as the creation of snow density maps from field snow density measurements in Reference [54] and described in Section 2.2, except that field measured snow depths, rather than field measured snow densities, are used as the predicted variable. Example showing the snow-on SfM point cloud in the center of the montane meadow plot (mostly white, green, and brown colored dots represent snow, evergreen canopy, and woody material) in relation to the lidar point cloud (red dots) for an area in the montane meadow plot. Note that in some areas, there is a difference between the lidar (from summer 2014) and SfM (from winter 2017) point clouds, due to post-wildfire effects (lost branches, fallen snags and trees), but for the snags and trees that did not change during this period, the SFM and lidar point clouds are very close to one another.
This processing workflow was automated using CloudCompare's command line mode (for the ICP algorithm), as well as the US Forest Service's FUSION/LDV [60] software and a python script to model and correct warping and tilting in the snow-on SfM point cloud. The ANN machine learning methodology is implemented in Matlab ® software, as described in Reference [54].

Using SfM Snowpack Data to Supplement Lidar Snowpack Data
We used these SfM snow depth maps to advance our understanding of snow distributions for our research sites in a number of ways. First, following the methodology in Reference [54], we combined them with maps of snow density generated using ANN modeling of snow density measurements to generate maps of SWE for our study plots. This enabled us to characterize differences in terms of snow depth, snow density, and SWE for different seasons. This is important because while the directly measured variable using lidar and SfM is snow depth, SWE is the snowpack variable that is of most interest in hydrological applications [61]. Although the distribution Example showing the snow-on SfM point cloud in the center of the montane meadow plot (mostly white, green, and brown colored dots represent snow, evergreen canopy, and woody material) in relation to the lidar point cloud (red dots) for an area in the montane meadow plot. Note that in some areas, there is a difference between the lidar (from summer 2014) and SfM (from winter 2017) point clouds, due to post-wildfire effects (lost branches, fallen snags and trees), but for the snags and trees that did not change during this period, the SFM and lidar point clouds are very close to one another. This processing workflow was automated using CloudCompare's command line mode (for the ICP algorithm), as well as the US Forest Service's FUSION/LDV [60] software and a python script to model and correct warping and tilting in the snow-on SfM point cloud. The ANN machine learning methodology is implemented in Matlab ® software, as described in Reference [54].

Using SfM Snowpack Data to Supplement Lidar Snowpack Data
We used these SfM snow depth maps to advance our understanding of snow distributions for our research sites in a number of ways. First, following the methodology in Reference [54], we combined them with maps of snow density generated using ANN modeling of snow density measurements to generate maps of SWE for our study plots. This enabled us to characterize differences in terms of snow depth, snow density, and SWE for different seasons. This is important because while the directly measured variable using lidar and SfM is snow depth, SWE is the snowpack variable that is of most interest in hydrological applications [61]. Although the distribution of SWE is broadly similar to that of snow depth, it is important to consider the effects of variable snow densities, which can vary substantially in both space and time [54,62].
Next, we investigated the relationships between snow depths for the larger mid-and high-elevation lidar footprints on different dates to see how well they could be used to extrapolate snow depths measured for the SfM study plots to the larger lidar footprints. In particular, we compared the relationships between snow depths from the mid-winter (on 1 February 2017) and late winter (on 7 March 2017 and 4 March 2019) lidar snow surveys with those from the mid-winter lidar snow survey and the late-winter SfM data. To ensure that only the best quality SfM data was used to construct these relationships, we only considered SfM data for open areas without an overhead canopy for these comparisons. Of particular interest was (1) the strength of the relationships between the lidar data on multiple dates (which indicated how well they could be used to predict snow depths for different times), and (2) how well the SfM-data over the small study plots could capture those relationships. Even though the late winter surveys occurred at the same time of year, the 7 March 2017 survey reflected more evolved snowpack conditions: For the mid-elevation lidar coverage, it reflected mid-ablation to nearly snow-free conditions (depending on elevation), and for the high-elevation lidar coverage, it reflected peak-SWE to mid-ablation conditions, depending on elevation. On the other hand, the 4 March 2019 survey reflected conditions just prior to Peak SWE for the high-elevation lidar coverage, and mid-ablation conditions for the mid-elevation lidar coverage. The 1 February survey reflected mid-winter snowpack conditions in both areas.
Finally, we tested how well snow depths could be estimated from the SfM data for areas outside the SfM plots (but within the lidar coverages) on these different dates using the relationships between the late-winter SfM data and the mid-winter (1 February 2017) lidar data. These relationships were modeled using a third order polynomial fit (which was chosen based on the data analysis described above), and the regression parameters were applied to the 1 February 2017 lidar data to generate the predictions. These simulated maps were then compared to the actual lidar observations on 7 March 2017 and 4 March 2019.

Comparison Between SfM, Field-Based, and Lidar Snow Depths
Overall, there is good agreement between the lidar and the SfM snow depth maps in the montane plots, which generally have sparse canopies. Figure 3a-i shows that most major snow depth features are consistently represented in both maps, such as the areas with deep snow (which usually correspond to north-facing aspects), and the extremely shallow snow on south-facing sides of large tree stands. Overall, there are high squared correlation coefficient (R 2 ) values (0.78 to 0.89), relatively low Root Mean Squared Error (RMSE) values (8.5 to 9.4 cm) and small average differences (SfM average depths are 1.0 to 3.5 cm shallower) between the SfM and lidar data ( Table 2). The coefficient of variation (CV) differences between the maps are also relatively low (0.01-0.04). At the same time, the agreement between the SfM and field measured snow depths is comparable to that between the lidar and field measured snow depths for the montane plots (Table 3, Figure S1). Although R 2 values between the SfM data and ground measurements (0.60-0.90) are a little lower than those for the lidar data (0.68-0.94), SfM RMSEs (8.6-13.4 cm) are similar to lidar RMSEs (8.5-12.4 cm). Both the SfM and lidar data are relatively unbiased with respect to ground observations for the montane plots.
There is generally a lower correspondence between the lidar and SfM snow depth maps for the dense forest plots. Major features can clearly be seen in both sets of maps (e.g., the road on the northern end of the plot as well as the snowbanks on either side of the road), though other details in the forest are less consistent (Figure 3j-r). For example, the SfM maps show more pockets of shallow snow (especially on 1 February 2017). Compared to the montane plots, the dense forest plot has lower R 2 values (0.23 to 0.57), higher RMSEs (12.4 to 24.0 cm), and higher CV differences (0.09 to 0.26; with the SfM data having higher CVs than the lidar data; Table 2). At the same time, though, average depth differences are only a little larger than those observed in the montane plots (3.8 to 6.3 cm). Compared to the ground observations, the lidar generally performs better than the SfM data, though the performance of the SfM data is also highly variable (Table 3, Figure S1) for the dense forest plots. RMSEs between SfM and field measured snow depths range from 8.8 and 17.2 cm (vs. 8.3 to 10.3 cm for the lidar data), and R 2 values are 0.26 to 0.29 (vs 0.33 to 0.62 for the lidar data). However, the SfM data has lower biases, compared to the ground observations, than the lidar data at the dense forest plot (−0.4 to 2.1 cm vs 0.8 to 5.8 cm). performance of the SfM data is also highly variable (Table 3, Figure S1) for the dense forest plots. RMSEs between SfM and field measured snow depths range from 8. 8  For the thinning comparison plot (Figure 3s-u), there is a higher correspondence between the lidar and SfM maps than for the dense forest plot. The thinning comparison plot also has very little snow (~10 cm on average) during the March 2019 survey (the only snow survey performed at this site), yet both technologies show similar snow patterns. Overall, RMSE (7.4 cm) between the lidar and SfM data is lower than for the other surveys, R 2 (0.41) is higher than some of the dense forest plots surveys but lower than those in the montane plots (with the opposite true for CV differences), and the average snow depth difference (1.2 cm) is small compared to the other surveys. Compared to the ground survey data, the performance of the lidar and SfM data is similar (with R 2 values of 0.37 and 0.46 and RMSEs of 7.1-7.9 cm), again with a relatively small bias (1.9 cm). Overall, the thinning comparison plot has lower vegetation density than the dense forest plot, but higher vegetation density than the montane plots (Table 1).     For the thinning comparison plot (Figure 3s-u), there is a higher correspondence between the lidar and SfM maps than for the dense forest plot. The thinning comparison plot also has very little snow (~10 cm on average) during the March 2019 survey (the only snow survey performed at this site), yet both technologies show similar snow patterns. Overall, RMSE (7.4 cm) between the lidar and SfM data is lower than for the other surveys, R 2 (0.41) is higher than some of the dense forest plots surveys but lower than those in the montane plots (with the opposite true for CV differences), and the average snow depth difference (1.2 cm) is small compared to the other surveys. Compared to the ground survey data, the performance of the lidar and SfM data is similar (with R 2 values of 0.37 and 0.46 and RMSEs of 7.1-7.9 cm), again with a relatively small bias (1.9 cm). Overall, the thinning comparison plot has lower vegetation density than the dense forest plot, but higher vegetation density than the montane plots ( Table 1).
The accuracy of these SfM data depends, in part, on the quality of the first guess snow depth maps (which are shown in Figure S2 Supplemental Materials). In general, these first guess maps are comparable, in magnitude to the SfM and lidar maps, shown in Figure 3, but they typically do not reflect as much variability as the SfM and lidar maps do. The overall agreement between theses maps is usually lower than that between the SfM and lidar maps (with the exception of the 1 February 2017 dense forest plot, R 2 values range from 0.09-0.33 lower, and RMSE values range from 0.4-3.6 cm higher). Note that for the forest thinning comparison plot, the snow depths were so shallow that a spatially variable first guess map was not needed. The first guess snow depth map showed virtually no spatial variability, due to the small snow depth signal compared to the uncertainty of these maps ( Figure S2g). Nevertheless, it was included here to maintain methodological consistency with all other surveys.
Finally, note that the same methodology used here can be used to generate SfM snow depth maps that based on snow-off SfM point clouds instead of snow-off lidar point clouds ( Figure S3). The agreement between the SfM and lidar snow depth maps is a little lower for these maps (compare Table 2 and Table S1). This discrepancy may be related to inconsistent ground filtering algorithms used for the lidar data (which was performed by QSI) and SfM data (see Section 2.3).

SWE Monitoring Using SfM
Due to the hydrological importance of quantifying water content of the snowpack, the snow density and SWE maps for each plot prepared using the methodology of Broxton et al. [54] (see Section 2.3) were useful for multi-temporal monitoring of snowpack. In this study, they were especially useful for the montane plots because this area is very important for predicting local water supplies. For this area, there were snow surveys for four subsequent water years near the time of peak SWE (in early March). Figure 4 shows spatial differences between snow depth, snow density, and SWE for this area from 2017-2020 (note that in March 2017, there was airborne lidar, but no SfM data in the area).
While the 2019 maps depict snowpack conditions just prior to peak SWE, the 2017 and 2020 maps reflect snowpack conditions that are further along in the melt season (note the relatively large areas with no snow on south-facing slopes), and the 2018 maps depict historically low snow conditions (as 2018 was a very dry and warm year). Despite the differences between snowpack conditions between 2018 and the other years, the 2018 data shows relatively deeper snowpack in the same areas as in the other sets of maps (e.g., on north-facing slopes in the montane valley plot and on the shaded sides of the montane meadow plot). There is a relatively high correlation between the 2017, 2019, and 2020 surveys (R 2 between the 2017 and 2019 survey data are 0.71 and 0.69 for depth and SWE, respectively, and the R 2 between the 2017 and 2020 surveys are 0.67 and 0.60 for depth and SWE, respectively). The R 2 between the 2017 and 2018 survey is lower, (R 2 = 0.36 for snow depth and 0.39 for SWE). Overall, the 2019 survey had the most snow, but snow variability (as measured by the CV of snow depth and SWE) was higher during 2017 (when there was more advanced snowmelt) and 2018 (when the snowpack was shallow; Table 4) surveys. Snow density was also lowest in 2018 (due to shallow snowpack) and highest in 2020 (as the survey occurred after a long period of warm, dry weather).
Remote Sens. 2020, 12, 2311 11 of 24 the snowpack was shallow; Table 4) surveys. Snow density was also lowest in 2018 (due to shallow snowpack) and highest in 2020 (as the survey occurred after a long period of warm, dry weather).   Figure 4 shows that while there are differences between the snowpack distributions for different years (e.g., between the peak-SWE conditions on 4 March 2019 vs the post-peak SWE conditions on 6 March 2017), there are also similarities (e.g., areas with deeper snow in one scene tend to have deeper snow in other scenes). To further understand the consistency between the snowpack at different times across much larger (~100 km 2 ) lidar coverages, we plotted the relationship between snow depth for three sets of lidar coverages (at different times in the snow season) for the mid-and high-elevation lidar domains ( Figure 5). The mid-winter (1 February 2017) lidar data is used as the abscissa, and the late-winter (7 March 2017 and 4 March 2019) lidar data are used as the ordinates in these plots. Figure 5 (red dots) shows that the relationships between the snow depth data change between different dates. For example, while the relationship between 1 February 2017 and 4 March    Figure 4 shows that while there are differences between the snowpack distributions for different years (e.g., between the peak-SWE conditions on 4 March 2019 vs the post-peak SWE conditions on 6 March 2017), there are also similarities (e.g., areas with deeper snow in one scene tend to have deeper snow in other scenes). To further understand the consistency between the snowpack at different times across much larger (~100 km 2 ) lidar coverages, we plotted the relationship between snow depth for three sets of lidar coverages (at different times in the snow season) for the mid-and high-elevation lidar domains ( Figure 5). The mid-winter (1 February 2017) lidar data is used as the abscissa, and the late-winter (7 March 2017 and 4 March 2019) lidar data are used as the ordinates in these plots. Figure 5 (red dots) shows that the relationships between the snow depth data change between different dates. For example, while the relationship between 1 February 2017 and 4 March 2019 lidar data for the high-elevation lidar coverage is linear, the relationship between the 1 February 2017 and 7 March 2017 lidar data for the high-elevation lidar box is nonlinear. The former represents a comparison between mid-winter and conditions just prior to peak SWE while the latter represents a comparison between mid-winter and primarily post-peak SWE conditions. Despite this changing shape, the strength of these relationships is consistent in both cases (with R 2 values of~0.9). For the mid-elevation lidar coverage, these relationships (which both represent a comparison of mid-winter conditions and mostly mid-ablation conditions) are nonlinear for both pairs of data, and are somewhat weaker, but still consistent (R 2 =~0.77 in both cases). This suggests there a high degree of intra-seasonal and inter-annual relatability between snow depth patterns for our study sites at different times. Note that these relationships are similar if all lidar pixels are considered (thick red dashed lines in Figure 5) vs. if only open pixels (where SfM performs well) are considered (thin red dashed lines in Figure 5). mid-elevation lidar coverage, these relationships (which both represent a comparison of mid-winter conditions and mostly mid-ablation conditions) are nonlinear for both pairs of data, and are somewhat weaker, but still consistent (R 2 = ~0.77 in both cases). This suggests there a high degree of intra-seasonal and inter-annual relatability between snow depth patterns for our study sites at different times. Note that these relationships are similar if all lidar pixels are considered (thick red dashed lines in Figure 5) vs. if only open pixels (where SfM performs well) are considered (thin red dashed lines in Figure 5).

Using Airborne Lidar to Extend the SfM Data
Next, we tested how well the SfM data within the small study plots can capture these relationships. Also plotted in Figure 5 are the relationships between the 1 February 2017 lidar snow depth and the 7 March 2017 and 4 March 2019 SfM data for the dense vegetation plot (for the midelevation panels) and for the montane plots (for the high-elevation panels), considering only open pixels (where the SfM data perform well). Note that there was no SfM survey for the montane plots on 7 March 2017. The data follow similar trends to those found for the larger lidar coverages (just a little flatter; Figure 5), though with lower R 2 values. This may be due to the lower quality of the SfM data in some environments (e.g., the dense forest plot), as well as the smaller dynamic range represented in these small plots compared with the overall lidar coverages (note that the dense forest plot, in particular, is flat with fairly uniform forest cover across most of the plot).  Next, we tested how well the SfM data within the small study plots can capture these relationships. Also plotted in Figure 5 are the relationships between the 1 February 2017 lidar snow depth and the 7 March 2017 and 4 March 2019 SfM data for the dense vegetation plot (for the mid-elevation panels) and for the montane plots (for the high-elevation panels), considering only open pixels (where the SfM data perform well). Note that there was no SfM survey for the montane plots on 7 March 2017. The data follow similar trends to those found for the larger lidar coverages (just a little flatter; Figure 5), though with lower R 2 values. This may be due to the lower quality of the SfM data in some environments (e.g., the dense forest plot), as well as the smaller dynamic range represented in these small plots compared with the overall lidar coverages (note that the dense forest plot, in particular, is flat with fairly uniform forest cover across most of the plot). simulated snow depth data for the mid-and high-elevation lidar coverages on 4 March 2019, illustrate that these extrapolations are fairly successful at estimating the observed lidar snow depth data. This can be seen for both the entire ~100 km 2 lidar coverages, as well as for the < 1 km 2 study plots. Note, however, that there can be significant differences between the observed and simulated snow depth patterns in areas that have undergone substantial vegetation changes (such as the thinning comparison plot in Figure 6g-i, which underwent mechanical thinning between the 2017 and 2019 snow surveys).  The strength of the relationships between the lidar snow depth data on different dates, as well as the ability of the SfM data to capture these trends suggests that the SfM data can be used to predict snow depth distributions over the larger airborne lidar coverages, given that a lidar snow depth map already exists. To evaluate how well this extrapolation works, we used the relationships between the 1 February 2017 lidar, 7 March 2017 and 4 March 2019 SfM data, shown in Figure 5, applied to the 1 February 2017 lidar map, as described in Section 2.4. Figures 6 and 7, which show observed and simulated snow depth data for the mid-and high-elevation lidar coverages on 4 March 2019, illustrate that these extrapolations are fairly successful at estimating the observed lidar snow depth data. This can be seen for both the entire~100 km 2 lidar coverages, as well as for the < 1 km 2 study plots. Note, however, that there can be significant differences between the observed and simulated snow depth patterns in areas that have undergone substantial vegetation changes (such as the thinning comparison plot in Figure 6g  In general, spatial statistics from these simulated snow depth maps are also fairly similar to those from the observed lidar data. Table 5 shows statistics related to the distribution of snow depths (mean, cv, skewness) as well as measures of autocorrelation (measured by the fractal dimension of snow depths for distances that are smaller or larger than a well-known scale break that occurs at ~25-30 m). Note that there was no SfM survey for the high-elevation coverage on 7 March 2017, so statistics for this date are given assuming a simulated map that would be obtained if the relationships between the lidar data (red lines in Figure 5) were captured. For the high-elevation coverage, the distribution parameters for the simulated maps are very close to those for the observed maps: The mean snow depth is within 5%, the CVs are ~10% lower, and the skewness is similar. Likewise, both the short-and long-range fractal dimensions are similar, with the latter above 2.9, suggesting spatially random variability for longer correlation distances, and the former ~2.7-2.8, indicating more spatial complexity at shorter correlation distances. The agreement between these statistics for the simulated and observed lidar maps is a little worse for the mid-elevation coverage, as the simulated maps had slightly lower (by ~10%) mean snow depths, lower (by ~25%) CV, but similar skewness. Compared with the high-elevation coverage, there were also slightly larger differences between the short-range In general, spatial statistics from these simulated snow depth maps are also fairly similar to those from the observed lidar data. Table 5 shows statistics related to the distribution of snow depths (mean, cv, skewness) as well as measures of autocorrelation (measured by the fractal dimension of snow depths for distances that are smaller or larger than a well-known scale break that occurs at~25-30 m). Note that there was no SfM survey for the high-elevation coverage on 7 March 2017, so statistics for this date are given assuming a simulated map that would be obtained if the relationships between the lidar data (red lines in Figure 5) were captured. For the high-elevation coverage, the distribution parameters for the simulated maps are very close to those for the observed maps: The mean snow depth is within 5%, the CVs are~10% lower, and the skewness is similar. Likewise, both the short-and long-range fractal dimensions are similar, with the latter above 2.9, suggesting spatially random variability for longer correlation distances, and the former~2.7-2.8, indicating more spatial complexity at shorter correlation distances. The agreement between these statistics for the simulated and observed lidar maps is a little worse for the mid-elevation coverage, as the simulated maps had slightly lower (by~10%) mean snow depths, lower (by~25%) CV, but similar skewness. Compared with the high-elevation coverage, there were also slightly larger differences between the short-range fractal dimensions for the observed and simulated maps, suggesting that the simulated maps for the mid-elevation coverage were not as good at capturing small scale snow depth variability as for the high-elevation coverage.
Because the fine scale (1 m) spatial variability depicted in these simulated maps can be somewhat different than in the observed maps, we compared the observed and simulated maps at a variety of spatial scales, by averaging the 1 m pixels to larger pixel sizes ( Figure 8). As expected, the simulated and observed maps have a closer fit for larger pixel sizes, with the largest increases occurring between 1 and 100 m, and smaller differences occurring after that. For the mid-elevation lidar coverage, the R 2 values increase from~0.7-0.8 for 1 m pixels to > 0.9 for 100 m pixels, while RMSEs decrease from 10-15 cm for 1 m pixels to~7-8 cm for 100 m pixels. For the high-elevation lidar coverage, the R 2 values increase from~0.9 for 1 m pixels to > 0.95 for 100 m pixels, while RMSEs again decrease from 10-12 cm for 1 m pixels to~6-7 cm for 100 m pixels. Figure 8 also shows the performance of simulated maps that would be obtained if the relationships between the lidar data (red lines in Figure 5) were captured perfectly. For these scenarios, there is some reduction in RMSE (by~20-30%), suggesting that some of the error between the simulated (using SfM data) and observed snow depth data is the result of this relationship not being fully captured by the SfM data. Similar results can be obtained for predictions of SWE values (compare Figures S5-S8 with Figures 5-8).

Discussion
In this study, we show that RPAS-SfM monitoring of snowpack can be a useful tool to provide temporal snowpack monitoring. This monitoring provides accurate snow depth maps over small study plots that do not have too much forest cover, and also offers the ability to capture relationships between current and past snow depths that can be used to estimate snow depths for larger areas. In particular, we show that for our study sites, limited spatial extent RPAS-SfM snow depth observations (when combined with past lidar snow depth observations) can be used to generate fairly

Discussion
In this study, we show that RPAS-SfM monitoring of snowpack can be a useful tool to provide temporal snowpack monitoring. This monitoring provides accurate snow depth maps over small study plots that do not have too much forest cover, and also offers the ability to capture relationships between current and past snow depths that can be used to estimate snow depths for larger areas. In particular, we show that for our study sites, limited spatial extent RPAS-SfM snow depth observations (when combined with past lidar snow depth observations) can be used to generate fairly accurate estimates of snow depth across the lidar domains (having R 2 and RMSE values of~0.75-0.9 and 10-15 cm, respectively, as compared with observed lidar snow depth distributions). These estimates improve further when aggregating to larger pixel sizes (Figure 8). Note that even at the 1 m pixel scale, these values are comparable with those of the SfM snow depth maps themselves and a little better (higher R2, lower RMSE) than those for the first guess snow depth maps (compare Figure 8 with Table 1 and Table S1). This implies that SfM should be useful for multitemporal snow monitoring that is applicable at basin-scales. This is important given the high costs of repeated airborne lidar acquisitions. With the exception of a very few high-budget operations (e.g., NASA's Airborne Snow Observatory [29]) this high cost makes it infeasible to collect many repeated snow-on lidar data for most sites. In addition, RPAS-SfM may be especially suitable for operational snowpack monitoring because the SfM data can be processed relatively quickly [39].
The success of the extrapolations shown here depends on two critical factors. First, they rely on the fact that there is a strong relationship between snow depth distributions at different times. We tested this by comparing the snow depth data across the mid-and high-elevation lidar coverages at three different times, representing a range of snowpack conditions, from mid-accumulation to mid-ablation ( Figure 5). Even though the shape of these relationships changed considerably at different times (i.e., it was more linear when comparing snow depths within the accumulation season and nonlinear when comparing accumulation season to ablation season snow depths), their strengths were consistent at different times (R 2~0 .77 and~0.90 for the mid-and high-elevation lidar coverage, respectively). This suggests that they should be useful for predicting snow depth distributions for other times as well.
The second critical factor affecting the success of these extrapolations is the ability of the distributed SfM data to capture these relationships. Our SfM data were generally able to capture these relationships, however, they did a better job for the montane plots than for the dense forest plot. This can be seen by the higher agreement statistics ( Figure 8) and smaller differences between statistics describing the spatial distribution and spatial autocorrelation (Table 5) of the observed and simulated maps for the high-elevation lidar coverage. This ability to reproduce these spatial statistics is particularly important for applications that are based on statistics from lidar snow data (such as snow-cover depletion models [66]). The lower agreement between the observed and simulated maps for the mid-elevation lidar coverage could be because the dense forest plot (from which the SfM data is extrapolated) has a flat topography and mostly uniform forest cover, and so contains a relatively small range of snow conditions. However, it could also be because SfM data quality is much poorer for densely forested conditions [20,67]. We tried to mitigate this by only considering pixels without an overhead canopy. While this improved the ability of the SfM data to capture these relationships, it also ignores under-canopy areas such as tree wells in their derivation. Nevertheless, this does not seem to have too much of an effect on the relationships themselves as they are only sensitive to snow depth changes between the survey dates, which occur similarly for both canopy-covered and open areas (note the similarity between the thick and thin red lines in Figure 8).
Although this method works well at our research sites, more research is needed to understand how well it works at other sites. It is well known that there is substantial similarity between snow distributions at similar points in the snow season (e.g., peak SWE) [68][69][70][71][72][73]. However, there are also substantial spatial differences between snow distributions during the accumulation and melt seasons [46,[74][75][76]. This does not necessarily mean, though, that there are not strong relationships relating melt-season snow depth to accumulation-season snow depth (for our study sites; there were strong relationships; Figure 5), especially since spatial patterns of snow disappearance are related to the distribution of snowpack at the start of the melt season [77].
One thing that is likely to lead to greater difficulties when using the pre-existing lidar observations of snow distributions is land cover change (such as forest thinning, wildfire, or normal forest growth) because these changes are known to alter the snow depth patterns [4,6,9,78]. We showed in Figure 6 that the observed and predicted snow depth distributions were different following forest thinning treatments at the forest thinning comparison site. Climate variability and climate change may also cause issues, since snow distributions for a particular year are influenced by differences in the sequence, timing and magnitudes of snowfall/accumulation events, redistribution (due to wind), interactions with canopy, and melt [50,79,80]. Nevertheless, the high correlations for the snow depth relationships, shown in Figure 5, and for the simulated snow depth maps resulting from these relationships (in Figures 6  and 7) suggest that these differences do not necessarily make it harder to simulate snow depth patterns for different years. That is, even though the shape of these relationships can change under changing climate conditions, it is unclear how their strength would change.
In this study, the SfM snow depth data, themselves are reasonably accurate, and under ideal conditions, they have comparable quality to airborne lidar. For example, at our montane plots, which have relatively gentle topography and sparse tree cover (average slope~10% and average canopy cover~15%), the SfM-generated snow depth maps generally have high correspondence with the airborne lidar-generated snow depth maps (with R 2~0 .75 to 0.85 between the two) and comparable agreement with ground observations (RMSE~10 cm for both technologies). For comparison, most studies that evaluate the performance of SfM for snow depth mapping generally find RMSEs of 5-15 cm in open settings [36,37,40,48,67,81]. SfM can even produce reasonable snowpack maps under very low snow conditions [37]. For example, at the thinning comparison plot in March 2019 (average snow depth~10 cm), there is good correspondence between the RPAS-SfM, airborne lidar, and ground snowpack measurements (Figure 3 and Tables 2 and 3). The performance of SfM (and, to a lesser extent, lidar) is markedly lower under dense forest conditions [20,67]). At our dense forest plot (average canopy cover~45%), the correspondence between the SfM data and the ground measurements is lower than between the lidar data and the ground measurements (RMSE~10 to 20 cm for the SfM maps vs 10 cm for the lidar maps). This performance difference is likely because SfM requires line of sight to collect visual data in order to reconstruct the ground surface, while lidar pulses can split and provide multiple returns, and thus, have a greater ability to penetrate dense tree canopies [27].
The quality of these SfM surveys depends on a variety of factors, including the conditions and parameters used during drone image acquisition, the type of environment being surveyed, as well as the success of each of the processing steps in Section 2.3. In terms of image acquisition, we found that our acquisition parameters (RPAS flight~75 m above ground at~5-7 m per second with~75-80% overlap for the images, and using the RPAS's automatic camera exposure adjustments) were sufficient to capture imagery sufficient to produce good quality point clouds for our study plots. The most significant issues that we experienced were due to deep tree shadows. Although these were not a problem under sparse canopy conditions (in the montane plots), they caused problems for point reconstruction under dense canopy conditions (in the dense vegetation plot). In heavily shaded areas, the dark and poorly resolved ground surface caused spurious noise in the point clouds, which in turn led to problems with the ground filtering (already a substantial source of uncertainty [82][83][84]). Ultimately, this led to spurious snow depth variability for the dense forest plot (such as the artifacts seen in Figure 3).
The automated post-processing steps described in Section 2.3 were also quite important for producing high-quality snow depth maps. Since many of our point clouds were directly georeferenced with onboard GPS data, we relied upon these steps to ensure accurate georegistration of the SfM point clouds, relative to the existing snow-off data. While this approach is capable of producing accurately georeferenced point clouds [85], it increases the reliance on these post-processing steps to produce accurate snow depth maps. In particular, it relies on sufficiently accurate first guess snow depth maps to add to the snow-off data so that they could be used to remove tilting and warping issues that sometimes affect SfM point clouds [86]. For our study sites, these maps ( Figure S2) were reasonably accurate (typically having RMSEs of 10-15 cm, a little higher than the final SfM maps; compare Table 2  and Table S1), leading to the ability to use them reliably in our SfM post-processing steps.
Although the high-resolution snow depth mapping afforded by lidar and SfM are already some of the best ways to monitor distributed snowpack properties in mountain forests, taking the extra step of estimating SWE is enormously valuable because SWE is the snowpack variable that is of most interest in hydrological applications [28,61]. Broxton et al. [54] showed that it is important to consider snow density differences across space at our sites, because even though snow depth variability accounts for a majority of SWE variability, snow density can still show substantial spatial variability on a given date (varying by more than 75% across these research sites). This is in line with other studies finding substantial spatial variability of snow density [62,[87][88][89][90]. Following [54], we use machine learning of how field snow density measurements relate to various lidar-derived physiographic variables to produce maps of snow density, which are then combined with the snow depth maps to produce maps of SWE (see Section 2.3). These maps validate fairly well with snow density / SWE observations at our research sites. The RMSEs of these snow density maps range from 0.024 to 0.033 g/cm 3 , and the RMSEs of the resulting SWE estimates range from~2 to 4 cm [54]. For comparison, the RMSEs of the extrapolated SWE maps range from~4-5 cm at the 1 m scale to~2-3 cm at the 1000 m scale ( Figure 8).
Overall, for multitemporal monitoring of snowpack where SfM is used as a tool to compliment airborne lidar, we would recommend: (1) Conducting an airborne lidar survey representing mid-winter to peak snowpack conditions, to provide a baseline coverage before too much melt so that there are no major snow-free areas over the study area (which would require an alternate technique to extrapolate snowpack conditions in these areas) (2) Planning SfM surveys in areas where SfM has a high chance of success (e.g., sparse canopy, gentle topography, not a whole lot of brush / low vegetation) and where a large range of snow depth variability can be observed (which can include multiple SfM plots), giving a greater ability for successful extrapolation to the larger airborne lidar domain (3) Designing complimentary ground snow surveys that capture snow conditions across a range of physiographic conditions to be able to generate a reliable first-guess snow depth map to aid in the creation of the SfM data (as well as allowing for good characterization of snow densities). These surveys should capture snow variability related to the important physiographic characteristics and snow processes for a given field site. For our area, snow depth variability was primarily influenced by the tradeoff between shading (from both terrain and vegetation and terrain influences) and interception [6]. Therefore, it was ideal to have transects spanning gradients of vegetation cover and tree shading (e.g., across forest gaps and into forest edges) and on opposing aspects.
The existence of snow-on lidar data makes these planning tasks easier as objective criteria (based on these data) can be used for site selection. The RPAS flight acquisition parameters (see Section 2.1) and processing steps (see Section 2.2) used in this study seemed to be adequate to produce good quality SfM snow maps, but it is also likely that additional improvements could be made. For example, using more exact positioning (e.g., using RTK GPS), which can provide cm level precision [91][92][93]), could help with model positioning and reducing the warping and tilting of the SfM models, thus reducing the importance of post-processing to correct the point cloud data. Also, the use of drone-based lidar would perform better than SfM in areas with thick canopy [20] (with the obvious drawback that it is much more expensive). Finally, using a fixed-wing aircraft that can cover larger areas may improve the ability to capture the snow spatial variability needed for basin-scale snow estimation, though this option may be limited in some denser forests, due to line-of-sight requirements in places such as the US.

Conclusions
Although airborne lidar has revolutionized our ability to make distributed measurements of snowpack in mountain forests, it is expensive and time-consuming to process, making it hard to use for real-time or multi-temporal snow monitoring. This study demonstrates that SfM based on multi-angle imagery from a multi-rotor RPAS and airborne lidar can provide complementary information for high-quality snowpack monitoring, as the use of RPAS allows for relatively low cost, on-demand snow monitoring over small plots, and airborne lidar provides monitoring over a much larger area and at a higher cost. The existence of snow-on lidar data may allow subsequent SfM data to be extrapolated to much larger areas by comparing the spatially distributed snow thicknesses from the SfM data with those of already-measured lidar snow depths. This would make SfM data much more useful for basin-scale hydrological applications such as Land Surface Model evaluation and water supply monitoring. At our study sites, using the relationships developed between the small plot SfM data and previously flown airborne lidar resulted in fairly accurate snowpack estimates for the larger domains (R 2 values between the extrapolated maps and observed lidar data that is used for verification are~0.7-0.9, and RMSE's between the two are~10-15 cm for a 1 m pixel scale, with even higher performance when aggregating the data to larger pixels). However, the methodology proposed here may not always be appropriate (e.g., after extensive land cover changes), and more research is needed to understand its effectiveness at other research sites. Ultimately, successful incorporation of SfM and airborne lidar will be important for cost-effective multitemporal monitoring of snowpack.