Mapping snow depth within a tundra ecosystem using multiscale observations and Bayesian methods

. This paper compares and integrates different strategies to characterize the variability of end-of-winter snow depth and its relationship to topography in ice-wedge polygon tundra of Arctic Alaska. Snow depth was measured using in situ snow depth probes and estimated using ground-penetrating radar (GPR) surveys and the photogrammetric detection and ranging (phodar) technique with an unmanned aerial system (UAS). We found that GPR data provided high-precision estimates of snow depth (RMSE = 2.9 cm), with a spatial sampling of 10 cm along transects. Phodar-based approaches provided snow depth estimates in a less laborious manner compared to GPR and probing, while yielding a high precision (RMSE = 6.0 cm) and a ﬁne spatial sampling (4cm × 4cm). We then investigated the spatial variability of snow depth and its correlation to micro-and macrotopogra-phy using the snow-free lidar digital elevation map (DEM) and the wavelet approach. We found that the end-of-winter snow depth was highly variable over short (several meter) distances, and the variability was correlated with microto-pography. Microtopographic lows (i.e., troughs and centers of low-centered polygons) were ﬁlled in with snow, which resulted in a smooth and even snow surface following macro-topography. We developed and implemented a Bayesian approach to integrate the snow-free lidar DEM and multiscale measurements (probe and GPR) as well as the topographic correlation for estimating snow depth over the landscape. Our approach led to high-precision estimates of snow depth (RMSE = 6.0 cm), at 0.5 m resolution and over the lidar domain (750m × 700m).


Introduction
Snow plays a critical role in ecosystem functioning of the Arctic tundra environment through its impacts on soil hydrothermal processes and energy exchange (e.g., Callaghan et al., 2011).Snow insulates the ground from intense cold during the Arctic winter, limiting the heat transfer between the air and the ground (Zhang, 2005).Snow depth affects active layer and permafrost temperatures throughout the year (Gamon et al., 2012;Stieglitz et al., 2003), and increased snow depth has resulted in permafrost degradation (Osterkamp, 2007).Snow's insulating capacity enhances conditions for active soil microbial processes and CO 2 /CH 4 production during the winter (Nobrega and Grogan, 2007;Schimel et al., 2004;Clein and Schimel, 1995;Jansson and Taş, 2014;Zona et al., 2016).In addition, snow serves as an important water source to tundra ecosystems during the growing season and therefore has a large impact on biological processes via hydrology.Snowmelt water can lead to extensive inundation of low-gradient tundra and large runoff events in early summer (Bowling et al., 2003;Kane et al., 1991;Liljedahl et al., 2016).Since soil biogeochemistry and vegetation are controlled by soil moisture (Sjögersten et al., 2006;Wainwright et al., 2015), the amount of snow affects ecosystem functioning throughout the season.
In order to investigate controls of snow on ecosystem properties, high-resolution estimates of snow are needed over large spatial regions.This is especially true in icewedge polygon tundra, which dominates a large portion of Published by Copernicus Publications on behalf of the European Geosciences Union.
the high Arctic (Zona et al., 2011).The ice wedges develop when frost cracks occur in the ground, and vertical ice wedges grow laterally over years (Leffingwell, 1915;MacKay, 2000).Soil movement associated with ice-wedge development creates small-scale topographic variations -microtopography -where the ground surface elevation can vary significantly over lateral length distances of several meters (e.g., Brown, 1967;MacKay, 2000;Engstrom et al., 2005;Zona et al., 2011).This microtopography leads to dramatically variable snow depth across short distances.Liljedahl et al. (2016) found that the differential snow distribution increased the partitioning of snowmelt water into runoff, leading to less water stored on the tundra landscape.Gamon et al. (2012) reported that snow depth heterogeneity results in differential thawing and active layer thickness variability.In addition, there is large-scale topographic variability at the scale of several hundred meters to kilometers -macrotopography -which is often associated with drained thaw lake basins or drainage features (Hinkel et al., 2003).Although the effect of macrotopography on snow depth has not been studied, Engstrom et al. (2005) quantified that both macrotopography and microtopography have a significant effect on soil moisture distribution.The snow representation of the Arctic tundra needs to be refined to account for the effect of such multiscale terrain heterogeneities on hydrology and ecosystem functioning by bridging between finer geographical scales (several meters) and large areal coverage (several hundred meters to kilometers).
Snow depth characterization in Arctic tundra environments has traditionally been performed using snow depth probes (Benson and Sturm, 1993;Hirashima et al., 2004;Derksen et al., 2009;Rees et al., 2014;Dvornikov et al., 2015) or modeled using terrain and vegetation information (Sturm and Wagner, 2010;Liston and Sturm, 1998;Pomeroy et al., 1997).Recently, there have been several new techniques for estimating snow depth in high resolution and in a noninvasive and spatially extensive manner.Groundpenetrating radar (GPR) has been widely used to characterize snow cover in alpine, arctic and glacier environments (e.g., Harper and Bradford, 2003;Machguth et al., 2006;Gusmeroli and Grosse, 2012;Gusmeroli et al., 2014).GPR measures the radar reflection from the snow-ground interface, which can be used to estimate snow depth.GPR can be collected by foot, snowmobile or airborne methods.In addition, light detection and ranging (lidar) and photogrammetric detection and ranging (phodar) airborne methods have recently been used to estimate snow depth at local and regional scales (e.g., Deems et al., 2013;Harpold et al., 2014;Nolan et al., 2015).Both techniques measure the snow surface elevation, using laser in lidar or a camera with a structurefrom-motion (SfM) algorithm in phodar.Both approaches allow us to estimate snow depth by subtracting the snowfree elevation from the snow surface elevation.While there is potential for providing detailed information about localscale snow variability using lidar and phodar snow depth es-timates, these techniques have not been extensively tested in ice-wedge polygonal tundra environments.
Such indirect geophysical methods are, however, known to have increased snow depth uncertainty relative to direct measurements (here ground-based snow depth probe measurements) (e.g., Hubbard and Rubin, 2005).The uncertainty of the snow depth probe measurements is sub-centimeter to several centimeters depending on the surface vegetation (Berezovskaya and Kane, 2007).In contrast, the snow depth estimates obtained using GPR can be affected by uncertainty associated with radar velocity, which depends on snow density (Harper and Bradford, 2003).In the environments with complex terrain such as ice-wedge polygonal tundra, GPRbased snow estimates could also be influenced by the errors stemming from radar positioning and ray path assumptions.The airborne lidar/phodar-based methods are subject to the errors associated with georeferencing, processing and calibration (e.g., Deems et al., 2013;Nolan et al., 2015).The accuracy of the airborne methods is usually several tens of centimeters, which is lower than the snow depth probe measurements.
Integrating different types of snow measurements can take advantage of the strengths of various techniques while minimizing the limitations stemming from using a single method.Bayesian approaches have proven to be useful for integrating multiscale, multi-type datasets to estimate spatially heterogeneous terrestrial system parameters in a manner that honors method-specific uncertainty (e.g., Wikle et al., 2001;Wainwright et al., 2014Wainwright et al., , 2016)).Bayesian methods also permit systematic incorporation of expert knowledge or process-specific information, such as the relationships between datasets and parameters.In particular, snow depth is known to be affected by topography and wind direction (e.g., Benson and Sturm, 1993;Anderson et al., 2014;Dvornikov et al., 2015).To our knowledge, such Bayesian data integration methods have never been applied to estimate end-ofwinter snow variability using multiple types of datasets.
The primary objectives of this study are to (1) compare point-scale snow depth probe, GPR and UAS-based phodar approaches for characterizing snow depth and the associated resolution and accuracy of the GPR and phodar methods; (2) quantify the spatial variability of end-of-winter snow depth in ice-wedge polygonal tundra landscape; (3) explore the relationship between snow depth and topography; and (4) develop a Bayesian method to integrate multiscale, multitype data to estimate snow depth over a lidar digital elevation model (DEM) covering an ice-wedge polygonal tundra landscape.In Sect.2, we describe our site and datasets, including snow depth probes, ground-based GPR and UASbased phodar.In Sect.3, we present the methodology to analyze the indirect snow depth measurements from GPR and phodar as well as to evaluate the heterogeneity of snow depth in relation to both microtopography (i.e., ice-wedge polygons) and macrotopography (i.e., large-scale gradient, drained thaw lake basins and interstitial upland tundra).We The Cryosphere, 11, 857-875, 2017 www.the-cryosphere.net/11/857/2017/then develop a Bayesian geostatistical approach to integrate the multiscale datasets to estimate snow depth over the lidar domain.The snow measurement and estimation results are presented in Sect. 4 and discussed in Sect. 5.
2 Data and site descriptions

