{'fraction_non_alphanumeric': 0.03203125, 'fraction_numerical': 0.019479166666666665, 'mean_word_length': 5.020171404682274, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
UAV-Based Estimate of Snow Cover Dynamics: Optimizing Semi-Arid Forest Structure for Snow Persistence

: Seasonal snow cover in the dry forests of the American West provides essential water resources to both human and natural systems. The structure of trees and their arrangement across the landscape are important drivers of snow cover distribution across these forests, varying widely in both space and time. We used unmanned aerial vehicle (UAV) multispectral imagery and Structure-from-Motion (SfM) models to quantify rapidly melting snow cover dynamics and examine the effects of forest structure shading on persistent snow cover in a recently thinned ponderosa pine forest. Using repeat UAV multispectral imagery (n = 11 dates) across the 76 ha forest, we ﬁrst developed a rapid and effective method for identifying persistent snow cover with 90.2% overall accuracy. The SfM model correctly identiﬁed 98% (n = 1280) of the trees, when compared with terrestrial laser scanner validation data. Using the SfM-derived forest structure variables, we then found that canopy shading associated with the vertical and horizontal metrics was a signiﬁcant driver of persistent snow cover patches ( R 2 = 0.70). The results indicate that UAV image-derived forest structure metrics can be used to accurately predict snow patch size and persistence. Our results provide insight into the importance of forest structure, speciﬁcally canopy shading, in the amount and distribution of persistent seasonal snow cover in a typical dry forest environment. An operational understanding of forest structure effects on snow cover will help drive forest management that can target snow cover dynamics in addition to forest health.


Introduction
Runoff from seasonal snowpack in semi-arid regions provides drinking water and agricultural irrigation for at least one-sixth of the world's population [1]. Snowpack also provides important ecosystem services, including water for vegetation, aquatic ecosystems, and shallow groundwater recharge [2][3][4][5][6]. Throughout the western United States, higher annual mean temperatures and changes to winter air humidity are contributing to an earlier and faster spring snowmelt [7]. Regional climate projections indicate that these effects will become more pronounced [8][9][10] and inconsistent snowpack will contribute to higher rates of drought-induced tree stress and mortality [11].
In the ponderosa pine (Pinus ponderosa) forests of northern and central Arizona, roughly half of the annual precipitation falls as rain during the summer, and the other half as snow during the winter [12]. For mature ponderosa pine trees, winter precipitation via snowpack is the dominant source of water throughout the year, underscoring the importance of snowpack to forest health [13,14]. In turn, forest cover influences the accumulation and ablation of snow on the ground surface, primarily by controlling canopy 1.
Quantify snow cover following three different winter storms and identify persistent snowpack across the study site; 2.
Examine forest structure shading effects on snowpack persistence; 3.
Model and predict persistent snowpack using the most influential forest structure metrics.
To assess the importance of specific characteristics and their spatial distribution, forest structure metrics were separated into groups that emphasize vertical and horizontal characteristics. We hypothesized that horizontal forest structure metrics would be more influential for predicting persistent snowpack than vertical forest structure metrics due to the variation in ground shading exhibited by trees, no matter their size or shape. In addition, we hypothesized that crown base height and crown volume would be important for predicting the size of persistent snow patches.

Study Area
Our study area is located in the Coconino National Forest in northern Arizona ( Figure 1). The study area includes 76 ha of forest, characterized by relatively flat topography ranging in elevation between 2200 and 2275 m, with slopes of 0-10% across most of the study area. The vegetation is characterized by ponderosa pine (Pinus ponderosa), which dominates both the region and the study area, and includes sporadic Gambel oak (Quercus gambelii) and Rocky Mountain juniper (Juniperus scopulorum). The climate is subhumid and characterized by distinct seasonal trends including early summer drought, with an average of 560 mm of precipitation (Western Regional Climate Center, 2020) falling half as snow during winter months and half as rain during the mid-summer monsoon season. The site had been undisturbed since its last naturally occurring fire in 1876, except for selective historical firewood harvesting [56]. However, the site was subjected to a prescribed fire in 1976, in which 63% of the smaller surface fuels and 69% of the woody surface fuels (up to 8 cm in diameter) were consumed [56]. and characterized by distinct seasonal trends including early summer drought, with an average of 560 mm of precipitation (Western Regional Climate Center, 2020) falling half as snow during winter months and half as rain during the mid-summer monsoon season. The site had been undisturbed since its last naturally occurring fire in 1876, except for selective historical firewood harvesting [56]. However, the site was subjected to a prescribed fire in 1976, in which 63% of the smaller surface fuels and 69% of the woody surface fuels (up to 8 cm in diameter) were consumed [56].
A mechanical thinning restoration treatment at the study area began during the fall of 2017 and was completed by the spring of 2018, yielding a 53-ha treated portion and 23ha untreated portion for this study (Figure 1). Similar to other regional restoration treatments that promote diversity in tree group and interspace size, shape, and spacing, this treatment aimed to reinstate pre-settlement forest conditions and included a range of thinning goals that would promote healthy overstory vegetation and the regeneration of understory vegetation [56][57][58][59]. Specifically, the restoration treatment prescription at our study area emphasized irregular tree group delineation, expansion of interspace, retention of all non-ponderosa pine species (e.g., Gambel oak and juniper), and significant reductions in smaller ponderosa pine trees within groups and interspaces.

Snow Covered Area and Persistent Snow Patches
Our examination of the relationship between persistent snow cover and forest structure began with quantifying snow covered area (SCA) across the entire study site throughout the melt-off period following three significant (>20 cm new snowfall depth) independent winter storm events during the winters of 2017-2018 and 2018-2019. For each storm, a set of 15-cm resolution multispectral UAV orthomosaic images were acquired from the Figure 1. An overview of the study site highlighting the locations and extents of the UAV Structure-from-Motion (SfM) data, terrestrial laser scanner (TLS) data, and field-based validation data. (A) shows an example of the field-measured and TLS validated plots. (B) shows the distribution of all field-measured and TLS validated plots within the SfM data extent, which includes both thinned and unthinned portions of the forest. A majority of our ground-based measurements are distributed across the larger, thinned portion of the study site. (C) shows the UAV and TLS instruments used in this study.
A mechanical thinning restoration treatment at the study area began during the fall of 2017 and was completed by the spring of 2018, yielding a 53-ha treated portion and 23-ha untreated portion for this study (Figure 1). Similar to other regional restoration treatments that promote diversity in tree group and interspace size, shape, and spacing, this treatment aimed to reinstate pre-settlement forest conditions and included a range of thinning goals that would promote healthy overstory vegetation and the regeneration of understory vegetation [56][57][58][59]. Specifically, the restoration treatment prescription at our study area emphasized irregular tree group delineation, expansion of interspace, retention of all non-ponderosa pine species (e.g., Gambel oak and juniper), and significant reductions in smaller ponderosa pine trees within groups and interspaces.

Snow Covered Area and Persistent Snow Patches
Our examination of the relationship between persistent snow cover and forest structure began with quantifying snow covered area (SCA) across the entire study site throughout the melt-off period following three significant (>20 cm new snowfall depth) independent winter storm events during the winters of 2017-2018 and 2018-2019. For each storm, a set of 15-cm resolution multispectral UAV orthomosaic images were acquired from the first day following the storm (i.e., maximum SCA) until the snow cover had fully Remote Sens. 2021, 13, 1036 4 of 20 disappeared (i.e., minimum SCA). To capture temporal changes in SCA following each storm, the interval at which images were acquired depended on the weather conditions in the preceding days, but we attempted to keep image date intervals consistent among independent storm events ( Table 1). The final SCA dataset consisted of three 'snow-series' (one snow-series per storm event), yielding a total of 11 different orthomosaic images. Table 1. Three storm events, subsequent image collection dates, and weather conditions. The total snowfall indicates the snow recorded for the storm event, while the mean temperatures and wind speeds are cumulative means between the image dates. Each storm event was imaged several times from full snow coverage on the first day until the snow cover nearly melted off. The individual image dates were selected considering daily weather conditions and study site access. Periods of near-freezing daily temperatures yielded longer image date intervals so that noticeable changes in snow coverage could be observed. In contrast, periods of higher mean temperature values required more frequent image date intervals as the snow melted rapidly. Weather data were recorded at the regional National Weather Service Forecast office, approximately 15 km from the study site and compiled from NOAA's Climate Data Online service. Data comprising the snow-series orthomosaic images were collected using a Sensefly eBee RTK fixed-wing UAV platform (SenseFly, Lausanne, Switzerland) fitted with a Parrot Sequoia multispectral sensor (Parrot Drones SAS, Paris, France). Each flight mission was pre-programmed using Sensefly's proprietary software, which controlled flight plans and customized all flight and data parameters. Flight specifications included high latitudinal and longitudinal overlaps (85% and 90%, respectively) in order to generate the SfM outputs, an average flight altitude of 120 m, and operation centered around solar noon on each image date. There was an average of~350 images captured for each flight based on the pre-determined flight path, with wind conditions largely responsible for the total number of images collected. During each flight, four georeferenced images were recorded at each photo location in the green (530-570 nm), red (640-680 nm), red edge (730-740 nm), and near infrared (770-810 nm) spectral bands. These images were post-processed using Sensefly's eMotion 3 software, which automatically excluded a few distorted images from each flight (SenseFly, Lausanne, Switzerland).

Storm
The images were then used to create the final datasets in Agisoft PhotoScan v1.4.0 photogrammetric processing software (Agisoft LLC, St. Petersburg, Russia) for each image date. All images from each flight were scanned for matching 'tie-points', oriented in three-dimensional space via bundle-adjustment, and then mosaicked together [39]. The general workflow in the software includes image alignment to create a sparse point cloud, incorporation of GCP locations, image alignment optimization, gradual filtering out of inaccurate and error-inducing points, and a full image realignment [39]. This workflow generated a final orthomosaicked image in 15 cm spatial resolution in all four spectral bands and dense 3-dimensional (3D) point cloud data photogrammetrically generated from the high resolution images using Structure-from-Motion (SfM) algorithms. We optimized the SfM algorithm parameters based on our previous results [39].
The snow-series images were used to first classify SCA across the study site in each image and subsequently delineate persistent snow cover patches. All images (n = 11) were classified into five classes: vegetation, sunlit bare ground, shaded bare ground, sunlit snow, and shaded snow, using a random forest supervised classification performed in Google Earth Engine [60,61]. The random forest classifier was parameterized based on previous remotely sensed image classification examples, with the number of trees set to 100 and the number of variables per split set to the square root of the number of variables [60,61]. The training dataset consisted of manually digitized polygons (n = 492 total) for each of the five classes. The image resolution (15 cm/pixel) allowed for the training dataset classes to be easily delineated. The final binary snow versus non-snow rasters were created by combining pixels classified as both sunlit and shaded snow into a single 'snow' class, while pixels in the remaining classes were combined into a single 'non-snow' class. The accuracy of this binary classification was assessed using a set of randomly generated snow/nonsnow samples (n = 500) from the image with roughly half its pixels snow covered, allowing for an unbiased accuracy assessment. Image classification relied on a multi-band image stack created from the four original spectral bands (green, red, red edge, near infrared) and six additional bands derived from the original four: normalized difference vegetation index (NDVI) [62], soil adjusted vegetation index (SAVI) [63], Gray Level Co-Occurrence Matrix (GLCM) [64] variance, GLCM homogeneity, GLCM contrast, and GLCM entropy. The GLCM texture bands were included to increase the effectiveness of the classification algorithm in discriminating the different combinations of bare ground, snow, and their shaded equivalents.
Finally, persistent snow cover was delineated using a single snow cover image composite. This image composite was created by stacking each binary classification (n = 11) and counting when each pixel was snow-covered out of the 11 total dates in the image composite. A simple post-processing procedure was used to eliminate spurious pixels and reduce noise along edges of snow pixels. Specifically, we conducted a progressive moving-window majority filtering using windows of 3 × 3 and 5 × 5 cells. From this image composite, persistent snow patches were identified by selecting isolated groups of adjacent pixels with snow cover in 10 or 11 out of the total 11 images ( Figure 2). Groups of pixels less than 200 m 2 were numerous and observed to melt completely between image dates, thus they were eliminated from the analysis. The resulting patches were refined using a general boundary cleaning procedure to eliminate spurious edges, then vectorized and attributed with an identification number and an area estimate (m 2 ).

Forest Structure Metrics
Forest structure metrics were quantified across the entire study site using the UAV SfM point cloud data and validated using both field measurements and a terrestrial laser scanner (TLS) point cloud dataset. The vertical forest structure metrics examined in our study are tree crown height (CH), crown diameter (CD), crown base height (CBH), the ratio of tree base height to crown height (CBH:Z), and crown volume (CV). The horizontal structure metrics in our study include trees per hectare (TPH), canopy cover (%) per area (CC), average of five nearest neighbor distances (KNN), Clark and Evans Index (CEI), a canopy height-based solar radiation footprint (SRF), and an average distance from the northern canopy edge (NCE).
The UAV SfM point cloud data were collected across the entire study site with the same platform and sensor used for snow-series data collection and following completion of the mechanical thinning restoration treatment [39]. The final orthomosaic has error estimates of 1.14 m and 1.80 in the X, Y, and Z dimensions, respectively [39]. In addition, the SfM point cloud was used to create a digital terrain model (DTM) using the Cloth Simulation Filter (CSF) tool in CloudCompare v2.9.1 [39]. The CSF tools was parameterized using 'relief' for Scenes, a Cloth Resolution of 1 m, and a Classification Threshold of 0.7 m. The accuracy of the DTM (R 2 = 0.95 and a RMSE = 2.98 m) was assessed by comparing points . UAV multispectral image analysis workflow. First, the original multispectral bands and the calculated indices were stacked into a single raster image, one for each image date (n = 11) and classified using a supervised random forest model to generate binary snow/non-snow classes. From the binary classification from each image date, the pixels with 10 or 11 snow-covered days (out of 11) were grouped to create the boundaries of individual persistent snow patches.

Forest Structure Metrics
Forest structure metrics were quantified across the entire study site using the UAV SfM point cloud data and validated using both field measurements and a terrestrial laser scanner (TLS) point cloud dataset. The vertical forest structure metrics examined in our study are tree crown height (CH), crown diameter (CD), crown base height (CBH), the ratio of tree base height to crown height (CBH:Z), and crown volume (CV). The horizontal structure metrics in our study include trees per hectare (TPH), canopy cover (%) per area (CC), average of five nearest neighbor distances (KNN), Clark and Evans Index (CEI), a canopy height-based solar radiation footprint (SRF), and an average distance from the northern canopy edge (NCE).
The UAV SfM point cloud data were collected across the entire study site with the same platform and sensor used for snow-series data collection and following completion of the mechanical thinning restoration treatment [39]. The final orthomosaic has error estimates of 1.14 m and 1.80 in the X, Y, and Z dimensions, respectively [39]. In addition, the SfM point cloud was used to create a digital terrain model (DTM) using the Cloth Simulation Filter (CSF) tool in CloudCompare v2.9.1 [39]. The CSF tools was parameterized using 'relief' for Scenes, a Cloth Resolution of 1 m, and a Classification Threshold of 0.7 m. The accuracy of the DTM (R 2 = 0.95 and a RMSE = 2.98 m) was assessed by comparing points extracted from the DTM to corresponding points from differentially corrected Trimble GeoXH GPS elevation values collected from the trees in the field-based validation dataset.
Individual tree segmentation was performed on both the SfM and TLS point clouds using the Li et al. [65] algorithm in the lidR package [66] in RStudio (R-Studio Team, 2015). The segmented point clouds were then used to estimate each tree's location (X, Y in UTM 12N m), maximum crown height (CH in m), and average of its widest and narrowest crown diameters (CD in m) using the rLiDAR package in R-Studio [67]. Each tree's crown base height (CBH in m) was estimated from the point cloud data using a novel multiple changepoint detection algorithm [68]. This algorithm used the area of two-dimensional Figure 2. UAV multispectral image analysis workflow. First, the original multispectral bands and the calculated indices were stacked into a single raster image, one for each image date (n = 11) and classified using a supervised random forest model to generate binary snow/non-snow classes. From the binary classification from each image date, the pixels with 10 or 11 snow-covered days (out of 11) were grouped to create the boundaries of individual persistent snow patches.
Individual tree segmentation was performed on both the SfM and TLS point clouds using the Li et al. [65] algorithm in the lidR package [66] in RStudio (R-Studio Team, 2015). The segmented point clouds were then used to estimate each tree's location (X, Y in UTM 12N m), maximum crown height (CH in m), and average of its widest and narrowest crown diameters (CD in m) using the rLiDAR package in R-Studio [67]. Each tree's crown base height (CBH in m) was estimated from the point cloud data using a novel multiple changepoint detection algorithm [68]. This algorithm used the area of two-dimensional convex hulls, which were calculated at every 50 cm along a tree's height to isolate the height at which trunk points transitioned into canopy points. Finally, each tree's CV (m 3 ) was then estimated by calculating the volume of the three-dimensional convex hull made from each tree's points above its CBH value.

Forest Structure Metrics Validation
We assessed the accuracy of the SfM point-cloud-derived tree metrics by comparing them to both the field-measured and terrestrial laser scanner (TLS) point cloud-derived tree metrics ( Figure 3). The TLS point cloud data provided accurate high-resolution estimates of all field-measured forest structure metrics and spatially extended the field-measured metrics while also alleviating known accuracy issues associated with UAV SfM-derived crown base height (CBH) and crown volume (CV) estimates [40]. Since we used a newer TLS instrument that has not been previously evaluated, an assessment of the omission and commission error rates was performed between the field-measured and TLS-derived datasets as well as between the TLS-and SfM-derived datasets. This facilitated direct comparison between trees from the TLS and SfM point clouds. tree metrics (Figure 3). The TLS point cloud data provided accurate high-resolution estimates of all field-measured forest structure metrics and spatially extended the field-measured metrics while also alleviating known accuracy issues associated with UAV SfM-derived crown base height (CBH) and crown volume (CV) estimates [40]. Since we used a newer TLS instrument that has not been previously evaluated, an assessment of the omission and commission error rates was performed between the field-measured and TLSderived datasets as well as between the TLS-and SfM-derived datasets. This facilitated direct comparison between trees from the TLS and SfM point clouds. The TLS point cloud data was acquired within 0.64 ha field plots (n=16) covering a total area of 5.17 ha and spread along a gradient of forest density present across the study site. The TLS data were collected using a Leica Geosystems BLK360 Imaging Laser Scanner (Leica Geosystems AG, 2020), which has a range of up to 60 m radius. The BLK360 captures 360,000 points per second at 830 nm wavelength with 300 and 360 degrees of vertical and horizontal field of view, respectively, and with 6mm-10m accuracies. Each TLS validation plot included at least three scans, which were tied together and georeferenced using four distinct reference targets. Reference targets were 50 × 50 cm reflective panels bolted to adjustable tripods at 1 m above ground, and their positions were recorded using a Trimble GeoXH GPS unit. We determined the locations and configuration of the The TLS point cloud data was acquired within 0.64 ha field plots (n = 16) covering a total area of 5.17 ha and spread along a gradient of forest density present across the study site. The TLS data were collected using a Leica Geosystems BLK360 Imaging Laser Scanner (Leica Geosystems AG, 2020), which has a range of up to 60 m radius. The BLK360 captures 360,000 points per second at 830 nm wavelength with 300 and 360 degrees of vertical and horizontal field of view, respectively, and with 6mm-10m accuracies. Each TLS validation plot included at least three scans, which were tied together and georeferenced using four distinct reference targets. Reference targets were 50 × 50 cm reflective panels bolted to adjustable tripods at 1 m above ground, and their positions were recorded using a Trimble GeoXH GPS unit. We determined the locations and configuration of the scanner and reference targets using guidelines developed in our previous study [34]. Using Cloud Compare software, the scans were co-registered using the target center points, then merged into a single point cloud and georeferenced using the GPS coordinates of the target; the overall X, Y, and Z positional accuracy of the 16 point clouds is 0.46 m. The average point spacing of all point clouds was 0.03 m with a range of 0.025-0.046 m, and each was subsampled to 0.01 m to reduce redundant points and ensure consistent point spacing across scans.
The TLS point cloud data were validated using field-based measurements taken from 137 trees within 30 × 30 m validation plots (n = 16), covering a total of 1.44 ha. The field-based validation plots were centered within the TLS validation plots, providing as much overlap between the datasets as possible. Using the Trimble GeoXH handheld GPS unit with an attached Trimble laser range finder module, each tree's geographic position (X,Y,Z coordinates), CH, CD, and CBH were measured and differentially corrected in GPS Pathfinder Office. The effect of forest canopy shading on persistent snow patches was quantified using distinct tree shading influence areas. Previous research shows that a tree's ground shading influence extends between 1 and 2 times its crown height during the winter season in forests of the southwestern U.S. [55,69,70]. To support and refine this, we used three different spatial extents of localized tree shading to assess the size of persistent snow patches. These tree shading influence areas (TSIAs) were established with respect to each individual persistent snow patch and included trees located within a distance of 1, 1.5, and 2 times their crown heights. To ensure that only the trees capable of influencing a persistent snow patch were included, trees were selected based on whether their shadows were within solar azimuth angle extents specific to study site location and snow-series image date. The resulting minimum bounding extents were termed TSIA1.0, TSIA1.5, and TSIA2.0 ( Figure 4). Once trees were selected in each TSIA, their vertical forest structure metrics were averaged to produce a final dataset.
The TLS point cloud data were validated using field-based measurements taken 137 trees within 30×30 m validation plots (n = 16), covering a total of 1.44 ha. The based validation plots were centered within the TLS validation plots, providing as overlap between the datasets as possible. Using the Trimble GeoXH handheld GPS with an attached Trimble laser range finder module, each tree's geographic po (X,Y,Z coordinates), CH, CD, and CBH were measured and differentially corrected in Pathfinder Office.

Forest Structure Predictor Variables
The effect of forest canopy shading on persistent snow patches was quantified distinct tree shading influence areas. Previous research shows that a tree's ground sh influence extends between 1 and 2 times its crown height during the winter seas forests of the southwestern U.S. [55,69,70]. To support and refine this, we used thre ferent spatial extents of localized tree shading to assess the size of persistent snow pa These tree shading influence areas (TSIAs) were established with respect to each ind ual persistent snow patch and included trees located within a distance of 1, 1.5, times their crown heights. To ensure that only the trees capable of influencing a pers snow patch were included, trees were selected based on whether their shadows within solar azimuth angle extents specific to study site location and snow-series i date. The resulting minimum bounding extents were termed TSIA1.0, TSIA1.5 TSIA2.0 ( Figure 4). Once trees were selected in each TSIA, their vertical forest stru metrics were averaged to produce a final dataset.  In addition to the vertical forest structure metric summaries, each TSIA was assigned a set of horizontal arrangement metrics to quantify the spatial distribution of trees. Forest density was estimated by calculating trees per hectare (TPH) and average canopy cover (CC) in percent, while spacing of trees was measured using the average nearest neighbor (KNN) distance in meters from each tree to its five nearest neighbors. The clustering of trees was measured using a unitless Clark and Evans Index (CEI) value, with values < 1 indicating ordering and values > 1 indicating clustering. The amount of ground shading was expressed using an average solar radiation footprint (SRF) value in w/m 2 for each TSIA polygon. The SRF was created using a sitewide canopy height model and the solar radiation toolset in ArcMap 10.8 to calculate the cumulative amount of incoming solar radiation (w/m 2 ) during each snowstorm event [71,72]. The ground surface pixel values were averaged across the three storms and then across each TSIA polygon. Finally, each pixel's distance to the northern canopy edge (DNCE) was calculated and averaged across each TSIA using a binary canopy cover raster based on the methodology developed by Mazzotti et al. [73]. The DNCE values provided a fine-scale metric that describes the directional within-stand differences of forest canopy spatial arrangement, and are well suited for assessing the effect of forest canopy radiative transfer.

Model Framework
An initial exploratory data analysis revealed sources of significant multicollinearity between forest structure metrics (independent variables) as well as numerous complex and non-linear relationships with persistent snow patch area (dependent variable). A variable selection process was used to eliminate the highly correlated and less descriptive forest structure metrics using variable inflation factor (VIF) scoring. To find a suitable model, we parameterized generalized linear models, random forest, support vector machine, and Multivariate Adaptive Regression Spline (MARS) algorithms for regression and tested their predictive accuracies using the Caret package [74] in RStudio. We selected the MARS model framework since it provided both the most accurate predictive accuracy as well as offering variable importance scores, which helped us interpret results into meaningful forest management recommendations.
The MARS model performs nonparametric multivariate regression without underlying assumptions of the data and is useful for regression problems with high-dimensional datasets and multiple predictor variable interactions. Our modeling framework assessed the relationship between persistent snow patch size (m 2 ) and forest structure metrics for each of the three different tree shadow influence areas (TSIAs). Initial attempts at model fitting and prediction presented inconsistent results, with model hyperparameters and predictive accuracies fluctuating based on different randomly selected training and testing datasets, which can be apparent when using machine learning methods on smaller datasets [75]. Due to our relatively small sample size (n = 99), we adopted a split-sample validation procedure that was repeated 500 times. This provided 500 different randomly selected testing and training data partitions, model construction, and predictive accuracy assessments. While the number of repetition is usually arbitrarily chosen, we found that 500 repetitions provided a large enough sample size to clearly discriminate between trends in variable importance. Each model iteration utilized an exhaustive grid search to select the optimal set of hyperparameters, allowing for the inclusion of all possible predictors as well as for two-way interactions between predictors. An optimal model was then selected from each iteration by minimizing the generalized cross-validation error estimate, from which variable importance scores ranging from 0 to 100 were calculated for each model term.
Model results were summarized based on the frequency of both the hyperparameter value and the variable importance scores. Separately, the optimal model framework for each iteration was used for prediction on the testing dataset. Instead of arriving at a single ideal model that definitively explains the relationship between persistence snow patch area and forest structure, this approach provides a robust conceptual understanding of the relationship.

Snow Cover Classification and Persistent Snow Patches
The snow cover classification accuracy assessment was performed using the third image from the second storm due to its mostly even distribution of snow/non-snow pixels. The snow classification performed well, with an overall accuracy of 90.2% and a kappa coefficient of 0.80. Relatively low rates of omission error were observed for the non-snow and snow classes at 13% and 6%, respectively, while similar rates of commission error were observed for the non-snow and snow classes at 6% and 14%, respectively (Table 2). Table 2. Accuracy of the snow/non-snow classification across the 76-ha study site. The image used to generate the error matrix included an even amount of snow and non-snow pixels, providing the least biased estimation of classification error. With the average daily high temperatures similar across all storm events (6.1 • C, 5.8 • C, and 5.7 • C for storms 1-3, respectively), storm 1 had the lowest initial snowfall amount (28.9 cm) and the lowest reduction in site-wide SCA (−51%) over 8 days. In contrast, storms 2 and 3 had greater initial snowfall amounts, 73.4 cm and 93.9 cm, respectively, and site-wide reductions in SCA were consistently higher and for longer time periods: −59% over 15 days and −70% over 26 days, respectively ( Table 1).

Class
The proportion of SCA in the treated portion of the study site on Day 1 was consistent across all storms, ranging from 91.2% to 89.4% to 90% for storms 1, 2, and 3, respectively. The proportion of SCA across the untreated portion of the study site on Day 1 of each storm was consistently lower, with values ranging from 43.4% to 62.3% to 58.1% for storms 1, 2, and 3, respectively. The effect of treatment condition on the magnitude of SCA reduction is evident both within and across storm events. The treated portion of the study site consistently exhibited a wider range of reductions in SCA within each storm compared to the narrower range in reductions observed in the untreated portion. The average of each storm's total SCA reduction in the treated portion was −76.5% compared to −38.6% in the untreated portion ( Figure 5).
Distinct persistent snow patches were identified from the final composite of classified snow-series images from grouped pixels that were covered in snow for 10 or 11 of the total days (n = 11) ( Figure 6A). A total of 99 individual snow patches were delineated across the entire study site, covering 8.4% (6.36 ha) of the total area and ranging in size from 203 to 2699 m 2 (SD = 469.7 m 2 ), with an average size of 646.9 m 2 and a majority (82%) being less than 1000 m 2 . The patches delineated across the treated portion of the study site (n = 67) had an average size of 722 m 2 and covered 9.2% (4.84 ha) of the treated area ( Figure 6B). In contrast, the patches in the untreated portion of the study site (n = 32) had an average size of 474 m 2 and covered 6.6% (1.52 ha) of that area ( Figure 6B).

Forest Structure Metrics Summary
The TLS-derived point cloud correctly identified 88% (n = 120) of trees from within the field-measured validation plots. The final TLS dataset consisted of 1293 trees identified within the TLS plot data extents (n = 16 validation plot, totaling 5.17 ha). The SfM-derived point cloud correctly identified 98% (n = 1280) trees within the TLS point cloud data extent, and the final dataset included a total of 8847 trees across the 76-ha study site. The individual tree metrics compared between all datasets are identical except for CV, which was calculated from TLS and SfM point clouds (Figure 2). Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 20 Figure 5. The percent of the treated and untreated regions of the study site that are classified SCA for each UAV snowseries image following a storm. This data is grouped by storm, and partitioned by forest treatment condition to illustrate the reductions in SCA within and across the storms. SCA patterns following storm events show a greater SCA reduction in the treated portion of the study site than in the untreated portion. Evident also are the relatively similar initial amounts of classified SCA on the first image date after a storm in the treated portion, contrasting the initial amounts in the untreated portion. Figure 6. Each UAV image (n = 11) was classified into snow/non-snow pixels (A) and the resulting rasters were assembled into a single composite, in which persistent snow patches were identified (B). The orthomosaic image used as the base imagery is from storm 2, day 3 (1 May 2019) and its binary snow classification (A) was used for the classification accuracy assessment because it contained a nearly even snow/non-snow (49/51%) area distribution. Figure 5. The percent of the treated and untreated regions of the study site that are classified SCA for each UAV snow-series image following a storm. This data is grouped by storm, and partitioned by forest treatment condition to illustrate the reductions in SCA within and across the storms. SCA patterns following storm events show a greater SCA reduction in the treated portion of the study site than in the untreated portion. Evident also are the relatively similar initial amounts of classified SCA on the first image date after a storm in the treated portion, contrasting the initial amounts in the untreated portion.
Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 20 Figure 5. The percent of the treated and untreated regions of the study site that are classified SCA for each UAV snowseries image following a storm. This data is grouped by storm, and partitioned by forest treatment condition to illustrate the reductions in SCA within and across the storms. SCA patterns following storm events show a greater SCA reduction in the treated portion of the study site than in the untreated portion. Evident also are the relatively similar initial amounts of classified SCA on the first image date after a storm in the treated portion, contrasting the initial amounts in the untreated portion. Figure 6. Each UAV image (n = 11) was classified into snow/non-snow pixels (A) and the resulting rasters were assembled into a single composite, in which persistent snow patches were identified (B). The orthomosaic image used as the base imagery is from storm 2, day 3 (1 May 2019) and its binary snow classification (A) was used for the classification accuracy assessment because it contained a nearly even snow/non-snow (49/51%) area distribution. Figure 6. Each UAV image (n = 11) was classified into snow/non-snow pixels (A) and the resulting rasters were assembled into a single composite, in which persistent snow patches were identified (B). The orthomosaic image used as the base imagery is from storm 2, day 3 (1 May 2019) and its binary snow classification (A) was used for the classification accuracy assessment because it contained a nearly even snow/non-snow (49/51%) area distribution.
The accuracy of the individual tree structure metrics used in the final model are a product of a multi-scale accuracy assessment, which compared the field-measured, TLSderived, and SfM-derived estimates. The relationships between the field-measured and TLS-derived tree metrics were generally strong, while the relationships between the TLSand SfM-derived tree metrics varied.
The final set of forest structure predictors used for modeling were selected to balance both accurate tree dimension estimation as well as to maximize the model's predictive capacity. The vertical forest structure metrics were selected by first examining the relationships between SfM-and TLS-derived estimates, then by assessing their multicollinearity and variance inflation factor (VIF) scores. For example, the weak relationship between SfMand TLS-derived CD (R 2 = 0.14; RMSE = 1.62 m) indicated that it should not be used as a predictor variable. In addition, CD was highly correlated with CV (R 2 = 0.72) and had a high VIF score (12.3), indicating that it would negatively impact model accuracy and interpretability. For the entire study site, the complete set of vertical and horizontal forest structure metrics included as model predictors are summarized by TSIA size (1, 1.5, and 2).

Relationship between Snow and Forest Structure
The most frequent number of model terms for TSIA 1.0 was 5 (50% of all model iterations), for TSIA 1.5 were 6 (35% of all model iterations), and for TSIA 2.0 were 7 (35% of all model iterations), while the most frequent interaction degree term was 1 for all TSIA groups. The average model error estimates calculated across all model iterations grouped by TSIA size were 0.14 for TSIA 1.0, 0.11 for TSIA 1.5, and 0.14 for TSIA 2.0. As the influence region size increased, more predictor variables were needed to stabilize the model and there was less consensus in the number of predictors selected.
For all TSIA sizes, the forest structure metrics with the highest variable importance scores were the mean canopy cover (CC), mean solar radiation footprint (SRF), and trees per hectare (TPH) values (Figure 7). While this remained relatively consistent overall, differences in variable importance scores existed across TSIA sizes as evidenced by the differences between variable importance scores from TSIA groups within CEI, mean SRF, and mean TPH metrics. Specifically, the greater importance were observed for CEI at the largest spatial scale and the lesser importance of mean SRF and mean TPH at the smallest spatial scale.
The models describing forest structure at TSIA 1.5 performed best overall with a mean prediction accuracy R 2 = 0.70 (RMSE = 267 m 2 ) (Figure 8), followed by those at TSIA 1.0 and TSIA 2.0 with mean prediction accuracies of R 2 = 0.66 (RMSE = 286 m 2 ) and R 2 = 0.61 (RMSE = 307 m 2 ), respectively. Of note is the increasing variability in model predictive accuracy on both ends of the persistent snow patch area, especially on the larger end where wide fanning indicates greater uncertainty. Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 20 Figure 7. Variable importance scores calculated for the forest structure metric predictor variables used in fitting the final MARS model framework. The variable importance score ranges from 0 to 100, with higher scores indicating that the predictor was more influential in model construction (n=500). The chart depicts the frequency with which each forest structure metric was within the respective range of the importance score. The data are partitioned by tree shading influence area (TSIA) size to better illustrate the interaction with forest structure metric. Overall, the most influential forest structure metric predictors are canopy cover (CC), mean solar radiation footprint (SRF) value, and mean trees per hectare (TPH).  Variable importance scores calculated for the forest structure metric predictor variables used in fitting the final MARS model framework. The variable importance score ranges from 0 to 100, with higher scores indicating that the predictor was more influential in model construction (n = 500). The chart depicts the frequency with which each forest structure metric was within the respective range of the importance score. The data are partitioned by tree shading influence area (TSIA) size to better illustrate the interaction with forest structure metric. Overall, the most influential forest structure metric predictors are canopy cover (CC), mean solar radiation footprint (SRF) value, and mean trees per hectare (TPH).
Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 20 Figure 7. Variable importance scores calculated for the forest structure metric predictor variables used in fitting the final MARS model framework. The variable importance score ranges from 0 to 100, with higher scores indicating that the predictor was more influential in model construction (n=500). The chart depicts the frequency with which each forest structure metric was within the respective range of the importance score. The data are partitioned by tree shading influence area (TSIA) size to better illustrate the interaction with forest structure metric. Overall, the most influential forest structure metric predictors are canopy cover (CC), mean solar radiation footprint (SRF) value, and mean trees per hectare (TPH).

Discussion
Seasonal snowpack provides vital water resources for both human-and ecosystemoriented services in semi-arid environments. Quantifying the spatially heterogenous and often ephemeral nature of this snow at the mid-(10-1000 ha) and landscape-scales (400+ ha) requires accurate spatially extensive and temporally dense datasets [37]. The structure and arrangement of trees in dry forests are central components influencing snow cover and persistent snowpack. In this study, we demonstrate an accurate method to provide near real-time estimates of mid-scale snow covered area (SCA) and assess how forest structure influences persistent snow patch size in a thinned ponderosa pine forest on the southern edge of the North American continental snow distribution. We found that forest structure metrics emphasizing the spatial arrangement of trees and tree groups were more influential on persistent snow cover (Figure 7), and that these effects were most pronounced when considering trees within 1.5 times their height to persistent snow cover (Figure 8).
We first quantified snow-covered area (SCA) across the discontinuous forest of our 76-ha study site using UAV snow-series datasets. In forested environments, SCA is often underestimated in remote sensing data propagating from lower resolution satellite and airborne imagery, with trees masking the ground surface and tree shadows being misclassified as 'non-snow' [76][77][78][79]. The overall accuracy of our high-resolution binary snow classification was 90.2%, indicating strategic UAV data acquisition, associated highresolution imagery, and a relatively straightforward classification process can be used to quantify SCA in forested areas [80]. Furthermore, by capturing inter-storm reductions in SCA ( Figure 5), we demonstrate this approach can effectively delineate regions harboring the most persistent SCA and do so at a fine spatial resolution ( Figure 6). This level of detail and accuracy provided the foundation for our assessment of forest structure impact on persistent snow cover patches.
Except for crown diameter (CD), we found good overall agreement between the UAV SfM estimates and our validation datasets. Horizontal forest structure metrics were more influential than vertical structure metrics in predicting persistent snow patch size (Figure 7). Specifically, tree canopy cover (CC), trees per ha (TPH), and solar radiation footprint (SRF) were the most explanatory variables. This result supports our hypotheses that snow cover and, therefore, subcanopy shortwave radiation, is moderated by the spatial arrangement of trees more than by the vertical structure of individual trees.
While the horizontal metrics of CC, TPH, and SRF were the most influential predictor variables, they were also the metrics with the greatest variability when grouped by persistent snow patch size. These are metrics of foliar cover and tree stem density that directly influence subcanopy shading, which subsequently impacts the relatively shallow, intermittent, and spatially heterogeneous snow cover found at our study site. In addition to the horizontal metrics, we expected the mean CEI and mean KNN distance values to be influential because they describe the location and spacing of trees. However, their lack of influence might indicate a need to adjust the scale and scope at which these metrics are calculated. For example, the spacing of groups of trees might be a more valuable metric than the spacing of individual trees when considering sub-canopy ground shading.
In contrast to the horizontal metrics, variables in the vertical metric group were only slightly or moderately useful predictors of persistent snow patch size (Figure 7). We anticipated that forest patches comprised of trees with noticeably different vertical structures would result in differences in the forest canopy shading and ultimately statistically significant impacts on persistent snow patch size. This is illustrated in Figure 9, where the amount of shading seemingly provided by trees with higher CBH:CH ratio and lower CV values (Panel A) is noticeably lower compared to the amount of shading provided by trees with overall larger crowns (Panel B). The post-treatment forest structure conditions in Figure 9A are more common than those in Figure 9B. A potential cause for this could be a combination of the relatively homogenous vertical structure metrics observed across the study site ( Figure 1B) and the relatively low agreement (R 2 = 0.50; RMSE = 3.23 m) between the TLS and field-measured CBH values. study site ( Figure 1B) and the relatively low agreement (R 2 = 0.50; RMSE = 3.23 m) between the TLS and field-measured CBH values. Mean canopy height (CH) had strong agreement between both the field-measured and TLS-estimated values and the TLS-and SfM-estimated values, as well as low variability when grouped by persistent snow patch size. Accurate tree height estimates are not surprising given the ability for SfM and TLS to accurately estimate tree height regardless of forest density [34,39,41] and the relatively homogenous post-settlement forest structure conditions present across the study site. Given the low variable importance scores, mean CH appeared unimportant when modeling persistent snow patch area. However, we believe the importance of mean CH is evidenced in the superior predictive ability of the TSIA 1.5 category. More specifically, using tree CH to define the north-south width of forest canopy gaps (in the northern hemisphere), specifically to 1.5 times the average crown height of adjacent trees, provides the most consistent estimates of both snow cover and snow persistence. Early research shows that snow persisted in 'zones of retention' corresponding to 1.5 to 2 times the height of adjacent tree stands in ponderosa pine forests [81]. While these early studies assumed that persistent snow retention zones were driven by forest structure, mainly tree CH, we use high-resolution UAV-based estimates of CH to confirm the 1.5 CH value as well as to identify horizontal metrics contributing to persistent snow. Our results indicate that mean tree CH and tree patch spacing at 1.5 times CH are important considerations in forest restoration, in addition to the horizontal metrics described above, for maximizing snow cover retention and water yield.
As thinning-based restoration is expanded throughout dry forests that harbor seasonal snowpack, it is crucial to understand how restoration-driven canopy reduction will impact forest water balance [82] and snowmelt water inputs into shallow groundwater Mean canopy height (CH) had strong agreement between both the field-measured and TLS-estimated values and the TLS-and SfM-estimated values, as well as low variability when grouped by persistent snow patch size. Accurate tree height estimates are not surprising given the ability for SfM and TLS to accurately estimate tree height regardless of forest density [34,39,41] and the relatively homogenous post-settlement forest structure conditions present across the study site. Given the low variable importance scores, mean CH appeared unimportant when modeling persistent snow patch area. However, we believe the importance of mean CH is evidenced in the superior predictive ability of the TSIA 1.5 category. More specifically, using tree CH to define the north-south width of forest canopy gaps (in the northern hemisphere), specifically to 1.5 times the average crown height of adjacent trees, provides the most consistent estimates of both snow cover and snow persistence. Early research shows that snow persisted in 'zones of retention' corresponding to 1.5 to 2 times the height of adjacent tree stands in ponderosa pine forests [81]. While these early studies assumed that persistent snow retention zones were driven by forest structure, mainly tree CH, we use high-resolution UAV-based estimates of CH to confirm the 1.5 CH value as well as to identify horizontal metrics contributing to persistent snow. Our results indicate that mean tree CH and tree patch spacing at 1.5 times CH are important considerations in forest restoration, in addition to the horizontal metrics described above, for maximizing snow cover retention and water yield.
As thinning-based restoration is expanded throughout dry forests that harbor seasonal snowpack, it is crucial to understand how restoration-driven canopy reduction will impact forest water balance [82] and snowmelt water inputs into shallow groundwater reservoirs. Previous research shows that snow accumulation in discontinuous or disturbed forests can be greater in less dense forests and within large canopy openings [37,73]. The restoration thinning at our study site significantly reduced both canopy cover (from about 40% to 10%) and tree density (from about 212 TPH to 65 TPH), while increasing the number of forest patches (from 39 to 133) and decreasing the mean forest patch area (from 0.68 ha to 0.13 ha) [39]. While this study focused on developing remote sensing techniques to quantify snowpack dynamics, contrasting the treated versus the untreated portions of our study area suggest important treatment impacts that should be confirmed by more replicated comparisons. Specifically, our results show there was more overall persistent snow cover in the thinned portion of the study site compared to the untreated portion (10% and 7%, respectively). Our results generally support and provide further insight into restoration-based reductions in forest density to promote snow cover. For example, the largest persistent snow patches were located adjacent to forest patches with 31-33% canopy cover (CC), while the smallest were located adjacent to patches with 50-52% CC, a trend which is supported by the optimal snow persistence CC value of 24% provided by regional satellite-derived estimates [39].
Our findings suggest that tailoring future forest restoration treatments to promote persistent snow cover in southwestern U.S. dry forest ecosystems should continue to focus on horizontal forest structure metrics like forest density and canopy cover at the landscapescale (>1000 ha). Our results also underscore and refine the importance of horizontal forest structure metrics like CC, TPH, and SRF, regardless of spatial scale. However, we propose these criteria be emphasized during the creation of forest patches at relatively fine patch scale (<4 ha). It is critical to regional restoration efforts to continue implementing the commonly accepted standards in southwestern dry forest restoration, like promoting variation in interspace size, tree group size and within-group tree spacing [80]. In addition, we propose landscape-scale restoration treatment goals should also be implemented at fine spatial scale, operating within and among individual forest patches. For example, assuming a restoration treatment reflects diversity in forest patch size, spacing, and density, an overarching goal may be to achieve 24-33% CC across a 1000-ha treatment area. While this level of CC could be measured across the entire treatment area, having individual forest patches at 24-33% CC and distributed with distances at 1.5 times the average tree height within the patch should also promote localized persistent snow cover.

Conclusions
This study considered the inclusion of detailed forest structure metrics in quantifying persistent snow cover in a dry southwestern U.S. forest. We found the size (m 2 ) of persistent snow patches can be effectively predicted using targeted forest structure metrics. Specifically, the most effective predictor metrics included tree shadows that are 1.5 times the tree heights, as well as tree density and canopy cover within this shaded area. While our findings underscore the importance of forest canopy shading on persistent snow cover, they also indicate the relationship between persistent snow cover and fine-scale forest structure is likely more complex, rooted in different variables, or present at different spatial scales. Maximizing persistent snow cover in dry forest environments can be achieved by controlling subcanopy shading at the ground surface through an optimal set of fine-scale forest structure and spacing metrics. Future research and restoration efforts can achieve this by coupling our UAV-based methodology for quantifying persistent snow cover with more descriptive measurements of snow dynamics, such as snow depth and snow water equivalent.
Our results support the utility of thinning-based forest restoration in dry southwestern forests to promote snow cover retention and forest health. We show there is a wide range of persistent snow patch sizes across thinned forest, and that differences in fine-scale forest structure are important for maximizing snow persistence. Adjusting existing restoration thinning prescriptions to reflect landscape-scale goals in fine-scale forest patches will help further this objective while promoting broader ecosystem health and resiliency.