Articles | Volume 16, issue 1
The Cryosphere, 16, 87–101, 2022
The Cryosphere, 16, 87–101, 2022

Research article 06 Jan 2022

Research article | 06 Jan 2022

Characterizing tundra snow sub-pixel variability to improve brightness temperature estimation in satellite SWE retrievals

Characterizing tundra snow sub-pixel variability to improve brightness temperature estimation in satellite SWE retrievals
Julien Meloche1,2, Alexandre Langlois1,2, Nick Rutter3, Alain Royer1,2, Josh King4, Branden Walker5, Philip Marsh5, and Evan J. Wilcox5 Julien Meloche et al.
  • 1Centre d'Applications et de Recherche en Télédétection, Université de Sherbrooke, Sherbrooke, J1K 2R1, Canada
  • 2Centre d'études Nordiques, Université Laval, Québec, G1V 0A6, Canada
  • 3Department of Geography and Environmental Sciences, Northumbria University, Newcastle upon Tyne, NE1 8ST, UK
  • 4Environment and Climate Change Canada, Climate Research Division, Toronto, M3H 5T4, Canada
  • 5Cold Regions Research Centre, Wilfrid Laurier University, Waterloo, N2L 3C5, Canada

Correspondence: Julien Meloche (


Topography and vegetation play a major role in sub-pixel variability of Arctic snowpack properties but are not considered in current passive microwave (PMW) satellite SWE retrievals. Simulation of sub-pixel variability of snow properties is also problematic when downscaling snow and climate models. In this study, we simplified observed variability of snowpack properties (depth, density, microstructure) in a two-layer model with mean values and distributions of two multi-year tundra dataset so they could be incorporated in SWE retrieval schemes. Spatial variation of snow depth was parameterized by a log-normal distribution with mean (μsd) values and coefficients of variation (CVsd). Snow depth variability (CVsd) was found to increase as a function of the area measured by a remotely piloted aircraft system (RPAS). Distributions of snow specific surface area (SSA) and density were found for the wind slab (WS) and depth hoar (DH) layers. The mean depth hoar fraction (DHF) was found to be higher in Trail Valley Creek (TVC) than in Cambridge Bay (CB), where TVC is at a lower latitude with a subarctic shrub tundra compared to CB, which is a graminoid tundra. DHFs were fitted with a Gaussian process and predicted from snow depth. Simulations of brightness temperatures using the Snow Microwave Radiative Transfer (SMRT) model incorporating snow depth and DHF variation were evaluated with measurements from the Special Sensor Microwave/Imager and Sounder (SSMIS) sensor. Variation in snow depth (CVsd) is proposed as an effective parameter to account for sub-pixel variability in PMW emission, improving simulation by 8 K. SMRT simulations using a CVsd of 0.9 best matched CVsd observations from spatial datasets for areas > 3 km2, which is comparable to the 3.125 km pixel size of the Equal-Area Scalable Earth (EASE)-Grid 2.0 enhanced resolution at 37 GHz.

1 Introduction

Snow cover is known to be highly variable at the local scale (10–1000 m) due to wind redistribution, sublimation (Liston and Sturm, 1998; Winstral et al., 2013) and vegetation trapping (Sturm et al., 2001). Physical properties of snow such as measurement of stratigraphy (Fierz et al., 2009) can be aggregated into layers, but their spatial distribution is highly variable given their dependence on total depth and surface roughness (Liljedahl et al., 2016; Rutter et al., 2014). Such variability leads to uncertainties in the retrievals of snow state variables such as snow water equivalent (SWE) using microwave remote sensing from local scales (King et al., 2018; Rutter et al., 2019) to global scales (Pulliainen et al., 2020). Improving our empirical understanding of the processes governing this variability would improve spaceborne snow monitoring, especially in Arctic regions where ground measurements and weather station networks are sparse.

Measurement of SWE using passive microwave satellite data (Larue et al., 2018; Pulliainen, 2006) is possible using a radiative transfer model to simulate snow emission at various frequencies, from which an inversion of the model can produce global estimates of snow depth (Takala et al., 2011). More specifically, passive microwave brightness temperatures (TB) are governed by radiometric properties of the layered snowpack. As such, each layer has its own absorption and scattering properties; the amount of scattering is proportional to snow total mass where the scattering and emission is frequency-dependent (Kelly et al., 2003). Scattering at higher frequencies such as 37 GHz, will lead to lower TB, so differences between TB at two frequencies (37–19 GHz) are related to snow mass (Chang et al., 1982). Arctic snowpack mainly consists of two distinct layers (wind slab and depth hoar), where each layer has unique scattering properties (Derksen et al., 2010). Complexity of the layered properties (density, temperature and microstructure) strongly influences radiative transfer modeling (King et al., 2015; Rutter et al., 2014). Furthermore, recent developments in radiative transfer modeling (SMRT: Picard et al., 2018; DMRT: Tsang et al., 2000; and MEMLS: Wiesmann and Mätzler, 1999), microstructure representation (Royer et al., 2017) and in situ measurement of snowpack properties (Gallet et al., 2009; Montpetit et al., 2012; Proksch et al., 2015) have provided significant agreement between models and in situ measurements. However, spatial distribution and heterogeneity of total snow depth and stratigraphy remain challenging to implement and are not considered for large-scale monitoring of SWE in tundra environments. Rutter et al. (2019) and Saberi et al. (2020), using three- and two-layer models respectively, demonstrated a relationship between the ratio of depth hoar and wind slab with respect to total depth, enabling the usage of the proportion of these two layers with total snow depth. Working with a simplified layer representation of a snowpack with well-defined physical properties may adequately characterize snowpack for large-scale SWE retrievals.

Two dominant processes governing snow depth variability in the Arctic are (1) wind redistribution with topography (Sturm and Wagner, 2010; Winstral et al., 2002) and (2) vegetation trapping (Domine et al., 2018; Sturm et al., 2001). Liston (2004) described snow depth heterogeneity using a log-normal distribution with a coefficient of variation of snow depth (CVsd), the ratio between standard deviation (σsd) and the mean of snow depth (μsd), indicating the extent and spread of a distribution (i.e., high variability over thin snow will lead to high values of CVsd). Also, Liston (2004) proposed nine categories of CVsd with values ranging from 0.9 to 0.06 for midlatitude treeless mountains to ephemeral snow, where arctic tundra type was 0.4. Snow depth variability is based on a parametrization of μsd and CVsd on the log-normal distribution scale parameters (λ, ζ). Gisnås et al. (2016) adapted that approach using scale parameters (α, β) of the gamma distribution. In all cases, CVsd is used to describe sub-grid variability (Clark et al., 2011), but its value remains challenging to quantify given that regional trends are linked to topography, vegetation and climate (Winstral and Marks, 2014). In this context, CVsd is used to quantify spatial heterogeneity of snow in climate modeling but so far has not been used in microwave SWE retrievals.

In SWE retrievals, snow depth is assumed to be uniform, and the mean depth is used to optimize brightness temperature and derive SWE from depth and assumed density (Kelly, 2009). There is potential for CVsd to be used as an effective parameter to estimate sub-pixel variability in brightness temperature. Bayesian frameworks are used in inversion schemes for SWE retrievals (Durand and Liu, 2012; Pan et al., 2017; Saberi et al., 2020) using a priori information (density, microstructure and temperature) from regional snowpack characteristics and inversion of radiative transfer models (Saberi et al., 2020). An iterative approach based on Bayesian theory is used (Takala et al., 2011) to match observed brightness temperature with modeled brightness temperature by iterating a priori information of the snowpack in order to derive snow depth and SWE. Saberi et al. (2020) conducted a case study for snow depth retrievals using a two-layer model from airborne microwave observations using a Bayesian framework (or Markov chain Monte Carlo) over tundra snow. However, high uncertainty (21.8 cm) in retrieved snow depth (via TB) resulted, which suggested the use of a term involving variation in snow depth and microstructure within the footprint instead of a uniform snow depth.

To address this research gap, we used a multi-year snow dataset from two Arctic locations to quantify sub-pixel variability of snow depth and microstructure and used CVsd as an effective parameter that controls snow sub-pixel variability. Firstly, we evaluate tundra snow depth spatial variability using probability density functions (log-normal and gamma) and its parameters, μsd and CVsd. Secondly, we present from in situ observations distinct snow microstructure and density values of both tundra main layers (depth hoar and wind slab), mean ratios of layer thickness, and the depth hoar fraction (DHF) relative to snow depth. Finally, we perform a Gaussian process fit to estimate DHF from snow depth and used probability density functions of snow depth to add variation of snow depth and microstructure within the footprint. Then we compare mean pixel snow properties simulations with simulations of sub-pixel variation in snow properties to evaluate biases between measured TB from a satellite sensor at 19 and 37 GHz and TB simulated by inversion of a radiative transfer model.

2 Methods

2.1 Study site

Data were collected in two regions of the Canadian Arctic, with different topography and vegetation yielding different snow depth distributions. Trail Valley Creek (TVC) research watershed, Northwest Territories (6844 N, 13333 W), located at the southern edge of arctic shrub tundra, is dominated by herbaceous tundra and dwarf shrubs and characterized by gently rolling hills with steep slopes. Greiner Lake watershed, Cambridge Bay (CB), Nunavut (6913 N, 10453 W), located within arctic tundra, is characterized by dwarf shrub and calcareous tills on upland sites with gently rolling hills and small ponds and lakes. TVC is considered to have more subarctic attributes with predominant vegetation than CB given its proximity to the northern edge of the boreal forest. Topographic maps (Fig. 1; ArcticDEM) show slightly higher variation in elevation at TVC with plateau and steep slopes compared to CB, which is dominated by ponds and small variation in topography.

Figure 1Locations of study areas in the Canadian Arctic, Cambridge Bay and Trail Valley Creek site. Grid shown is the enhanced 3.125 km EASE-Grid 2.0 (Brodzik et al., 2018) used for satellite data. The ArcticDEM used to show elevation has a 2 m resolution (Porter et al., 2018) derived from stereo high-resolution visible imagery for the entire Arctic domain, which is freely available.

2.2 Data

Snow pits (315) at each site (TVC: 68, CB: 248) provided information on snow layering, vertical profiles of geophysical properties (includes temperature, grain type classification, hardness, density, microstructure and depth). Measurements of visual stratigraphy and grain type classification were conducted following Fierz et al. (2009). Density was measured using 100 cm3 density cutters and digital scales. Snow specific surface area (SSA) was measured using an InfraRed Integrating Sphere (IRIS) (Montpetit et al., 2012) in Cambridge Bay and an A2 Photonic Sensors IceCube in TVC, both based on 1300 nm laser reflectometry (Gallet et al., 2009). Snow depth measurements, linear transects and circular transect around snow pits used a magnaprobe from SnowHydro LLC (Sturm and Holmgren, 2018), which is equipped with a standard GPS unit. Measured snow depth distributions were used to identify subsequent pit locations (on site) from a predefined transect across CB watershed in order to ensure the snow pit locations were representative of wider spatial variability (Table 1). For TVC, pit locations were chosen based on previous snow depth distribution (2016), slope and elevation. Multiple snow depth maps at 1 m resolution from remotely piloted aircraft system (RPAS) surveys derived from photogrammetry conducted in March 2018 (Walker et al., 2020a) were used to estimate snow depth distribution in TVC with total spatial coverage of 5.3 km2. An airborne lidar dataset of TVC snow depths (93 km2 at 10 m resolution) collected by an aircraft in April 2013 (Rutter et al., 2019) was also used. Monte Carlo simulations of both the μsd and CVsd were performed on each snow depth map. Simulations randomly selected pixels as the center of a circular mask with a random radius. The mask was used to select all pixels within the circle so the statistical parameters (μsd and CVsd) could be calculated. Also, a small RPAS map was available for CB with spatial coverage of 0.2 km2 at 1 m resolution. Maps of normalized difference vegetation index (NDVI) were created from Sentinel-2 (10 m resolution) images from late summer (1 September 2019 for TVC and 8 September 2019 for CB).

Table 1Summary of the number of snow depth measurements (magnaprobe and RPAS) and snow pit sites per year. The availability of SSA and density measurements across sites and years is also noted (×). See Table 2 for full dates.

Download Print Version | Download XLSX

2.3 Measured brightness temperatures and Snow Microwave Radiative Transfer (SMRT)

Microwave TB values were used to evaluate simulations from SMRT at 37 and 19 GHz from the Special Sensor Microwave/Imager and Sounder (SSMIS) sensor, with the Equal-Area Scalable Earth (EASE)-Grid 2.0 resampled at 3.125 km (6.25 km for 19 GHz) resolution (Brodzik et al., 2018), for both TVC and CB regions. TB values were estimated by averaging all pixels within snow pit area (CB: 24 pixels, TVC: 8 pixels for 37 GHz). Each pixel with at least one snow pit inside was used. Since all snow pits were aggregated to obtain mean value and distribution of snow properties for SMRT, averaged TB covering the snow pits area was used.

The area was also filtered to remove any contribution from sea or deep lakes, as pixels with liquid water exhibit large biases even if the signal at 37 GHz is mostly sensitive to snow (Derksen et al., 2012). For CB, an area with the same spatial coverage but a slightly different location was used since the snow pit area was within 25 km (full resolution of SSMIS) from the ocean. TB values were temporally averaged to match times of field measurements, representing peak winter snow accumulation (Table 2). Also, TB values were corrected for atmospheric contributions using the linear relation with precipitable water from the 29 atmospheric NARR layers (Vargel et al., 2020; Roy et al., 2013).

A multi-layered snowpack radiative transfer model (SMRT, Picard et al., 2018) was used to simulate snow emission at 19 and 37 GHz. Model inputs are snow temperature, density and microstructure of each snow layer. Correlation length of snow microstructure in each layer was estimated from mean density and SSA measurements of each layer when available (wind slab: WS; depth hoar: DH) using Debye's equation scaled by a factor (κ37=1.39,κ19=1.71) for arctic snow as suggested by Eqs. (3b) and (4) in Vargel et al. (2020) with the improved Born approximation (IBA-Exp) configuration. Soil emission was simulated using the Wegmüller and Mätzler (1999) model with permittivity and roughness values from a field study of frozen soil emission based in CB (Meloche et al., 2020). The soil parameters from CB (Meloche et al., 2020) closely match values from a study in TVC (King et al., 2018) and were used for both site simulations. The lakes in CB shown in Fig. 1 were not considered in the soil emission contribution because most of the water was frozen (ε 4–6) (Mironov et al., 2010), which had a similar permittivity to frozen soil (ε 2–4) (Mavrovic et al., 2021) than liquid water (ε>25).

The basal layer temperature was set to the mean soil–DH interface measurements from snow pits of each site. The temperature of the WS layer was estimated from the North American Regional Reanalysis (NARR) air surface temperature, which closely matched snow pit surface layer temperature. NARR air surface temperatures were used because they provide a global estimate that matches spatial coverage of the EASE-Grid 2.0, which is continuous (spatially and temporally) compared to the sparse snow pit observations.

Table 2Summary of mean basal and air surface temperatures for SMRT simulations, precipitable water (PWAT) used for atmospheric correction, and measured (corrected) TB at both vertical (V) and horizontal (H) polarization by the SSMIS sensor (platform F18).

Download Print Version | Download XLSX

2.4 Gaussian processes

Gaussian processes (GP) are a non-parametric Bayesian method used in regression models. These processes are effective and flexible tools to fit complex functions with small training datasets (Quiñonero-Candela and Rasmussen, 2005). Gaussian processes provide uncertainties on predictions, using training data and prior distributions to produce posterior distributions for predictions. Mean (m(x)) and covariance (k(xx)) functions from the multi-variate Gaussian distribution are used to fit data (x: snow depth, y: ratio of layer, DHF). The m(x) function describes the expected value of the distribution, and the k(xx) describes the shape of the correlation between data points (xi). Different mean and covariance kernels can be chosen to fit the data. From the Bayes rule in Eq. (1), where y (DHF) and X (snow depth) are observed data and f the GP function, posterior predictions of DHF can be produced. Posterior predictions were calculated using the standard method of Markov chain Monte Carlo (MCMC) sampling using PyMC3 (Salvatier et al., 2016).


Equation (2) defined f as a function of mx,kx,x. A mean function m(x), following an inverse logic function (ϕ) (Eq. 3), was chosen due to the close fit with observations. The covariance function kx,x is a Gaussian white noise covariance function and is defined with noise (σ) and the Kronecker delta function (δx,x) (Eq. 4), to best fit the observations. By using a scaling function (ϕ), the covariance function (uniform noise in this case) can be modified as a function of x. The scaling function used is also an inverse logic function (ϕ) that takes the same form as Eq. (3). Finally, a deterministic transformation is applied to the prior (GP) to constrain values to a ratio (0, 1). The likelihood of DHF observations is defined by a Beta distribution (0, 1).

3 Results

3.1 Snow depth distribution

Distributions of snow depth are needed when integrating over large areas to calculate sub-grid snow variability for distributed models (Clark et al., 2011; Liston, 2004). The μsd and the CVsd of snow depth are used as parameters in probability density functions to estimate the shape of the log-normal and gamma distributions. To find which distribution best fits the depth observations, we tested the log-normal and gamma distributions using the Kolmogorov–Smirnov two-sample test with snow depth observations (shown in blue in Fig. 2). The statistical fits for each distribution are shown in Table 3. For both the log-normal and gamma distributions the null hypothesis is validated at the 5 % significance level from p value > 0.05 (i.e., the two samples were drawn from the same distribution), which agrees with previous assessments of Arctic snow (Clark et al., 2011; Gisnås et al., 2016).

Table 3Kolmogorov–Smirnov (KS) test for two samples of probability distribution function (PDF).

Download Print Version | Download XLSX

Figure 2Log-normal and gamma distribution fit to the measured snow depths.


Distributions with parameterization using measured μsd and CVsd (Fig. 2) differ from the best fit with regular parameters, especially compared with the log-normal distribution in CB (black dashed line in Fig. 2b). Liston (2004) reported a CVsd of 0.4 for Arctic tundra snow, which is in close agreement with the values of 0.43 for TVC and 0.56 for CB. These values were obtained from spatially distributed snow depth measurements around snow pits. For comparison, maps of snow depth from RPAS for TVC (n= 6 325 365 with total spatial coverage of 5.3 km2) showed a much larger CVsd= 0.78 than magnaprobe data (n= 15 731) with CVsd= 0.43 (Table 4). A RPAS dataset was also available for CB but with a much smaller spatial coverage (0.2 km2) showing a CVsd of 0.49. In Figs. 3 and 4, we investigated the relationship between spatial coverage of sampling and the CVsd parameter. Datasets include RPAS-derived data at TVC (TVC18-RPAS) containing six areas with various size from 1–3 km2, a CB18-RPAS map of 0.2 km2 and the larger lidar-derived snow map in TVC (TVC13-Lidar). Figure 3a shows snow accumulation of TVC13-Lidar and TVC18-RPAS with snow drift visible in dark blue, and a sub-grid of 1 km2 showed areas with high CVsd (Fig. 3b) containing more drift. For both areas, 500 Monte Carlo simulations were performed by randomly selecting subregions within each domain (Fig. 4) so the mean and variability as a function of coverage could be investigated. Simulations showed subsampling of μsd and CVsd converged to the values of the full area. The mean of each area was similar in value with less variation in the simulations compared to CVsd. A difference of 0.2 between the full CVsd of the RPAS (5 km2) and lidar (93 km2) maps (Fig. 4) was found. Also included in Fig. 4 are in situ CVsd estimates with variable high-density sampling (magnaprobe) over different spatial extents at Daring Lake, NWT (Derksen et al., 2009; Rees et al., 2014); Puvirnituq, QC (Derksen et al., 2010); and at Eureka, NU (Saberi et al., 2017). The two points at the limit coverage scale correspond to areas of respectively 625 km2 (CVsd= 1, Daring Lake site; Chris Derksen, personal communication, 2009) and 198 km2 (CVsd= 0.89, Eureka site; Saberi et al., 2007).

Table 4Statistical parameters of snow depth distributions. Bold font represents multi-year average.

Download Print Version | Download XLSX

Figure 3RPAS and lidar dataset of snow depth at TVC (TVC13-Lidar and TVC18-RPAS). TVC13-Lidar is the largest dataset covering 93 km2. TVC18-RPAS is a smaller dataset within the area of TVC13-Lidar. Panel (a) shows the snow depth map at 10 m resolution from 2013. Panel (b) shows a sub-grid of 1 km with CVsd for each cell; (c) same as panel (b) but for μsd.


Figure 4Snow depth mean (μsd(a, b) and variability (CVsd) (c, d) as a function of coverage for sampling area: (a, c) small area; (b, d) large area. Monte Carlo simulations were done using the two datasets in TVC. CB18-RPAS was also added in panels (a) and (c) because of the similar coverage. The μsd and CVsd of both full areas are shown by the black dotted and dashed line.


3.2 Analysis of SSA and density per layer

After combining measurements from all snow pits at TVC and CB (n= 315) the mean proportion of DH layer thickness was 46 % and WS was 54 %. The goal was to classify DH as large grain snow (large facets, depth hoar cups and chains) and then all other snow layers above the DH as wind slab (WS). Some layers were more difficult to classify as they contained mixed crystals or were a transitional slab-to-hoar layer (also referred to as indurated hoar) (Derksen et al., 2009). Slab that contained small faceted crystals (< 2 mm) were classified as WS. Indurated hoar, a wind slab metamorphosed into depth hoar, was classified into DH with a typical density ∼300 kg m−3. For this reason, the peak of each distribution appeared close to each other in Fig. 5c and d. For retrieval of snow properties using satellite remote sensing, a two-layer model using WS and DH can be used to simplify much of the layer complexity found in arctic snowpacks (Rutter et al., 2019; Saberi et al., 2017). A small amount of surface fresh snow (SS) was present in some pits but was not included in this calculation as this type of snow was a short-lived layer, combining fresh precipitation that rapidly transformed into rounded grains due to destructive metamorphism and defragmentation by wind. Distributions of SSA are more distinct between layers than density (Fig. 5a and b); see Rutter et al. (2019). Figure 5c and d show that the mean values for density of WS (335 kg m−3) and DH (266 kg m−3) were closer together. SSA distributions also showed a gap between both mean values (WS: 19.7 m2 kg−1; DH: 11.1 m2 kg−1) (Fig. 5, Table 5). Even if snow properties can show high heterogeneity at local scales, simple distributions approximate this variability well. Temporal (year) and spatial (regional between site) variation is low, and snow properties (density and SSA) can be approximated by a distribution for each distinct layer, WS and DH as in Fig. 5. Therefore, snow properties were simplified in distributions for each layer (WS and DH) representing tundra snow.

Figure 5SSA and density variability of surface snow (SS), wind slab (WS) and depth hoar (DH) for the two studied sites (TVC and CB) and different dates (see Table 5). In panel (c) the best log-normal fit distribution is shown in black; (d) same as panel (c) but for the normal fit distribution. In panels (c) and (d), the kernel density estimates (KDEs) of the histogram of each layer are also shown (in color).


Table 5Parameters for best fitting distribution of SSA and density for layers of DH and WS.

Download Print Version | Download XLSX

Parr et al. (2020) found a key threshold of μsd+1σsd to define snow drifts in tundra environments. This threshold (> 0.6–0.8 m), based on data presented in Table 4, is an important metric in Fig. 6 since above this depth the variability and the mean DHF are greatly reduced as the snowpack is dominated by wind slab for larger depth (drift). As defined in Parr et al. (2020), the transported snow from wind accumulates at these particular locations (drift) where it was scoured from wind affected area yielding lower depth with high DHF.

Figure 6(a) Depth hoar fraction (DHF) as a function of total depth for snow pit data from 2015–2019 in Cambridge Bay and 2018–2019 for Trail Valley Creek. Both datasets were separated in equal bins (10 cm) to estimate the mean value shown with the dashed line. The black line represents the mean for both sites with the 95 % confidence interval. (b) DHF is shown as a function of NDVI from the snow pit area with the mean DHF and NDVI per site shown by dashed lines and the Gaussian distributions of DHF by the solid lines.


Vegetation also strongly influences variability of DHF in shallower snowpacks, where arctic shrubs and tussocks promote depth hoar formation (Domine et al., 2016; Royer et al., 2021; Sturm et al., 2001). However, there is no clear link between DHF and NDVI (a proxy for vegetation type) at local scales (Fig. 6b). Since shrubs provide shelter to snow up to their own height (Gouttevin et al., 2018), vegetation height rather than type would be required. However, at the regional scale, differences were evident between both regions, where mean NDVI and DHF are greater at TVC (NDVI = 0.5, DHF = 0.54) than CB (NDVI = 0.27, DHF = 0.38). This finding is in agreement with Royer et al. (2021) over a northeastern latitudinal gradient, showing that sites with shrubs and tussocks have a greater DHF than those without.

3.3 DHF predictions using snow depth with Gaussian processes

The impact on microwave scattering of variability of layer microstructures with snow depth was previously accounted for in Saberi et al. (2020) by defining two categories, a high scattering thin snow layer (high DHF) and a thicker self-emitting layer (low DHF). Instead, using GP, DHFs were fitted and predicted based on snow depth values (Fig. 7). In order to use GP, the mean function m(x), following an inverse logic function (ϕ1: Eq. 3), was chosen with parameters a=-5, x0=0.6, b=0.35 and c=0.2 to best match the mean line observation for both sites in Fig. 6. The mean function set the mean value across the snow depth range. The correlation function was set to a uniform noise, but this noise was reduced from depth > 40 cm by using a scaling function (ϕ2:a=-5, x0=0.6, b=1.5 and c= 0.25). An inverse logic function (ϕ1, ϕ2) was used twice in the fitting (1) for the mean value and (2) to reduce the variability (noise) as snow depth increased. The snow pit dataset (n= 315, Fig. 6) was used to build posterior predictions using MCMC sampling.

For prediction of DHF, any number of snow depths can feed into the posterior prediction or GP fit. Snow depths were generated from a log-normal distribution with parameters (μsd, CVsd) from the previous section in Table 4. Posterior predictions of DHF were similar to observed data (Fig. 7) and followed closely the posterior probability representation in red (GP fit). Again, higher variability in DHF was reproduced for depths < 0.5 m, which was then reduced for depths > 0.5 m following the red posterior prediction representation in Fig. 7.

Figure 7Prediction on DHF (cyan) using a GP fit trained on observed data (black). Snow depth is from samples from a log-normal distribution with parameters from Table 4. The GP fit is illustrated in red, where darker red represents a high posterior probability that follows the mean function.


3.4 SMRT simulation of sub-grid variability within sensor footprint

SMRT simulations using measured snowpack properties were compared with the satellite measurements of TB. Two simulations were evaluated using (1) mean measured depth, each layer's density and SSA, and DHF and (2) a log-normal distribution of snow depth and the GP fit (predicted DHF). We hypothesized that the EASE 2.0 grid pixels can be separated into n smaller sub-grid pixels. Sub-grid pixels (n= 500) represent the observed snow variability, where n snow depths will follow a log-normal distribution with parameters μsd and CVsd. The ratio of each layer is predicted using the GP fit with depth as input from the log-normal distribution. Mean SSA (DH: 11 m2 kg−1; WS: 20 m2 kg−1) and density (DH: 266 kg m−3; WS: 335 kg m−3) per layer were determined from measurements (Fig. 5).

For one standard EASE-Grid 2.0 pixel, a distribution of sub-grid TB was simulated to reproduce a realistic distribution of TB within the radiometer footprint. This variability was derived from spatially distributed observations from snow pits and snow depths observation. Snow depths followed a log-normal distribution with the mean measured depth (μsd) of each region (Table 4) and a depth variability (CVsd) that was evaluated from a range of 0.1 to 1. The GP mean function from Fig. 7 was used to predict the DHF for each region. When using CVsd= 0.7 , the simulated distribution showed a wide sub-pixel variability (±40 K) with a mean value of TB37V= 194.7 K (blue line in Fig. 8a), which is very close to the satellite-measured TB37V of 196.5 K (green dotted line in Fig. 8a). In this case, the TB value simulated from the mean measured snow depth and mean DHF was slightly lower (190.7 K, i.e., a bias of 5.8 K) (black dotted line in Fig. 8a). To represent the signal measured by the sensor, the mean of the simulated TB was chosen, and it was assumed that the sub-pixel effect combined linearly at this scale. Because the simulated TB37V distribution was not exactly a normal distribution, it appeared that the mean TB of this distribution increased when CVsd increased (Fig. 8b). This meant that snow depth variability (CVsd) must be accounted for when estimating the average TB at 37 GHz, in addition to the mean snow depth values. The influence of the GP simulation on the mean simulated TB37V was approximately 10 K (Fig. 8b) as CVsd varies from 0.1 to 1. The addition of snow variability in simulation (Fig. 8c and d) of 19 GHz has negligible effect on TB19 and shows a constant simulation across the CVsd range of 0.1 to 1. Simulation of TB19 shows higher biases at horizontal polarization than at vertical polarization.

Figure 8Brightness temperature variability simulation distribution (a, c) of simulated TB within a pixel, where vertical lines represent the mean of this distribution for V pol (blue), measured by satellite (green), and TB value simulated from the mean measured snow depth and mean DHF (black). In panels (b) and (d), the mean of the simulated TB for H pol (red) and V pol (blue) as a function of CVsd with mean values (dotted black lines). The CVsd that minimized biases is located at the red/blue–green intersection. Shaded blue and red areas correspond to a 2σ range representing uncertainty inherent from our Bayesian simulations in estimating the mean of simulated TB for the pixel.


GP simulation reduced biases by 5 K with a higher optimized CVsd (intersection of red/blue–green line, Fig. 8b). A similar pattern was observed for CB (not shown here), but the measured TB at CB was much higher than the GP simulation, resulting in a large bias for CB ( 20 K) compared to TVC (Table 6). Both sites suggested a larger CVsd, which agreed with a CVsd for larger spatial coverage measured in Fig. 4. Observed large biases at CB vary over the years from 5 to 29 K. The total RMSE of both sites and years linearly decreased as a function of CVsd (Fig. 9). Total RMSE is minimized with higher CVsd (0.8–0.9) typical of a large sampling scale (over 4 km2) as shown in Fig. 4.

Table 6Bias between SMRT simulated and measured TB from the SSMIS sensor at each site.

Download Print Version | Download XLSX

Figure 9Overall RMSE (year and site) with the mean simulation (dashed black line) and the GP simulation in blue as a function of the coefficient of variation of snow depth.


4 Discussion

We show that as spatial coverage increased, the CVsd parameter converged to the full area values (Fig. 4). Monte Carlo simulations of snow depth distribution and coverage show high variation in CVsd (from 0.1 to 2) for areas < 10 km2. Snow accumulation varied at the mesoscale (100 m to 10 km) due to topography and vegetation (Pomeroy et al., 2002) by varying wind-flow direction (Liston and Sturm, 1998). At the mesoscale, variability in CVsd was high due to topographic differences; plateau, slope and valley create favorable conditions from wind-flow direction to promote snowdrift, scour and sublimation processes (Parr et al., 2020; Rutter et al., 2019). Vegetation facilitates snow holding capacities by decreasing wind speed near the ground within and downwind of shrub (Marsh et al., 2010; Sturm et al., 2001). Some areas include both extreme drifts and thin snow, resulting in high CVsd (dark green areas in Fig. 3b), which are commonly found in TVC (Walker et al., 2020a). CVsd was lower for areas without drifts (light green areas in Fig. 3b). In areas >10 km2  (Fig. 4d), variation in CVsd was reduced and yielded higher values > 0.6.

Convergence to higher CVsd as spatial coverage increased matched the PMW optimized values found in this study using GP simulation (0.8–1.0). Our analysis in Fig. 4d showed that CVsd of TVC13-Lidar converged to 0.6 at 93 km2, but two in situ points from other studies at 625 km2 had higher CVsd (0.9–1) due to larger coverage or different site characteristics. This indicates that a CVsd between 0.6–1.0 is more suitable to represent snow depth variability in SWE retrievals for PMW SWE products at 25 km for the EASE-Grid 2.0 and 625 km2 for GlobSnow 3.0 (Pulliainen et al., 2020). For active sensors (resolution < 1 km), the high variability in CVsd under 1 km2 due to high variation in snow depth (Fig. 4b) can affect back scattering since the active sensor at the Ku band is also sensitive to volume scattering (King et al., 2018). The need for prediction of μsd and CVsd based on topography could become essential at these scales not only for microwave remote sensing but also to improve snow modeling or land data assimilation (Kim et al., 2021).

Spatial complexities of Arctic snowpacks can be adequately characterized with distributions of snow depth (Fig. 2) and simplified by considering density and SSA of two main layers (Fig. 5). Such simplifications could be potentially useful for satellite SWE retrievals across Arctic tundra regions. Since Bayesian SWE optimization needs a strong first guess from regional a priori information, multiple distributions of snow depth, density and SSA presented here can be used for tundra type snow in MCMC sampling (Pan et al., 2017; Saberi et al., 2020). Additionally, a similar approach to our GP simulation can be added so the CVsd parameter can also be used as a priori information with a distribution from 0.8 to 1, since it improved TB RMSE by ∼8 K (Fig. 9). This approach improved TB simulation compared to using only mean values of snowpack properties by adding variability within the footprint. The CVsd parameter (describing variation in snow depth) has a considerable effect on brightness temperature (10 K) when used as an effective parameter to account for sub-pixel variability of snow depth. The amount of scatterers (snow grain and structure) within the radiometer's footprint is adjusted via the DHF predicted from snow depth (CVsd). The relationship found in Fig. 6 used to predict DHF (Fig. 7) could also be used deterministically with the mean function (ϕ1) or with a linear relation of DHF decreasing from 50 % to 20 %. However, the Bayesian Gaussian process was used because SWE retrievals are currently implemented in a Bayesian framework (Takala et al., 2011).

Considering that the difference between 19 and 37 GHz is used in SWE retrievals (Takala et al., 2011), using the CVsd to account for variability of scatterers only affected simulation of 37 GHz with a weak effect on 19 GHz (Fig. 8). If standard deviation of snow increases (more drift), then relatively fewer large scatterers from depth hoar are present within the footprint due to a low DHF generally observed in large drifts. The net result is then an increase in TB at 37 GHz resulting from an increase in CVsd (Fig. 8).

This idea of modulating the amount of scatterers based of DHF prediction and a distribution of snow depth (μsd and CVsd) can be extended to a future active Ku-band mission (Garnaud et al., 2019; King et al., 2018) as it is known that microwave spatial variability affects backscatter signal (King et al., 2015) and SWE retrievals (Vander Jagt et al., 2013). The CVsd parameter is proposed as an effective parameter to account for variability inside the grid cell, while the mean depth (μsd) is assimilated by in situ measurements at weather stations in data assimilation schemes (Takala et al., 2011) or by the physical snow model (Larue et al., 2018). The CVsd could be optimized or predicted using relationships with spatial coverage (Fig. 4) and statistical topographic regression (Grünewald et al., 2013). Future works would need a dataset covering a large area where μsd and CVsd could be investigated with topography in smaller sub-areas.

5 Conclusion

This study evaluated the use of parameters controlling Arctic snow depth distributions to improve passive microwave SWE retrievals by characterizing tundra snow sub-pixel variability. In shrub and graminoid tundra environments, mean values of snow depths (μsd= 0.33–0.44 m) and coefficient of variations (CVsd= 0.4–0.8) were similar to those previously reported in Arctic tundra (Derksen et al., 2014; Liston, 2004; Derksen et al., 2009). Monte Carlo simulations were applied to investigate μsd and CVsd as a function of spatial coverage. An increase in CVsd matched increased spatial coverage of snow depth sampling, indicating that a higher CVsd (0.6–0.9) is more suited to estimate snow depth variation at the 3.125 km resolution EASE-Grid 2.0. Also, simulations show high variation in CVsd (> 0.9) for areas < 10 km2, suggesting a need for topography-based prediction of μsd and CVsd at this scale. The CVsd is shown to be an effective parameter to account for snow depth variability in simulation of snow TB. A two-layer snowpack model (depth hoar and surface wind slab), which simplifies snowpack properties into distributions, was used to initialize the SMRT model via a GP fit of the DHF related to snow depth. DHF is fitted to snow depth using a Bayesian Gaussian process, which accounts for variation in snow scattering using CVsd. SMRT simulation was used successfully to simulate satellite TB, but there are still substantial uncertainties in the simulated values which are likely linked to microstructural properties not captured by SSA (Krol and Löwe, 2016). SMRT simulations of TB were reduced by 8 K after optimizing CVsd to higher values (0.8–1.0), thereby matching CVsd of spatially distributed snow depth from TVC18-RPAS accounting for variation in snow properties inside the footprint of satellite sensor. The CVsd parameter is proposed as an effective parameter to account for variability inside the footprint to minimize the difference between microwave measurements and simulations in SWE retrievals algorithm. This would be beneficial to the data assimilation scheme of the European Space Agency GlobSnow product (Takala et al., 2011) and to model the large-scale climate trend of tundra snow (Mortimer et al., 2020; Pulliainen et al., 2020).

Code and data availability

Data and code for the Gaussian process fit and GP simulation are available at (last access: 27 December 2021;, Meloche, 2021). The RPAS map and magnaprobe from TVC are available at (Walker et al., 2020b).

Author contributions

JM led the formal analysis, performed the investigation and wrote the original draft. AL participated in the writing, review, and editing; participated in the supervision; participated in the investigation; and led the funding and resource acquisition. NR participated in the writing, review, and editing; supervision; and investigation. AR participated in the writing, review, and editing; supervision; and investigation. JK and BW participated in the writing, review, and editing as well as the data acquisition. PM and EW participated in the data acquisition.

Competing interests

At least one of the (co-)authors is a member of the editorial board of The Cryosphere. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors would like to thank the entire GRIMP research team from Université de Sherbrooke and CHARS staff team for field work assistance from 2015–2019 in Cambridge Bay. We thank Chris Derksen and Peter Toose (Environment and Climate Change Canada) for logistics and field work at Trail Valley Creek Research Station. We also thank Victoria Dutch from Northumbria University for her help with dataset management of TVC. The authors would like to thank the reviewers and the editor, Carrie Vuyovich, for their help in improving this paper.

Financial support

This research has been supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), Polar Knowledge Canada, the Canadian Foundation for Innovation (CFI), Environment and Climate Change Canada (ECCC), Fonds de recherche du Québec – Nature et technologies (FRQNT), the Northern Scientific Training Program (NSTP), and research funding from Northumbria University, UK.

Review statement

This paper was edited by Carrie Vuyovich and reviewed by Matthew Sturm and one anonymous referee.


Brodzik, M. J., Long, D. G., and Hardman, M. A.: Best practices in crafting the calibrated, Enhanced-Resolution passive-microwave EASE-Grid 2.0 brightness temperature Earth System Data Record, Remote Sens., 10, 1793,, 2018. 

Chang, A. T. C., Foster, J. L., Hall, D. K., Rango, A., and Hartline, B. K.: Snow water equivalent estimation by microwave radiometry, Cold Reg. Sci. Technol., 5, 259–267,, 1982. 

Clark, M. P., Hendrikx, J., Slater, A. G., Kavetski, D., Anderson, B., Cullen, N. J., Kerr, T., Örn Hreinsson, E., and Woods, R. A.: Representing spatial variability of snow water equivalent in hydrologic and land-surface models: A review, Water Resour. Res., 47, W07539,, 2011. 

Derksen, C., Sturm, M., Liston, G. E., Holmgren, J., Huntington, H., Silis, A., and Solie, D.: Northwest Territories and Nunavut snow characteristics from a subarctic traverse: Implications for passive microwave remote sensing, J. Hydrometeorol., 10, 448–463,, 2009. 

Derksen, C., Toose, P., Rees, A., Wang, L., English, M., Walker, A., and Sturm, M.: Development of a tundra-specific snow water equivalent retrieval algorithm for satellite passive microwave data, Remote Sens. Environ., 114, 1699–1709,, 2010. 

Derksen, C., Toose, P., Lemmetyinen, J., Pulliainen, J., Langlois, A., Rutter, N., and Fuller, M. C.: Evaluation of passive microwave brightness temperature simulations and snow water equivalent retrievals through a winter season, Remote Sens. Environ., 117, 236–248,, 2012. 

Derksen, C., Lemmetyinen, J., Toose, P., Silis, A., Pulliainen, J., and Sturm, M.: Physical properties of Arctic versus subarctic snow: Implications for high latitude passive microwave snow water equivalent retrievals, J. Geophys. Res.-Atmos., 119, 7254–7270,, 2014. 

Domine, F., Barrere, M., and Morin, S.: The growth of shrubs on high Arctic tundra at Bylot Island: impact on snow physical properties and permafrost thermal regime, Biogeosciences, 13, 6471–6486,, 2016. 

Domine, F., Picard, G., Morin, S., Barrere, M., Madore, J.-B., and Langlois, A.: Major Issues in Simulating some Arctic Snowpack Properties Using Current Detailed Snow Physics Models. Consequences for the Thermal Regime and Water Budget of Permafrost, J. Adv. Model. Earth Sy., 11, 34–44,, 2018. 

Durand, M. and Liu, D.: The need for prior information in characterizing snow water equivalent from microwave brightness temperatures, Remote Sens. Environ., 126, 248–257,, 2012. 

Fierz, C., Armstrong, R. L., Durand, Y., Etchevers, P., Greene, E., McClung, D. M., Nishimura, K., Satyawali, P. K., and Sokratov, S. A.: The international classification for seasonal snow on the ground, UNESCO, IHP–VII, Tech. Doc. Hydrol. No. 83, IACS Contrib. No. 1 80,, 2009. 

Gallet, J.-C., Domine, F., Zender, C. S., and Picard, G.: Measurement of the specific surface area of snow using infrared reflectance in an integrating sphere at 1310 and 1550 nm, The Cryosphere, 3, 167–182,, 2009. 

Garnaud, C., Bélair, S., Carrera, M. L., Derksen, C., Bilodeau, B., Abrahamowicz, M., Gauthier, N., and Vionnet, V.: Quantifying Snow Mass Mission Concept Trade-Offs Using an Observing System Simulation Experiment, J. Hydrometeorol., 20, 155–173,, 2019. 

Gisnås, K., Westermann, S., Schuler, T. V., Melvold, K., and Etzelmüller, B.: Small-scale variation of snow in a regional permafrost model, The Cryosphere, 10, 1201–1215,, 2016. 

Gouttevin, I., Langer, M., Löwe, H., Boike, J., Proksch, M., and Schneebeli, M.: Observation and modelling of snow at a polygonal tundra permafrost site: spatial variability and thermal implications, The Cryosphere, 12, 3693–3717,, 2018. 

Grünewald, T., Stötter, J., Pomeroy, J. W., Dadic, R., Moreno Baños, I., Marturià, J., Spross, M., Hopkinson, C., Burlando, P., and Lehning, M.: Statistical modelling of the snow depth distribution in open alpine terrain, Hydrol. Earth Syst. Sci., 17, 3005–3021,, 2013. 

Kelly, R.: The AMSR-E Snow Depth Algorithm: Description and Initial Results, J. Remote Sens. Soc. Japan 29, 307–317,, 2009. 

Kelly, R. E., Chang, A. T., Tsang, L., and Foster, J. L.: A prototype AMSR-E global snow area and snow depth algorithm, IEEE T. Geosci. Remote Sens., 41, 230–242,, 2003. 

Kim, R. S., Kumar, S., Vuyovich, C., Houser, P., Lundquist, J., Mudryk, L., Durand, M., Barros, A., Kim, E. J., Forman, B. A., Gutmann, E. D., Wrzesien, M. L., Garnaud, C., Sandells, M., Marshall, H.-P., Cristea, N., Pflug, J. M., Johnston, J., Cao, Y., Mocko, D., and Wang, S.: Snow Ensemble Uncertainty Project (SEUP): quantification of snow water equivalent uncertainty across North America via ensemble land surface modeling, The Cryosphere, 15, 771–791,, 2021. 

King, J., Kelly, R., Kasurak, A., Duguay, C., Gunn, G., Rutter, N., Watts, T., and Derksen, C.: Spatio-temporal influence of tundra snow properties on Ku-band (17.2 GHz) backscatter, J. Glaciol., 61, 267–279,, 2015. 

King, J., Derksen, C., Toose, P., Langlois, A., Larsen, C., Lemmetyinen, J., Marsh, P., Montpetit, B., Roy, A., Rutter, N., and Sturm, M.: The influence of snow microstructure on dual-frequency radar measurements in a tundra environment, Remote Sens. Environ., 215, 242–254,, 2018. 

Krol, Q. and Löwe, H.: Relating optical and microwave grain metrics of snow: the relevance of grain shape, The Cryosphere, 10, 2847–2863,, 2016. 

Larue, F., Royer, A., De Sève, D., Roy, A., and Cosme, E.: Assimilation of passive microwave AMSR-2 satellite observations in a snowpack evolution model over northeastern Canada, Hydrol. Earth Syst. Sci., 22, 5711–5734,, 2018. 

Liljedahl, A. K., Boike, J., Daanen, R. P., Fedorov, A. N., Frost, G. V., Grosse, G., Hinzman, L. D., Iijma, Y., Jorgenson, J. C., Matveyeva, N., Necsoiu, M., Raynolds, M. K., Romanovsky, V. E., Schulla, J., Tape, K. D., Walker, D. A., Wilson, C. J., Yabuki, H., and Zona, D.: Pan-Arctic ice-wedge degradation in warming permafrost and its influence on tundra hydrology, Nat. Geosci., 9, 312–318,, 2016. 

Liston, G. E.: Representing subgrid snow cover heterogeneities in regional and global models, J. Climate, 17, 1381–1397,<1381:RSSCHI>2.0.CO;2, 2004. 

Liston, G. E. and Sturm, M.: A snow-transport model for complex terrain, J. Glaciol., 44, 498–516,, 1998. 

Marsh, P., Bartlett, P., MacKay, M., Pohl, S., and Lantz, T.: Snowmelt energetics at a shrub tundra site in the western Canadian Arctic, Hydrol. Process., 24, 3603–3620,, 2010. 

Mavrovic, A., Pardo Lara, R., Berg, A., Demontoux, F., Royer, A., and Roy, A.: Soil dielectric characterization during freeze–thaw transitions using L-band coaxial and soil moisture probes, Hydrol. Earth Syst. Sci., 25, 1117–1131,, 2021. 

Meloche, J.: JulienMeloche/Gaussian_process_smrt_simulation, Release publication, Zenodo [code] [data set],, 2021. 

Meloche, J., Royer, A., Langlois, A., Rutter, N., and Sasseville, V.: Improvement of microwave emissivity parameterization of frozen Arctic soils using roughness measurements derived from photogrammetry, Int. J. Digit. Earth, 14, 1380–1396,, 2020. 

Mironov, V. L., De Roo, R. D., and Savin, I. V.: Temperature-dependable microwave dielectric model for an arctic soil, IEEE T. Geosci. Remote, 48, 2544–2556,, 2010. 

Montpetit, B., Royer, A., Langlois, A., Cliché, P., Roy, A., Champollion, N., Picard, G., Domine, F., and Obbard, R.: New shortwave infrared albedo measurements for snow specific surface area retrieval, J. Glaciol., 58, 941–952,, 2012. 

Mortimer, C., Mudryk, L., Derksen, C., Luojus, K., Brown, R., Kelly, R., and Tedesco, M.: Evaluation of long-term Northern Hemisphere snow water equivalent products, The Cryosphere, 14, 1579–1594,, 2020. 

Pan, J., Durand, M. T., Vander Jagt, B. J., and Liu, D.: Application of a Markov Chain Monte Carlo algorithm for snow water equivalent retrieval from passive microwave measurements, Remote Sens. Environ., 192, 150–165,, 2017. 

Parr, C., Sturm, M., and Larsen, C.: Snowdrift Landscape Patterns: An Arctic Investigation, Water Resour. Res., 56, e2020WR027823,, 2020. 

Picard, G., Sandells, M., and Löwe, H.: SMRT: an active–passive microwave radiative transfer model for snow with multiple microstructure and scattering formulations (v1.0), Geosci. Model Dev., 11, 2763–2788,, 2018. 

Pomeroy, J. W., Gray, D. M., Hedstrom, N. R., and Janowicz, J. R.: Prediction of seasonal snow accumulation in cold climate forests, Hydrol. Process., 16, 3543–3558,, 2002. 

Porter, C., Morin, P., Howat, I., Noh, M.-J., Bates, B., Peterman, K., Keesey, S., Schlenk, M., Gardiner, J., Tomko, K., Willis, M., Kelleher, C., Cloutier, M., Husby, E., Foga, S., Nakamura, H., Platson, M., Wethington, M., Williamson, C., Bauer, G., Enos, J., Arnold, G., Kramer, W., Becker, P., Doshi, A., D'Souza, C., Cummens, P., Laurier, F., and Bojesen, M.: ArcticDEM, Harvard Dataverse [data set], V1,, 2018. 

Proksch, M., Löwe, H., and Schneebeli, M.: Density, specific surface area, and correlation length of snow measured by high-resolution penetrometry, J. Geophys. Res.-Earth, 120, 346–362,, 2015. 

Pulliainen, J.: Mapping of snow water equivalent and snow depth in boreal and sub-arctic zones by assimilating space-borne microwave radiometer data and ground-based observations, Remote Sens. Environ., 101, 257–269,, 2006. 

Pulliainen, J., Luojus, K., Derksen, C., Mudryk, L., Lemmetyinen, J., Salminen, M., Ikonen, J., Takala, M., Cohen, J., Smolander, T., and Norberg, J.: Patterns and trends of Northern Hemisphere snow mass from 1980 to 2018, Nature, 581, 294–298,, 2020. 

Quiñonero-Candela, J. and Rasmussen, C. E.: A unifying view of sparse approximate Gaussian process regression, J. Mach. Learn. Res., 6, 1939–1959, 2005. 

Rees, A., English, M., Derksen, C., Toose, P., and Silis, A.: Observations of late winter Canadian tundra snow cover properties, Hydrol. Process., 28, 3962–3977,, 2014. 

Roy, A., Picard, G., Royer, A., Montpetit, B., Dupont, F., Langlois, A., Derksen, C., and Champollion, N.: Brightness Temperature Simulations of the Canadian Seasonal Snowpack Driven by Measurements of the Snow Specific Surface Area, IEEE T. Geosci. Remote, 51, 4692–4704,, 2013. 

Royer, A., Roy, A., Montpetit, B., Saint-Jean-Rondeau, O., Picard, G., Brucker, L., and Langlois, A.: Comparison of commonly-used microwave radiative transfer models for snow remote sensing, Remote Sens. Environ., 190, 247–259,, 2017. 

Royer, A., Domine, F., Roy, A., Langlois, A., Marchand, N., and Davesne, G.: New northern snowpack classification linked to vegetation cover on a latitudinal mega-transect across northeastern Canada, Écoscience, 28, 225–242,, 2021. 

Rutter, N., Sandells, M., Derksen, C., Toose, P., Royer, A., Montpetit, B., Langlois, A., Lemmetyinen, J., and Pulliainen, J.: Snow stratigraphic heterogeneity within ground-based passive microwave radiometer footprints: Implications for emission modeling, J. Geophys. Res.-Earth, 119, 550–565,, 2014. 

Rutter, N., Sandells, M. J., Derksen, C., King, J., Toose, P., Wake, L., Watts, T., Essery, R., Roy, A., Royer, A., Marsh, P., Larsen, C., and Sturm, M.: Effect of snow microstructure variability on Ku-band radar snow water equivalent retrievals, The Cryosphere, 13, 3045–3059,, 2019. 

Saberi, N., Kelly, R., Toose, P., Roy, A., and Derksen, C.: Modeling the observed microwave emission from shallow multi-layer Tundra Snow using DMRT-ML, Remote Sens., 9, 1327,, 2017. 

Saberi, N., Kelly, R., Pan, J., Durand, M., Goh, J., and Scott, K. A.: The Use of a Monte Carlo Markov Chain Method for Snow-Depth Retrievals: A Case Study Based on Airborne Microwave Observations and Emission Modeling Experiments of Tundra Snow, IEEE T. Geosci. Remote, 59, 1876–1889,, 2020. 

Salvatier, J., Wiecki, T. V., and Fonnesbeck, C.: Probabilistic programming in Python using PyMC3, PeerJ Comput. Sci., 2, e55,, 2016. 

Sturm, M. and Holmgren, J.: An Automatic Snow Depth Probe for Field Validation Campaigns, Water Resour. Res., 54, 9695–9701,, 2018. 

Sturm, M. and Wagner, A. M.: Using repeated patterns in snow distribution modeling: An Arctic example, Water Resour. Res., 46, 1–15,, 2010. 

Sturm, M., McFadden, J. P., Liston, G. E., Stuart Chapin, F., Racine, C. H., and Holmgren, J.: Snow-shrub interactions in Arctic Tundra: A hypothesis with climatic implications, J. Climate, 14, 336–344,<0336:SSIIAT>2.0.CO;2, 2001. 

Takala, M., Luojus, K., Pulliainen, J., Derksen, C., Lemmetyinen, J., Kärnä, J. P., Koskinen, J., and Bojkov, B.: Estimating northern hemisphere snow water equivalent for climate research through assimilation of space-borne radiometer data and ground-based measurements, Remote Sens. Environ., 115, 3517–3529,, 2011. 

Tsang, L., Chen, C. Te, Chang, A. T. C., Guo, J., and Ding, K. H.: Dense media radiative transfer theory based on quasicrystalline approximation with applications to passive microwave remote sensing of snow, Radio Sci., 35, 731–749,, 2000. 

Vander Jagt, B. J., Durand, M. T., Margulis, S. A., Kim, E. J., and Molotch, N. P.: The effect of spatial variability on the sensitivity of passive microwave measurements to snow water equivalent, Remote Sens. Environ., 136, 163–179,, 2013. 

Vargel, C., Royer, A., St-jean-rondeau, O., Picard, G., Roy, A., Sasseville, V., and Langlois, A.: Remote Sensing of Environment Arctic and subarctic snow microstructure analysis for microwave brightness temperature simulations, Remote Sens. Environ., 242, 111754,, 2020. 

Walker, B., Wilcox, E. J., and Marsh, P.: Accuracy assessment of late winter snow depth mapping for tundra environments using Structure-from-Motion photogrammetry, Antarct. Sci., 17, 1–17,, 2020a. 

Walker, B., Wilcox, E., and Marsh, P.: Structure-from-Motion Snow Depth Products and In Situ Observations for Late-Winter Tundra Mapping Project, Trail Valley Creek Research Station, Spring 2018, Scholars Portal Dataverse, V1 [data set],, 2020b. 

Wegmüller, U. and Mätzler, C.: Rough bare soil reflectivity model, IEEE T. Geosci. Remote, 37, 1391–1395,, 1999. 

Wiesmann, A. and Mätzler, C.: Microwave emission model of layered snowpacks, Remote Sens. Environ., 70, 307–316,, 1999.  

Winstral, A. and Marks, D.: Long-term snow distribution observations in a mountain catchment: Assessing variability, time stability, and the representativeness of an index site, Water Resour. Res., 50, 293–305,, 2014. 

Winstral, A., Elder, K., and Davis, R. E.: Spatial Snow Modeling of Wind-Redistributed Snow Using Terrain-Based Parameters, J. Hydrometeorol., 3, 524–538,<0524:SSMOWR>2.0.CO;2, 2002. 

Winstral, A., Marks, D., Gurney, R.: Simulating wind-affected snow accumulations at catchment to basin scales. Adv. Water Resour. 55, 64–79., 2013. 

Short summary
To estimate snow water equivalent from space, model predictions of the satellite measurement (brightness temperature in our case) have to be used. These models allow us to estimate snow properties from the brightness temperature by inverting the model. To improve SWE estimate, we proposed incorporating the variability of snow in these model as it has not been taken into account yet. A new parameter (coefficient of variation) is proposed because it improved simulation of brightness temperature.