Study site
Snow survey data were collected within a study site (approximately 750 m × 700 m) located on the Barrow Environmental Observatory near Barrow, Alaska, as part of the Department of Energy's Next-Generation Ecosystem Experiment (NGEE) Arctic project (Fig. 1).This study domain has been characterized intensively in the NGEE Arctic project, leading to various ecosystem and subsurface datasets, including snow depth measurements (Wainwright et al., 2015;Dafflon et al., 2016).Mean annual air temperature at the Barrow site is −11.3 • C and mean annual precipitation is 173 mm (Liljedahl et al., 2011).Snowmelt usually ends in early to mid-June.The wind direction is predominantly from east to west throughout the year.Ice-wedge polygons are prevalent in the region, including low-centered polygons in drained thaw lake basins and highcentered polygons with well-developed troughs in the upland tundra (Hinkel et al., 2003;Wainwright et al., 2015).The dominant plants are mosses (Dicranum elongatum, Sphagnum), lichens and vascular plants (such as Carex aquatilis); plant distribution at the site is governed by surface moisture variability (e.g., Hinkel et al., 2003;Zona et al., 2011).There are currently no tall shrubs or woody plants established within the study site, and therefore complex topography is most likely to control the snow depth distribution within the study domain (Sturm et al., 2005;Dvornikov et al., 2015).
Three long transects and four representative plots were chosen within the study site to explore snow variability and its relationship to topography (Fig. 1).Typical for lowgradient tundra terrain, ice-wedge polygons and microtopographic variations are superimposed on macrotopographic trends at the study site.The elevation is higher in the center of the domain (interstitial upland tundra) and lower near the drainage features in the south.The elevation is also relatively lower in the drained thaw lake basins (DTLB) region, which is located in the northeastern and northwestern edges of the study site.The four intensive plots (A-D), each 160 m × 160 m, were chosen to represent specific polygon types or macrotopographic positions within the study area.The three parallel transects, each ∼ 500m long, were designed to traverse multiple polygon types in a continuous fashion (Hubbard et al., 2013;Wainwright et al., 2015).We refer to those transects by "the 500 m transects".

