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

Abstract

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 (${\mathit{\mu }}_{\mathrm{sd}}$) values and coefficients of variation (CV${}_{\mathrm{sd}}$). Snow depth variability (CV${}_{\mathrm{sd}}$) 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 (CV${}_{\mathrm{sd}}$) is proposed as an effective parameter to account for sub-pixel variability in PMW emission, improving simulation by 8 K. SMRT simulations using a CV${}_{\mathrm{sd}}$ of 0.9 best matched CV${}_{\mathrm{sd}}$ observations from spatial datasets for areas $\mathit{>}$ 3 km${}^{\mathrm{2}}$, 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 ($T_{\mathrm{B}}$) 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 $T_{\mathrm{B}}$, so differences between $T_{\mathrm{B}}$ 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 (CV${}_{\mathrm{sd}}$), the ratio between standard deviation (${\mathit{\sigma }}_{\mathrm{sd}}$) and the mean of snow depth (${\mathit{\mu }}_{\mathrm{sd}}$), indicating the extent and spread of a distribution (i.e., high variability over thin snow will lead to high values of CV${}_{\mathrm{sd}}$). Also, Liston (2004) proposed nine categories of CV${}_{\mathrm{sd}}$ 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 ${\mathit{\mu }}_{\mathrm{sd}}$ and CV${}_{\mathrm{sd}}$ on the log-normal distribution scale parameters ($\mathit{\lambda }$, $\mathit{\zeta }$). Gisnås et al. (2016) adapted that approach using scale parameters ($\mathit{\alpha }$, $\mathit{\beta }$) of the gamma distribution. In all cases, CV${}_{\mathrm{sd}}$ 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, CV${}_{\mathrm{sd}}$ 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 CV${}_{\mathrm{sd}}$ 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 $T_{\mathrm{B}}$) 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 CV${}_{\mathrm{sd}}$ 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, ${\mathit{\mu }}_{\mathrm{sd}}$ and CV${}_{\mathrm{sd}}$. 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 $T_{\mathrm{B}}$ from a satellite sensor at 19 and 37 GHz and $T_{\mathrm{B}}$ 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 (68${}^{\circ }$44${}^{\prime }$ N, 133${}^{\circ }$33${}^{\prime }$ 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 (69${}^{\circ }$13${}^{\prime }$ N, 104${}^{\circ }$53${}^{\prime }$ 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 1

Locations 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.

[Figure omitted. See PDF]

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 cm${}^{\mathrm{3}}$ 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 km${}^{\mathrm{2}}$. An airborne lidar dataset of TVC snow depths (93 km${}^{\mathrm{2}}$ at 10 m resolution) collected by an aircraft in April 2013 (Rutter et al., 2019) was also used. Monte Carlo simulations of both the ${\mathit{\mu }}_{\mathrm{sd}}$ and ${\mathrm{CV}}_{\mathrm{sd}}$ 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 (${\mathit{\mu }}_{\mathrm{sd}}$ and ${\mathrm{CV}}_{\mathrm{sd}}$) could be calculated. Also, a small RPAS map was available for CB with spatial coverage of 0.2 km${}^{\mathrm{2}}$ 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 1

Summary 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 ($\times$). See Table 2 for full dates.

Site	Date	Depth measurement	Snow pit	SSA	Density
TVC	15–25 March 2019	8541	32	$\times$	$\times$
	15–23 March 2018	7190	36	$\times$	$\times$
TVC18-RPAS	12 March–22 April 2018	Pixels: 6 325 365; resolution: 1 m
TVC13-Lidar	April 2013	Pixels: 969 168; resolution: 10 m
CB18-RPAS	15 April 2018	Pixels: 72 902; resolution: 1 m
CB	15–29 April 2019	982	64	$\times$	$\times$
	12–24 April 2018	–	50	$\times$	$\times$
	1–8 May 2017	4045	51		$\times$
	2–10 April 2016	3403	35		$\times$
	9–16 April 2015	12 282	48		$\times$

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

Microwave $T_{\mathrm{B}}$ 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. $T_{\mathrm{B}}$ 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 $T_{\mathrm{B}}$ 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. $T_{\mathrm{B}}$ values were temporally averaged to match times of field measurements, representing peak winter snow accumulation (Table 2). Also, $T_{\mathrm{B}}$ 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 (${\mathit{\kappa }}_{\mathrm{37}}=\mathrm{1.39},{\mathit{\kappa }}_{\mathrm{19}}=\mathrm{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 (${\mathit{\epsilon }}^{\prime }\approx$ 4–6) (Mironov et al., 2010), which had a similar permittivity to frozen soil (${\mathit{\epsilon }}^{\prime }\approx$ 2–4) (Mavrovic et al., 2021) than liquid water (${\mathit{\epsilon }}^{\prime }>$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 2

Summary of mean basal and air surface temperatures for SMRT simulations, precipitable water (PWAT) used for atmospheric correction, and measured (corrected) $T_{\mathrm{B}}$ at both vertical (V) and horizontal (H) polarization by the SSMIS sensor (platform F18).

Sites	$T_{\mathrm{base}}$	$T_{\mathrm{surface}}$	PWAT	$T_{\mathrm{B}}\mathrm{37}$	$T_{\mathrm{B}}\mathrm{37}$
	(K)	NARR (K)	(mm)	H (K)	V (K)
CB (15–29 April 2019)	257	261.5	3.61	195.3	211.0
CB (12–24 April 2018)	257	260.1	3.72	179.3	195.7
CB (1–8 May 2017)	263	261.3	3.33	187.1	205.0
CB (2–10 April 2016)	256	258.8	2.80	190.1	215.4
CB (9–16 April 2015)	254	256.2	2.34	193.0	215.9
TVC (15–25 March 2019)	266	261.8	7.04	177.0	199.5
TVC (15–23 March 2018)	264	261.8	4.21	176.6	197.6

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\left(x\right)$) and covariance ($k\left(x{x}^{\prime }\right)$) functions from the multi-variate Gaussian distribution are used to fit data ($x$: snow depth, $y$: ratio of layer, DHF). The $m\left(x\right)$ function describes the expected value of the distribution, and the $k(x{x}^{\prime }$) describes the shape of the correlation between data points ($x_i$). 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). $$\begin{array}{}\text{1}& {\displaystyle }\begin{array}{rl}\mathrm{Posterior}& ={\displaystyle \frac{\mathrm{likelihood}\cdot \mathrm{prior}}{\mathrm{marginal}\mathrm{likelihood}}}\\ & =p\left(f\mathrm{|}y,X\right)={\displaystyle \frac{p\left(y\mathrm{|}X,f\right)\cdot p\left(f\right)}{p\left(y\mathrm{|}X\right)}},\end{array}\text{2}& {\displaystyle }f\left(x\right)\sim \mathrm{GP}\left(m\left(x\right),k\left(x,x^{\prime }\right),\mathit{\varphi }(x\right)).\end{array}$$ Equation (2) defined $f$ as a function of $m\left(x\right),k\left(x,x^{\prime }\right)$. A mean function $m\left(x\right)$, following an inverse logic function ($\mathit{\varphi }$) (Eq. 3), was chosen due to the close fit with observations. The covariance function $k\left(x,x^{\prime }\right)$ is a Gaussian white noise covariance function and is defined with noise ($\mathit{\sigma }$) and the Kronecker delta function (${\mathit{\delta }}_{x,x}$) (Eq. 4), to best fit the observations. By using a scaling function ($\mathit{\varphi }$), 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 ($\mathit{\varphi }$) 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). $$\begin{array}{}\text{3}& {\displaystyle }m\left(x\right)=\mathit{\varphi }\left(x\right)=c+b\left[{\displaystyle \frac{e^{a\left(x-x_{\mathrm{0}}\right)}}{\mathrm{1}+e^{a\left(x-x_{\mathrm{0}}\right)}}}\right]\text{4}& {\displaystyle }k\left(x{x}^{\prime }\right)={\mathit{\sigma }}^{\mathrm{2}}{\mathit{\delta }}_{x,x^{\prime }}\mathit{\varphi }\left(x\right)\end{array}$$

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 ${\mathit{\mu }}_{\mathrm{sd}}$ and the CV${}_{\mathrm{sd}}$ 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 $\mathit{>}$ 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 3

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

Site	PDF	KS stats	$p$ value
TVC	log-normal	0.029	0.41
	gamma	0.039	0.11
CB	log-normal	0.024	0.63
	gamma	0.017	0.95

Figure 2

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

[Figure omitted. See PDF]

Distributions with parameterization using measured ${\mathit{\mu }}_{\mathrm{sd}}$ and CV${}_{\mathrm{sd}}$ (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 CV${}_{\mathrm{sd}}$ 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 km${}^{\mathrm{2}}$) showed a much larger CV${}_{\mathrm{sd}}=$ 0.78 than magnaprobe data ($n=$ 15 731) with CV${}_{\mathrm{sd}}=$ 0.43 (Table 4). A RPAS dataset was also available for CB but with a much smaller spatial coverage (0.2 km${}^{\mathrm{2}}$) showing a CV${}_{\mathrm{sd}}$ of 0.49. In Figs. 3 and 4, we investigated the relationship between spatial coverage of sampling and the CV${}_{\mathrm{sd}}$ parameter. Datasets include RPAS-derived data at TVC (TVC18-RPAS) containing six areas with various size from 1–3 km${}^{\mathrm{2}}$, a CB18-RPAS map of 0.2 km${}^{\mathrm{2}}$ 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 km${}^{\mathrm{2}}$ showed areas with high ${\mathrm{CV}}_{\mathrm{sd}}$ (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 ${\mathit{\mu }}_{\mathrm{sd}}$ and ${\mathrm{CV}}_{\mathrm{sd}}$ 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 ${\mathrm{CV}}_{\mathrm{sd}}$. A difference of 0.2 between the full ${\mathrm{CV}}_{\mathrm{sd}}$ of the RPAS (5 km${}^{\mathrm{2}}$) and lidar (93 km${}^{\mathrm{2}}$) maps (Fig. 4) was found. Also included in Fig. 4 are in situ ${\mathrm{CV}}_{\mathrm{sd}}$ 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 km${}^{\mathrm{2}}$ (${\mathrm{CV}}_{\mathrm{sd}}=$ 1, Daring Lake site; Chris Derksen, personal communication, 2009) and 198 km${}^{\mathrm{2}}$ (${\mathrm{CV}}_{\mathrm{sd}}=$ 0.89, Eureka site; Saberi et al., 2007).

Table 4

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

Site	$n$	$\mathit{\mu }$ (m)	$\mathit{\sigma }$ (m)	CV${}_{\mathrm{sd}}$
TVC19	8541	0.44	0.14	0.33
TVC18	7190	0.39	0.21	0.54
TVC	15 731	0.42	0.19	0.43
TVC18-RPAS	6 325 365	0.46	0.36	0.78
TVC13-Lidar	969 168	0.40	0.23	0.58
CB19	982	0.42	0.17	0.40
CB18	577	0.34	0.18	0.53
CB18-RPAS	7290	0.39	0.19	0.49
CB17	4045	0.42	0.19	0.46
CB16	3403	0.28	0.16	0.61
CB15	12 282	0.32	0.18	0.57
CB	20 712	0.36	0.18	0.52

Figure 3

RPAS and lidar dataset of snow depth at TVC (TVC13-Lidar and TVC18-RPAS). TVC13-Lidar is the largest dataset covering 93 km${}^{\mathrm{2}}$. 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 ${\mathrm{CV}}_{\mathrm{sd}}$ for each cell; (c) same as panel (b) but for ${\mathit{\mu }}_{\mathrm{sd}}$.

[Figure omitted. See PDF]

Figure 4

Snow depth mean (${\mathit{\mu }}_{\mathrm{sd}}$) (a, b) and variability (${\mathrm{CV}}_{\mathrm{sd}}$) (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 ${\mathit{\mu }}_{\mathrm{sd}}$ and ${\mathrm{CV}}_{\mathrm{sd}}$ of both full areas are shown by the black dotted and dashed line.

[Figure omitted. See PDF]

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 ($\mathit{<}$ 2 mm) were classified as WS. Indurated hoar, a wind slab metamorphosed into depth hoar, was classified into DH with a typical density $\sim \mathrm{300}$ kg m${}^{-\mathrm{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${}^{-\mathrm{3}}$) and DH (266 kg m${}^{-\mathrm{3}}$) were closer together. SSA distributions also showed a gap between both mean values (WS: 19.7 m${}^{\mathrm{2}}$ kg${}^{-\mathrm{1}}$; DH: 11.1 m${}^{\mathrm{2}}$ kg${}^{-\mathrm{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 5

SSA 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).

[Figure omitted. See PDF]

Table 5

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

Snow property	Best fit PDF		$\mathit{\mu }$	$\mathit{\sigma }$
SSA (m${}^{\mathrm{2}}$ kg${}^{-\mathrm{1}}$)	log-normal	DH	11.1	3.8
		WS	19.7	7.8
			$\mathit{\mu }$	$\mathit{\sigma }$
Density (kg m${}^{-\mathrm{3}}$)	normal	DH	266.3	48.9
		WS	335.2	57.1

Parr et al. (2020) found a key threshold of ${\mathit{\mu }}_{\mathrm{sd}}+\mathrm{1}{\mathit{\sigma }}_{\mathrm{sd}}$ to define snow drifts in tundra environments. This threshold ($\mathit{>}$ 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.

[Figure omitted. See PDF]

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\left(x\right)$, following an inverse logic function (${\mathit{\varphi }}_{\mathrm{1}}:$ Eq. 3), was chosen with parameters $a=-\mathrm{5}$, $x_{\mathrm{0}}=\mathrm{0.6}$, $b=\mathrm{0.35}$ and $c=\mathrm{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 $\mathit{>}$ 40 cm by using a scaling function (${\mathit{\varphi }}_{\mathrm{2}}:a=-\mathrm{5}$, $x_{\mathrm{0}}=\mathrm{0.6}$, $b=\mathrm{1.5}$ and $c=$ 0.25). An inverse logic function (${\mathit{\varphi }}_{\mathrm{1}},{\mathit{\varphi }}_{\mathrm{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 (${\mathit{\mu }}_{\mathrm{sd}}$, CV${}_{\mathrm{sd}}$) 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 $\mathit{<}$ 0.5 m, which was then reduced for depths $\mathit{>}$ 0.5 m following the red posterior prediction representation in Fig. 7.

Figure 7

Prediction 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.

[Figure omitted. See PDF]

3.4 SMRT simulation of sub-grid variability within sensor footprint

SMRT simulations using measured snowpack properties were compared with the satellite measurements of $T_{\mathrm{B}}$. 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 ${\mathit{\mu }}_{\mathrm{sd}}$ and CV${}_{\mathrm{sd}}$. The ratio of each layer is predicted using the GP fit with depth as input from the log-normal distribution. Mean SSA (DH: 11 m${}^{\mathrm{2}}$ kg${}^{-\mathrm{1}}$; WS: 20 m${}^{\mathrm{2}}$ kg${}^{-\mathrm{1}}$) and density (DH: 266 kg m${}^{-\mathrm{3}}$; WS: 335 kg m${}^{-\mathrm{3}}$) per layer were determined from measurements (Fig. 5).

For one standard EASE-Grid 2.0 pixel, a distribution of sub-grid $T_{\mathrm{B}}$ was simulated to reproduce a realistic distribution of $T_{\mathrm{B}}$ 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 (${\mathit{\mu }}_{\mathrm{sd}}$) of each region (Table 4) and a depth variability (CV${}_{\mathrm{sd}}$) 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 CV${}_{\mathrm{sd}}=$ 0.7 , the simulated distribution showed a wide sub-pixel variability ($\pm \mathrm{40}$ K) with a mean value of $T_{\mathrm{B}\mathrm{37}\mathrm{V}}=$ 194.7 K (blue line in Fig. 8a), which is very close to the satellite-measured $T_{\mathrm{B}\mathrm{37}\mathrm{V}}$ of 196.5 K (green dotted line in Fig. 8a). In this case, the $T_{\mathrm{B}}$ 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 $T_{\mathrm{B}}$ was chosen, and it was assumed that the sub-pixel effect combined linearly at this scale. Because the simulated $T_{\mathrm{B}\mathrm{37}\mathrm{V}}$ distribution was not exactly a normal distribution, it appeared that the mean $T_{\mathrm{B}}$ of this distribution increased when CV${}_{\mathrm{sd}}$ increased (Fig. 8b). This meant that snow depth variability (CV${}_{\mathrm{sd}}$) must be accounted for when estimating the average $T_{\mathrm{B}}$ at 37 GHz, in addition to the mean snow depth values. The influence of the GP simulation on the mean simulated $T_{\mathrm{B}\mathrm{37}\mathrm{V}}$ was approximately 10 K (Fig. 8b) as CV${}_{\mathrm{sd}}$ varies from 0.1 to 1. The addition of snow variability in simulation (Fig. 8c and d) of 19 GHz has negligible effect on $T_{\mathrm{B}\mathrm{19}}$ and shows a constant simulation across the ${\mathrm{CV}}_{\mathrm{sd}}$ range of 0.1 to 1. Simulation of $T_{\mathrm{B}\mathrm{19}}$ shows higher biases at horizontal polarization than at vertical polarization.

Figure 8

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

[Figure omitted. See PDF]

GP simulation reduced biases by 5 K with a higher optimized CV${}_{\mathrm{sd}}$ (intersection of red/blue–green line, Fig. 8b). A similar pattern was observed for CB (not shown here), but the measured $T_{\mathrm{B}}$ at CB was much higher than the GP simulation, resulting in a large bias for CB ($\sim$ 20 K) compared to TVC (Table 6). Both sites suggested a larger CV${}_{\mathrm{sd}}$, which agreed with a CV${}_{\mathrm{sd}}$ 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 CV${}_{\mathrm{sd}}$ (Fig. 9). Total RMSE is minimized with higher CV${}_{\mathrm{sd}}$ (0.8–0.9) typical of a large sampling scale (over 4 km${}^{\mathrm{2}}$) as shown in Fig. 4.

Table 6

Bias between SMRT simulated and measured $T_{\mathrm{B}}$ from the SSMIS sensor at each site.

		Bias (K)
		CB		TVC		RMSE (K)
SMRT simulation type	Year	H pol	V pol	H pol	V pol	H pol	V pol
Mean depth and DHF	2019	28.2	25.9	6.9	10.3	17.8	19.1
	2018	8.0	5.3	5.1	6.8
	2017	19.9	18.9	–	–
	2016	16.9	23.2	–	–
	2015	24.7	29.1	–	–
GP simulation CV $=$ 0.9	2019	18.6	15.7	$-$4.4	$-$1.2	9.7	10.4
	2018	$-$3.7	$-$6.2	$-$4.9	$-$3.2
	2017	10.4	9.3	–	–
	2016	7.1	13.5	–	–
	2015	10.0	13.9	–	–

Figure 9

Overall 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.

[Figure omitted. See PDF]

4 Discussion

We show that as spatial coverage increased, the ${\mathrm{CV}}_{\mathrm{sd}}$ parameter converged to the full area values (Fig. 4). Monte Carlo simulations of snow depth distribution and coverage show high variation in ${\mathrm{CV}}_{\mathrm{sd}}$ (from 0.1 to 2) for areas $\mathit{<}$ 10 km${}^{\mathrm{2}}$. 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 ${\mathrm{CV}}_{\mathrm{sd}}$ 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 ${\mathrm{CV}}_{\mathrm{sd}}$ (dark green areas in Fig. 3b), which are commonly found in TVC (Walker et al., 2020a). ${\mathrm{CV}}_{\mathrm{sd}}$ was lower for areas without drifts (light green areas in Fig. 3b). In areas >10 km${}^{\mathrm{2}}$ (Fig. 4d), variation in ${\mathrm{CV}}_{\mathrm{sd}}$ was reduced and yielded higher values $\mathit{>}$ 0.6.

Convergence to higher ${\mathrm{CV}}_{\mathrm{sd}}$ 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 ${\mathrm{CV}}_{\mathrm{sd}}$ of TVC13-Lidar converged to 0.6 at 93 km${}^{\mathrm{2}}$, but two in situ points from other studies at 625 km${}^{\mathrm{2}}$ had higher ${\mathrm{CV}}_{\mathrm{sd}}$ (0.9–1) due to larger coverage or different site characteristics. This indicates that a ${\mathrm{CV}}_{\mathrm{sd}}$ 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 km${}^{\mathrm{2}}$ for GlobSnow 3.0 (Pulliainen et al., 2020). For active sensors (resolution $\mathit{<}$ 1 km), the high variability in ${\mathrm{CV}}_{\mathrm{sd}}$ under 1 km${}^{\mathrm{2}}$ 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 ${\mathit{\mu }}_{\mathrm{sd}}$ and ${\mathrm{CV}}_{\mathrm{sd}}$ 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 CV${}_{\mathrm{sd}}$ parameter can also be used as a priori information with a distribution from 0.8 to 1, since it improved $T_{\mathrm{B}}$ RMSE by $\sim \mathrm{8}$ K (Fig. 9). This approach improved $T_{\mathrm{B}}$ simulation compared to using only mean values of snowpack properties by adding variability within the footprint. The CV${}_{\mathrm{sd}}$ 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 (CV${}_{\mathrm{sd}}$). The relationship found in Fig. 6 used to predict DHF (Fig. 7) could also be used deterministically with the mean function (${\mathit{\varphi }}_{\mathrm{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 ${\mathrm{CV}}_{\mathrm{sd}}$ 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 $T_{\mathrm{B}}$ at 37 GHz resulting from an increase in ${\mathrm{CV}}_{\mathrm{sd}}$ (Fig. 8).

This idea of modulating the amount of scatterers based of DHF prediction and a distribution of snow depth (${\mathit{\mu }}_{\mathrm{sd}}$ and CV${}_{\mathrm{sd}}$) 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 CV${}_{\mathrm{sd}}$ parameter is proposed as an effective parameter to account for variability inside the grid cell, while the mean depth (${\mathit{\mu }}_{\mathrm{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 CV${}_{\mathrm{sd}}$ 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 ${\mathit{\mu }}_{\mathrm{sd}}$ and CV${}_{\mathrm{sd}}$ 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 (${\mathit{\mu }}_{\mathrm{sd}}=$ 0.33–0.44 m) and coefficient of variations (CV${}_{\mathrm{sd}}=$ 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 ${\mathit{\mu }}_{\mathrm{sd}}$ and ${\mathrm{CV}}_{\mathrm{sd}}$ as a function of spatial coverage. An increase in ${\mathrm{CV}}_{\mathrm{sd}}$ matched increased spatial coverage of snow depth sampling, indicating that a higher ${\mathrm{CV}}_{\mathrm{sd}}$ (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 ${\mathrm{CV}}_{\mathrm{sd}}$ ($\mathit{>}$ 0.9) for areas $\mathit{<}$ 10 km${}^{\mathrm{2}}$, suggesting a need for topography-based prediction of ${\mathit{\mu }}_{\mathrm{sd}}$ and ${\mathrm{CV}}_{\mathrm{sd}}$ at this scale. The CV${}_{\mathrm{sd}}$ is shown to be an effective parameter to account for snow depth variability in simulation of snow $T_{\mathrm{B}}$. 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 CV${}_{\mathrm{sd}}$. SMRT simulation was used successfully to simulate satellite $T_{\mathrm{B}}$, 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 $T_{\mathrm{B}}$ were reduced by 8 K after optimizing CV${}_{\mathrm{sd}}$ to higher values (0.8–1.0), thereby matching CV${}_{\mathrm{sd}}$ of spatially distributed snow depth from TVC18-RPAS accounting for variation in snow properties inside the footprint of satellite sensor. The CV${}_{\mathrm{sd}}$ 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 https://github.com/JulienMeloche/Gaussian_process_smrt_simulation (last access: 27 December 2021; 10.5281/zenodo.5806672, Meloche, 2021). The RPAS map and magnaprobe from TVC are available at 10.5683/SP2/PWSKKG (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.

Disclaimer

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

Acknowledgements

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.