Snowpack relative permittivity and density derived from near‐coincident lidar and ground‐penetrating radar

Depth‐based and radar‐based remote sensing methods (e.g., lidar, synthetic aperture radar) are promising approaches for remotely measuring snow water equivalent (SWE) at high spatial resolution. These approaches require snow density estimates, obtained from in‐situ measurements or density models, to calculate SWE. However, in‐situ measurements are operationally limited, and few density models have seen extensive evaluation. Here, we combine near‐coincident, lidar‐measured snow depths with ground‐penetrating radar (GPR) two‐way travel times (twt) of snowpack thickness to derive >20 km of relative permittivity estimates from nine dry and two wet snow surveys at Grand Mesa, Cameron Pass, and Ranch Creek, Colorado. We tested three equations for converting dry snow relative permittivity to snow density and found the Kovacs et al. (1995) equation to yield the best comparison with in‐situ measurements (RMSE = 54 kg m−3). Variogram analyses revealed a 19 m median correlation length for relative permittivity and snow density in dry snow, which increased to >30 m in wet conditions. We compared derived densities with estimated densities from several empirical models, the Snow Data Assimilation System (SNODAS), and the physically based iSnobal model. Estimated and derived densities were combined with snow depths and twt to evaluate density model performance within SWE remote sensing methods. The Jonas et al. (2009) empirical model yielded the most accurate SWE from lidar snow depths (RMSE = 51 mm), whereas SNODAS yielded the most accurate SWE from GPR twt (RMSE = 41 mm). Densities from both models generated SWE estimates within ±10% of derived SWE when SWE averaged >400 mm, however, model uncertainty increased to >20% when SWE averaged <300 mm. The development and refinement of density models, particularly in lower SWE conditions, is a high priority to fully realize the potential of SWE remote sensing methods.

tested three equations for converting dry snow relative permittivity to snow density and found the Kovacs et al. (1995) equation to yield the best comparison with in-situ measurements (RMSE = 54 kg m À3 ).Variogram analyses revealed a 19 m median correlation length for relative permittivity and snow density in dry snow, which increased to >30 m in wet conditions.We compared derived densities with estimated densities from several empirical models, the Snow Data Assimilation System (SNODAS), and the physically based iSnobal model.Estimated and derived densities were combined with snow depths and twt to evaluate density model performance within SWE remote sensing methods.The Jonas et al. (2009) empirical model yielded the most accurate SWE from lidar snow depths (RMSE = 51 mm), whereas SNODAS yielded the most accurate SWE from GPR twt (RMSE = 41 mm).Densities from both models generated SWE estimates within ±10% of derived SWE when SWE averaged >400 mm, however, model uncertainty increased to >20% when SWE averaged <300 mm.The development and refinement of density models, particularly in lower SWE conditions, is a high priority to fully realize the potential of SWE remote sensing methods.

K E Y W O R D S
ground-penetrating radar, lidar, remote sensing, snow density, snow modeling, snow water equivalent