Datasets
Airborne lidar data were collected at the site on 4 October 2005 and used to provide a high-resolution DEM of the snow-free ground at 0.5 m × 0.5 m resolution (Hubbard et al., 2013).The DEM effectively resolves both microand macrotopography at the study site (Fig. 1).The original reported accuracy is 0.3 m in the horizontal direction and 0.15 m in the vertical direction.To further evaluate the accuracy of the airborne DEM, we measured the ground surface elevation in September 2011 at 1286 points around the 500 m transects, using a high-precision centimeter-grade RTK differential GPS (DGPS) system (the reported precision about 2 cm in the horizontal direction and 3 cm in the vertical direction).The root mean square error (RMSE) of the lidar DEM compared to the DGPS data was 6.08 cm.
The majority of the snow depth data were collected on 6-12 May 2012, during which no snowfall occurred and little change in snow depth was observed.Snow depth was measured in the four intensive study plots and along three transect lines (Fig. 1).Two sets of snow depth measurewww.the-cryosphere.net/11/857/2017/The Cryosphere, 11, 857-875, 2017 ments using a snow depth probe were collected.The "finegrid" dataset was aimed to characterize the fine-scale heterogeneity by ∼ 7200 snow depth point measurements (every ∼ 0.3 m along transects with a 4 m spacing) across a small domain (∼ 50 m × 50 m) within plots A-D (Liljedahl et al., 2014).This was done using a GPS snow depth probe (Magnaprobe by Snow-Hydro) which had a reported vertical precision of < 0.01 m and horizontal precision of 2-10 m.The corner coordinates within each grid were surveyed with the RTK DGPS, while each snow depth point measurement was associated with latitude and longitude positional information recorded by the Magnaprobe's built-in GPS receiver.All the snow depth point measurements were made along regularly spaced transects.Comparisons between coordinates surveyed with both the RTK DGPS and the Magnaprobe's built-in GPS confirmed constant biases in the horizontal directions, which allowed a constant bias adjustment for all the GPS surveyed snow depth point measurements.
A second "coarse-grid" set of snow depth measurements covered the entire area in plots A-D (∼ 160 m × 160 m) with lower sampling density (Peterson et al., 2015a).The coarsegrid snow data were collected using a tile probe, which had a precision of approximately 0.01 m.Snow depth was measured every 8 m along a measurement tape on the five parallel transects in the coarse grid, which were spaced 40 m apart.The total number of data points was 380 (95 points in each plot).Along the 500 m transects, we used the tile probe along with a measurement tape and measured eight points along each of the three lines.The start and end coordinates of each transect were surveyed with a RTK DGPS and used to georeference the measurement locations.
Ground-based GPR data were acquired over the four study plots and along the three 500 m transects (Peterson et al., 2015b).The instrument (Mala ProEx with 500 MHz antenna) was pulled on a sled.In each plot, we acquired the GPR data at 0.1 m intervals (triggered by an odometer wheel) along 37 lines of 4 m spacing.The start and end coordinates of each transect were surveyed with a RTK DGPS and used to georeference the measurement locations.We compared the distance from the wheel with the distance on tape and confirmed that the difference is generally very small at this site.The error of horizontal positioning is estimated to be about 0.1 m.Several of the GPR lines were co-located with the "coarse-grid" snow depth probe measurements.The GPR technique allowed for denser sampling within the plot relative to the snow depth probe, with more than 50 000 points in each plot.Due to the microtopography at this site, the positioning errors between in situ measurements and GPR data could lead to an error in the radar velocity and snow depth estimation.We evaluate the effect of such positioning errors extensively, as described in Sect.3.1.
The GPR reflection signal from the bottom of snowpack (i.e., the ground surface) was clear, which allowed us to measure the travel time between the top and bottom of snowpack.The GPR processing routine consisted of (1) zero-time ad-justment, (2) average tracer removal, (3) picking the travel time (manually with automated snapping in the ProMAX ® software) of the reflected GPR signal that traveled from the snow surface to the snow-ground interface and back to the snow surface and (4) dividing this travel time by two to obtain a one-way travel time between the snow surface and ground surface.We processed the GPR data including traveltime picking before accounting for topography.More details on GPR processing and theory can be found in Annan (2005) and Jol ( 2009), while more detailed explanation on the use of GPR in the tundra can be found in Hubbard et al. (2013).Differing from previous studies (e.g., Harper and Bradford, 2003), we did not observe echoes from snow layering.This is possibly because of the low antenna frequency (500 MHz), relatively thin snow layers (if present), and the low contrast between various snow layers.In addition, hoar layers or ice layers were not visible in our data or sensed using the probe.Although ice may form at the ground surface, causing the uncertainty of a few centimeters, we did not consider this effect in this study.
Additional campaigns were carried out in 2013-2015 along the 500 m transects only.UAS-based phodar data were collected in July 2013 and 2014 to estimate snow-free ground surface elevation and in May 2015 for estimating snow depth along the transects (Dafflon et al., 2015).To make these measurements, we lifted a consumer-grade digital camera (Sony Nex-5R) to about 40 m above ground level using a kite and acquired downward-looking red-green-blue landscape images, as well as collected some surface elevation data (method described in Smith et al., 2009).The reconstruction procedure was performed using a commercial computer vision software package (PhotoScan from Agisoft LLC).Reconstruction involved automatic image feature detection/matching, structure-from-motion and multiview stereo techniques for 3-D point-cloud generation, and georeferenced mosaic reconstruction (Nolan et al., 2015).High-accuracy georeferencing was enabled by using a network of ground control points placed on the ground (in summer) and on the snow (in winter) that were surveyed with a high-precision centimeter-grade RTK DGPS system.The reconstructed phodar surface elevation models at this site show a resolution of 4 cm × 4 cm.We investigated the accuracy in detail as described in Sect.3.2.
The snow-free ground surface elevation measurements were then subtracted from the snow surface data to estimate the snow depth over the area.The snow depth probe measurements were taken at 183 locations along one of the 500 m transects to validate the phodar-based snow depth estimates.The locations were marked on a measurement tape, the start and end coordinates of which were surveyed with a RTK DGPS and used to georeference the measurement locations.

GPR snow depth analysis
Snow depth can be inferred by multiplying GPR one-way travel time by radar velocity.The radar velocity is determined with the dielectric constant, which depends on snow density in dry snow (Tiuri, et al., 1984;Harper and Bradford, 2003).Depending on site conditions, the snow density can vary in both vertical and horizontal directions (Proksch et al., 2015).In this study, we assume that the depth-averaged radar velocity -which is a function of depth-averaged snow density -is sufficient for estimating snow depth.Thus, we compute the radar velocity based on the known snow depth from co-located snow depth probe measurements as (radar velocity) = (probe-based snow depth) / (GPR one-way travel time).In addition, we investigate whether the lateral variations in snow density are significant at our site.
Identifying co-located points between the GPR and snow depth probe measurements, however, is not a trivial task in polygonal ground, since the topography and snow depth can vary significantly within a meter.To address these issues, we investigate the correlations between the radar velocity and the submeter-scale variability of topography.To link the DEM elevation data to the snow depth probe and GPR data, we selected the DEM elevation (0.5 m×0.5 m resolution) and GPR measurement at the nearest locations to the tile probe measurements.We assume that the effect of positioning errors is larger near the edge of polygons or in the region where the submeter-scale topographic variability is high.We consider that the uncertainty of radar velocity can be reduced by not using the co-located snow depth probe measurements in regions of high submeter-scale variability.To define the submeter-scale variability, we compute the elevation difference within a 1 m radius of each snow depth probe measurement.In addition, the reflections from the troughs could originate from the edge of polygons rather than the location right below the GPR instrument.Such an "edge reflection" effect can lead to overestimation of the radar velocity.We assume that we could detect the presence of the edge reflection by evaluating the systematic bias (i.e., overestimation) in the radar velocity in relation to the submeter-scale topographic variability.

UAS-based phodar snow depth analysis
We first evaluate the accuracy of the phodar-derived digital surface model (DSM) by comparing it to the RTK GPS elevation measurements along the 500 m transects acquired in 2011.Since the phodar-derived DSM was obtained at very high lateral resolution (4 cm × 4 cm), it was more prone to noise or small-scale variability (Nolan et al., 2015).As such, we test three schemes to explore the vertical agreement between the two datasets: (1) nearest points, (2) average elevation within the 0.5 m radius and (3) minimum elevation within the 0.5 m radius.We use the same scheme (the best scheme among the three) for determining the snow-free and snow surface elevation at the co-located points.We then compare the snow depth estimates from phodar and snow depth probe measurements at co-located points (the May 2015 snow data).Since we assume that the phodar snow depth estimates would suffer from the same positioning errors associated with the snow depth probe data as GPR, we eliminate the snow depth probe measurements in the regions where the submeter-scale topographic variability is high.

