Effect of snow microstructure variability on Ku-band radar snow water equivalent retrievals

Spatial variability in snowpack properties negatively impacts our capacity to make direct measurements of snow water equivalent (SWE) using satellites. A comprehensive data set of snow microstructure (94 profiles at 36 sites) and snow layer thickness (9000 vertical profiles across nine trenches) collected over two winters at Trail Valley Creek, NWT, Canada, was applied in synthetic radiative transfer experiments. This allowed for robust assessment of the impact of estimation accuracy of unknown snow microstructural characteristics on the viability of SWE retrievals. Depth hoar layer thickness varied over the shortest horizontal distances, controlled by subnivean vegetation and topography, while variability in total snowpack thickness approximated that of wind slab layers. Mean horizontal correlation lengths of layer thickness were less than a metre for all layers. Depth hoar was consistently ∼30 % of total depth, and with increasing total depth the proportion of wind slab increased at the expense of the decreasing surface snow layer. Distinct differences were evident between distributions of layer properties; a single median value represented density and specific surface area (SSA) of each layer well. Spatial variability in microstructure of depth hoar layers dominated SWE retrieval errors. A depth hoar SSA estimate of around 7 % under the median value was needed to accurately retrieve SWE. In shallow snowpacks <0.6 m, depth hoar SSA estimates of ±5 %–10 % around the optimal retrieval SSA allowed SWE retrievals within a tolerance of ±30 mm. Where snowpacks were deeper than ∼30 cm, accurate values of representative SSA for depth hoar became critical as retrieval errors were exceeded if the median depth hoar SSA was applied.

1 Introduction

Seasonally snow-covered non-glaciated Arctic terrestrial environments north of tree lines cover approximately 5.05×10⁶ km² (Walker et al., 2005). Layering of snow, where distinct differences in snow properties exist between vertically adjacent strata (Fierz et al., 2009), is spatially heterogeneous in Arctic regions with a dense wind slab layer overlaying less-dense depth hoar (Benson and Sturm, 1993; see Fig. 1; Derksen et al., 2009). As depth hoar and wind slab have strongly diverging microwave scattering properties (Hall et al., 1991), the relative proportion of each strongly influences Ku-band radar backscatter (Yueh et al., 2009; King et al., 2015, 2018). Knowledge of how layers vary within Arctic snowpacks is therefore critical to the assessment of uncertainty in radar-based retrievals of snow water equivalent (SWE) and forward models of snow radiative transfer.

Figure 1 Locations of snow trenches excavated in April 2013 are indicated by the pink squares. Locations of snow pits excavated in April 2013 and March 2018 are indicated by the blue crosses. Brown, yellow, green, and blue colours on the map denote the different terrain types identified by the topographic position index. In 2018, five additional pits were located on Husky Lakes (not on map), approximately 5 km to the east of the presented domain.

Subnivean topography, wind redistribution, and vertical thermal gradients dominate the formation of layers in Arctic snowpacks (Benson and Sturm, 1993; Sturm and Benson, 2004). Grassy tussocks are common in tundra environments, allowing early-winter snowfall to collect in wind-protected hollows between tussock mounds. Strong thermal gradients in shallow early-season snowpacks cause extreme thermal metamorphism required to initiate growth of large depth hoar crystals and chains (Sturm et al., 1997). Growth of depth hoar in this lowest layer of the snowpack then continues throughout the winter, usually more so than in any other snowpack layer as it is subject to strong thermal gradients for the longest period of time. Subsequent snowfall events throughout the winter contribute to the development of high-density layers of wind-compacted snow crystals through aeolian-driven redistribution (Derksen et al., 2014). These wind-compacted layers become prominent when the total snow depth exceeds the height of the tussocks and the snow surface is fully exposed to winds. Tundra shrubs have a similar influence on snow catchment and metamorphism: reducing local wind velocities and providing shelter for early-season snow deposition, which favours development of depth hoar (Sturm et al., 2001).

Highly heterogeneous microstructural properties of snowpacks are problematic for retrieval of snow water equivalent (SWE) via remote sensing. While passive microwave remote sensing historically has provided estimates of global SWE distributions (Kelly, 2009), much uncertainty exists, largely due to the impact of seasonal snow microstructural evolution on microwave scattering in SWE retrieval models (e.g. Derksen et al., 2014). However, recent experimental work using Ku-band radar suggests using two frequencies, each with different sensitivities to wind slab and depth hoar, may mitigate this primary source of uncertainty (King et al., 2018; Lemmetyinen et al., 2018). In order to apply this two-frequency approach in a distributed manner, we need an understanding of layer length scales: horizontal distances over which the physical properties of each layer decouple and become statistically uncorrelated to each other. The spatial scales of interest depend on the application. For SWE distribution within tundra catchments, it is critical to understand layer variability at the landscape scale (10–1000 m), where different topographic units (e.g. plateau, slope and valley) are subject to different snowdrift, scour, and sublimation processes. Understanding this resolution of landscape-scale variability in snow properties (layers, density, and microstructure) is particularly relevant to future active microwave satellite mission concepts, as well as distributed hydrological modelling which uses meteorological inputs from high-resolution numerical weather prediction models. While the range of spatial variability in snowpack layering and snowpack properties can be estimated from multiple one-dimensional profiles distributed throughout a catchment, this alone will not adequately describe variability in horizontal correlation length scales of snowpack properties. Knowing how these correlation length scales change between different landscape topographic units can help upscale centimetre-scale field measurements of snow properties characterising spatial uncertainty in the statistical distributions of parameters with strongly different microwave-scattering capabilities. Such statistical distributions could then be used in a radiative transfer model, such as the recently developed Snow Microwave Radiative Transfer model (SMRT) (Picard et al., 2018), to address how accurately the information on snowpack properties needs to be known to inform a viable SWE retrieval. Consequently, for the first time, this opens up the potential to explore how variability in snow layers might impact radar backscatter and retrievals of SWE from backscatter. To do so, this study will do the following:

1. quantify the spatial variability in snow depth and layer thickness for surface snow, wind slab, and depth hoar within the Trail Valley Creek catchment;

2. determine representative snow microstructural parameters, specifically specific surface area (SSA) and density, and their associated variability for individual tundra snowpack layers;

3. use these relationships to construct a series of synthetic snowpacks with a realistic range of parameters to quantify the impact of spatial uncertainty in snow microstructural parameters on SWE retrieval accuracy from radar backscatter at two frequencies (13.4 and 17.2 GHz) with SMRT.

2 Methods

2.1 Field data

Field data presented in this study were collected between 4–9 April 2013 and 14–22 March 2018 within the research basin of Trail Valley Creek (TVC), NWT, Canada (68^∘44^′ N, 133^∘33^′ W), located at the southern edge of the Arctic tundra. Measurements of snow microstructure were made mainly on graminoid tundra, which dominates the land cover, as well as in patches of taller shrubs (willow or alder) found on south-facing slopes and in proximity to drainage channels and water features (Marsh et al., 2010). Snow pit and snow trench locations (Fig. 1) were focussed on gently undulating (<5^∘ slope angle) upper tundra areas that were exposed to high winds while also incorporating some slope and valley bottom areas, representative of more sheltered areas in the catchment. Hourly air temperatures were measured throughout both winters to provide temporal context of freeze-up, mid-winter melt events, and snowmelt onset.

2.1.1 Snowpack properties

To investigate the spatial variability in snowpack layering within the snowpack at TVC, nine snow trenches each ranging in length from 5 to 50 m were excavated in 2013. Of the nine snow trenches, one 50 m (trench 4) and six 5 m trenches were located in gently undulating upper tundra plateau areas, while two 5 m trenches (trench 6 and 7) were located within valley bottoms (Fig. 1). Following the methods of Tape et al. (2010), at each trench, near-infrared (NIR) photography was used to quantify two-dimensional changes in snowpack layering. The positions of all layer boundaries were subsequently geolocated at 1 cm resolution throughout the length of each trench (Watts, 2015); a 5 m trench therefore provides 500 vertical profiles of snowpack layering. Within each 5 m trench section (including ten 5 m sections of the 50 m trench), measurements of density and specific surface area (SSA) were made in a single vertical profile and subsequently assigned to individual layers, which were assessed through visual inspection and hardness. SSA was measured using both an InfraRed Integrating Sphere (IRIS) (Montpetit et al., 2012) for the trenches and pits in 2013 and 2018 and an A2 Photonic Sensors IceCube measurement system for the pits in 2018; both measurement systems followed principles presented in Gallet et al. (2009). Layers were often discontinuous due to the spatial variability in subnivean topography and aeolian processes. One-dimensional vertical profiles of snow properties, in combination with inspection of two-dimensional NIR imagery, allowed individual layers to be manually classified into one of three microstructure types: surface snow (SS), wind slab (WS), or depth hoar (DH) (Fig. 2). While the surface snow layer is often of low density (<100 kg m⁻³) and dominated by decomposing and fragmented precipitation particles, surface snow may have been subject to metamorphism or melt, creating rounded grains and melt forms which can increase the layer density (100–300 kg m⁻³). Depth hoar and wind slab layer properties followed the classifications of Fierz et al. (2009). Spatial variability in layer thicknesses was assessed using unidirectional semivariogram to estimate correlation length scales of layer thicknesses: horizontal distances over which the thickness of a layer along a trench becomes statistically uncorrelated. The range to sill of a semivariogram, using a stable bounded fitting model, was used to quantify this distance for each layer (SS, WS, DH) in all trenches.

Figure 2 (a) Stitched NIR images of the trench face (trench 10). (b) Layer boundaries derived from NIR imagery; blue lines highlight boundaries between snow and air, surface snow and wind slab, and wind slab and depth hoar. The brown area is subnivean soil or vegetation. Symbols describe snow type following the classification by Fierz et al. (2009).

In addition to trenches, measurements of the same snow properties were made in 85 snow pits at 54 locations throughout TVC, where each in situ snow measurement was attributed to one of three layers. Whilst only trench data were used in the semivariogram analysis, data from the trenches and pits in 2013 and the pits in 2018 were combined to determine layer thickness as a function of snow depth and the microstructural properties of the snow within those layers.

2.1.2 Land surface slope classification