| INTRODUCTION
Seasonal snow covers up to 60% of the Northern Hemisphere land area (Hammond et al., 2018;Kim, 2018) and serves as a vital water resource for ecosystems spanning prairies, mountains, tundra, and boreal forests (Sturm et al., 1995).Mountains tend to accumulate deep snowpacks that provide water resources for one sixth of the world's population (Barnett et al., 2005;Mankin et al., 2015).In North America, mountains comprise 25% of the land area but store 60% of the total snow water equivalent (SWE; Wrzesien et al., 2018).Warming in the mountainous western United States (U.S.) has caused SWE losses of 15%-30% over the last 70 years (Mote et al., 2018), while a further 25% loss in SWE is predicted by 2050 (Siirila-Woodburn et al., 2021).These changes, compounded with human dependence upon snow water resources, make the pursuit of global SWE estimates a highly prioritized, trillion-dollar endeavour (National Academies of Sciences, Engineering, and Medicine, 2018;Sturm, 2015).Currently, no single method or ensemble of methods has proven capable of measuring SWE to the high standard of accuracy established for global monitoring (Dozier et al., 2016).Recent campaigns, such as the U.S.-based National Aeronautics and Space Agency (NASA) Snow Experiment (SnowEx; Durand et al., 2018) and the Europe-based NoSREx and APRESS (Tsang et al., 2022), evaluated a suite of remote sensing approaches (e.g., lidar, radar) for SWE-mapping applications.At the watershed scale, light detection and ranging (lidar) operations, such as the Airborne Snow Observatory (ASO), have demonstrated operational feasibility (Deems et al., 2013;Painter et al., 2016), but can be cost prohibitive.Satellite remote sensing methods for SWE-mapping have been under development for decades (Dietz et al., 2012;Nolin, 2010;Shi et al., 2016) and see continued interest, particularly in data assimilation applications, wherein derived snow products are integrated within physically based snow models (Largeron et al., 2020).For simplification, we discuss two major approaches of snow remote sensing: depth-based and radar-based methods.Both approaches have high spatial resolution and accuracy and are thus forerunners in the development of satellite-based SWE retrievals in mountains (Dozier et al., 2016).
Measuring SWE from depth-based remote sensing approaches relies on differencing repeat surface elevation measurements (i.e., snow-free from snow-on elevation surfaces) to derive snow depth (Currier et al., 2019;Deems et al., 2013).For this approach, SWE is the product of the derived snow depth (d s ) and snow density (ρ s ): Depth-based methods are primarily limited to the optical to near infrared portion of the electromagnetic spectrum, where atmospheric transmission is high, and wavelengths are significantly smaller than snow grains.Depth-based methods, including photogrammetry and lidar, have demonstrated potential for satellite platforms (Enderlin et al., 2022;Shean et al., 2016).
Radar-based methods further require an estimate of relative permittivity (ε s ) to characterize the electromagnetic wave velocity of the snowpack to derive snow depth, and thereby SWE, from the signal path length (l p ; e.g., two-way travel time; Marshall et al., 2005).Thus, radar-based SWE is generally expressed as a function, where the relative permittivity, a measure of the ability of a material to store charge relative to free-space (Daniels, 2004), is controlled by the snow density and volumetric liquid water content (LWC; θ LWC ).In dry snow, dielectric permittivity is primarily determined by snow density.However, wet snow permittivity exhibits a large dependence on LWC because the relative permittivity of water is $60 times that of snow (Bonnell et al., 2021).LWC induces the imaginary component of relative permittivity and acts as a frequency-dependent attenuator of signal strength (Koch et al., 2014), an attribute which was leveraged by Bradford et al. (2009) as the first study to derive spatially distributed LWC along a ground-penetrating radar (GPR) transect.At least six equations have been published for the relative permittivity of wet snow, whereas there are >19 published equations relating the relative permittivity of dry snow to its density, effectively increasing the uncertainty of radar-based SWE retrievals due to a large range in permittivity for a given snow density (Di Paolo et al., 2020).Not all radar SWE retrieval approaches utilize the signal path length through the snow.Some approaches, such as the Ka-band interferometry for wet snow conditions (Moller et al., 2017), may be better described as depth-based approaches, whereas Synthetic Aperture Radar (SAR) backscatter approaches rely on empirical models that are at least partly dependent on the snowpack relative permittivity (Tsang et al., 2022).
Density, a required input for both SWE approaches, can be either modelled or measured in snow pits, along snow courses, or by automated weather stations with a depth sensor co-located above a snow pillow.Density varies at the hillslope scale (Alford, 1967) due to differences in overburden pressure, radiation inputs, and wind and precipitation patterns (Winkler et al., 2021).However, density tends to vary less than snow depth, and because manual measurements are timeconsuming, studies generally assume limited spatial snow density variability (L opez-Moreno et al., 2013;Shook & Gray, 1996;Sturm et al., 2010).Thus, current density sampling approaches may acquire measurements too sparsely to capture its inherent spatial variability (Meehan, 2022), making snow density models an appealing alternative for remote sensing approaches.
Selecting an appropriate density model can be difficult.There are three broad categories: empirical models (Avanzi & De Michele, 2015), physically based models (e.g., Havens et al., 2020;Lehning, Bartelt, Brown, & Fierz, 2002;Lehning, Bartelt, Brown, Fierz, & Satyawali, 2002;Marks et al., 1999), and data assimilation systems (Largeron et al., 2020), and within each category, numerous models exist.Additionally, validation efforts outside of the original publications are limited.Empirical models are developed using a statistical relation between snow density, a time parameter (e.g., month), and, in more complicated equations, snow depth and geographic metrics.Physically based models can be powerful predictive tools but tend to be computationally expensive, require extensive meteorological forcing datasets, and often use simplified densification formulas (Hedrick et al., 2018).Empirical models can be simpler to implement and produce statistically robust estimates at the interannual scale, however few accurately capture short-term variability in snow density (McCreight & Small, 2014) Largeron et al., 2020).Recently published density models include a semi-empirical model that balances computational efficiency by implementing a simplified physically based snow compaction equation within an empirical framework (Winkler et al., 2021), and empirical models that estimate snow density from a suite of lidar-derived parameters (Bisset et al., 2022;Meehan, 2022).Both model types are optimized for depth-based remote sensing methods: the semiempirical model ingests repeat-measured snow depths acquired with low temporal baselines (≤7 days) and may be implemented for any given year/location, whereas lidar-derived models are designed and implemented for a single lidar survey (Meehan, 2022).
GPR measures the two-way travel time (twt) of the radar wave through the snowpack, which can be combined with lidar-measured snow depths to estimate radar velocity and relative permittivity.Several studies have established this method and converted relative permittivity to LWC by constraining snow density using snow pit measurements (Bonnell et al., 2021;Heilig et al., 2015;Webb et al., 2018;Webb, Wigmore, et al., 2020).More recent studies have used this technique to derive snow density by coupling uncrewed aerial vehicle (UAV) Structure from Motion (SfM) measurements of snow depth with GPR-measured twt and identified spatial variabilities that were larger than variabilities mapped by previous in-situ studies (McGrath et al., 2022;Yildiz et al., 2021).We expanded on previous work by leveraging a time series of GPR and lidar datasets from Grand Mesa and Cameron Pass, Colorado, and a one-off survey conducted at Ranch Creek, Colorado.We combined lidar-measured snow depths with GPR-measured twt of the snowpack thickness to derive spatially distributed relative permittivity in both dry and wet conditions.For dry-snow surveys only, relative permittivity is converted to density using three different equations to illustrate the variability in published equations.We then compared the derived densities with in-situ measurements and selected the most representative equation.Given the upcoming launches of L-band (1-2 GHz) SAR satellites (e.g., ALOS-4, NISAR, ROSE-L, and TanDEM-L) and their potential for global SWEmonitoring (Deeb et al., 2011;Guneriussen et al., 2001;Marshall et al., 2021;Tarricone et al., 2023), we estimated the spatial variability of derived relative permittivity and snow density at the approximate scale of SAR satellite platforms by conducting a variogram analysis.