Spatial variability analysis of topography and snow depth
To quantify the topographic effects in a complex terrain of ice-wedge polygons and to partition micro-and macrotopography, we apply the wavelet transform method to the airborne lidar DEM, which is commonly used for 2-D image processing.The wavelet approach has been applied to DEM in geomorphic studies, including terrain analysis and landslide analysis (Bjørke and Nilsen, 2003;Kalbermatten, 2010;Kalbermatten et al., 2012).In this transform, a high-pass filter (a mother wavelet) and a low-pass filter (a father wavelet) are applied to decompose the DEM into four images at each scale: low-pass, high-pass horizontal, high-pass vertical and high-pass diagonal images.The scale is a parameter in the wavelet transform, representing the width of the filter and the scale of topographic variability (Kalbermatten et al., 2012).Depending on the scale of the wavelet transform, the method yields different images, corresponding to different scales of topographic features.We define this wavelet scale as a topography separation scale.We consider the low-pass image as macrotopographic elevation (i.e., the smoothed version of the original DEM) and the high-pass diagonal image as microtopographic elevation (i.e., the topographic variability associated with ice-wedge polygon development).Removing the large-scale topography has been done in the previous studies in order to capture or quantify the effect of microtopography on carbon fluxes (Wainwright et al., 2015) or soil properties (Gillin et al., 2015).
Correlations between the topographic metrics and snow depth are identified using the Pearson product-moment correlation coefficient (Anderson et al., 2014).At each spatial scale, we can compute micro-and macrotopographic metrics such as slope and curvature as well as their correlations with corresponding probe-measured snow depth.The curvature is of particular interest, since Dvornikov et al. (2015) reported strong correlations between snow surface curvature and snow depth as well as a dependency of this correlation on the DEM resolution (the lower resolution led to lower correlation coefficients).Note that the DEM resolution (0.5 m) in this study is much higher than the one (25 m) in Dvornikov et al. (2015).We compute a wind factor in a similar manner as Dvornikov et al. (2015), with a slight modification.Here we define the wind factor as the inner product of the slope direction and www.the-cryosphere.net/11/857/2017/The Cryosphere, 11, 857-875, 2017 predominant wind direction.With this calculation, the wind factor is smallest in the slope against the wind direction and largest in the slope in line with the wind, which is reasonable and also consistent with visual observations at the site.When the correlation is statistically significant, the metrics are included in a regression analysis (Davison, 2003) to represent the snow depth as a function of the topographic metrics.
A geostatistical approach has been used to investigate the spatial variability of snow depth as well as the scales of variability (Anderson et al., 2014).The standard geostatistical analysis starts with creating an empirical variogram, followed by estimating the spatial correlation parameters (Diggle and Ribeiro Jr., 2007).The spatial correlation parameters include (1) magnitude of variability (or spatial heterogeneity) as variance, (2) fraction of correlated and uncorrelated variability (nugget ratio), (3) spatial correlation length (range) and (4) covariance model (i.e., the shape of decay in the spatial correlation as a function of distance), such as exponential and spherical models.The covariance models (equivalent to variogram models) can be selected to minimize the weighted sum of squares during variogram fitting.
Such spatial variability and correlation are particularly important for interpolating the sparse in situ snow depth measurements.The interpolation can be applied not only for snow depth itself but also for snow surface (snow depth plus elevation) or residual snow depth after removing topographic correlations in the regression analysis.The same geostatistical analysis method is therefore performed for snow surface and residual snow depth.We used the geoR package in statistical software R (Ribeiro Jr. and Diggle, 2001; https: //www.r-project.org/).

Bayesian geostatistical estimation method
We first define that the snow depth at each pixel y i (i = 1, . .., n) is a hidden variable which can be observed only with an added measurement error.In this study, we set the pixel size to 0.5 m × 0.5 m, which corresponded to the lidar DEM resolution.The snow depth distribution (or field) is defined by a vector y = {y i | i = 1, . .., n}.We integrate three datasets: snow depth probe data z p , GPR data z g and lidar DEM z d .The goal of the estimation is to determine the posterior distribution of snow depth conditioned on all the given datasets, p(y | z p , z g , z d ).Following a Bayesian hierarchical approach, we divide this posterior distribution into three sets of statistical submodels (Wikle et al., 2001;Wainwright et al., 2014Wainwright et al., , 2016)).First, data models represent each data value as a function of snow depth at each pixel, depending on different data types.Second, process models describe the spatial distribution of snow depth (i.e., snow depth field) as function of topography and correlation parameters.Finally, prior models define the prior information of parameters.The hierarchical approach breaks down a complex posterior distribution into a series of simple models and hence enables us to capture complex relationships easily.In addition to the snow field vector and data vectors, two parameter vectors are defined: the process-model parameter vector a to represent the heterogeneous pattern of snow depth and the data-model parameter vector b to describe the correlations between the snow depth and the GPR travel time.
We assume a linear model to describe the snow depth field, where A is the design matrix as a function of the topographic metrics as explanatory variables (and hence a function of DEM z d ).The process-model parameter vector a describes the correlation between the topographic metrics and the snow depth field.We assume that the residual of this correlation τ represents the unexplained variability by the topographic metrics and that τ is spatially correlated.The residual term τ is described by a multivariate normal distribution with a covariance , which is determined by a geostatistical analysis (Diggle and Ribeiro Jr., 2007).Although we may include the uncertainty of those geostatistical parameters in the Bayesian estimation (Diggle and Ribeiro Jr., 2007;Lavigne et al., 2017), we assume that those parameters are fixed during the Bayesian estimation process in this study.This is because we have a large amount of point measurements (snow depth probe data), and also it is known that indirect information (such as geophysics) does not significantly improve the estimation of geostatistical parameters (Day-Lewis and Lane Jr., 2004; Murakami et al., 2010).
The data model for the snow depth probe measurements defines the snow depth probe data z p as a function of snow depth y: (2) We assume that the vector ε p is an uncorrelated normally distributed measurement error at each data location with the standard deviation (SD) of σ p .We determine the error based on the precision estimate of each snow depth probe.The snow depth probe data vector z p follows a multivariate normal distribution with the mean vector y and the covariance matrix D p , which is a diagonal matrix with diagonal elements of σ 2 p .Although it is not considered this study, we could include a systematic bias of snow probe measurements as an added shift (Berezovskaya and Kane, 2007).
The data model for the GPR data describes the GPR data z g as a function of the snow depth y at the GPR locations.The GPR data model can be represented by a linear model: where B is a matrix, the diagonal elements of which is b 1 .The error vector ε g is an uncorrelated normally distributed measurement error with the SD of σ g .The SD is computed from comparing the GPR-based snow depth to the probebased one.At the same time, the GPR data model can be written as a function of the parameter vector b such that The Cryosphere, 11, 857-875, 2017 www.the-cryosphere.net/11/857/2017/where Y is the design matrix with the first column being y, and the second column being all one.The parameter vector b = {b 1 , b 0 } represents the linear correlations between the GPR data and snow depth.This alternative model is useful during the estimation procedure described below.The GPR data vector z g follows a multivariate normal distribution with the mean vector y and the covariance matrix D g that is a diagonal matrix with diagonal elements of σ 2 g .The posterior distribution of the snow depth conditioned on the datasets p(y|z d , z p , z g ) is a marginal distribution of p(y, a, b|z d , z p , z g ).By applying Bayes' rule and following the conditional dependencies defined above, we can decompose this posterior distribution as (5) Table 1 defines all the distributions on the right-hand side of Eq. ( 5) based on the models defined in Eqs. ( 1)-(4).We also assume multivariate normal distributions for the prior distributions of the parameter vectors a and b.The posterior distribution in Eq. ( 5) can be computed using the Markovchain Monte Carlo (MCMC) method (Gamerman and Lopes, 2006).Since all the distributions are defined as multivariate normal distributions, it is possible to use efficient Gibbs' algorithm.The MCMC procedure is described in Appendix A.
The convergence can be confirmed by the Geweke convergence diagnostic (Geweke, 1992).The entire workflow is included in Appendix B.