To classify the surface topography of TVC by slope position and landform category (Fig. 1), a topographic position index (TPI) was calculated following the methods of Weiss (2001). A 1 m resolution digital elevation model (DEM) was used, created during snow-free conditions in August 2008 (Hopkinson et al., 2011), with an Optech ALTM 3100 lidar. The absolute vertical accuracy of the 1 m DEM was at best ±25 cm, and the horizontal positional accuracy was ±50 cm. The TPI values were combined with slope information to classify the study domain into the following landforms: flat upland plateau ($<\mathrm{5}{}^{\circ }$), flat valley bottom ($<\mathrm{5}{}^{\circ }$), slopes ($>\mathrm{5}{}^{\circ }$), lakes (lake extent was extracted from 1:50 000 Canadian topographic maps).

2.1.3 Airborne lidar snow depths

A spatially continuous surface of snow heights in TVC was measured in April 2013 using an airborne laser-scanning Riegl LMS-Q240in lidar (Johnson et al., 2013). The lidar, flown at 500 m above ground level, had a laser shot footprint diameter of ∼20 cm, which was aggregated to a 1 m×1 m surface height product across a swath width of 500 m. Typical measurement errors associated with this system were up to ±20 cm (Johnson et al., 2013). Airborne lidar snow depth was then estimated based on elevation differences between these snow surface heights and the snow-free summer 2008 DEM data (Hopkinson et al., 2011); see Sect. 2.1.2. The 1 m×1 m airborne lidar snow depth raster was then resampled to 10 m×10 m resolution using a cubic interpolation.

In a subset of all lidar snow depths, King et al. (2018) demonstrated that lidar and in situ measured snow depths closely agreed in relatively flat, open environments. Lidar and in situ measured mean snow depths in upland tundra areas of TVC were shown to have a RMSE of 8.5 cm (King et al., 2018), while additional unpublished analysis showed a ∼14 cm positive bias of lidar snow depths over the whole TVC domain.

2.2 SWE retrieval errors using SMRT

The Snow Microwave Radiative Transfer (SMRT) backscatter and emission model (Picard et al., 2018) was used to illustrate potential retrieval error from a dual Ku-band radar system (cf. King et al., 2018; Lemmetyinen et al., 2018). Three layers were assumed within the snowpack for snow up to 0.7 m in depth: depth hoar (DH), wind slab (WS), and surface snow (SS), with layer thickness dependent on total snow depth and based on relationships derived from snow trenches. For deeper snow, generally only wind slab and depth hoar layers were found (see Sect. 3.1), so only two layers were simulated for snow depths greater than 0.7 m. While ice lenses were occasionally present in the 2018 snowpack and volumetric field samples of snow density and SSA contained sections of ice lenses, their impact on backscatter was not explicitly modelled by SMRT. There is a need for detailed field measurements of ice lenses coincident with radar measurements as a priority for future SMRT evaluations. Ice lenses may be simulated in SMRT in terms of dielectric discontinuities, although coherent effects as a result of ice lens thickness are yet to be included.

A set of three experiments was performed that considered spatial variability in SSA in each of the DH, WS, and SS layers in turn. Synthetic scenarios were used to simulate “truth” backscatter of the scene, with information from observed spatial variability in microstructural properties. Whilst spatial variation in SSA was represented in one layer in the truth simulations, horizontally homogeneous snow was assumed for all layers in the retrieval, with snow depth retrieved as a function of the estimated snow microstructure. Davenport et al. (2008) used a similar synthetic scene methodology to quantify soil moisture error from passive microwave observations.

For simulation of the truth backscatter, density was held constant for each layer. SSA observations were classified by five equally sized intervals across the observed SSA range. Five intervals were chosen to describe the SSA distribution adequately whilst maintaining simplicity. Truth scenes were constructed from five backscatter simulations and aggregated with weights taken from the histogram frequency:

$$\begin{array}{}\text{(1)}& {\mathit{\phi }}_k^{\mathrm{tru}}={\sum }_{n=\mathrm{1}}^{\mathrm{5}}{w}_n{\mathit{\phi }}_{n,k},\end{array}$$

where ${\mathit{\phi }}_k^{\mathrm{tru}}$ is the aggregated truth backscatter at frequency k and φ_n,k is the backscatter simulated for a three-layer snowpack with nth SSA value. To isolate the impact of individual layers, SSA in two layers was kept spatially constant whilst being allowed to vary in the third layer, as illustrated in Fig. 3. The weights, w_n, derived from the histogram frequencies f_n for each SSA interval ($\mathrm{1}\le n\le \mathrm{5}$) are

$$\begin{array}{}\text{(2)}& w_n={\displaystyle \frac{f_n}{\mathrm{\Sigma }{f}_n}}.\end{array}$$

Aggregated together, the set of five simulations represents spatial variability in that the frequency and weights represent the proportion of a scene that might take each SSA value.

Figure 3 Illustration of three-layer (SS: surface snow; WS: wind slab; DH: depth hoar) truth scene with spatial variability in specific surface area of WS (SSA_WSa to SSA_WSe) and the retrieval scene with horizontally uniform SSA. Spatial distribution in truth SSA given by observed spatial distribution at TVC over two winters (2013 and 2018).