| FIELD SITES
We used 1.  (Williams, 2021).This site included repeat 0.5 km GPR transects in 2020 and repeat 0.9 km transects in 2021 (Bonnell et al., 2022;McGrath et al., 2021).The nearby Joe Wright SNOTEL station (Figure S1b,e) recorded snow depth changes of +0.97 m (+280 mm SWE) and À0.17 m (+40 mm SWE) between the three WY 2020 surveys, and +0.10 m (+59 mm SWE), +0.41 m (+119 mm SWE), and À0.84 m (À104 mm SWE) between the four surveys conducted in WY 2021.Ranch Creek was surveyed on 7 April 2021 using a lidar system borne by a UAV (Bauer et al., 2023), while 2.9 km of GPR data were collected in a spiral survey design (Bonnell & McGrath, 2023).The Ranch Creek survey was conducted 10 days after peak SWE was observed by the nearby USGS weather station  Note: Water Year is abbreviated to WY. Grand Mesa vegetation summary sourced from Webb, Raleigh, et al. (2020).Cameron Pass and Ranch Creek vegetation notes were taken in the field and verified with Huckaby and Moir (1998) and Fassnacht et al. (2018).

| Calculating relative permittivity and snow density
Lidar-measured snow depths were collected from one of three A common practice of GPR in snow applications is to convert the twt of the snowpack thickness to snow depth (d s ; e.g., Lundberg et al., 2006;Marshall et al., 2005), using an estimate of the snowpack radar velocity (v s ; Daniels, 2004): Radar velocity is estimated from the relative permittivity of the snowpack, where, c is the velocity of electromagnetic radiation in a vacuum.
Here, we constrain the snow depth using lidar and directly estimate relative permittivity: In dry snow, snow density can be estimated directly from a relative permittivity equation.Because the choice of an equation is not straightforward (Di Paolo et al., 2020), we calculated snow density from three permittivity equations.Of the published equations and for a given relative permittivity, the Kovacs et al. (1995) equation (Equation 6) estimates the median density, the Kuroiwa (1954) equation (Equation 7) estimates the minimum density, and the Webb et al.
(2021) equation (Equation 8) estimates the maximum density (Di Paolo et al., 2020).The three equations are: The equations are written such that units for density (ρ s ; kg m À3 ) are consistent.
Measurements of twt were binned within lidar grid cells (3 m Â 3 m) by calculating the median twt value per cell.Grid cells that did not meet a minimum threshold of 15 twt measurements within the cell were removed from the analysis.Relative permittivity was calculated from coincident snow depth and twt cells (Equation 5).Previous studies have established a large randomly distributed error in the relative permittivity estimates that results from uncertainties in snow depths and twt, but with sufficient sampling and filtering, a robust estimate can be established (Bonnell et al., 2021;McGrath et al., 2022;Meehan, 2022).
Erroneous relative permittivity values (e.g., ε s < 1) were reduced by removing all values outside of the inter-quartile range.

| Uncertainty analysis
We estimated the uncertainty in relative permittivities through Monte Carlo simulations for each survey date.The mean snow depth and twt were calculated from the 3 m rasters.Grand Mesa lidar snow depth uncertainty was estimated from the extensive comparisons between the Grand Mesa airborne lidar and snow depth probes established by Currier et al. (2019).For Cameron Pass, the lidar snow depth uncertainty was estimated using the standard deviation of elevational differences between the bare-earth and snow-on DEMs along the CO-14 highway surface.The Ranch Creek lidar snow depth uncertainty was estimated from comparison with surveyed ground control points (Bauer et al., 2023).Uncertainties in twt were estimated from the mean withinpixel twt standard deviation for each survey date.Then, Monte Carlo simulations with 100 000 realizations from Equations ( 4) and (5) were performed using a random normal distribution, where the uncertainty estimates were considered representative of the standard deviation around the mean snow depth and twt.This established estimates for the mean and standard deviations of derived relative permittivity and snow density for each survey date, which is estimated as the uncertainty range in our derived relative permittivity and snow density datasets.All Monte Carlo simulation parameters and estimates are listed in Table S1.