GPR radar velocity analysis
Our results (based on the GPR data and tile probe data collected in May 2012) indicate that the estimated radar velocity itself does not have a systematic dependency on (or trend with) the snow depth or submeter-scale variability of topography in May 2012 (Fig. 2a and b).The correlation coefficient between the radar velocity and snow depth is 0.11 and between the radar velocity and submeter-scale variability is 0.15.The variability of the radar velocity, however, depends on those two factors (i.e., the variability of snow depth and topography).Hence, the variability is higher in areas with shallower snow depths (Fig. 2a).The SD of the radar velocity is 0.039 m ns −1 at the snow depth smaller than 1 SD minus the median snow depth and is 0.019 m ns −1 at snow depth larger than 1 SD plus the median.The radar velocity variability is higher also in localized regions of large submeter-scale topographic variability (Fig. 2b).The SD of the radar velocity is 0.015 m ns −1 at the submeter-scale topographic variability (i.e., elevation difference within a 1 m radius) smaller than 0.05 m, and 0.036 m ns −1 at the one larger than 0.05 m.By selecting the points with the submeter-scale topographic variability < 0.05 m, we obtained a mean radar velocity of 0.25 m ns −1 , which was used for subsequent analysis.
Using the mean velocity value in May 2012, the calculated GPR-based snow depth estimates were compared with www.the-cryosphere.net/11/857/2017/The Cryosphere, 11, 857-875, 2017 the snow depth probe measurements (Fig. 2c).The correlation between the measured and estimated snow depth is high (the correlation coefficient is 0.88), with the RMSE being 5.4 cm and no significant under-or overestimation (the mean bias error −0.16 cm).The selected points in the regions of low submeter-scale topographic variability (red circles) are more tightly distributed around the one-to-one line.In these regions, the RMSE of GPR-based snow depth improved to 2.9 cm with a increased correlation coefficient between the GPR-based and probe-based snow depth to 0.94.These results confirm that snow density variations are limited, and using a constant mean GPR velocity is acceptable.

Snow depth measurements in different polygon types
Figure 3 shows the lidar DEM as well as snow depth probe measurements and GPR estimates in plots A-D (May 2012).
The lidar DEM (Fig. 3a) illustrates the difference among four plots in terms of both macro-and microtopography.For example, plot A has better defined polygon rims and troughs than plot D, although plots A and D are both low-centered polygons.Plot B has round-shaped high-centered polygons, while plot C has flat-centered polygons with well-defined troughs.The average size of polygons is also different, with smaller polygons in plot B and larger polygons in plots A, C and D. In addition, these figures illustrate some macrotopo-graphic trends.Plot C is gradually sloping down towards the east, and plot D has a depression (i.e., DTLB) in the northeastern half.In Fig. 3b shows the snow depth probe data collected using the fine-grid and coarse-grid scheme collected in May 2012.The fine-grid data reveal the detailed heterogeneity of snow depth around a single polygon.For example, the finegrid data in plot A show the snow depth distribution in a lowcentered polygon, including thin snow along the polygon rim and thick snow at the polygon center and trough.Comparison of the fine-grid snow data with the DEM reveals the microtopographic effect such that the troughs and center of the polygon have larger snow depth.The coarse-grid dataset covers the entire plot, although it is much more difficult to ascertain the relationship between the snow depth and microtopography.The snow depth probe data show that the snow depth is highly variable, ranging from 0.2 to 0.8 m in a single plot.
In Fig. 3c, the May 2012 snow depth was estimated from GPR using a fixed radar velocity 0.25 m ns −1 along the lines within the plots, and then interpolated with a simple linear interpolation in between the lines.The high-resolution GPR snow depth estimates are useful for determining if microtopographic features can influence the distribution of snow depths across each study plot.The high-resolution snow estimates over the large area allow us to visually identify the macrotopographic control on snow depth.In plot C, for example, the snow depth does not have an increasing or decreasing trend, even though the elevation gradually decreases towards east.Plot D, in contrast, has more snow accumulation in the eastern part of the domain, which is in the depression associated with DTLB.