In retrieval snowpack scenes, SSA was assumed constant within layers, as shown by Fig. 3. For the spatially constant layers, the same SSA was used in both truth and retrieval simulations. For the layer with spatially varying SSA, a range of values was used to represent an unknown homogeneous SSA. As such, the retrieval does not attempt to retrieve SSA directly, but comparisons of the SWE retrieval errors for a given true snow depth allow for selection of the SSA that best represents the spatially variable truth. Simulations with a range of estimates of unknown SSA were used to determine the backscatter of homogeneous scenes and retrieved depth (subsequently converted to SWE) identified as the minimum of a cost function. Median observed layer densities used for the truth simulations were also used for the retrieval backscatter and SWE calculations. For these simulations, a simplified version of the cost function given in Lemmetyinen et al. (2018) was used (Eq. 3). From the optimal dual-frequency approach to SWE retrieval determined by Lemmetyinen et al. (2018), the retrieval algorithm in this study was based on the backscatter difference at two frequencies (13.4 and 17.2 GHz at VV polarisation) and a fixed incidence angle of 35^∘. Retrievals followed an equivalent SMRT model configuration to the forward modelling of Lemmetyinen et al. (2018) and King et al. (2018), i.e. an exponential snow microstructure model and the improved Born approximation electromagnetic theory. The model approaches differ in the solution of the radiative transfer equations: a multiflux solver (Picard et al., 2004, 2014) was used in the SMRT simulations, whereas a six-flux solver was used in King et al. (2018).

Perfect radar observations were assumed (i.e. no noise added to truth backscatter), and an isothermal temperature of 265 K and constant soil backscatter contribution of −13 dB were assumed in both truth and retrieval simulations. The assumed soil properties are for illustration purposes only, given the lack of bare frozen soil backscatter observations at this frequency. This appears to be a plausible value from the snow to snow-free transition period shown in Scipal et al. (2002), and an alternative soil backscatter of −10 dB was used to test the robustness of conclusions from the −13 dB simulations; however, a full error budget study should consider a range of values. The simplified cost function was

$$\begin{array}{}\text{(3)}& \mathrm{CF}(d,\mathrm{SSA})={\left[{\mathit{\phi }}_{\mathrm{diff}}^{\mathrm{sim}}(d,\mathrm{SSA},x_{\mathrm{1}},\mathrm{\dots },x_m)-{\mathit{\phi }}_{\mathrm{diff}}^{\mathrm{tru}}\right]}^{\mathrm{2}},\end{array}$$

where φ_diff is the backscatter difference φ_17.2 VV− φ_13.4 VV, d is snow depth, the estimated microstructure is determined from the SSA, and x_m represents other parameters assumed to be known exactly, i.e. the same as used in the truth simulations. For these simulations the exponential autocorrelation length used in the SMRT simulations (both truth and retrieval) is calculated with Eqs. (5) and (10) of Mätzler (2002).

In summary, the methodology was used to isolate the impact of the spatial variability of SSA in a single layer on retrieval accuracy from all other sources of retrieval error. Therefore, the following was assumed: (1) perfect radar observations, (2) perfect knowledge of snow layer densities, (3) perfect knowledge of the soil, (4) correct transformation of SSA into exponential autocorrelation length, and (5) perfect knowledge of SSA in all but one layer in the snowpack. Single (homogenous) values were used to represent the unknown SSA of the remaining layer. Only depth is estimated by the retrieval, which is then transformed to SWE via the known densities to compute the SWE retrieval error.

3 Results

3.1 Snowpack properties

The winter of 2017–2018 was warmer than 2012–2013 (Fig. 4). Similar air temperatures were observed in both winters from October through November, but December through March in 2017–2018 was on average 9 ^∘C warmer. Importantly, during 2017–2018 there were three short (<1 d) periods where mid-winter air temperatures increased above −5 ^∘C and approached melting point. The warm period on 15 January 2018 coincided with reports of light freezing rain at Inuvik airport, the nearest weather observing station 49 km away. Therefore, while similar meteorological conditions influenced early-season winter snowpack accumulation and depth hoar metamorphism, a spatially extensive ice lens (>1 mm thick) was formed in the 2017–2018 snowpack because of the January warm event. Ice lenses provide potential to restrict vertical vapour diffusion, although such impacts may be limited due to the additional presence of a dense wind slab layer of already tightly packed snow grains. Ice lens formation can also restrict the flux of blowing snow, reducing the potential for subsequent drifting. Consequently, through the aggregation of snowpack measurements over two winters with strongly differing meteorological conditions, the resulting data set provides a highly valuable description of tundra snowpack properties in a warming Arctic.

Figure 4 Air temperature at Trail Valley Creek during winter 2012–2013 and 2017–2018.