| Variogram analysis
We conducted a variogram analysis using a lag spacing of 10 m, which approximates the spatial resolution of SAR satellites.Variance at the lag spacing, γ(h), for experimental variograms is given as, where, N(h) is the number of point-pairs at the given lag spacing and x is the variable of interest (Anderson et al., 2014;Schwanghart, 2022a;Webster & Oliver, 2001).Omni-directional experimental variograms were computed from the datasets for snow depth, twt, relative permittivity, snow density, and SWE (achieved by multiplying the derived snow density by its corresponding snow depth).Experimental variograms were computed using the same grid cells across all variables for the given survey date.Based on the shape of calculated experimental variograms, we identified the exponential model as the most representative variogram model for each of the variables.The exponential variogram model with a nugget effect is given as where, a is the correlation distance, s 0 is the nugget effect, and s is the exponential model contribution to the sill (Isaaks & Srivastava, 1989).
Variogram model parameters were estimated by least-squares fitting (Schwanghart, 2022b).For the Grand Mesa datasets, sufficient GPR

| Modeling snow density
We tested density estimates from four empirical models, SNODAS, and iSnobal against our derived density dataset.was chosen to represent data assimilation methods, given its use as a benchmark for evaluating larger scaled models (Broxton et al., 2016) and history of validation efforts (Clow et al., 2012;Hedrick et al., 2015;Lv et al., 2019).The iSnobal model (Havens et al., 2020;Marks et al., 1999) was chosen because of its incorporation into lidar SWE products geared toward operational water supply applications (Painter et al., 2016).iSnobal was run over Cameron Pass for both unscaled (iSn un ) and rescaled (iSn re ; e.g., Hale et al., 2023;Kiewiet et al., 2022;Voegeli et al., 2016) S3).Both the Kovacs et al. (1995) and Kuroiwa (1954) equations yielded densities with an overall negative bias, whereas the Webb et al. ( 2021) equation yielded densities with a positive bias (Table S3).However, the mean residual for the Kuroiwa (1954) S3), but frequently yielded physically unrealistic snow densities (e.g., <50 kg m À3 on 10 February 2021).Therefore, we selected the Kovacs et al. (1995) equation to use at both Grand Mesa and Cameron Pass because of its overall lower RMSE and more physically realistic snow densities.
We established survey-dependent uncertainty ranges for the derived relative permittivity and snow density datasets (Table S1) and found that the uncertainties in derived relative permittivity are more sensitive to snow depth than twt.At lower mean snow depths (<1 m), the effect of the snow depth uncertainty is increased, and large surveys.An overview of the mean and standard deviations for each of the datasets can be found in Table S4.
The lowest mean snow depth/twt at Grand Mesa was observed on forested transects on 8 February (1.19 m, 9.90 ns; Table S4 S5. Variogram analysis reveals longer correlation lengths for snow depth and twt at Grand Mesa in the forests than in the open, with forests having a median difference of +6 m for snow depth and +11 m for twt (Figure 4a-e; Table S5).However, relative permittivity, snow St 10 consistently overestimated snow density (mean difference = +37 kg m À3 ).S&F 14 performed comparably to J 09 at Cameron Pass.We found that the variability from empirically estimated snow densities was more limited compared to derived (Figure 5a), and that empirical model accuracy decreased in lower snow depth conditions.SNODAS densities We then used the derived, modelled, and the mean of in-situ den- The combination of GPR and lidar for the derivation of snow densities or LWC is a promising method that could be employed to supplement in-situ methods for both operational surveys and large-scale snow campaigns (e.g., NASA SnowEx).From our study, we can compare the combination of different lidar platforms with different GPR survey designs (e.g., TLS with a transect GPR survey) to suggest 'best practices' for future studies.Combined lidar-GPR surveys for deriving relative permittivity need to consider the large uncertainty ranges that result from low snow depths (<1 m) and high snow depth uncertainty (>0.10 m; McGrath et al., 2022).We found that a spatial filter reduced the uncertainty range and yielded physically realistic densities, particularly when surveys were performed as grids or spirals.Spiral and grid surveys have the added potential of spatial modeling to estimate relatively complete maps of relative permittivity and snow density in the surveyed region (e.g., Meehan, 2022).Airborne lidar platforms yielded lower snow depth uncertainties (±0.05 to ±0.08 m; Table S1) than the TLS (±0.07 to ±0.22 m), but the TLS snow depth uncertainty could be reduced substantially by surveying from a stable surface, rather than from a potentially shifting platform on the snow surface.Regardless of platform, we recommend that surveyors take careful notes regarding the GPR sled compression and the timing of the lidar collection relative to GPR surveys.

| Comparison with previous studies
The local-scale spatial variability of snow density is time consuming to measure using traditional in-situ sampling methods (L opez-Moreno et al., 2020).Therefore, snow density measurements tend to be sparser than other snow metrics, such as depth (Sturm et al., 2010), and thus density's spatial variability has seen few in-situ studies.We found that the average range in derived densities at Grand Mesa (135 kg m À3 ) and Cameron Pass (179 kg m À3 ; Figure 3, Figures S3-S5) compare favourably to the ranges in snow densities measured by insitu studies conducted in the Spanish Pyrenees, where snow depths had a similar range to those measured in our study, but warmer temperatures led to rainfall during the winter of one study year (Fassnacht et al., 2010;L opez-Moreno et al., 2013).Density ranges up to 300 kg m À3 were observed from four snow density measurement Pyrenees, densities ranged from 350 to 450 kg m À3 in January 2009 (Fassnacht et al., 2010).The ranges reported by L opez-Moreno et al.
(2013) fully encompass the maximum range observed in our derived density datasets, whereas Fassnacht et al. (2010) observed ranges more comparable to those we derived at Grand Mesa.It can be noted that, depending on the dominant processes, snowpack density can vary spatially at the 0.10 m scale (Fassnacht, 2021).
Several GPR-based methodologies have been developed to better understand the spatial distribution of snow densities and have increased the number of density estimates by two to three orders of magnitude compared to in-situ studies.However, GPR-based methodologies have generated snow density ranges that are physically unrealistic, especially when compared to in-situ measurements.Velocity migration analysis has derived densities with a large range (100-500 kg m À3 ) along relatively short ($100 m) transects in Wyoming's Medicine Bow Mountains (St. Clair & Holbrook, 2017), whereas multi-antenna pair GPR systems have derived densities as low as 150 kg m À3 during the melt season near Davos, Switzerland (Griessinger et al., 2018).Still, the multi-antenna pair methodology is promising; densities tended to range <80 kg m À3 for single $100 m transects for Griessinger et al. (2018), whereas another multi-antenna pair GPR system was used to derive a narrow range (370-420 kg m À3 ) along a 78 km transect across the Greenland Ice Sheet (Meehan et al., 2021).Our approach differs in that it integrates independently  S4) by a factor of two, relative to Yildiz et al. (2021) and McGrath et al. (2022).Although the overall RMSE for our study is still high, the RMSE at Grand Mesa (45 kg m À3 ) is lower than either of the UAV SfM/GPR studies (68-74 kg m À3 ) and derived density ranges are physically realistic.Thus, as this method is further developed and tested, it may be used to provide unprecedented detail on the spatiotemporal distribution of snow densities.

| Spatial variability of snowpack parameters
Previous studies have shown that increases in average LWC are associated with larger LWC spatial variability and that radar-based methods for measuring SWE will have increased uncertainty when LWC is present in the snowpack (Bonnell et al., 2021;Webb et al., 2022).Although we do not estimate LWC, we identify a large shift between dry and wet relative permittivity sills (Table S5), indicating that relative permittivity variance is higher for wet snow surveys.forests exhibited cyclicity at the sill level (Figure S7), possibly due to the snow accumulation dynamics in the forest islands of Grand Mesa (Webb, Raleigh, et al., 2020).Sills for derived relative permittivity and density at Cameron Pass (Figure 4, Figures S8 and S9) tended to be high in December and February and decrease with accumulation, where higher sills may result from snowpack variability induced by topography and vegetation before the features are blanketed by snow.

| Model performance and selection
Of the empirical models, J 09 estimated densities and depth-based SWE with the lowest overall RMSE.St 10 showed a consistent, positive bias (Figure 5), which was also noted by McCreight and Small (2014), largely due to the DOY parameter, which is less sensitive to shortterm fluctuations in snow density than the month parameter used by J 09 .We found that the DOY-based M&P 08 systematically underestimated snow density and SWE at our sites, which could be improved for sites where calibration is an option (Pistocchi, 2016).Each evaluated model has known limitations.Empirical models are limited to a narrow, unrealistic range of density estimates (Raleigh & Small, 2017) and have been shown to decrease in accuracy with increased elevation (Avanzi & De Michele, 2015).Evaluations of SNODAS indicate that modelled SWE and snow depth have lower errors in subalpine regions (Clow et al., 2012;Dozier et al., 2016), but SNODAS tends to underestimate accumulation in deeper snowpacks (Hedrick et al., 2015) and may be unsuitable in alpine regions, where snow deposition and redistribution processes are more complex (Dozier et al., 2016)  Pass, correlation lengths decreased by 10 m over the span of 1.5 months.We observed that relative permittivity has a higher variance in wet snow than in dry snow, but correlation lengths were longer for the wet snow surveys.We compared derived densities with densities estimated from four empirical models, SNODAS, and iSnobal and found that J 09 was the most accurate density model for depth-based SWE retrievals (±10% for seven of nine surveys), whereas SNODAS densities yielded the most accurate SWE retrievals for radar-based methods (±10% for eight of nine surveys).However, for lower SWE environments (<300 mm) uncertainty increased to >20% for all models, which points to a potential issue in density modeling.Regardless, selecting an appropriate density model is difficult, and new models continue to be developed (e.g., Bisset et al., 2022;Meehan, 2022;Winkler et al., 2021).Empirical models evaluated in this study are directly applicable for depth-based SWE remote sensing methods and the Sentinel-1 co/cross-polarization backscatter method (Lievens et al., 2019(Lievens et al., , 2022)), but site-specific calibration in both open and forested environments may improve results.
Physically based models may be particularly powerful for InSAR methods of SWE retrieval, which can be applied over expansive regions (e.g., Oveisgharan et al., 2023) and require a surface density estimate to derive a change in SWE (Marshall et al., 2021), but these

(
29 March; Figure S1c,f), which had lost 0.31 m snow depth (135 mm SWE) over that same time period.F I G U R E 1 (a) Location of field sites within Colorado and 1 January to 1 July snow persistence (Moore et al., 2015).(b) Digital elevation model (DEM) and (c) National Agricultural Imagery Program (NAIP imagery of Cameron Pass, with snow pit locations.(d) DEM (Bauer et al., 2023) and (e) NAIP imagery of Ranch Creek, with weather station and snow pit locations.(f) DEM (Painter & Bormann, 2020) and NAIP imagery of Grand Mesa, with locations of SNOTEL stations and snow pits.Figure scales and colour ramps differ by field site.The vertical datums are the WGS84 Ellipsoid for Grand Mesa and the NAVD88 Geoid 18 for Cameron Pass and Ranch Creek.NAIP imagery acquired from USGS Earth Explorer (https://earthexplorer.usgs.gov/,accessed 10 October 2022).T A B L E 1 Technical details of the instrumentation and datasets for each of the field sites in Colorado.
platforms: airborne (Grand Mesa), terrestrial (Cameron Pass), and UAV-borne (Ranch Creek).GPR radargrams were collected as common-offset surveys via a sled coupled to the snow surface and utilized L-band centre-frequency.Processing of lidar point clouds collected by terrestrial and UAV platforms generally followed the workflow outlined by Currier et al. (2019), whereas radargram processing and picking of snowpack twt thickness followed McGrath et al. (2019).Appendix A.1 provides further details of the data processing.
Uncertainties were further reduced by smoothing the remaining relative permittivity estimates with a 21 m Â 21 m moving window median filter.A filter of this size was chosen to retain spatial variability along the transect-oriented surveys.We categorized surveys as dry or wet based on the presence of any LWC noted in snow pits and corroborated by pit temperatures (Appendix A.2; FigureS2).Relative permittivity estimates obtained in drysnow conditions (all surveys, except the 7 April and 27 May 2021 surveys) were then converted to density using theKovacs et al. (1995),Kuroiwa (1954), andWebb et al. (2021) equations and compared to in-situ density measurements to calculate the RMSE and Pearson correlation coefficient to determine the most representative equation.
observations were collected in forests to enable variogram analyses of snowpacks in both forests and open environments.Following McGrath et al. (2019), we chose a 2 m vegetation height metric and used the ASO vegetation heights dataset collected in Fall 2016 (Painter & Bormann, 2020) to generate a binary forest/open mask.Variograms were not calculated for snow density or SWE when LWC was present in the snowpack (i.e., 7 April and 27 May 2021).

|
Overview of the derived relative permittivity and snow density datasetsWe evaluated the lidar-GPR derived snow densities from the three relative permittivity equations using 40 snow pits from Grand Mesa (median = 322 kg m À3 , range = 93 kg m À3 ) and six snow pits from Cameron Pass (median = 255 kg m À3 , range = 92 kg m À3 ).Ranch Creek was excluded from this analysis because of wet snow conditions.Grand Mesa snow pits had an average of 19 derived snow densities within 30 m of each pit, whereas Cameron Pass averaged 11.When compared to snow pit measurements, theKovacs et al. (1995),Kuroiwa (1954), andWebb et al. (2021) derived densities yielded overall RMSEs of 54 kg m À3 (r = 0.09), 97 kg m À3 (r = 0.07), and 83 kg m À3 (r = 0.08), respectively (Figure2; Table derived densities was three times the mean residual of theKovacs et al. (1995) derived densities.At Cameron Pass, theKuroiwa (1954) RMSE was 25% lower thanKovacs et al. (1995; Table relative permittivity uncertainties result, leading to unrealistic estimates (e.g., ε s < 1 and ε s > 88).The average uncertainty for densities derived from TLS (196 kg m À3 ) is approximately twice the average for those derived from airborne lidar (88 kg m À3 ).This is primarily due to the lower average snow depths observed at Cameron Pass and the larger snow depth uncertainty range for the TLS platform.Differences between GPR systems did not have a substantial effect on estimated twt uncertainty, likely because of the similar vertical resolutions (calculated as one fourth of the wavelength;Daniels, 2004) of the 1.6 GHz ($0.04 m) and 1.0 GHz ($0.06 m) systems.A subset of lidar snow depths, twt, derived relative permittivities, and derived snow densities are shown for Grand Mesa, Cameron Pass, and Ranch Creek in Figure 3.Although the full spatial extent of lidar and twt datasets are displayed, only snow depths that were used to derive relative permittivity values are discussed.All snow depth, twt, derived relative permittivity, and derived density datasets are plotted in Figures S3-S6.Results for Cameron Pass 2020 and 2021 are discussed as a time series.Although the Grand Mesa lidar was collected in sequential dates, Grand Mesa results are not treated as a time series because each week had a different region of interest for GPR et al. (1995), and Webb et al. (2021).(b) Median derived snow density (considering all values within 30 m of the snow pit) versus measured snow density.Snow densities were derived using the three relations (Equations 6-8) presented in panel (a).F I G U R E 3 Examples from each field site and study year of the datasets used in this study.Rows are organized by field site, with (a-d) the 16 February 2017 Grand Mesa survey (GM 2017), (e-h) the 12 March 2020 Cameron Pass survey (CP 2020), (i-l) the 22 March 2021 Cameron Pass survey (CP 2021), and (m-o) the 7 April 2021 Ranch Creek survey (RC 2021).Columns are organized from left to right as snow depth, twt, relative permittivity, and snow density.
density, and SWE have longer correlation lengths in the open than in the forests, with open environments having a median difference of +7 m for relative permittivity and snow density and +3 m for SWE.At Cameron Pass, no parameters revealed an obvious relation with time for WY 2020, but relative permittivity and snow density correlation lengths for WY 2021 decreased from 24 m on 10 February to 10 m on 22 March.Subsequently, the correlation length for relative permittivity increased substantially to 55 m during the wet snow survey on 27 May.A longer correlation length was also observed for relative permittivity on 7 April at Ranch Creek (34 m).These are two of the longest relative permittivity correlation lengths found in the study, indicating that wet snow may have longer correlation lengths for relative permittivity than dry snow.4.3 | Evaluation of density modelsSpatially distributed snow densities at Grand Mesa and Cameron Pass enable an evaluation of snow density models across a range of dates and snow conditions (Figure5a-c).Here, we use the derived densities as a benchmark for comparison with modelled densities.A comparison between modelled densities and in-situ densities is available in Section-A.4.Of the empirical models, J 09 yielded the lowest RMSE at both Grand Mesa (RMSE = 16 kg m À3 ) and Cameron Pass (RMSE = 43 kg m À3 ), whereas M&P 08 yielded large RMSEs at Cameron Pass (RMSE = 87 kg m À3 ) and Grand Mesa (RMSE = 71 kg m À3 ), and systematically underestimated snow density (mean difference = À75 kg m À3 ).
U R E 4 Exponential variogram models calculated from the best fit of the experimental variograms.Rows represent different field sites and are organized from top to bottom as (a-e) Grand Mesa, (f-j) Cameron Pass, and (k-m) Ranch Creek.Columns represent different variables and are organized from left to right as snow depth, twt, relative permittivity, snow density, and SWE.U R E 5 (a-c) Comparison between mean snow densities, (d-f) mean depth-based SWE, and (g-i) mean radar-based SWE estimated from the derived, in-situ, and modelled densities.Depth-based SWE was calculated using lidar snow depths, whereas radar-based SWE was calculated using twt.Columns are organized by field site, from left to right: Grand Mesa surveys from WY 2017, Cameron Pass surveys from WY 2020, and Cameron Pass surveys from WY 2021.Error bars show ±1 standard deviation.Depth-based empirical density models cannot be incorporated into radar-based SWE calculations and are thus absent in (g-i).at Grand Mesa yielded an RMSE (20 kg m À3 ) that was comparable to J 09 , whereas SNODAS densities yielded the lowest RMSE (32 kg m À3 ) at Cameron Pass.Unscaled iSnobal densities for Cameron Pass in WY 2020 yielded the highest RMSE (119 kg m À3 ), whereas rescaled iSnobal densities improved the RMSE by 20 kg m À3 .iSnobal density accuracy improved for the March surveys, where RMSEs were 20%-60% lower than the overall RMSE.For both unscaled and rescaled iSnobal densities, estimates exhibited a substantial negative bias (mean residuals: iSn un = À85 kg m À3 , iSn re = À67 kg m À3 ).
sities to calculate depth-based SWE from the lidar-measured snow depths (Figure 5d-f) and radar-based SWE from the GPR-measured twt (Figure 5g-i).We considered the SWE calculated from derived densities as the benchmark for this analysis.J 09 yielded the best performance for depth-based SWE estimates, estimating SWE within 10% for seven of nine surveys, whereas unscaled iSnobal yielded the largest errors, estimating depth-based and radar-based SWE within 20% for only two surveys.M&P 08 and rescaled iSnobal performed comparably: M&P 08 yielded two depth-based SWE surveys within 20% and five radar-based surveys within 20%, whereas rescaled iSnobal estimated depth-based and radar-based SWE within 20% for four surveys.M&P 08 did not produce SWE estimates within 10% for any survey, but rescaled iSnobal improved density estimates by an average of 7% and estimated depth-based and radar-based SWE within 10% for two surveys.St 10 estimated depth-based SWE within 10% for three of nine surveys, but yielded estimates within 20% for an additional five surveys.S&F 14 performed comparably to J 09 at Cameron Pass and estimated depth-based SWE within 10% for four surveys.SNODAS yielded the best overall performance for radar-based SWE, estimating SWE within 10% for eight of nine surveys, although its performance for depth-based SWE was reduced to four surveys within 10%.We found that the average in-situ density measurement yielded SWE within 10% for four of nine depth-based surveys and five of nine radar-based surveys.5 | DISCUSSION 5.1 | Considerations for lidar-GPR surveys deriving relative permittivity and snow density campaigns conducted during accumulation (February 2010 and 2011) and around peak SWE (April 2010 and 2011) across three different locations ranging in elevation from 1517 to 3015 m (L opez-Moreno et al., 2013).Along a 5.4 km section of the Rio Esera in the Spanish measured snow depths with twt to derive permittivity and density, similar to Yildiz et al. (2021), who derived densities by combining UAV SfM snow depths with GPR in the Ilgaz Mountains, Turkey, and McGrath et al. (2022), who conducted a time series of UAV SfM snow depth/GPR surveys along a 150 m transect at a field site adjacent to Cameron Pass during the 2021 study year.However, our use of a spatial filter reduced the derived snow density standard deviation (21-88 kg m À3 ; Table However, we found increased correlation lengths during wet snow surveys, which seems counterintuitive, given the complicated spatial and temporal variability of LWC.The increased correlation length indicates that modeling bulk volumetric LWC may be possible for field sites that have little topographic complexity, such as Cameron Pass, and could reduce the uncertainty of SWE retrievals from radar-based methods.Median correlation lengths were nearly identical for twt (27 m) and snow depth (30 m) and were lower than probed snow depth correlation lengths measured by Anderson et al. (2014) in southwest Idaho (median correlation length = 46 m).At Grand Mesa, median correlation lengths for snow depth and twt were longer in the forests than in the open, corroborating the findings of McGrath et al. (2019), who also found longer correlation lengths for probed depths and GPR-derived snow depths in the forest, but our mean correlation lengths in the open were nearly four times longer.Snow density and relative permittivity variograms yielded a median correlation length of 19 m, which is about twice that observed by Yildiz et al. (2021) but is much smaller than what L opez-Moreno et al. (2013) observed (<150 m).Experimental variograms of derived densities in Grand Mesa J 09 estimated depth-based SWE within 10% for seven of the nine surveys and may represent a viable option for depth-based SWE remote sensing methods.Of the more complex models, SNODAS had a small negative bias across all surveys at Grand Mesa (overall mean residual = À14 kg m À3 ) but exhibited no consistent bias for the Cameron Pass surveys, whereas unscaled and rescaled iSnobal densities were negatively biased for all Cameron Pass surveys (Figure5).However, rescaling precipitation yielded improved iSnobal density estimates that exhibited lower RMSE (<50 kg m À3 ) for WY 2020 surveys with deeper snowpacks (>1.30m).Of all models, SNODAS yielded the most consistently accurate radar-based SWE (eight of nine surveys within 10%) and the second most consistently accurate depth-based SWE (nine of nine surveys within 15%), thus SNODAS may offer a suitable solution for snow density estimates in North America for both depth-based and radar-based methods.
By combining GPR-measured twt with lidar-measured snow depths, we derived relative permittivity and snow density estimates along GPR transects at three sites during 11 surveys.For dry-snow surveys, this method yielded a mean RMSE of 54 kg m À3 (RMSE minimum = 45 kg m À3 , maximum = 92 kg m À3 ) compared with nearby in-situ density measurements.At Grand Mesa, we observed that median correlation lengths for relative permittivity and density were 6 m longer in the open than in the forests, whereas at Cameron models are inherently limited by the chosen forcing data.With the increased availability of depth-based and radar-based remote sensing datasets, accurate density estimates are required to advance SWE remote sensing techniques.The accuracy and spatial coverage capabilities of the lidar-GPR snow density method makes it a promising avenue for model evaluation, training, and development.
(Meyer et al., 2023)ios because the chosen atmospheric model used for meteorological forcing within iSnobal, the High Resolution Rapid Refresh (HRRR) model, has been shown to underreport SWE by up to 25%(Meyer et al., 2023).Only surveys with dry snow conditions were modelled.An overview of the model parameters is given in TableS2 and furtherdetails regarding model runs are provided in Appendix A.3.All modelled densities were compared with derived densities and the average of the in-situ measurements.Modelled densities were incorporated with lidar snow depths and GPR twt to calculate depthbased SWE and radar-based SWE.In-situ measurements were averaged for each date and used to calculate depth-based SWE from the lidar snow depths and radar-based SWE from twt to evaluate the efficacy of snow models compared to the 'gold standard'.The RMSE was then calculated for modelled versus derived densities, modelled and in-situ versus derived depth-based SWE, and modelled and in-situ versus derived radar-based SWE.
. Accurate forcing data are necessary for robust iSnobal model output, but the HRRR precipitation forcing data used in our study was consistently lower than Joe Wright SNOTEL SWE for WY 2020 and 2021 (Figure S11).For WY 2021, HRRR precipitation underestimated SWE by >50%, which is a larger error than what Meyer et al. (2023) reported, and caused the substantially negative bias observed for estimated densities.Rescaling HRRR precipitation based on the 19 March 2021 Cameron Pass airborne lidar survey yielded more realistic spatial patterns for snow deposition but failed to generate accurate precipitation values at our field site, hence our introduced rescaling factor.Despite the higher RMSE values and biases, iSnobal remains a promising model for estimating snow densities, particularly where weather station coverage is sufficient or whereHRRR precipitation is either accurate or can be corrected.As a final note, we found that lower SWE conditions (<300 mm) resulted in larger percent errors from models, increasing the J 09 and SNODAS SWE errors to 18% and 13%, respectively.St 10 and S&F 14 yielded SWE errors as high as 23%, whereas rescaled iSnobal SWE errors exceeded 30%.Thus, further evaluation of modelled densities in lower snow conditions is warranted.