Phodar-based snow depth measurements
In the region of the 500 m transects, the phodar-derived snow-free DSMs (Fig. 4a) collected in July 2013 and August 2014 were first compared with the RTK DGPS data (acquired in 2011) in Table 2, using the different schemes to identify co-location.We included the results of both years to confirm the consistency between the two snow-free DSM products at the same terrain.Although all the schemes yielded an excellent accuracy (the RMSE less than 7.0 cm), taking the average provides the lowest RMSE in both years (6.41 cm in 2013 and 6.19 cm in 2014), which is approximately the same as the lidar data (RMSE = 6.08 cm).The phodar-derived snow depth estimates in May 2015 were obtained by differencing the snow surface and snow-free DSM (Fig. 4b).The comparison between the phodar-based snow estimates and the snow depth probe data is favorable (Fig. 4c), with a RMSE of 6.0 cm.When we removed the points that had a large submeter-scale topographic variability in the vicinity (in the same way and the same cutoff values as the GPR snow depth analysis), the RMSE improved to 4.6 cm (Fig. 4c).
The phodar-derived snow depth (Fig. 4b) around the 500 m transects in May 2015 reveals a similar pattern of snow dis- tribution as the GPR data in Fig. 3, having deeper snow in the troughs and the centers of low-centered polygons.The high-resolution image of the phodar data reveals more detail of the microtopographic effect than the interpolated image of the GPR data, particularly in the narrow troughs.The large aerial coverage also shows the minimal effect of macrotopography: while the elevation decreases towards south, the snow depth does not have a large-scale trend.

Variability among different polygon types
Figure 5 shows the box plots of the snow depth, elevation and microtopographic elevation ( elevation) in each plot measured in May 2012.We used the coarse-grid snow depth probe measurements, since the samples are uniformly distributed over each plot.The median snow depth (Fig. 5a) is fairly similar among four plots, even though they have different geomorphologic features and polygon types.Tukey's pairwise comparison test (Table 3) shows that only plot B (small high-centered polygons) is significantly different from the other plots.
The absolute elevation distribution varies among the four plots (Fig. 5b), although the snow depth for each of the plots has similar median values and distributions.Plot A (welldefined low-centered polygons), for example, is at a higher elevation than plots C (flat-centered polygons) and D (lowcentered polygons in DTLB), but the difference in the average snow depth is not statistically significant (Table 3).The microtopographic elevation is computed based on the wavelet transform with the scale of 32 m as described in Sect.3.3 (Fig. 5b).The scale of 32 m was selected to yield the best correlation between snow depth and microtopographic elevation.Plot D (low-centered polygons in DTLB), for example, has less variability in both elevation and snow depth, because plot D has less distinct microtopography than others.In contrast, plot B has the largest variability in both microtopography and snow depth.

Plots Snow depth
A-B 6.34

Correlations between snow depth and topographic indices in May 2012
Among the topographic indices of macro-and microtopography, the snow depth in May 2012 (measured by the snow depth probe) was significantly correlated only to the microtopographic elevation for all plots (Fig. 6a).The correlation coefficient changes with the scale of the wavelet transform that separates micro-and macrotopography.The correlation coefficient is up to −0.8 at plot B (small high-centered polygons) and up to −0.7 at all the data points.The correlation coefficient is different among different plots (i.e., dif- ferent polygon types); the correlation is less significant at plot D (low-centered polygons in DTLB) than other plots.The best correlation (i.e., the largest absolute value) can be achieved at a different scale in each plot (plot B < plot A and plot C < plot D).A significant correlation between snow depth and wind factor of macrotopography was identified only in plot D (low-centered polygons in DTLB; Fig. 6b).The correlation coefficient is up to 0.41 at the scale of 38 m.Other topographic indices (i.e., the slope and curvature of both microand macrotopography, the wind factor of microtopography) are not shown here, since we did not find any significant correlation.Although Dvornikov et al. (2015) reported a strong correlation between snow depth and curvature (snow-free DEM), we did not find any significant correlation in our data.This is possibly because the microtopography at our site was completely filled by snow, and the overall elevation gradient at our site (the elevation difference in our domain is 3.1 m) is much smaller than the one that Dvornikov et al. (2015) The Cryosphere, 11, 857-875, 2017 www.the-cryosphere.net/11/857/2017/reported (the elevation difference in their domain was more than 60 m).

Geostatistical analysis of snow depth
Spatial correlation exists for all three variables in May 2012: snow depth, snow surface and residual snow depth after removing the correlation to the microtopographic elevation (Table 4).The correlation range is less than 20 m for the snow depth, which is consistent with the large variability in a short distance.The snow surface, in contrast, has a larger correlation range (253 m).The estimation of a snow surface height (elevation + snow depth) effectively removes the influence of microtopography, resulting in much a larger correlation range.The variance is comparable between the snow depth and snow surface, while the variance is much lower in the residual snow depth, since the topographic correlation explains a large portion of the snow depth variability.

Snow depth estimation based on lidar DEM
Based on the snow-topography analysis in Sect.4.2, we included the linear correlation between snow and microtopographic elevation in Eq. ( 1) to describe the snow variability in May 2012.We used the Shapiro-Wilk normality test to confirm that the residual of the linear correlation, defined by τ in Eq. ( 1), follows a normal distribution (the p value of rejecting this hypothesis was 0.21).The first column of the design matrix A is the microtopographic elevation at all the pixels, and the second one is a vector of all 1s.The parameter vector a is a 2-by-1 vector with the linear correlation parameters (slope and intercept).The Bayesian method (Sect.3.4) yielded 10 000 equally likely fields of the snow depth from the posterior distribution in Eq. ( 5).The Bayesian estimated mean snow depth field over the full study domain in May 2012 (Fig. 7a) captures the effects of microtopography, such as more snow accumulation in polygon troughs and centers of low-centered polygons.The snow depth does not have any large-scale trends over the full study domain, which is different from the lidar DEM in Fig. 1b, but consistent with the interpolated GPR snow depths depicted in Fig. 3c and the measured UAS snow depth measurements depicted in Fig. 4b.The variability is larger in the southern region where there are high-centered polygons with deep troughs.
In addition, we compared this result (Fig. 7a) with the mean field from estimating the snow surface elevation and subtracting the ground surface elevation (Fig. 7b).In this estimation, we used the same Bayesian algorithm described in Sect.3.4, except that we removed the topographic correlations and assumed a standard geostatistical model for snow surface (Diggle and Ribeiro Jr., 2007).In other words, we had the same algorithm except that we modified Eq. ( 1) to y = c − z + τ , where y + z represents the surface elevation and c is a constant.Although the two mean fields (Fig. 7) are similar in the central regions that have many measurements, the regions without any measurements have a significant deviation.This is because the snow surface estimation did not capture the change in macrotopography (e.g., the drainage feature in the southern part of the domain).
The estimated SD of the Bayesian-derived snow depth over the study domain (Fig. 8a) also shows a significant difference from the one based on the snow surface interpolation (Fig. 8b).This SD represents the uncertainty in the estimation.In both cases, the SD is smaller near the measurement locations along the transects and within the four plots.However, when the topographic correlation is included (Fig. 8a), the SD increases more rapidly as the pixel is farther away from the data points.This is due to the fact that the spatial correlation range is small for the residual snow depth after removing the topographic correlation (Table 4).
Validation of the snow depth estimates over the study area (plots A-D and the 500 m transects) was performed by comparing the estimates with the snow depth probe data (May 2012) not used in the Bayesian snow depth estimation.We selected 100 points randomly from the snow depth probe data (all the locations in plots A-D and the 500 m transects), using a uniform distribution.The validation results (Fig. 9) show that the estimated confidence interval captures the probe-measured snow depth.The estimated snow depth is distributed along with the one-to-one line without any significant bias.The estimation, including the topographic correlation (Fig. 9a), has a tighter confidence interval and better estimation results than the one from interpolating the snow surface (Fig. 9b).The RMSE for the Bayesian method of estimating snow depth including the topographic correlation is 6.0 cm, while the RMSE for the interpolated snow surface is 8.8 cm.