Figure 5 compares statistical distributions of snow depths from lidar in three topographically delimited subdomains (Fig. 1: flat upper plateau, slope, and flat valley bottom) of the TVC catchment. Median snow depths were very similar (0.60–0.65 m) across all subdomains. The interquartile range of snow depths on slopes was largest, reflecting enhanced drift processes that preferentially redistributed SWE to wind-sheltered areas of sloping terrain. However, similarity in frequencies of snow depths greater than 0.9 m between flat upper plateau and sloped subdomains suggested preferential SWE deposition also occurred over relatively flat (<5^∘) terrain, in addition to drifts in leeward slopes aligned perpendicular to prevailing northwesterly winds (King et al., 2018). This is of importance as the flat upland plateau, where the majority of trenches and pits were located, was the areally dominant subdomain (66 % of the total TVC lidar coverage in contrast to 19 % slopes and 8 % flat valley bottom). Consequently, field measurements well represented both the range and frequency of the total TVC lidar snow distribution. This type of exposed, flat, largely unforested terrain is representative of pan-Arctic tundra environments, allowing potential for implications to be drawn beyond TVC.

Figure 5 Snow depths (limited to 2 m) in TVC from airborne lidar: (a) to (c) histograms of three topographically delimited subdomains of the TVC catchment (red) overlaid on the histogram of all snow depths (blue); (d) distributions of snow depth by land surface type: blue box (interquartile range), red line (median), and whiskers (dashed lines) extend from the end of each box to 1.5 times the interquartile range; outliers beyond this range are omitted.

Thickness of tundra snowpack layers was spatially variable and frequently laterally discontinuous. Layers expand and contract horizontally, often creating several different discrete entities of the same layer across a trench face (Fig. 6). The number of layer entities in snow trenches ranged from 5 to 14 in the 5 m trenches and up to 36 layer entities in the 50 m trench (Table 1). Mean snow depth of trenches ranged between 26 and 53 cm for all but trench 6 (79 cm), located in a valley bottom. Consequently, even when considering a positive lidar measurement bias of ∼14 cm, trenches were generally located between the median and lower quartile of the total TVC lidar snow depth distribution (Fig. 5). A coherent layer of surface snow was evident in four trenches (Table 1), consisting of up to 26 % of the mean layer thickness. The proportions of wind slab layers (35 % to 80 %) and depth hoar layers (20 % to 46 %) exhibited a larger range. Figure 7 shows that the mean proportion of depth hoar was consistently just under 30 % of total snow depth. The mean proportion of wind slab was consistently greater than 50 % and showed an increasing trend with increasing total snow depth, indicating that (other factors being equal) where wind slab thickness was greater so was the total depth of the snowpack. A decreasing trend in the mean proportion of surface snow (approximately 25 % to 0 %) with increasing total depth was most likely a result of greater wind erosion and redistribution from the surface where the snowpack was deeper and more wind affected. While interquartile ranges around these trend lines express the natural variability in measured proportional thickness, where total snow depth is known, trend lines made it possible to estimate the percentage of wind slab and depth hoar in a snowpack of unknown microstructure, thereby allowing potential for application of these relationships over larger spatial scales.

Figure 6 Cross section with vertical exaggeration of layer boundaries of individual stratigraphic layers in trench 4 (50 m). Blue lines highlight boundaries between snow and air, surface snow and wind slab, and wind slab and depth hoar. The brown area is subnivean soil or vegetation (two tussocks are labelled). Black lines show boundaries of individual layers aggregated within the wind slab and depth hoar layers.

Figure 7 Relative change in thickness of snowpack layers with total depth: median (solid) and interquartile range (shaded) of surface snow (red), wind slab (blue), and depth hoar (black) layers. Dotted line describes the dependence of layer thickness on depth based on linear trend line fits. Layer thickness data from 2013 trenches only.

Table 1 Descriptive statistics of snow trenches excavated in April 2013. See text for explanation of subnivean roughness metric. For thickness statistics, individual layers were aggregated into three common types: SS represents surface snow, WS represents wind slab, and DH represents depth hoar (nota bene trenches 1–3 were discarded from analysis due to measurement error).

Centimetre-scale variability in tundra snowpack layer thickness was quantified from trench measurements in a spatially distributed manner throughout a catchment for the first time. The range to sill (i.e. horizontal correlation length) of semivariograms (Fig. 8) was used to quantify spatial variability, which varied for all layer thicknesses between 16 and 158 cm (Table 2). While the semivariance of layer thickness in trenches (Fig. 8) changed as a function of absolute layer thickness (Table 1), the mean range to sill of layer thickness increased from depth hoar (45 cm) via wind slab (59 cm) to surface snow (81 cm). The mean range to sill of total snow depths (61 cm) was only slightly greater than that of wind slab, suggesting horizontal variability in wind slab thickness has a strong control over total snowpack thickness. Subnivean roughness, the boundary between snow and the underlying substrate, is likely to have a strong influence on depth hoar thickness. To estimate the importance of this influence, roughness was quantified as twice the value of the root mean square of residuals between the snow–substrate boundary and a linear best fit line to that boundary (Table 1). This provides an estimate of the peak-to-trough amplitude between adjacent subnivean topographic features, often controlled in tundra environments by tussock grasses (e.g. Fig. 6). The roughness metric was calculated across 2 m moving windows. The starting position of each window moved in 1 cm horizontal increments along each trench, and then the roughness of all moving windows was averaged per trench. The 2 m distance was chosen to exceed the maximum measured horizontal correlation length of layer thicknesses while also being broadly representative from visual field inspection of the spacing between tussocks. Across all trenches, roughness of the subnivean boundary ranged from 9 % to 32 % of total trench depth, with a mean of 18 % for all trenches. This is consistent with the premise that depth hoar layer thickness (∼30 % of total snowpack depth) is strongly influenced, but not exclusively controlled, by subnivean roughness.