Different observational platforms
Our analysis showed that GPR data provided the endof-winter snow depth distribution with high accuracy (RMSE = 2.9 cm) and resolution (10 cm along each line).The GPR-based estimation, however, requires care, particularly regarding the estimation of radar velocity and associated possible errors, such as those due to positioning.Although the radar velocity is known to depend on the snow density, we attribute the variability of radar velocity at our site to random or positioning errors.Three results support this claim.First, the variability of radar velocity is smaller in a thicker snow pack, suggesting the small contribution of the error relative to the overall snow depth.The relatively low topographic variability over the site (compared to mountainous terrains) would have contributed to this fairly uniform radar velocity.Second, the radar velocity variability depends on the submeter-scale variability of the topography in the vicinwww.the-cryosphere.net/11/857/2017/The Cryosphere, 11, 857-875, 2017  ity of the calibration points, suggesting the impact of positioning errors.Third, there was no systematic trend in the radar velocity as a function of the snow depth or topographic positions.We developed a simple methodology (described in Sect.3.1) to select co-located calibration points based on the submeter-scale variability of topography, which proved to be useful to compute accurate velocity.We note thateven though the depth-averaged radar velocity and hence the depth-averaged snow density have little variability over the space -the snow density could be variable vertically along the depth.From snow coring, we indeed found some layers of ice created by winter rain events that were not detected by the GPR or with probe; it is possible that there might be a difference in the depth-averaged density and radar velocity at a later time, when the snow pack starts to melt in a heterogeneous manner.UAS-based phodar provided an attractive alternative for estimating snow depth at high resolution over a large area.With much less labor and time, UAS-based phodar can provide many more sample points than GPR.The phodar-based snow depth, however, was less accurate than ground-based GPR or snow depth probe measurements (RMSE = 6.0 cm).The main contribution of this error resulted from the snowfree elevation, since RMSE for the surface DSM is around 6 cm.We note that the RMSE of 6.0 cm is still significantly more accurate than the previous lidar and other airborne surveys (e.g., Deems et al., 2013;Harpold et al. 2014;Nolan et al., 2015).
The phodar-based approach is expected to continue its trajectory of continuous improvements in terms of technical as-pects, ease of use and accuracy.At the time of our campaign, we were allowed to use only a kite due to regulations, which led to a limited number of pictures that could be used to reconstruct the DSM.The accuracy will significantly improve with the use of a light unmanned aerial vehicle (UAV).Although UAS-based lidar acquisition technology continues to improve (e.g., Anderson and Gaston, 2013), as is expected to be a powerful alternative to characterize snow, the lidar device is still significantly more expensive than a conventional camera (roughly by factor of 100).Given that the vegetation height is fairly small in the Arctic tundra, the phodar technique is an affordable alternative.
For all the types of measurements, accurate positioning was critical in the polygonal tundra due to microtopography.The GPS snow depth probe (Snow-Hydro), for example, had a positioning error larger than several meters and required extra post-processing to correct the locations.However, measuring the RTK DGPS at all the snow depth measurement locations would not be realistic since it would take time.We found that having a measurement tape and measuring the start and end points by the DGPS was a reasonable approach, when the snow surface is smooth and hard.In this study, we used the snow depth probe data as the true snow depth to compare with other measurements (i.e., GPR, phodar, and Bayesian estimation).To improve the accuracy further, it would be necessary to quantify the uncertainty in the snow depth probe associated with the vegetation and other issues (Berezovskaya and Kane, 2007).

Snow depth variability
The end-of-year snow depth distribution at the ice-wedge polygons was highly variable over a short distance in May 2012.The snow depth was, however, significantly correlated with the microtopographic elevation, suggesting that the snow depth could be described by microtopography.The wind-blown snow transport leads to significant snow redistribution and fills microtopographic lows (i.e., troughs and centers of low-centered polygons) with thicker snow pack (e.g., Pomeroy et al., 1993).The redistribution also results in the smooth snow surface, following the macrotopography.The exception was observed at the edge of the DTLB, where the abrupt change in macrotopography led to increased accumulation in the depression.This is a similar effect to that observed along the riverbanks by Benson and Sturm (1993).Although the tundra ecosystem studies have focused on the effect of microtopography (e.g., Zona et al., 2011), the macrotopography also may be important when we characterize snow distribution over a larger area.
The "average" (or median) snow depth over a hundredmeter scale (i.e., the size of plots A-D), however, was fairly uniform across the site despite the different polygon types in May 2012.Plots A (well-defined low-centered polygons) and C (flat-centered polygons), for example, have different polygon types, but they have a similar median snow depth.This is because microtopography and microtopographic features (i.e., polygon troughs, rims) mainly control the snow distribution.Plot B (small high-centered polygons) is an exception, having smaller median snow depth than the other plots.Plot B has the largest variability in microtopography, characterized by the small round high-centered polygons, like numerous small mounds (Fig. 3).Such mounds are prone to erosion by the wind, and hence lead to less snow trapping and accumulation.
Identifying such correlations between snow depth and topography requires an effective approach to separate microand macrotopography.Our wavelet analysis revealed that the separation scale depends on the polygon sizes; for example, the larger polygons in plots A (well-defined low-centered polygons) and C (flat-centered polygons) lead to a larger separation scale than the smaller polygons in plot B (small high-centered polygons).It is a challenge to map macrotopography accurately over a larger area, particularly at the www.the-cryosphere.net/11/857/2017/The Cryosphere, 11, 857-875, 2017 present site, where different types and sizes of polygons mix together.Although we used the same scale for the estimation, an improved polygon delineation algorithm will possibly enable us to separate micro-and macrotopography in the future (e.g., Wainwright et al., 2015).