Table 2 Mean range to sill (cm) for thicknesses of each layer (surface snow, wind slab, and depth hoar) and total snow depth using the stable variogram model.

Figure 8 Semivariograms in individual trenches for (a) total snow depth, and layer thickness for (b) surface snow, (c) wind slab, and (d) depth hoar.

Figure 9 demonstrates relationships between mean percentage thickness of different layers, density, and SSA from a combination of snow microstructural profiles from trenches and pits. Layers were primarily delimited in the field by vertical profiling (visual inspection and hardness). Figure 9a and b show median densities of 104 kg m⁻³ (surface snow), 253 kg m⁻³ (depth hoar), and 316 kg m⁻³ (wind slab); differences in median densities between layers were greater than the 5 %–9 % sampling error associated with gravimetric cutters (Proksch et al., 2016). The interquartile ranges of each layer density did not overlap, indicating clear differences between layer densities, even though there was overlap between full measurement ranges. The upper quartile of surface snow densities overlapped densities of both wind slab and depth hoar, as although it is structurally distinct from the lower wind slab layer, surface snow may have been subject to decomposition, melt, or some wind-packing effects. Additionally, overlap between the lower quartile of wind slab densities and the upper quartile of depth hoar densities resulted from densities in lower sections of wind slab which exhibited the hardness of wind slab yet also microstructural similarity to depth hoar. This is often reported in Arctic tundra snowpacks that undergo strong temperature gradient metamorphism, and it has previously been classified as a unique layer type such as indurated hoar (Sturm et al., 1997; Fierz et al., 2009; Derksen et al., 2014; Domine et al., 2016). Differences between SSA of different layers were more distinct than for densities (Fig. 9c and d); interquartile ranges did not overlap as median SSA increased from depth hoar (11 m² kg⁻¹) via wind slab (24 m² kg⁻¹) to surface snow (45 m² kg⁻¹). Differences in median SSA between layers were much greater than the 10 % measurement uncertainties in SSA (for snow <60 m² kg⁻¹) from IR hemispherical reflectance techniques at 1310 nm (Gallet et al., 2009). Consequently, as the density and SSA values of each layer are nicely separate from each other, it is reasonable to expect that with respect to radar backscatter it is the relative proportions of these snow components, and their attributes, which will drive the radar results. In the next section we explore this using a three-layer radiative transfer model.

Figure 9 Change in snow layer (a) density and (c) SSA with relative thickness of snowpack layers: surface snow (red circle), wind slab (blue cross), and depth hoar (black star). Distributions of snow layer (b) density and (d) SSA: blue box (interquartile range), red line (median), and whiskers (dashed lines) extend from the end of each box to 1.5 times the interquartile range; outliers beyond this range are omitted.

3.2 SWE retrieval accuracy from radar backscatter

Notable differences between distributions of SSA and density in surface snow, wind slab, and depth hoar layers allowed parameterisation of snowpack microstructure in the SMRT model (Table 3). Coupled with fitted relationships between total snow depth and layer thickness as a proportion of total depth (Table 3 and dotted trend lines in Fig. 7), SWE retrieval errors from SMRT simulations were calculated as a function of measured variability in SSA in different layers. SSA measurements of snowpack layers show positively skewed distributions (Fig. 10 a, c, and e); surface snow, wind slab, and depth hoar distributions comprised 64, 77, and 85 averaged layer measurements, respectively. SWE retrieval errors were calculated as truth SWE minus retrieved SWE. SWE retrieval error due to heterogeneity in SSA for SS, WS, and DH layers is presented in Fig. 10b, d, and f for snowpacks of depths between 0.2 and 1 m, and for estimates of layer SSA within 20 % of the known median value.

Table 3 SMRT parameters derived from TVC data. Thickness relationship derived from NIR-derived stratigraphy observed in 2013. Density and fixed SSA taken from median of all observations for that layer type in 2013 and 2018. Fixed SSA values were used in place of the SSA distribution when assessing the impact of spatial variability in other layers.

Figure 10 (a, c, e) Histograms of SSA within each layer. (b, d, f) Synthetic error budget study for a range of snow depths assuming homogeneity in the retrieval. Contour lines show acceptable SSA error as a function of snow depth to remain within a ±30 mm SWE retrieval accuracy limit. White contour line shows perfect retrieval.

Spatial variability in surface snow SSA has a negligible effect on the retrieval error (Fig. 10b). Instead, the retrieval is most sensitive to spatial variability in depth hoar (Fig. 10f), despite only a small range in measured SSA values. A 20 % underestimation in SSA leads to an underestimation in retrieved SWE by as much as 100 mm SWE in 0.8 m deep snowpacks: an underestimation of between one- and two-thirds of total SWE assuming typical tundra snow densities. Overestimation of SSA leads to nearly as large an error in the other direction. Errors in assigning SSA values to wind slab also affect the SWE retrieval, but to a much lesser extent. While a perfect retrieval is possible at all depths (shown by the white contour line in Fig. 10f), an estimated SSA lower than the median is required. For depth hoar, approximately 7 % less than the median is required. For wind slab, closer to 10 % under the median is required.