Snow depth estimation
The developed Bayesian approach enabled us to estimate the snow depth distribution over a large area based on the lidar DEM and the correlation between the snow depth and topography.Although this paper only used the ground-based GPR and snow depth probe measurements collected at the same time, phodar could be easily included in the same framework.The Bayesian method allowed us to integrate three types of datasets (lidar DEM, snow depth probe and GPR) in a consistent manner and also provided the uncertainty estimate for the estimated snow depth.Taking into account the topographic correlation explicitly improved the accuracy of estimation significantly (RMSE = 6.0 cm), compared to interpolating the snow surface and subtracting the DEM (RMSE = 8.8 cm).
Our approach can be extended to snow estimates over both time and space.The correlations between snow depth and topography may change over time.In early and later winter, for example, the snow depth would be more affected by curvature and slope of microtopography, since the microtopographic lows (troughs and centers of the low-centered polygons) are not filled by snow.It would be possible to quantify the seasonal changes in the topography-snow correlations by designing a full season ground-based measurement campaign and acquisition of remote sensing snow depth measurements (by phodar or lidar) that monitored the same site over several years to account for interannual variability.The Bayesian method presented here is flexible enough to account for changes in parameters over time for the spatiotemporal data integration (e.g., Wikle et al., 2001).Although physically based snow distribution models can be used for the same purposes (e.g., Pomeroy et al., 1993;Liston andSturm, 1998, 2002), it is difficult to parameterize all the processes, such as sublimation and turbulent transport.Our datadriven approach provides a powerful alternative to distribute snow depth based on various datasets.

Summary
In this study, we explored various strategies to estimate the end-of-year snow depth distribution over an Arctic icewedge polygon tundra region.We first developed an effective methodology to calibrate GPR and phodar in the presence of submeter-scale variability of topography.We then investigated the characteristics and accuracy of three observational platforms: snow depth probe, GPR and phodar.The phodarderived snow depth estimates have great potential for accurately characterizing snow depth over larger regions (with an RMSE of 4.6 cm), relative to the in situ snow depth measurements.The GPR snow depth estimates were slightly more accurate (with an RMSE of 2.9 cm), but they required considerable more effort to obtain and require complex postprocessing to minimize errors associated with radar positioning.
We investigated the spatial variability of the snow depth and its dependency on the topographic metrics.At the peak snow depth during our data acquisition, the snow depth was highly correlated with microtopographic elevation (the correlation coefficient of up to −0.8), although it was highly variable over short distances (the correlation range of 12.3 m).It is considered that the wind redistribution filled the microtopography by snow and created a snow surface following macrotopography at the site.The challenge was to separate macro-and microtopography, since the separation scale was not arbitrary and depended on the polygon size.The wavelet analysis provided an effective approach to identify this separation scale.
The Bayesian method was effective at integrating different measurements to estimate snow depth distribution over the site.Although our estimation is based on the data collected from a one-time campaign, and the correlations to topography may change over time, the approach developed here is expected to be applicable for estimating both spatial and temporal variability of snow depth at other sites and in other landscapes.

Figure 1 .
Figure 1.(a) Location of Barrow, Alaska, USA, and Barrow Environmental Observatory (BEO) from Hubbard et al. (2013).(b) NGEE Arctic site with the digital elevation map from the airborne lidar (in meters).The black boxes are the intensive sampling plots (plots A, B, C and D).The white rectangles are the fine-grid snow depth measurements by a snow depth probe.The three black lines represent the 500 m transects.

Figure 2 .
Figure 2. Radar velocity as a function of (a) co-located snow depth measured by a snow depth probe and (b) elevation difference (i.e., topographic variability) within 1 m.(c) Comparison between the probe-derived and GPR-derived snow depth at all the co-located locations (blue circles) and at selected locations (red circles) where topographic variability is low.In (a), the black vertical line is the median snow depth, and the dotted lines are ±1 SD from the median snow depth.In (b), the black line is the cutoff elevation difference of 0.05 m.
p y, a, b|z d , z p , z g ∝ p z g |y, b p z p |y p (y|a, z d ) p(a)p(b).

Figure 3 .
Figure 3. Elevation and snow depth in plots A, B, C and D. Panel (a) is lidar DEM (in meters), panel (b) is the probe measured snow depth (in meters), and panel (c) is the interpolated snow depth estimated using GPR (in meters).

Figure 4 .
Figure 4. (a) Phodar-derived DSM in meters (August 2014), (b) phodar-derived snow depth in meters (May 2015) and (c) comparison between the phodar-based and probe-based snow depth at all the locations (blue circles) and at selected locations (red circles) having low topographic variability (the sub-meter elevation variability less than 0.05 m).The black line in (b) represents the snow depth probe measurements every 3 m along the 500 m transect.

Figure 5 .
Figure 5. Box plots of (a) snow depth, (b) elevation and (c) microtopographic elevation in plots A-D.

Figure 6 .
Figure 6.Correlation coefficients between snow depth and topographic metrics as a function of the wavelet scale: (a) the microtopographic elevation and (b) the wind factor of macrotopography.The different colors represent different plots (plots A-D) or all the data (All).Each dash line represents the scale that maximize the magnitude of the correlation coefficient.

Figure 7 .
Figure 7.The estimated mean snow depth across the site (in meters) based on (a) the proposed Bayesian method including the correlation to microtopography and (b) the kriging-based interpolation of the snow surface.The spatial extent is the same as Fig. 1b.

Figure 8 .
Figure 8.The estimated standard deviation of snow depth across the site (in meters) based on (a) the proposed Bayesian method, including the correlation to microtopography, and (b) the kriging-based interpolation of the snow surface.The spatial extent is the same as Fig. 1b.

Figure 9 .
Figure9.Estimated mean and confidence intervals from the Bayesian method, compared to the probe-measured snow depth by (a) using the correlation to microtopography and (b) interpolating the snow surface.The red circles represent the snow depth at the validation locations (the snow depth probe measurements not used in the estimation), the blue lines are the confidence intervals based on the standard deviation (SD) multiplied by 1.9 (94 % confidence intervals), and the black lines are the one-to-one line.

Table 1 .
Multivariate normal distribution defined for each variable.

Table 2 .
Root mean square error (RMSE) between the phodarderived DSM and RTK DGPS elevation measurements based on the three schemes: nearest neighbor, average, and minimum elevation within the 0.5 m radius.

Table 4 .
Estimated geostatistical parameters and covariance models for snow depth, snow surface and residual snow depth.