A radar SWE retrieval algorithm will have constraints on the allowable accuracy. Here, we have set a limit of ±30 mm SWE, indicated as black contour lines in Fig. 10, as the acceptable limit on SWE. This error constrains the setting of SSA values (black contour lines in Fig. 10f), which become increasingly stringent for deeper snow. The constraints in SSA are considerably more liberal for wind slab (Fig. 10d) than for depth hoar as the retrievals are less sensitive to the spatial variability in SSA. Fig. 10f (DH) shows that as snow depth exceeds 0.6 m, the need for accurate values of SSA for depth hoar becomes critical. If the median SSA is used for depth hoar, the maximum SWE error is 40 % of the truth SWE. For depth hoar SSA 20 % above the median, the SWE retrieval error can exceed the actual SWE, particularly for shallow snowpacks. A higher soil backscatter contribution of −10 dB does not change the optimal SSA estimates but does result in more stringent retrieval requirements. For example, at 60 cm snow depth the SSA should be between 89 % and 97 % of the median value to remain within error budget for a soil backscatter of −13 dB. For a soil backscatter of −10 dB, this range is reduced to between 90 % and 96 %. This highlights the need for field observations of backscatter from bare frozen soil to better constrain this value.

Estimates of the single SSA of layers in retrieval scenes will always be lower than the median because of nonlinearity between scattering and SSA. Snow layers with smaller SSA will have a disproportionately larger effect in the truth scene as these microstructure elements will scatter more (optical grain diameter is inversely proportional to SSA). SMRT simulations can compensate for the heterogeneity by assuming slightly smaller SSA in the homogeneous scene. The amount of compensation required could differ as a function of ground contributions, but the representative SSA will always be smaller.

4 Discussion and conclusions

Comparison between the mean snow depth of all trenches (40 cm) and distributions of lidar-derived snow depths demonstrated that trenches were highly representative of snow depths across the whole TVC catchment, incorporating a range of different topographic landscape units. Snow depths in trenches were also consistent with snowpacks found over much wider spatial scales; mean snow depths of 39 and 23 cm for sub-Arctic and Arctic snowpacks, respectively, were reported by Derksen et al. (2014). In addition, trenches at TVC were highly representative of the previously documented complex and often discontinuous layering of tundra and sub-Arctic snowpacks (Benson and Sturm, 1993; Sturm and Benson, 2004; Domine et al., 2012; Rutter et al., 2014, 2016). Consequently, through the coincident combination of spatially distributed measurements of snow microstructure (94 profiles at 36 sites across two winters) and snow layer thickness (9000 vertical profiles across 9 trenches during one winter), we present a unique and robust data set that is likely to be representative of tundra snowpacks in general. The application of these data in a Snow Microwave Radiative Transfer model (Fig. 10), therefore, is of high potential relevance to remote sensing of seasonal snow in the Northern Hemisphere.

Horizontal correlation lengths of layer thicknesses (i.e. distances over which variability in the thickness of individual layers decouples and becomes independent) have not previously been reported. High-resolution layer boundary identification using stitched and georeferenced NIR imagery (Tape et al., 2010; Watts, 2015) now allows the spatial variability in layer thickness to be quantified using semivariograms in analytical approaches similar to assessment of snow depth distributions (e.g. Deems, 2008; Trujillo, 2015). Horizontal correlation lengths (i.e. range to sill of semivariograms) and variance in major layer thicknesses at TVC suggest depth hoar layers vary the most over the shortest distances because of the subnivean vegetation and topography. In tundra environments, especially on topographically exposed plateaus, undulating subnivean topography traps early-winter snowfall in troughs between tussocks. Discontinuities in snowpack layering occur between adjacent ridges and troughs (often caused by tussock grasses), and trapped snow is then subject to strong temperature gradient metamorphism (Colbeck, 1983, 1987), creating large depth hoar crystals. The decreasing spatial variability (increasing distance in horizontal correlation length) of wind slab and surface snow layer thickness reflects the decreasing influence of topography once snow has filled local hollows and the increasing predominance of spatially smoother wind-driven processes such as snow redistribution and compaction. For comparison, the horizontal correlation length of total snowpack thickness, a consequence of variability in all three layers, approximates that of wind slab. Quantification of layer thickness variability using trench measurements suggests a relative hierarchy of variability in commonly occurring tundra snowpack layer types, which all vary significantly at the submetre scale.

Correlation lengths highlight minimum measurement distances over which sampling of snowpack properties should take place. While the intensity of such sampling may have practical limitations, targeted submetre-scale sampling as part of wider ground-based catchment-scale snow measurement campaigns is highly relevant for evaluation of current and future satellite sensors that operate over resolutions on the scale of metres to tens of metres (e.g. Cline et al., 2009; Yueh et al., 2009; King et al., 2018). Rapid acquisition of vertical profiles of snowpack properties using snow micropenetrometers (Schneebeli et al., 1999; Proksch et al., 2015) may increasingly provide the enhanced field measurement capacity required to achieve this spatial sampling resolution.

Correlation lengths can be applied in distributed modelling of snowpack properties using kriging techniques to enable spatial interpolation. Application of the length scales of major snowpack layers, as well as variability of properties within each layer, has the potential for use in catchment-scale models along with knowledge of other landscape elements such as snow drifts. While accurate snowpack modelling at the tens-of-metres scale in tundra environments is challenging (Essery et al., 1999; Clark et al., 2011), snowpack layer correlation lengths could be used to parameterise models so variability in snowpack properties is reliably accounted for when modelling at coarser resolutions. This will become an important parameterisation for future linkages between physical snowpack models and one-dimensional radiative transfer models (Sandells et al., 2017), to better model surface and volume scattering in snow (Zhu et al., 2018), and in observing system simulation experiments (e.g. Garnaud et al., 2019), to assess the performance of future Ku-band radar satellite sensor configurations.

Relationships between layer thicknesses and total snow depth are presented for purposes of simplifying snowpack layer complexity for applications over large scales. Linear trend lines provide a guide to changes in layer proportions for snowpacks up to 1 m deep. Interestingly, the proportion of depth hoar is consistently ∼30 % irrespective of total depth, a trend which continued in the few snow pits in 2018 that were deeper than 1 m. This is somewhat surprising, due to the previously described strong controls of ground surface roughness over depth hoar creation and the potential for large wind slabs to then dominate as the total snow depth increases above the top of ridges. However, the proportion of wind slab layer instead increases at the expense of the decreasing surface snow layer. Consequently, this proportional approach to layer partitioning of snowpacks provides an alternative to applying a maximum depth hoar thickness of 25 cm (King et al., 2018); therefore it may be more appropriate for wider applications across pan-Arctic tundra over multiple winter seasons.

Although rates of compaction, and hence density, of a dry snow layer under its own mass will have a proportional relationship to its thickness (Colbeck, 1972; Anderson, 1976), the weak relationship for fresh surface tundra snow is expected due to other influential aeolian and metamorphic processes. In addition, there is unlikely to be a direct relationship between density and proportional layer thickness for either wind slab or depth hoar, due to wind and thermal gradients dominating processes which control the density of these respective layers. Such processes are currently very challenging to simulate accurately in Arctic snowpacks (Barrere et al., 2017; Domine et al., 2019). So, as distinct differences are evident between distributions of layer properties, using a single median value to represent the density for each layer is appropriate. Differences are even more profound between distributions of SSA in different layers which exhibit low interquartile ranges, especially for depth hoar layers which dominate Ku-band microwave scattering (King et al., 2015). This is particularly important for depth hoar as nonlinear microwave scattering with respect to SSA (SSA is inversely proportional to exponential correlation length or optical grain diameter) means that the smaller SSA will have a disproportionate scattering effect compared to larger SSA. This is demonstrated by synthetic experiments that evaluate the impact of spatial uncertainty in snow microstructural parameters on SWE retrievals, where a single SSA describing a homogeneous scene must be smaller than the median of a heterogeneous scene to compensate for the nonlinear scattering response. Small SSA corresponds to high scattering, so it follows that an underestimation in the SSA estimate (too much scattering) results in less scattering material needed (an underestimation in the retrieved SWE). Whilst the effects of spatial variability in SSA can be countered, the value of representative SSA is critical. Although the SWE retrieval accuracy requirements for depth hoar are more stringent than for wind slab, it is known more precisely than wind slab SSA for the TVC data (Fig. 9). The insensitivity of the retrieval to the spatial variability in surface snow SSA suggests that a median value of SSA will suffice for SWE retrievals. Future work will consider whether a two layer retrieval system could be used to represent the scattering of a three-layer snowpack.

By applying a comprehensive data set of snow microstructural properties collected over two winters at TVC in synthetic radiative transfer experiments, the impact of estimation of unknown snow microstructural characteristics on the viability of SWE retrieval strategies relevant for satellite mission design has been robustly assessed. Spatial variability in microstructure of depth hoar layers dominates SWE retrieval errors. A depth hoar SSA estimate of around 7 % under the median value is needed to accurately retrieve SWE for this snow. In shallow snowpacks <0.6 m, depth hoar SSA estimates of ±5 %–10 % around this value allow retrievals within a tolerance of ±30 mm SWE. Where snowpacks are deeper than around 30 cm, accurate values of representative SSA for depth hoar become critical as the retrieval error will be exceeded if the median depth hoar SSA is applied. Importantly, these experimental results allow potential for uncertainty in SSA in Arctic tundra snow to be used to produce future SWE retrieval quality flags in remotely sensed products and also provide benchmark accuracies for physical snowpack models to deliver SSA estimates. As modelling microwave scattering in large-scale applications has currently been limited to single-layer snowpacks (Takala et al., 2011), the potential for using multilayer information is exciting and progressive, especially when considering the strongly different scattering properties of wind slab and depth hoar (King et al., 2018), the combined influence of which dominates control over microwave scattering in Arctic tundra snowpacks.