**Research article**| 27 Apr 2022

# Snow water equivalent change mapping from slope-correlated synthetic aperture radar interferometry (InSAR) phase variations

Jayson Eppler Bernhard Rabus and Peter Morse

^{1},

^{1},

^{2}

**Jayson Eppler et al.**Jayson Eppler Bernhard Rabus and Peter Morse

^{1},

^{1},

^{2}

^{1}School of Engineering Science, Simon Fraser University, Burnaby, BC V5A 1S6, Canada^{2}Geological Survey of Canada, Natural Resources Canada, Ottawa, ON, K1A 0E8, Canada

^{1}School of Engineering Science, Simon Fraser University, Burnaby, BC V5A 1S6, Canada^{2}Geological Survey of Canada, Natural Resources Canada, Ottawa, ON, K1A 0E8, Canada

**Correspondence**: Jayson Eppler (jaysone@sfu.ca)

**Correspondence**: Jayson Eppler (jaysone@sfu.ca)

Received: 20 Nov 2021 – Discussion started: 15 Dec 2021 – Revised: 23 Mar 2022 – Accepted: 23 Mar 2022 – Published: 27 Apr 2022

Area-based measurements of snow water equivalent (SWE) are important for understanding earth system processes such as glacier mass balance, winter hydrological storage in drainage basins, and ground thermal regimes. Remote sensing techniques are ideally suited for wide-scale area-based mapping with the most commonly used technique to measure SWE being passive microwave, which is limited to coarse spatial resolutions of 25 km or greater and to areas without significant topographic variation. Passive microwave also has a negative bias for large SWE. Another method is repeat-pass synthetic aperture radar interferometry (InSAR) that allows measurement of SWE change at much higher spatial resolution. However, it has not been widely adopted because (1) the phase unwrapping problem has not been robustly addressed, especially for interferograms with poor coherence, and (2) SWE change maps scaled directly from repeat-pass interferograms are not an absolute measurement but contain unknown offsets for each contiguous coherent area. We develop and test a novel method for repeat-pass InSAR-based dry-snow SWE estimation that exploits the sensitivity of the dry-snow refraction-induced InSAR phase to topographic variations. The method robustly estimates absolute SWE change at spatial resolutions of < 1 km without the need for phase unwrapping. We derive a quantitative signal model for this new SWE change estimator and identify the relevant sources of bias. The method is demonstrated using both simulated SWE distributions and a 9-year RADARSAT-2 (C-band, 5.405 GHz) spotlight-mode dataset near Inuvik, Northwest Territories (NWT), Canada. SWE results are compared to in situ snow survey measurements and estimates from ERA5 reanalysis. Our method performs well in high-relief areas, thus providing complementary coverage to passive-microwave-based SWE estimation. Further, our method has the advantage of requiring only a single wavelength band and thus can utilize existing spaceborne synthetic aperture radar systems.

- Article
(13490 KB) -
Supplement
(7883 KB) - BibTeX
- EndNote

The works published in this journal are distributed under the Creative Commons Attribution 4.0 License. This licence does not affect the Crown copyright work, which is re-usable under the Open Government Licence (OGL). The Creative Commons Attribution 4.0 License and the OGL are interoperable and do not conflict with, reduce or limit each other.

© Crown copyright 2022

Snowmelt is one of the greatest sources of water in snow-affected areas (Barnett et al., 2005), snow accumulation on glaciers is critical to their mass balance and longevity (Zemp et al., 2009), and snow cover variation is a dominant control on ground thermal regimes in cold regions (Mackay and MacKay, 1974); thus quantification of snowpack conditions is essential to understanding the role of snow on earth system processes. Snow water equivalent (SWE) is the amount of liquid water in a vertical column of snow and is a key parameter of snowpack conditions. The SWE of a snow layer is equal to the product of its depth and mean density. Traditionally, SWE variation over space and time is interpolated from sparse point-based measurements (Grünewald et al., 2010), which generally show significant spatial variation due to factors such as topography and vegetation (Anderton et al., 2004; Jost et al., 2007), as well as temporal bias if the acquisition of each set of points is distributed over a range of dates (Smyth et al., 2020). Ideally the SWE measurement is carried out over a short time span to limit within-set variation due to changing meteorological conditions. Further bias may be introduced if different point-based methods are integrated together to estimate SWE over large areas. To avoid biases arising from the interpolation of point-based measurements of typical snow packs, several remote sensing techniques have been developed for wide-scale, area-based determination of SWE. With precipitation regimes changing globally under persistent climate change (Pachauri et al., 2015), remote sensing techniques offer the potential to monitor hydrological conditions over vast areas of cold regions where point-based measurements are typically sparse and collected intermittently with a range of techniques.

Snow depth has been determined using surface elevation models developed using lidar, as well as structure-from-motion (optical photogrammetry) systems, either by differential repeat-pass measurement or comparison with a pre-existing snow-free reference elevation model (Deems et al., 2013; Nolan et al., 2015). However, as these systems are presently only feasible on airborne platforms they have not yet been used for repeated, wide-scale SWE mapping, which is required to substantively improve hydrological monitoring. Alternatively, spaceborne systems are designed for repeated wide-scale monitoring, and both passive- and active-microwave methods to determine SWE have been demonstrated (Saberi et al., 2020; Lemmetyinen et al., 2018). However, despite the potential for spaceborne microwave determination of SWE, there are key limitations inherent within each of the current approaches that affect the sensitivity and uncertainty of the measurements and prevent adoption of the methods for wide-scale SWE monitoring.

Passive-microwave SWE estimation is based on the measurement of the attenuation of Earth-emitted thermal microwave radiation by an overlaying snowpack (Kunzi et al., 1982). It has been implemented on several spaceborne systems, generating global SWE maps with a revisit frequency of up to daily (e.g., Takala et al., 2011). However, the passive-microwave method has severe limitations (Takala et al., 2011). These are as follows: (1) increased uncertainty for larger SWE (> 150 mm) due to attenuation saturation, (2) dependence on properties of snow microstructure, (3) signal contamination from topographic variation, which restricts the method to low-relief regions, and (4) coarse spatial resolution, e.g., > 25 km in the case of Advanced Microwave Scanning Radiometer (AMSR)–E daily products.

In contrast to passive systems, active-microwave SWE estimation methods are based on interpreting variation in backscattered radiation following interaction with the snowpack. One method that promises high spatial resolution is dual-frequency active microwave estimation, which measures total SWE in a single pass of a synthetic aperture radar (SAR) based on the difference between the volume backscatter power returned by the snowpack by two radio frequency bands (e.g., X-band vs. Ku-band). The dual-band requirement of this method currently limits it to ground-based and airborne implementations, although method-capable spaceborne missions have been proposed (e.g., Rott et al., 2012; Derksen et al., 2019). Limitations of this method include a dependence on snow microstructure and an applicability restriction to predominately dry-snow and vegetation-free areas (Oveisgharan et al., 2020).

Active-microwave-based estimation of dry-snow SWE change can also be achieved through repeat-pass SAR interferometry (InSAR) by exploiting the phase
contribution that results from refraction in dry snow (Guneriussen et al., 2001; Deeb et al., 2011; Leinss et al., 2015). This method measures
temporal SWE changes between repeat acquisitions rather than total SWE; the latter must be inferred later through integrating the SWE change maps over
time. Performance depends on sensing frequency: higher frequencies (e.g., X-band compared to C- or L-band) have better phase sensitivity to SWE change
but are more prone to temporal decorrelation and to interference from volume scattering. This method promises high spatial resolution, but it is
constrained by the following requirements: (1) observability is limited to conditions in space and time with low decorrelation, (2) a spatial
reference with known SWE change since the InSAR phase carries an unknown offset, and (3) phase unwrapping to resolve 2*π* ambiguities, which occur
in significant number for typical snowpacks (Rott et al., 2003; Leinss et al., 2015). Decorrelation increases with liquid water content, changes in
snow distribution, and volume scattering. Particularly, due to the need and difficulty of phase unwrapping, the method has not yet been operationally
adopted for wide-scale, repeat monitoring of SWE despite being proposed and demonstrated over 20 years ago (Guneriussen et al., 2001).

Furthermore, polarimetric refraction-based methods have been proposed to exploit the structural anisotropy of snow to provide additional information about SWE change within a snowpack (Leinss et al., 2016, 2020), although in this article we focus on single-polarimetric methods.

The Delta-K method for SWE estimation (Engen et al., 2004) is a variant of the InSAR method that differences the InSAR dry-snow phase between two radar carrier frequency sub-bands. This approach has the benefits of not requiring a spatial phase reference and eliminates the need for phase unwrapping. However, these are at the cost of a much-reduced sensitivity that requires substantial spatial smoothing to suppress the noise, which results in only coarse spatial resolution SWE estimates. For example, Engen et al. (2004) report 100 mm accuracy at 10 km resolution when applying the method to European Remote-Sensing Satellite (ERS-1) (C-band) data.

We hypothesize that it is possible to overcome some of the key limitations of repeat-pass InSAR for dry-snow SWE change estimation, at the cost of only moderately reduced spatial resolution, by exploiting the sensitivity of the dry-snow refractive InSAR phase signal to topographic slope. Our objective is to develop a spatial analysis method for repeat-pass InSAR that allows for unambiguous SWE change estimation directly from the wrapped InSAR phase without the need for a spatial phase reference and phase unwrapping and that is insensitive to other InSAR phase components not correlated with topographic slope at the scale of the output resolution (size of SWE estimation window). Our method exploits the dependency of the refractive phase delay within the snowpack with respect to the local terrain slope. This novel method is conceptually similar to established methods for estimating temporal changes in the layered tropospheric delay in repeat-pass InSAR (e.g., Lin et al., 2010) and offers the same benefits as the Delta-K method but with substantially better sensitivity and spatial resolution for most terrains.

The proposed method was tested with a multi-year RADARSAT-2 dataset situated in the western arctic region of Canada. The study site was chosen for three primary reasons: (1) relatively long dry-snow season, (2) moderate topographic variation, e.g., more than flat prairie-like conditions but less than in an alpine region, and (3) dataset availability since an archived stack of 120 Spotlight images over 9 years is somewhat unusual and allows for detailed temporal analysis of results.

Our paper is organized as follows: first we describe the study site, SAR dataset, and ancillary data that we use to demonstrate and test the novel method. We then develop a quantitative model for the repeat-pass InSAR phase resulting from dry-snow refraction that explicitly considers topographic slope variation. Subsequently, we introduce a spatial-correlation-based method for estimating SWE change from a repeat-pass interferogram, first from a locally unwrapped spatial phase signal and then directly from the wrapped phase. We compare the expected performance and precision of our method to that of the Delta-K method. Next, we show results from applying our novel SWE retrieval method to the full RADARSAT-2 InSAR dataset. We compare these results to several in situ SWE transect surveys located within the imaged footprint, as well as the SWE temporal history estimated by the ERA5 climate reanalysis model. We then discuss sources of estimation error in our new approach and establish their relative significance. Sources considered include the effects of (1) violated model assumptions, (2) other InSAR phase components, and (3) digital elevation model (DEM) errors. Finally, we present our conclusions with a focus on the feasibility of our method for SWE change mapping in the context of existing methods.

## 2.1 SAR dataset and study site

The InSAR dataset used for this study is a time series of 120 RADARSAT-2 Spotlight-A (RS2-SLA) descending-pass single-look complex (SLC) images
spanning the period 22 February 2012 to 29 January 2021 at 24 d intervals, and it is nearly continuous with only a few time gaps
(Table 1). The study area is situated in the Northwest Territories, Canada, at the eastern margin of the Mackenzie Delta, about
120 km south of the Beaufort Sea coast. It covers 238 km^{2} defined by the SAR imaging footprint and includes the town of Inuvik. The
mean backscatter amplitude image for the dataset is shown in Fig. 1a.

Our proposed SWE change estimation method requires knowledge of the variation in topographic slope. To derive the slope, we used a 12 m spatial resolution TanDEM-X digital elevation model (DEM) provided by the German Aerospace Centre (DLR). This DEM was derived from source images acquired between 13 January 2011 and 8 August 2014 including snow-season acquisitions, so some snow-related bias may affect the estimated elevations. Regarding DEM accuracy, our method described in Sect. 3.3 depends on local slope variations and is therefore sensitive only to local relative rather than absolute DEM errors. Wessel et al. (2018) report a root-mean-square error (RMSE) of 1.8 m for TanDEM-X elevations over forested terrain. We also used this DEM for InSAR topographic phase correction and as an input for snow transport modelling, discussed in Sect. 5.2. The DEM shows that the Mackenzie Delta (near sea level) in the eastern portion of the study area has little topographic variation, but to the east of the delta's margin, elevation increases to a maximum of approximately 200 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{m}.\mathrm{s}.\mathrm{l}.$ (Fig. 1b). Several riparian zones drain this upland area towards the delta and create substantial local-scale topographic variation.

The study area is south of the treeline and ecologically is broadly categorized as subarctic boreal forest. However, according to the North American
Land Change Monitoring System (NALCMS) 2015 dataset, which is a 30 m resolution land classification based on Landsat and RapidEye imagery using the
methodology described by Latifovic et al. (2017), there is significant variation in vegetation classification within the study area
(Fig. 1c). The upland area east of the delta features forest, shrubland, and grassland, whereas the delta proper is predominantly forested
but is interspersed with areas of shrubs and sedge wetlands (Burn and Kokelj, 2009). There are 6.5 km^{2} of developed lands in the study
area that include the centrally located town of Inuvik, segments of the Inuvik–Tuktoyaktuk highway (ITH) and the Dempster highway that extend, respectively,
northward and westward from the town site, and the Inuvik airport in the southeastern corner of the study footprint.

The study area lies within the (spatially) continuous permafrost region (Heginbottom, 1995; Nguyen et al., 2009). Permafrost
is ground that remains at or below 0 ^{∘}C for 2 or more years that is overlain by an active layer that thaws seasonally. The
permafrost in the study area is commonly ice-rich (Burn et al., 2009) and occurs within surficial deposits of several origins that are
predominantly unconsolidated (Fig. 1d). As a result of climate warming in the region, permafrost temperatures and active layer thickness
are increasing and producing ground surface settlement where there is net thaw of near-surface, ice-rich permafrost (Kokelj et al., 2017; O'Neill et al., 2019).

As shown in Fig. 2, Inuvik normally receives some level of snowfall for all months of the year except July (ECCC, 2021a). The mean annual snowfall is 1586 mm. For dry-snow conditions, temperatures must remain below freezing; the daily maximum temperatures in Inuvik normally remain below zero for the months of October to April. At upland locations, dry snow in the area is redistributed by wind, and large drifts develop in relation to topographic relief and shrub growth, while in the delta proper, forest and shrub vegetation combine with the lack of relief to minimize wind effects on the snow pack (Burn et al., 2009; Mackay and MacKay, 1974; Palmer, 2007; Morse et al., 2012).

## 2.2 In situ measurements

A number of mid-to-late winter snow transect surveys were compiled (Eppler, 2022) to facilitate method validation which is described in Sect. 4. These are plotted in Fig. 3 and summarized in Table 2, and their positions A–E are labelled in Fig. 1a. These surveys consist of (1) three transects at sites A, B and C that are immediately north of Inuvik, all surveyed by the authors on the same day, (2) one transect, located at site D, provided by Tim Ensom (Wilfred Laurier University) in a fully tree-covered area adjacent to a small creek, and (3) a series of four transect measurements made at the same location, site E, by the Government of Northwest Territories – Department of Environment and Natural Resources (GNWT-ENR) at 2-week intervals at a site with mixed trees and shrubs that is adjacent to the Inuvik Satellite Station Facility. The Tim Ensom transect was surveyed with a 6 cm ESC-30 snow sampler, and all other transects were surveyed with a 6 cm Metric Prairie snow sampler.

## 3.1 Background

Guneriussen et al. (2001) first described the method of SAR interferometric repeat-pass SWE change estimation, which we summarize as follows. In
the case of a dry homogenous layer of snow over ground, a microwave radar signal will penetrate the snow layer, be scattered at the ground surface, and
then return through the snow layer. Surface and volume scattering occur at the air–snow interface and within the snowpack, respectively, but for
sufficiently dry snow, it can be expected that these contributions will be much less that the primary ground-scattered return
(Leinss et al., 2015). The dry-snow layer has a higher real permittivity than air and therefore causes refraction at the air–snow interface, which
corresponds to a reduction in the propagation speed of the signal within the snowpack and a corresponding geometric path-length increase. This is
shown in Fig. 4, which defines the local 3D geometry and depicts both snow-free and snow-covered cases which diverge at the point where the
wavefront reaches the snow surface. The unwrapped phase of the SAR signal due to the dry-snow layer (Φ_{s}) is given by

and it increases with both snow density *ρ* (defined dimensionless as the ratio of the gravimetric densities of snow and water) and snow layer
thickness *D*_{s} (projected along the slope normal direction, represented by direction vector ** n**). Functional dependency on

*ρ*is via the real part of the relative dielectric permittivity

*ϵ*, with other parameters in Eq. (1) being the local incidence angle

*θ*and the radar wavelength

*λ*. Note that

*θ*is defined as the angle between

**and −**

*n***(reversed SAR line-of-sight direction vector) and therefore depends on both the magnitude and aspect of the local slope. As such, this geometry can be defined everywhere within a SAR scene given maps of spatially varying**

*l***, expressed in the local (east, north, vertical) coordinate system and DEM-derived slope magnitude and aspect maps.**

*l*In the case of a stratified snowpack, Leinss et al. (2015) have shown that Φ_{s} is well approximated by a single layer with mean
vertical *ρ*,

where *S* is the SWE of the snowpack with depth *Z*_{s} (measured vertically). Note that ${Z}_{\mathrm{s}}={D}_{\mathrm{s}}/\mathrm{cos}\mathit{\alpha}$, where
*α* is the local topographic slope angle, defined as the angle between ** n** and

**. Equations (1) and (2) can be combined to provide the explicit linear relation between dry-snow phase and SWE:**

*z*The relation holds for the general case of repeat-pass interferometric dry-snow phase between any two dry-snow states by replacing total SWE, *S*,
with SWE change, Δ*S*, and interpreting *ρ* as the magnitude of the mean *ρ* of the removed or added snow layer. As such,
Δ*S* can be either positively or negatively valued and *ρ* ≥ 0.

## 3.2 Dry-snow phase dependence on local topographic slope

SAR images are formed by focusing complex-valued pulse returns obtained sequentially along a 1D flight track (e.g., satellite orbit
segment in the case of spaceborne systems). Each image resolution cell is formed by coherently combining signals received over a finite
aperture. Φ_{s} depends on the local incidence angle, which varies along the aperture. The net effect on the focused
image Φ_{s}, however, is well approximated by considering only the Doppler centroid (i.e., beam center) incident angle at the target
(Eppler and Rabus, 2021).

For a focused image, in a local spatial region with constant topographic slope where incidence angle is uniform (gradual across-swath variations can
be neglected), a spatially uniform snow layer results in a uniform Φ_{s}. In contrast, a local spatial region with topographic slope
variations will generate Φ_{s} variations corresponding to variations in the local incidence angle and projection of the snow depth onto
the slope normal. The validity of assuming a uniform snow layer under a local window and selecting an appropriate spatial scale are discussed in
Sect. 5.2 and 5.5.

Figure 5 depicts the topographic modulation of dry-snow phase explicitly at foreslope, horizontal, and backslope locations. The local incidence angle is less on the foreslope and becomes greater along the transition towards the backslope. The dry-snow phase increases with incidence angle, and therefore this effect causes the phase contribution to increase along the transition from foreslope to backslope.

The projection of the snow depth onto the slope normal scales with the cosine of the slope angle and is therefore greatest in horizontal areas and
less on both the foreslope and backslope. Considering Eq. (3), if the topographic variation within a SAR scene is known, then the spatially
variable Φ_{s} sensitivity to a uniform SWE layer (*ξ*) can be computed as

where we assume that *ρ* is known. A spatially varying change in dry SWE during a time interval spanned by two SAR acquisitions will cause a
corresponding spatially variable Φ_{s}, modulated by the topographic sensitivity:

This effect is depicted in Fig. 6, which shows *ξ* computed over the RS2-SLA spatial footprint and a late snow-season interferogram which
shows spatial phase variations that are correlated with *ξ*.

## 3.3 Estimating absolute ΔSWE by exploiting topographic variation

According to Eq. (5), if Φ_{s} can be recovered, then the spatially varying SWE change can be directly estimated at the same
spatial resolution as the interferogram. However, this requires phase unwrapping which can be difficult for snow-covered interferograms due to both
aliasing and local incoherence of the InSAR phase. Furthermore, even if the total unwrapped phase (Φ) can be determined, the phase offset is
unknown. One or more reference targets with known SWE change (ΔSWE) are required to estimate the offset. For these reasons, direct
estimation of ΔSWE from interferograms is not straightforward. However, we suggest that the effect of topographic modulation
of Φ_{s}, described in Sect. 3.2, allows an absolute estimate without requiring
phase unwrapping or knowing the phase offset.

Consider a local estimation window, e.g., 1 km × 1 km. The phase within this window consists of the superposition
of Φ_{s} and other phase components (e.g., due to decorrelation, surface displacement, imperfect topographic correction, soil moisture,
or tropospheric delay). For the purposes of estimating ΔSWE, we treat the other phase components as a single error
term, Φ_{e}. Due to the difficulty of phase unwrapping and offset estimation we assume that only the locally de-meaned superposition
phase ($\stackrel{\mathrm{\u0303}}{\mathrm{\Phi}}$) is available:

where 〈⋅〉 denotes the mean over the window. Initially, for simplicity we assume also that phase variation within the window is
limited to a 2*π* ambiguity interval so that (unwrapped) $\stackrel{\mathrm{\u0303}}{\mathrm{\Phi}}$ is easily obtained from the corresponding wrapped-phase, *ϕ*, within the
window:

where ^{∗} denotes complex conjugation.

The unknown spatially variable SWE change within the window can then be described by

where $\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}S}$ corresponds to the horizontal SWE change variation. *ξ* can be similarly decomposed as

where $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$ is the variation in *ξ*. Substituting Eqs. (8) and (9) into Eq. (5) and expanding yields

Substituting Eq. (10) into Eq. (6) and noting that $\langle \stackrel{\mathrm{\u0303}}{\mathrm{\Delta}S}\rangle =\langle \stackrel{\mathrm{\u0303}}{\mathit{\xi}}\rangle =\mathrm{0}$ gives

Assuming $\stackrel{\mathrm{\u0303}}{\mathrm{\Phi}}$ correlates with $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$ with the proportionality 〈Δ*S*〉, we introduce the following
correlation-based estimator for 〈Δ*S*〉:

This estimator correlates $\stackrel{\mathrm{\u0303}}{\mathrm{\Phi}}$ with $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$ and relies on the assumption that the horizontal SWE change distribution is uniform within
the window and that Φ_{e} components are uncorrelated with $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$. The bias of $\widehat{\mathrm{\Delta}S}$ with respect to 〈Δ*S*〉 is given by

where *E*(⋅) denotes the expectation. Therefore, $\widehat{\mathrm{\Delta}S}$ is biased by the components of $\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}S}$ that are
systematically correlated with either $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$ or ${\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}^{\mathrm{2}}$, as well as by Φ_{e} components systematically correlated with
$\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$. The significance of the different bias terms is examined in Sects. 5.2
and 5.3.

## 3.4 Wrapped-phase (“SlopeVar”) estimator

The estimator defined in Eq. (12) requires $\stackrel{\mathrm{\u0303}}{\mathrm{\Phi}}$ within the spatial window of consideration, and we have assumed so far that this
can be obtained from *ϕ* through Eq. (7) implicitly assuming phase variations in the window are limited to a maximum of 2*π*. As shown
by Just and Bamler (1994), decorrelation to any degree results in a per-target phase distribution that is not bound by the ±*π* interval, and so
this condition is never strictly met. Except for high-coherence areas with limited short-scale phase variation, either the phase within the window must
be unwrapped or an alternative wrapped-phase version of the estimator is required.

Here we propose such an alternative estimator (which we denote the “SlopeVar” estimator) that uses periodogram optimization to iteratively estimate
〈Δ*S*〉 using the set of wrapped phase *ϕ* under the window

Considering Eq. (11) and the fact that *e*^{jϕ}=*e*^{jΦ}, Eq. (14) can be expanded to show the components of
*ϕ*,

In the case where Φ_{e}=0 and $\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}S}=\mathrm{0}$, this yields an unbiased estimate, i.e., ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}=\langle \mathrm{\Delta}S\rangle $, otherwise the $\langle \mathit{\xi}\rangle \stackrel{\mathrm{\u0303}}{\mathrm{\Delta}S}$, $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}S}$, and $\u2329\stackrel{\mathrm{\u0303}}{\mathit{\xi}}\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}S}\u232a$ terms in Eq. (15) will bias ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$.

Figure 7 shows an example of the SlopeVar estimator being applied to a 1 km×1 km
(ground range × azimuth) window for the inset area labelled in Fig. 6. $\stackrel{\mathrm{\u0303}}{\mathrm{\Phi}}$ is moderately correlated with $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$
over the window as seen by comparing Fig. 7a and b and by examining the 2D histogram in
Fig. 7d. Figure 7e shows the periodogram with a distinct peak at ΔSWE of 28 mm. Figure 7c shows
the effect of removing the Φ_{s} component assuming constant ΔSWE of 28 mm under the window, which results in a
much more uniform phase, although residual variations can still be seen.

## 3.5 Estimator implementation

The estimator was implemented as a Python software module (Eppler, 2022). Input interferograms were generated by co-registering the set of RS2-SLA SLCs to a common master SLC and computing
sequential topographic phase-corrected 3 × 12 (range pixels × azimuth pixels) multi-looked interferograms using the GAMMA software
package (GAMMA Remote Sensing AG, Switzerland). Multi-looking was performed to better match the sample spacing of the interferograms to that of the
DEM-derived *ξ* map. This improved estimator computation runtime and was seen to have a negligible effect on the resulting ΔSWE
estimates. All interferograms were visually inspected, and those with sparse coherence, typically those spanning the spring melt and early winter
periods, were discarded from further analysis.

For simplicity, we only processed sequential interferograms since these provide the best coherence during the dry-snow season. However, for distributed scatterers which predominate in natural terrain, additional interferograms (e.g., 48, 72 d, …) are known to allow for a reduction in the variance in estimated sequential phases using a phase-linking estimator (e.g., Guarnieri and Tebaldini, 2008). This is because such scatterers exhibit complex Gaussian scattering statistics which are not fully characterized by the set of sequential interferograms alone.

We found the estimator to be sensitive to errors in *ξ*, especially those resulting from DEM errors and interpolation errors. High-spatial-frequency
errors in the *ξ* map tend to bias the estimated ΔSWE magnitude towards zero (discussed in Sect. 5.4). We initially used
the raw TanDEM-X DEM, which we cubic-resampled to the multi-looked SAR geometry, to compute the slope and aspect angle maps required to
compute *ξ*. The resulting *ξ* map contained high-frequency artifacts from the cubic interpolation, as well as substantial errors due to stitching
artifacts in the raw DEM. We did not thoroughly investigate the issue of the most appropriate interpolator, and so it may be possible to reduce these
errors by using a different interpolator. Estimator performance improved significantly after we (1) manually identified and masked local DEM stitching
artifacts and interpolated over them, (2) smoothed the DEM using a 2D Gaussian filter with a standard deviation of three DEM pixels (resulting in
approximately 90 m spatial resolution) to reduce high-frequency DEM errors, and (3) computed the slope angle maps prior to resampling to the
SAR geometry. Note that to compute *ξ* according to Eq. (4), we computed *ϵ* of dry snow as $\mathit{\u03f5}\left(\mathit{\rho}\right)=\mathrm{1}+\mathrm{1.5995}\mathit{\rho}+\mathrm{1.861}{\mathit{\rho}}^{\mathrm{3}}$, according to Mätzler (1987), and for simplicity assumed *ρ*=0.3 for all estimates (discussed further in
Sect. 5.1).

We then implemented the SlopeVar estimator described in Sect. 3.4. For each candidate value of
Δ*S*_{w}, the periodogram magnitude, $\mathrm{|}\langle {e}^{j(\mathit{\varphi}-\mathrm{\Delta}{S}_{\mathrm{w}}\stackrel{\mathrm{\u0303}}{\mathit{\xi}})}\rangle \mathrm{|}$, can be computed
efficiently using a 2D fast-Fourier-transform-based rectangular smooth filter over the complex-valued interferogram, allowing for the generation of
spatially oversampled maps at the interferogram sample spacing instead of the much coarser estimation window size. We evaluated the periodogram over
the ΔSWE range of [−50 mm, 80 mm] with an interval spacing of 2 mm. Peak locations for each spatial point
were then refined by fitting a quadratic by least squares using three points centered on the maximum value found over the search grid.

Solutions where the maximum was within two grid intervals from the edge of the search range were labelled as invalid. This provides a means for excluding poor results from low-coherence areas since for these areas the periodogram analysis typically does not result in a peak within the search range. We considered adding an additional exclusion threshold based on coherence magnitude but found it unnecessary. Water-body areas were also labelled as invalid.

Figure 8a shows the resulting ΔSWE map after applying the estimator to all areas of the example interferogram depicted in Fig. 6 using a 500 m square estimation window.

## 3.6 Monte Carlo estimation of ΔSWE estimation variance

The amount of decorrelation can vary significantly both spatially within each interferogram and temporally depending on changing surface conditions and pass-to-pass variations in the SAR orbit. It is therefore important to estimate the ΔSWE variance contribution as a space- and time-varying quantity.

Specifying the transfer function between phase noise and ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ variance is made difficult by the fact that the
distribution of $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$ is generally unique for each estimation window. We decided to apply the Monte Carlo approach by simulating a
sufficiently large ensemble of zero-signal partially decorrelated interferograms, applying the ΔSWE estimation on each instance, and then
computing statistics over the ensemble of estimates. We approximate the phase noise as a circular complex stationary white process with zero mean
phase as described by Just and Bamler (1994) and which is uniquely specified by the coherence magnitude, *γ*. We estimate *γ* as the
ensemble coherence of the residual phase after removing Φ_{s} estimated using ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$:

Figure 8b–d, respectively, show the residual phase coherence map, the standard deviation
of *ξ* under the window (i.e., $\langle {\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}^{\mathrm{2}}{\rangle}^{\mathrm{1}/\mathrm{2}})$, and the resulting Monte-Carlo-based error estimation for the same
interferogram depicted in Fig. 6. Local factors affecting the ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ variance include (1) the residual
coherence magnitude, (2) the diversity of *ξ* within the estimation window, and (3) the spatial support under the window (areas adjacent to water
bodies have a smaller ensemble of valid phase samples and hence a larger variance).

## 3.7 Precision comparison with Delta-K ΔSWE estimator

Our method has some similarity with the Delta-K method for estimating ΔSWE in that both methods exploit a diversity in the dry-snow
phase sensitivity so that absolute ΔSWE can be estimated without the need for phase unwrapping. Furthermore, both methods operate over
an estimation window, yielding results at a resolution lower than the single-look resolution. In order to compare the two methods, we consider their
relative precision estimated over a spatial window with *N* multi-looked-resolution cells and assume Gaussian phase noise, which is reasonable
when multi-looked interferograms with sufficient looks are used (i.e., > 10 as noted by Hagberg et al., 1995). Note that for the
Delta-K method, multi-looking of sub-band interferograms prior to computing the band difference phases has been suggested by
Bamler and Eineder (2005) to preserve the Gaussian phase noise assumption and shown by De Zan et al. (2015) to yield better precision
than the alternative of first computing the sub-band difference and then multi-looking.

The Delta-K method exploits the *ξ* diversity that exists between separate radar-frequency sub-bands corresponding to $\mathit{\xi}\pm \mathrm{\Delta}\mathit{\xi}/\mathrm{2}$, where
*ξ* in this case is the value at the central carrier frequency, and $\mathrm{\Delta}\mathit{\xi}=B(\mathrm{1}-b)\mathit{\xi}$, where $B\in [\mathrm{0},\mathrm{1}]$ is the bandwidth of
the SAR normalized with respect to carrier frequency and $b\in [\mathrm{0},\mathrm{1}]$ is the sub-bandwidth fraction used. Assuming multi-looked
phase noise with variance ${\mathit{\sigma}}_{\mathrm{\Phi}}^{\mathrm{2}}$, the mean sub-band interferometric phase will have the variance ${\mathit{\sigma}}_{\mathrm{\Phi}}^{\mathrm{2}}/bN$, and therefore
the Delta-K difference phase will have the variance ${\mathit{\sigma}}_{\mathrm{\Phi}\mathit{\_}DK}^{\mathrm{2}}=\mathrm{2}{\mathit{\sigma}}_{\mathrm{\Phi}}^{\mathrm{2}}/bN$. Delta-K estimates ΔSWE according
to ${\widehat{\mathrm{\Delta}S}}_{\text{DK}}={\mathrm{\Phi}}_{\text{DK}}/\mathrm{\Delta}\mathit{\xi}$, and therefore the ${\widehat{\mathrm{\Delta}S}}_{\text{DK}}$ standard deviation
(*σ*_{ΔS_DK}) is

For the SlopeVar method, the estimation of ΔSWE corresponds to a linear regression over *N* samples, and therefore, assuming
$\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$ is Gaussian-distributed under the estimation window with variance $\langle {\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}^{\mathrm{2}}\rangle $, the estimation
standard deviation (*σ*_{ΔS_SV}) is

The relative estimation precision ratio between the two methods is therefore given by

where we have replaced *ξ* in Eq. (17) with the mean under the estimation window, 〈*ξ*〉. We use $b=\mathrm{1}/\mathrm{3}$, which,
neglecting system noise, is the optimum value for minimizing *σ*_{ΔS_DK}, and *B*=0.0185, which is the normalized range bandwidth
for the RADARSAT-2 Spotlight-A mode. Figure 9 shows the distribution of ${\mathit{\sigma}}_{\mathrm{\Delta}S\mathit{\_}\text{DK}}/{\mathit{\sigma}}_{\mathrm{\Delta}S\mathit{\_}\text{SV}}$
computed for all 500 m square estimation windows in the RS2-SLA study footprint, excluding water bodies. The distribution has a median value
of 3.8 and is > 1 for 97.5 % of the area, indicating that the SlopeVar estimator can be expected to provide substantially better
ΔSWE precision compared to Delta-K for the great majority of areas within the study footprint.

## 4.1 ΔSWE estimation maps

We computed 95 coherent sequential interferograms from the assembled set of 120 RS2-SLA SLCs, and these were used with the SlopeVar estimator to generate ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ maps using an estimation window of 500 m × 500 m (ground range × azimuth). This corresponds to 46 snow-season interferograms over two partial and eight complete snow seasons and 49 non-snow interferograms over nine snow-free seasons. The estimator was run on the snow-free-season interferograms to serve as control cases since these are known to have zero dry-snow ΔSWE. Additionally, the Monte-Carlo-based variance estimator described in Sect. 3.6 was run on all interferograms. Water bodies and areas with no valid solution from the periodogram peak-finding analysis were labelled as invalid and masked out. Figure 10 shows an example set of five such ΔSWE maps from the 2017–2018 snow season (all snow-season ΔSWE maps from the 2012–2021 dataset are shown in Fig. S1 in the Supplement).

## 4.2 Comparison of ΔSWE estimates with in situ measurements

To compare the results with the in situ transect measurements presented in Sect. 2.2, we generated snow-season cumulative SWE maps by integrating the sequential ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ maps through time, starting with the first dry-snow interferogram identified for each snow season. As discussed in Sect. 3.5, some ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ values are labelled as invalid due to a lack of a distinct periodogram peak. On average these correspond to about 10 % of non-water spatial samples, and therefore, in order to preserve temporal continuity for the integration, these were replaced with values obtained by smoothing all valid estimates within a 2 km square window. The cumulative SWE maps do not represent estimates of true absolute SWE because early snow-season interferograms are generally not sufficiently coherent to allow ΔSWE estimation, and therefore some amount of accumulated SWE is not accounted for. Therefore, for the comparison, the in situ measurements correspond to an upper limit for the expected SlopeVar-estimated total SWE maps, neglecting any biases. However, given the error sources described in Sect. 5, actual estimates may exceed this expected upper limit.

We also estimated variances for the cumulative SWE maps by summing the Monte-Carlo-derived per-interval variance estimates, assuming they are uncorrelated through time.

For each transect dataset, all SWE measurements were averaged to generate a mean value for the local area and the sample standard deviation was used to compute a corresponding standard error. This was done because the transect lengths, listed in Table 2, were 200 m or less, being well within the 500 m SlopeVar estimator window size. The in situ measurements were each made on a particular date, whereas the InSAR cumulative SWE maps are produced at 24 d intervals. The cumulative SWE maps temporally bracketing each transect date were therefore sampled at the spatial mean position of the transect and then linearly interpolated in time to generate a cumulative SWE estimate for that point in space and time. Likewise, the error variance maps were also interpolated and converted to standard deviations.

Figure 11 shows a scatterplot comparison of the InSAR-derived cumulative SWE estimates and the transect mean SWE values. Vertical and
horizontal error bars correspond to 1 standard deviation and 1 standard error, respectively. Seven of the eight transects are within 1 standard
deviation of the 1:1 line and demonstrate good agreement between the measured and estimated values. Treating the transect mean values as truth and
neglecting the unaccounted early snow-season SWE, the RMSE for all transect comparisons is 14.8 mm and the bias is −6.6 mm. It is
noteworthy that the author-surveyed transects (locations “A”, “B”, and “C”) correspond to areas of low topographic variation as indicated by the
two rightmost columns of Table 2, and, as shown in Fig. 11, locations “A” and “C” have larger predicted total SWE standard
deviations (i.e., > 25 mm). This is consistent with the fact that estimation variance is expected to increase with decreasing *ξ*
diversity.

## 4.3 Comparison of ΔSWE estimates with predicted temporal change

The estimator results were examined temporally by comparing the spatial mean of each ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ map, after projection to ground-range geometry, with the corresponding SWE change predicted by the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 reanalysis model (Hersbach et al., 2020) over the same time interval. We chose ERA5 for comparison based on Mortimer et al. (2020), who report that it compares favourably with independent snow course data and more so than currently available global passive microwave products do. The gridded ERA5 results were spatially interpolated to the RS2-SLA scene center. Note that all but five of the analyzed interferograms span 24 d intervals. The additional five intervals each span 48 d due to missed SAR acquisitions. For the temporal comparison, each of these five intervals were split into two 24 d intervals with half the ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ arbitrarily assigned to each interval.

Figure 12a compares the continuous ERA5 24 d interval ΔSWE history time series with the SlopeVar estimator results. Figure 12b shows the same data as a scatterplot comparison. Both plots are labelled according to the time of year as four annual quarters to highlight seasonal effects. Treating the ERA5 estimates as truth, Table 3 shows error statistics (bias, RMSE, and correlation coefficient) for a number of temporal subsets that include the four seasonal quarters, all intervals, the nominal dry-snow period (October–March), and the “no snow” intervals (all intervals with zero ERA5 SWE at both interval endpoints).

Regarding correspondence with ERA5, we are primarily concerned with the October–December and January–March estimates because these correspond to the dry-snow period when the estimator is expected to provide useful results. Figure 12b shows that these do follow the 1:1 correspondence line. However, the October–December points appear to have a significant positive bias, whereas the January–March points do not. This is quantified in Table 3, which shows that, in aggregate, the October–March time period has a correlation of 0.57 with respect to ERA5, which is moderately strong. There is an 8.6 mm per-interval ΔSWE bias for October–December but only 1.7 mm for January–March. This is consistent with the expected impact of early freeze-season heave, which is discussed further in Sect. 5.3.2. The no-snow intervals are of interest because they serve as a control set with known zero ΔSWE. These show a bias of +1.7 mm, which is discussed further in Sect. 5.3.2. The RMSE of the no-snow set is 4.2 mm, which serves as a measure of the estimation uncertainty due to factors other than ΔSWE inhomogeneity. The April–June intervals are included for completeness even though these are not expected to agree because this corresponds to the snowmelt period. This is especially apparent for the two intervals with large negative ERA5 ΔSWE values corresponding to periods of intense snowmelt. The SlopeVar estimates for these periods are near zero likely because only snow-free areas (snow already melted) are coherent and hence contribute to the estimates, and seasonal active-layer thaw has not yet commenced.

## 4.4 Distribution of the estimated ΔSWE

Spatial variations in the ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ maps are due to a combination of true ΔSWE spatial variations, stochastic
error due to decorrelation, and spatially variable biases as described in Sect. 5. An independent set of densely sampled
ΔSWE measurements over our study area does not exist, which precludes direct estimation of the errors. Instead, we compare temporally
averaged seasonal subsets of the ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ maps since we expect the bias to differ between seasons. Figure 13a, c,
and e show the temporal mean ΔSWE maps for the early snow season (October–December), late snow season (January–March), and all
snow-free intervals. These show a moderate degree of negative correlation between the October–December and snow-free maps, especially in the upland
area east of the delta margin. The Pearson correlation coefficients between the three maps are
*R*_{October–December/January–March} = 0.28; *R*_{October–December/Snow-free} = −0.52; and
*R*_{January–March/Snow-free} = −0.17. These are consistent with a dominant bias contribution from surface normal heave and subsidence
contributing to the October–December and snow-free subsets but with opposite polarity. The weaker correlations between these and the January–March
subset suggest a weaker surface-normal bias contribution during this period which is consistent with the expected annual cycle of active-layer heave
and subsidence discussed further in Sect. 5.3.2.

Figure 13b, d, and f show the normalized distributions for all ΔSWE
estimates partitioned by season. Since the snow-free intervals have a known ΔSWE of zero, a global RMSE value can be computed. This value,
computed over all non-water areas and all snow-free maps, is RMSE_{snow-free} = 21 mm. We know from Eqs. (13)
and (30) that both additive and multiplicative bias contributions are modulated by the local *ξ* variance inverse. We repeated the RMSE calculation
for the top *ξ* variance quartile of spatial samples (corresponding to threshold of $\langle {\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}^{\mathrm{2}}{\rangle}^{\mathrm{1}/\mathrm{2}}$
> 7.6 rad mm^{−1}), resulting in RMSE_{snow-free} = 15 mm, which shows that the impact of bias can be reduced
by restricting the estimator to areas with greater *ξ* diversity, albeit with a loss of spatial coverage. Normalized distributions for this
restricted subset are also shown in Fig. 13b, d, and f.

## 5.1 Snow density misspecification

Our proposed method relies on knowing *ξ*, which, as shown in Eq. (3), depends on *ρ*. Therefore *ρ* must either be known a
priori or an assumed value must be specified. Here we assess the effect of misspecifying *ρ*. Since the method estimates ΔSWE from the
interferometric phase, the required *ρ* corresponds to the snowpack change rather than the actual bulk *ρ* of the snowpack. For the case of
fresh snowfall on a dense late season snowpack, the *ρ* of the added layer will be less than the mean *ρ* before and after the change.

A misspecification of *ρ* will bias the resulting ΔSWE estimate. This effect can be quantified by parameterizing *ξ* with respect
to *ρ* so that Φ_{s}(*ρ*)=*ξ*(*ρ*)Δ*S*. The fractional ΔSWE bias caused by specifying *ρ*^{′} rather than
the true value, *ρ*, is

This fractional bias is shown in Fig. 14a for true and assumed *ρ* over the interval [0, 0.5], and it indicates that
assumed *ρ* greater or less than the true *ρ* results in a positive or negative bias, respectively. Figure 14a also shows that the
bias magnitude is < 5 % of the true SWE change for the considered *ρ* range.

In the absence of any ΔSWE, snowpack evolution may result in a bulk *ρ* change (e.g., due to settling) which will cause a
Φ_{s} contribution. This contribution will be spuriously attributed to a SWE change. Consider the case of a snowpack with total SWE of
*S* that undergoes a *ρ* change from initial *ρ*_{1} to final *ρ*_{2}. According to Eq. (5), this will cause a Φ_{s}
change:

Assuming that the temporal mean *ρ* of $({\mathit{\rho}}_{\mathrm{1}}+{\mathit{\rho}}_{\mathrm{2}})/\mathrm{2}$ is somehow known and used for estimation, the resulting ΔSWE
error, *S*_{e}, as a fraction of total SWE will be

This fractional error is shown in Fig. 14b for initial and final *ρ* over the interval [0, 0.5]. This shows that a
*ρ* increase or decrease will, respectively, cause a negative or positive ΔSWE error. Since snowpack evolution generally results in
densification, this effect will tend to negatively bias estimated ΔSWE. However, for the considered *ρ* range, the bias
is < 5 % of the total SWE.

For simplicity we assumed a value of *ρ* = 0.3 for all ΔSWE estimates. This value is likely too large for the dry-snow season in
Inuvik. For example, the mean snow density measured across all transects summarized in Table 2 is *ρ* = 0.17, and according to
Eq. (20), the use of the assumed value results in a 2.1 % bias which is small compared to the other errors discussed in
Sect. 5. This small bias could be mostly mitigated by assuming a more appropriate *ρ* according to the location and time of year
(e.g., Sturm and Wagner, 2010).

Regarding the effect of vertical density layering, our study area is prone to wind slab formation in which the late season snowpack can consist of a dense
wind slab overlaying a low-density hoar layer (Rutter et al., 2019; King et al., 2018). Considering the extreme case of near-zero density hoar
overlain by *ρ* = 0.5 wind slab, assuming a uniform *ρ* = 0.3 results in a +2.5 % estimation bias which is still a small error
compared to the other bias sources considered in our analysis.

## 5.2 Violation of ΔSWE horizontal homogeneity assumption

Temporal changes in SWE can be caused by a number of different processes including deposition of new snow, as well as redistribution by wind, melt, and sublimation. Each of these processes is spatially modulated in natural terrain by the configuration of topography and vegetation (Morse et al., 2012; Anderton et al., 2004; Palmer et al., 2012). Furthermore, vegetation distribution is itself modulated by topography (Walker, 2000). For these reasons, spatial variation in ΔSWE over a particular time period is likely to show some correlation with topography, either directly with elevation or with derived indices such as slope, aspect, or measures of curvature.

As shown in Sect. 3.3, the estimation of ΔSWE within a spatial window through the correlation with $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$ is subject to bias when there are horizontal ΔSWE variations within the window that are correlated with either $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$ itself or with ${\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}^{\mathrm{2}}$. To investigate the significance of this effect, we simulated the seasonal evolution of spatially variable ΔSWE distributions due to variations in topography and vegetation using SnowModel (Liston and Elder, 2006), implemented as a Fortran software package that includes sub-models for metrological forcing conditions, surface energy exchanges, snowpack evolution, and 3D wind transport. We simulated an evolving snowpack over several snow seasons for two cases: (1) bare-earth case with no snow-holding vegetation to examine the impact of topographic variation only, and (2) spatially variable vegetation heights (i.e., variable snow-holding capacity) to consider the additional influence of vegetation. The spatially varying vegetation was modelled using the NALCMS 2015 land-classification data for our study area.

We chose a study area and time period for the simulation that correspond to the spatial footprint and a similar time period as the RS2-SLA dataset we
used for the method demonstration described in Sect. 2.1. SnowModel also takes as input the digital elevation model of the
study area and one or more time series of meteorological driving data in the form of regularly sampled temperature, wind, and precipitation data. We
ran the model using the same 12 m TanDEM-X DEM used for InSAR topographic phase correction and SWE change estimation. Therefore, the model did
not simulate variations in SWE at scales finer than 12 m. Note that spatial SWE variations finer than the resolution of the SlopeVAR input
data, *ξ* and *ϕ*, will not bias the estimates but instead will contribute to phase decorrelation and therefore affect the estimation variance,
as discussed in Sect. 3.6.

Meteorological forcing data was input from three locations: (1) Environment and Climate Change Canada (ECCC) Inuvik climate station
(68.32^{∘} N, 133.53^{∘} W), adjacent to position “E” in Fig. 1a and 6 km from the scene center
(ECCC, 2021b); (2) ECCC Trail Valley Creek station (68.74^{∘} N, 133.50^{∘} W), 43 km north of the scene
center; and (3) the Government of Northwest Territories – Department of Environment and Natural Resources (GNWT-ENR) southern ITH station (68.54^{∘} N,
133.77^{∘} W), 21 km north of the scene center (dataset provided directly from GNWT-ENR). SnowModel weights the station inputs
according to an inverse distance squared parameter, and therefore the Inuvik Climate Station dominates the model input. The model was run with default
settings at a daily time step interval with results output every 24 d corresponding to the imaged dates of the RS2-SLA dataset. The resulting
absolute SWE maps, produced in the original DEM spatial reference, were resampled to the RS2-SLA range–Doppler geometry and sequentially differenced to
produce 24 d ΔSWE maps.

These ΔSWE maps were then used along with the *ξ* map, computed from the DEM with a nominal value of *ρ*=0.3, to determine the
expected ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ error according to the first two terms in Eq. (13), i.e., those related to
$\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}S}$. The modelled ΔSWE error maps were computed with rectangular estimation windows with
500 m × 500 m ground footprint.

We analyzed our simulations for all 24 d intervals from 17 July 2011 to 29 April 2020 falling within the expected dry-snow season, nominally chosen to be 1 October to 31 March of each year based on the climate normals for Inuvik (see Fig. 2). The results, summarized in Fig. 15, show the spatial distribution of the ΔSWE error averaged over all modelled time intervals and histograms of all expected ΔSWE errors and the per-interval spatially averaged mean ΔSWE error.

For the bare-earth case (Fig. 15a–c), a single outlier interval (5 January 2012 to 29 January 2012) was omitted from the summary analysis because the errors for that interval were very large compared to all other intervals. This appears to be due to a strong correlation of the simulated spatial ΔSWE distribution with $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$ for that 24 d interval which can occur when a strong wind-transport event directionally aligns with the horizontal projection of the SAR line-of-sight. Hence, the simulated topography-modulated scouring and deposition along the wind-direction axis can, by chance, align with the SAR geometry ground-range direction which results in a very large estimation bias. It may be possible to detect these events by analysis of the wind history (e.g., thresholding of a suitable blown-snow index while also considering the wind direction with respect to the SAR ground-range direction). We investigated this approach by computing the time-varying blowing snow probability as modelled by Li and Pomeroy (1997) which provides the probability that blowing snow conditions will occur as a function of wind speed and temperature. We computed these at daily time steps, assigned the mean wind direction to each time step, and then integrated these 2D vectors over the 24 d simulation intervals. The result is a directional cumulative blowing snow index in units of hours (of blowing snow) and is shown in Fig. 16 for all simulated 24 d intervals. This shows that the SnowModel outlier case is also an outlier with respect to the blowing snow index magnitude (and also when projected to the SAR look direction axis) suggesting that wind-driven redistribution is responsible for the predicted large ΔSWE bias.

Another potential mitigation is to use SAR images from the opposite pass direction. The bias from the near-opposite horizontal direction should have the opposite sign, and therefore the estimated ΔSWE should differ significantly between pass directions, allowing for the effect to be detected and perhaps mitigated by averaging the results.

For the bare-ground case, the temporal mean ΔSWE error map shows mean errors in the ± 2 mm range with spatial features (Fig. 15a) that are generally aligned with topographic variations (Fig. 1b). The distribution of all ΔSWE errors (Fig. 15b) has a mean of just −0.2 mm which corresponds to a small negative overall estimation bias and an RMS value of 2.0 mm which corresponds to approximately 20 % of the simulated mean per-interval ΔSWE value. The histograms of the per-interval spatially averaged errors (Fig. 15c) show that the majority of 24 d intervals have a mean error near zero, suggesting that averaging estimates over wider areas reduces the bias effect introduced by spatially variable ΔSWE distribution.

For the NALCMS 2015 vegetation case, the temporal mean ΔSWE error map (Fig. 15d) shows mean errors with magnitudes that are much larger than for the bare-earth case (up to 100 mm in some areas) corresponding to areas with heterogenous vegetation classifications (see Fig. 1c). The distribution of ΔSWE errors (Fig. 15e and f) has a mean and RMS value of 4.3 mm and 39.6 mm which are very substantial. This suggests that spatially varying vegetation snow-holding conditions may substantially bias the ΔSWE estimates. However, it should be noted that results from real data, presented in Sect. 4, do not show such large variations as would be expected by the error magnitudes in Fig. 15d–f. The reason could be that the simulation driven by the NALCMS vegetation classification significantly over-represents the effect of spatially variable vegetation. Nevertheless, these results highlight the relative importance of vegetation as a source of method bias and suggest that some form of mitigation may be required to yield useful estimates from the proposed method in areas with variable vegetation snow-holding heights.

We compared the relative values of the $\langle \mathit{\xi}\rangle \langle \stackrel{\mathrm{\u0303}}{\mathit{\xi}}\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}S}\rangle $ and $\langle {\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}^{\mathrm{2}}\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}S}\rangle $ contributions to the ΔSWE error which showed that for all 24 d intervals, and both bare-earth and NALCMS 2015 cases, the RMS of the $\langle \mathit{\xi}\rangle \langle \stackrel{\mathrm{\u0303}}{\mathit{\xi}}\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}S}\rangle $ component was much greater than that of the $\langle {\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}^{\mathrm{2}}\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}S}\rangle $ component so that with respect to the influence of $\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}S}$, the spatial correlation with respect to ${\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}^{\mathrm{2}}$ can be neglected.

## 5.3 Correlated phase components

As shown in Eq. (12), additional phase components, other than Φ_{s}, will contribute to ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$
error according to their spatial correlation with $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$: $\langle \stackrel{\mathrm{\u0303}}{\mathit{\xi}}{\mathrm{\Phi}}_{\mathrm{e}}\rangle /\langle {\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}^{\mathrm{2}}\rangle $. We distinguish between (1) systematic spatial correlation between *ξ* and the Φ_{e} components that have a topographic
dependence and (2) spurious correlations with Φ_{e} components that are not systematically dependent on topography. In this section we
consider the former by examining several Φ_{e} components, their spatial characteristics with respect to topographic slope, and their
expected contribution to ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ error.

Φ_{e} components not systematically correlated with $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$ will contribute to the ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ error due
to finite sampling over the estimation window. These include decorrelation phase noise and the non-$\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$-correlated parts of other phase
components. In the case of these spurious correlations, the resulting estimation error magnitude can be expected to diminish with larger estimation
window sizes and, in the global sense, will not systematically bias the estimates but instead contribute to the variance.

### 5.3.1 Atmospheric delay

The troposphere contributes substantial phase delay during SAR imaging (Hanssen, 2001). This delay can be decomposed into static and
dynamic components. Temporal differences in the static component between image acquisitions results in a Φ_{e} component that is
modulated by topography. Neglecting any height dependence of the delay over the height range within the estimation window, the static atmospheric
phase contribution (Φ_{sa}) can be modelled as a simple linear function (Lin et al., 2010),

where Δ*K* is the linear transfer function coefficient corresponding to the volumetric mean delay difference between SAR acquisitions, *h* is
the topographic height, and *β* is an offset constant which we set to zero without loss of generality since we consider only the correlation
of Φ_{sa} with zero-mean $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$. Note that the horizontal variation in Δ*K* is gradual, and therefore Δ*K* can be
assumed constant within the estimation window. Therefore, the spatial variation in Φ_{sa} within the ΔSWE window (assumed to
be ≤ 1 km as discussed in Sect. 5.5) can be attributed solely to variation in *h*.

The dynamic component is driven by turbulent mixing, and therefore its spatial power spectrum follows a power-law relationship, which limits the
amplitude at spatial scales below a few kilometres (Hanssen, 2001). Therefore, for estimation windows < 1 km, the dynamic component can
reasonably be neglected. For larger estimation windows (not explicitly considered in this paper), the influence of the dynamic component can be
mitigated by spatial high-pass filtering both the InSAR phase and *ξ* prior to ΔSWE estimation, similar to the method employed
by Lin et al. (2010) for estimating Δ*K* from interferograms.

Figure 17a and b show the expected ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ error for the case
of Δ*K* = 10^{−5} (i.e., 10 mm km^{−1}). For this case the RMS error is 1.1 mm. However, the distribution is centered
around zero so that the mean error is only 0.04 mm. Therefore, spatial averaging of the estimates over a sufficient area will tend to mitigate
the error due to static atmospheric delay.

Regarding the temporal distribution of Δ*K*, given that it is a pairwise differential quantity, it has zero expectation on an annualized
timescale. However, there is some seasonal dependence on tropospheric delay (Cong et al., 2012), and therefore Δ*K*, for sequential
dry-snow-season interferograms, should not be considered temporally random. For a particular spatial location, the errors may accumulate
constructively over a snow season to yield a consistent, albeit small, bias effect for that location. However, methods exist for the estimation of
Δ*K* either by correlation with topographic height over estimation windows corresponding to the full-scene or tens of kilometres
(Lin et al., 2010; Bekaert et al., 2015) or from climate reanalysis data (Cong et al., 2012). Therefore, interference from
Φ_{sa} is not considered a limiting factor for the proposed method since the effect can be largely mitigated.

### 5.3.2 Surface displacement

Surface displacements may result from freezing that causes upward heave and thawing that causes downward subsidence. Surface displacements associated
with freezing or thawing are predominantly in the direction of the local surface normal, and therefore spatial variation in the InSAR phase signal
induced by such displacement (Φ_{he}) is approximated by

where *d*_{h} is the displacement magnitude away from the reference surface, with positive and negative displacements, respectively,
corresponding to freezing (heave) and thawing (subsidence).

In equilibrium climate conditions, seasonal active-layer displacement tends to follow a pattern of heave in early to mid-winter related to freeze-back to the top of permafrost, which is followed by static conditions until the spring. Then the active layer thaws gradually throughout the summer. Therefore, the dry-snow season can be divided into an initial period of positively valued displacement followed by a static period. Longer-term positively or negatively valued displacement may be superimposed onto this seasonal pattern due to changes in permafrost conditions. Annual displacement amplitudes for seasonally frozen terrain typically occur in the range of 20 to 120 mm (Gruber, 2020).

Areas with heavily organic soil (e.g., peat) can have significantly higher seasonal displacements. In the vicinity of the study area, we have
observed that these occur in localized areas without significant topographic variation, i.e., areas already not well suited for the method because of
limited *ξ* diversity. It may be possible to exclude these areas by masking based on local *ξ* diversity or analysis of summer interferograms (we
have observed that these areas tend to have low temporal coherence due to their significant displacement phase).

To assess the significance of heave as an interfering factor we consider a reference case in which *d*_{h} = 10 mm over a
24 d interval. Figure 17c and d show the expected
${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ error for this reference case. The error is always positively valued and corresponds to an approximately
+9 mm ΔSWE bias over the 24 d estimation interval which is very substantial, considering that 24 d
ΔSWE values are typically < 20 mm for the study area (see Fig. 12). Hence, there is potential for displacement
processes, if present, to significantly interfere with the proposed estimation method. Regarding whether we can assume that *d*_{h} is constant
over the estimation window, surface conditions affecting heave can certainly vary at scales shorter than the 500 m window size, e.g., as
observed by Liu et al. (2012), and this spatial variability will contribute to estimation error to some degree, although we have not
attempted to quantify this effect.

Regarding the results presented in Sect. 4.3, the ΔSWE values are expected to yield a positive bias until freeze-back of the active layer is completed. Considering the reference case presented here in which 10 mm of heave results in a mean ΔSWE bias of 9 mm, the October–December bias can be explained by $(\mathrm{8.6}\phantom{\rule{0.125em}{0ex}}\mathrm{mm}/\mathrm{0.9})(\mathrm{92}\phantom{\rule{0.125em}{0ex}}\mathrm{d}/\mathrm{24}\phantom{\rule{0.125em}{0ex}}\mathrm{d})$ = 37 mm of mean upward displacement amplitude, which is a plausible value for the area during active-layer freeze-back (Gruber, 2020).

The +1.7 mm bias for the snow-free periods presented in Sect. 4.3 is unexpected since the snow-free period is expected to correspond to subsidence and therefore a negative bias with magnitude similar to but less than October–December because the subsidence is expected to occur over a longer time interval than the re-freeze. This disagreement suggests other factors that are difficult to separate from one another; perhaps the thaw is directionally asymmetric with the heave (i.e., contribution from solifluction) as described in Sect. 5.3.3.

One mitigation approach is to limit use of the estimator spatially and temporally to correspond to known areas and periods of low displacement activity. For example, temporal analysis of summer interferograms can be used to identify areas with significant subsidence, and these can be masked out during ΔSWE estimation under the assumption that thaw-period subsidence is spatially correlated with freeze-period heave. Of course, this is not possible for the case of widespread seasonal surface displacement as is common in periglacial regions. Such areas are relatively common (e.g., Obu et al., 2019, report that 22 % of the exposed land area of the Northern Hemisphere potentially contains permafrost). It may also be possible to estimate and correct for the early winter heave component by estimating the thaw subsidence from summer-season interferograms and then applying a model for the cyclical deformation.

### 5.3.3 Solifluction

We consider slow-moving, gravity-driven surface soil flux, i.e., solifluction, as a potential Φ_{e} component since it is correlated
with topographic slope. We do not expect solifluction to be a significant bias source during the dry-snow season but include it because it may affect
snow-free interferograms analyzed for the purpose of method validation.

Solifluction is driven by cyclical surface displacement (diurnal or seasonal) coupled with sufficient topographic slope magnitude and surface
conditions (Matsuoka, 2001), and it is common in periglacial regions. The flow direction is principally in the local downslope
direction, and therefore the resulting InSAR phase (Φ_{sf}) is approximated by

where *d*_{sf} is the downslope surface displacement magnitude due to solifluction, and ** l** and

**are the unit magnitude SAR line-of-sight and downslope direction vectors. Spatial variation in**

*d**d*

_{sf}can be modelled as a function of

*d*

_{h}and

*α*only. An upper bound for

*d*

_{sf}suggested by Matsuoka (2001) on an annualized basis is

which attributes downslope soil motion to the heave process minus any retrograde motion due to soil cohesion. In order to assess the significance of
solifluction as an interfering factor we consider a reference case in which *d*_{sf} = 10 mm ⋅tan *α*. Figure 17e
and f show the expected ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ error for this reference case. The
error is always positively valued and corresponds to a 0.6 mm bias per 10 mm heave on an annualize basis, which is significantly less
than that due directly to heave discussed in Sect. 5.3.2.

### 5.3.4 Soil moisture

Near-surface soil moisture conditions are a known contributor to the InSAR phase signal (Nolan et al., 2003), and these are likely modulated by
topographic slope both in terms of gravity-driven effects on near-surface hydrology and the influence of slope aspect on solar heating. As such, there
may be some correlation with *ξ* leading to estimation bias. As with solifluction, we do not expect soil moisture to significantly affect dry-snow
interferograms but consider its effect on snow-free interferograms used as method validation controls.

Although the correlation of soil moisture phase with topographic slope is difficult to model, it is possible to establish a bound on its impact on
${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ by considering that the soil moisture phase magnitude, Φ_{sm}, is limited to some maximum value,
Φ_{sm_max}, according to sensor wavelength and soil type (De Zan et al., 2014). Considering a case where $\mathit{\xi}\in [{\mathit{\xi}}_{min},{\mathit{\xi}}_{max}]$, the spatial distribution of Φ_{sm} producing maximum possible ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ bias is

which corresponds to perfect correlation between Φ_{sm} and *ξ*. According to Eq. (12) this yields a ΔSWE bias
limit of ${\mathrm{\Phi}}_{\text{sm\_max}}/({\mathit{\xi}}_{max}-{\mathit{\xi}}_{min})$. Considering our study case, *ξ*∈ [0.22 rad mm^{−1},
0.28 rad mm^{−1} ] (see Fig. 6a), and assuming Φ_{sm_max} = 90^{∘}, which is reasonable considering the model
of De Zan et al. (2014) and the simulation results from Rabus et al. (2010), this corresponds to a bias magnitude limit of 1.3 mm
which is relatively small especially when considering that any actual spatial correlation is likely much less than this ideal upper limit.

## 5.4 DEM error

Errors in the DEM used for analysis contribute to the ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ error due to both the residual topographic phase, as a
component of Φ_{e}, and error in $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$ used by the estimator.

The topographic phase error (Φ_{topo}) scales linearly with both perpendicular orbit baseline, *B*_{⊥}, and DEM error,
Δ*h*. Considering the Φ_{e} term in Eq. (13) and the fact that Δ*h* varies spatially, the resulting
ΔSWE bias scales with $E\left({B}_{\perp}\right)E(\langle \stackrel{\mathrm{\u0303}}{\mathit{\xi}}\mathrm{\Delta}h\rangle )$ with respect to all interferograms. We assume that
*E*(*B*_{⊥}) is zero for a sensor with a maintained orbit, and therefore this effect is likely unbiased with respect to a temporal set of
ΔSWE maps.

Next, we consider the effect of errors in $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$ caused by DEM errors. The estimator defined in Eq. (12) can be better expressed using the DEM-derived ${\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}_{D}=\stackrel{\mathrm{\u0303}}{\mathit{\xi}}+{\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}_{\mathrm{\Delta}h}$ instead of true $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$, in which ${\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}_{\mathrm{\Delta}h}$ is the error in ${\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}_{D}$. As a first-order analysis of this error, we neglect the second-order terms involving both ${\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}_{\mathrm{\Delta}h}$ and the bias terms in Eq. (12). Therefore,

If we assume that ${\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}_{\mathrm{\Delta}h}$ is uncorrelated with $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$,

which corresponds to a multiplicative bias towards lower amplitudes that is more pronounced in areas with low *ξ* variance. We have not attempted to
quantify $\langle {\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}_{\mathrm{\Delta}h}^{\mathrm{2}}\rangle $ for our TanDEM-X DEM-derived *ξ* map, but we made efforts to reduce it by smoothing the DEM
(see Sect. 3.5), noting that this approach has its limits since over-smoothing will tend to increase $\langle {\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}_{\mathrm{\Delta}h}^{\mathrm{2}}\rangle $. Note that smoothing the DEM has the benefit of reducing potential bias introduced by short-scale horizontal
heterogeneity in the snowpack (e.g., those noted by Sturm and Benson, 2004). Instead, the effect of these variations will be limited to an
increase in the uncorrelated phase components discussed in Sect. 3.6.

Regarding requirements for DEM generation frequency to accommodate landform changes, at the scale applied in our study (DEM smoothed to 90 m resolution), DEM generation frequency can likely be > 10 years. However, if applying the method using more finely scaled data, a more recent DEM may be required since fine-scale landform changes occur more frequently.

## 5.5 Spatial resolution and geographic suitability

Regarding the choice of estimator window size, larger windows improve precision by mitigating the effect of spatially random noise. However, as discussed in Sect. 5.3, increasing window size does not necessarily mitigate deterministic bias contributions. Therefore, in terms of the trade-off between estimator error and spatial resolution, there is an upper limit for which increasing window size will significantly reduce the expected error. For our study dataset we empirically determined 500 m to be best choice for the window size since larger windows did not yield significant reductions in the observed spatial variance of ΔSWE estimates.

As discussed in Sect. 3.6, the precision of the SlopeVar estimator depends on both the coherence level and the local
*ξ* variance. These factors are spatially variable, and therefore in order to generalize the achievable spatial resolution to other locations, one
needs to consider both the degree of topographic slope variation and conditions affecting coherence (e.g., vegetation). Within a single dataset, this
heteroskedasticity of the estimator could be addressed by applying an adaptively sized estimation window during the estimation to achieve
approximately uniform variance but variable spatial resolution, in other words, growing the window size to compensate in areas of low
*ξ* diversity. Another approach would be to simply mask out areas expected to have high variance based on coherence or *ξ* diversity.

These considerations affect the geographic suitability of the method; it requires dry-snow conditions and at least moderate topographic variation. Such conditions are present at least part of the time in many geographic regions. The method is less suitable in areas with frequent wet-snow conditions such as coastal areas and in regions that are mostly flat, such as prairie.

## 5.6 Correcting for the dry-snow phase component

The main focus of this article is the estimation of ΔSWE. However, for repeat-pass InSAR applications other than ΔSWE
estimation (e.g., surface deformation estimation), Φ_{s} can be considered an interfering factor. In this case, knowing the spatial
ΔSWE distribution for a given interferogram allows for removal of Φ_{s} according to Eq. (5). Furthermore,
ΔSWE can, in principle, be estimated at scales finer than the SlopeVar estimation window size by (1) first using SlopeVar to generate a
coarse-scale estimate, (2) removing Φ_{s} corresponding to this initial estimate, and then (3) attributing the remaining phase to
additional (high-spatial-frequency) ΔSWE which can then be estimated per resolution cell according to Eq. (5). We performed
a preliminary investigation of this by scaling the ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ maps by the *ξ* map to generate spatial maps of estimated
Φ_{s} and then demodulating the original interferograms. The result was, in general, a significant increase in observed spatial phase
gradients and short-scale phase artifacts likely because of short-scale errors in the ${\widehat{\mathrm{\Delta}S}}_{\mathrm{w}}$ maps. This approach
therefore requires further investigation because, in order to maintain spatial phase smoothness required for subsequent phase unwrapping, the
ΔSWE maps likely require imposition of smoothness constraints.

We have introduced a novel spatial analysis method, “SlopeVar”, that estimates dry SWE change from wrapped repeat-pass interferograms by spatially correlating the InSAR phase to a DEM-derived, slope-dependent dry-snow phase sensitivity map over a suitably sized estimation window. The method does not require phase unwrapping or any spatial reference with known SWE change and therefore addresses some of the key challenges posed by InSAR-based SWE change estimation.

We tested our method using a 9-year stack of RADARSAT-2 Spotlight-A mode images over an area surrounding the town of Inuvik, Northwest Territories, Canada, and corresponding TanDEM-X DEM for the area. We generated SlopeVar ΔSWE maps from the dataset for both snow-covered and snow-free periods and compared these with a number of in situ SWE transect measurements made over two winters, and the results showed good agreement with an absolute SWE RMSE of 15 mm. We compared our spatially averaged results with the ERA5 ΔSWE time history, and the results showed a moderately strong correlation and 7.9 mm RMSE (per 24 d period) over the October–March period. Seasonal analysis revealed a substantial early winter (October–December) bias (8.5 mm per 24 d period) that is largely absent in late winter (January–March), which we suspect is due to seasonal active-layer heave.

We empirically determined that a 500 m estimation window size provides the best combination of error mitigation and spatial resolution for our dataset. However, the achievable spatial resolution depends both on terrain factors affecting the coherence and on the degree of slope variation, both of which are location-dependent.

The method is subject to bias from factors that spatially correlate with the dry-snow phase sensitivity map. We investigated the role of spatial SWE change inhomogeneity as a possible bias factor and determined that vegetation homogeneity should be controlled for during estimation, e.g., by limiting the estimator to homogenous areas. We investigated the role of static atmospheric delay and surface deformation as potential bias factors and also investigated solifluction and soil moisture phase, two bias sources expected to only affect snow-free interferograms used as zero ΔSWE controls for method validation. We found that static atmospheric delay, soil moisture, and solifluction are likely not significant sources of bias and that surface deformation due to heave, if present, may lead to substantial estimation bias. We also considered the effect of DEM error and found that DEM errors lead to damping multiplicative ΔSWE bias, and therefore DEM filtering, pre-processing, or both should be performed to minimize this effect.

We compared the SlopeVar and Delta-K methods in terms of their relative achievable precision and found that, for our test area, SlopeVar is typically about 4× more precise than Delta-K for a given estimation window size. Regarding the benefits of the method compared to dual-frequency volume-scattering-based estimation, our method is less suitable in areas with relatively low topographic variation, although it may be possible to exploit even sparse areas of variation such as along waterways to sparsely map SWE change in otherwise flat areas. In other areas, the constraints are similar in that they both require dry-snow conditions, and there is opportunity for cross-validation in these areas. Our method has the advantage of only requiring a single band and is therefore feasible with existing spaceborne systems, especially L-band. In comparison to passive microwave-based estimation, we expect the SlopeVar method to perform better in high-relief areas which potentially offers complementary coverage in alpine regions that have previously been difficult areas to evaluate.

The code used for this paper has been deposited on the Zenodo repository and is publicly accessible via https://doi.org/10.5281/zenodo.6480208 (Eppler, 2022). Third-party software used consists of the GAMMA software tools (GAMMA Remote Sensing AG).

The SAR imagery and DEM used are uniquely identified by the description in the paper. They have been obtained under non-transferable user licenses and cannot be redistributed. The in situ dataset is included in the same Zenodo repository as the code (https://doi.org/10.5281/zenodo.6480208, Eppler, 2022).

The supplement related to this article is available online at: https://doi.org/10.5194/tc-16-1497-2022-supplement.

JE and BR conceptualized the proposed method. JE performed analysis, produced figures, and wrote the original draft. JE and PM conducted fieldwork to collect validation data. All authors contributed to the manuscript and editing.

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.

Glen Liston provided the SnowModel software used to generate simulated time-evolving SWE distributions. Tim Ensom provided snow survey data. Francesco De Zan provided useful discussions regarding spatial distribution of the dry-snow phase and the possible role of solifluction as an interfering factor. DLR provided the TanDEM-X DEM under science proposal DEM_CALVAL2850. MDA provided the RADARSAT-2 imagery used for the study. The Government of Northwest Territories – Department of Environment and Natural Resources provided metrological station data and SWE survey data. This work was funded in part by the NSERC MDA CSA Industrial Research Chair in SAR Technologies, Techniques, and Applications at SFU. This paper is a contribution to the Geological Survey of Canada's Climate Change Geoscience program and is NRCan contribution number 20210443.

This research has been supported by the Natural Sciences and Engineering Research Council of Canada (NSERC MDA CSA Industrial Research Chair in SAR Technologies, Techniques, and Applications).

This paper was edited by Alexandre Langlois and reviewed by Silvan Leinss and two anonymous referees.

Anderton, S. P., White, S. M., and Alvera, B.: Evaluation of spatial variability in snow water equivalent for a high mountain catchment, Hydrol. Process., 18, 435–453, https://doi.org/10.1002/hyp.1319, 2004. a, b

Bamler, R. and Eineder, M.: Accuracy of Differential Shift Estimation by Correlation and Split-Bandwidth Interferometry for Wideband and Delta-k SAR Systems, IEEE Geosci. Remote S., 2, 151–155, https://doi.org/10.1109/LGRS.2004.843203, 2005. a

Barnett, T. P., Adam, J. C., and Lettenmaier, D. P.: Potential impacts of a warming climate on water availability in snow-dominated regions, Nature, 438, 303–309, https://doi.org/10.1038/nature04141, 2005. a

Bekaert, D. P. S., Hooper, A., and Wright, T. J.: A spatially variable power law tropospheric correction technique for InSAR data, J. Geophys. Res.-Sol. Ea., 120, 1345–1356, https://doi.org/10.1002/2014JB011558, 2015. a

Burn, C. R. and Kokelj, S. V.: The environment and permafrost of the Mackenzie Delta area, Permafrost Periglac., 20, 83–105, https://doi.org/10.1002/ppp.655, 2009. a

Burn, C. R., Mackay, J. R., and Kokelj, S. V.: The thermal regime of permafrost and its susceptibility to degradation in upland terrain near Inuvik, N.W.T., Permafrost Periglac., 20, 221–227, https://doi.org/10.1002/ppp.649, 2009. a, b

Cong, X., Balss, U., Eineder, M., and Fritz, T.: Imaging Geodesy–Centimeter-Level Ranging Accuracy With TerraSAR-X: An Update, IEEE Geosci. Remote S., 9, 948–952, https://doi.org/10.1109/LGRS.2012.2187042, 2012. a, b

De Zan, F., Parizzi, A., Prats-Iraola, P., and Lopez-Dekker, P.: A SAR Interferometric Model for Soil Moisture, IEEE T. Geosci. Remote, 52, 418–425, https://doi.org/10.1109/TGRS.2013.2241069, 2014. a, b

De Zan, F., Prats-Iraola, P., and Rodriguez-Cassola, M.: On the Dependence of Delta-k Efficiency on Multilooking, IEEE Geosci. Remote S., 12, 1745–1749, https://doi.org/10.1109/LGRS.2015.2424272, 2015. a

Deeb, E. J., Forster, R. R., and Kane, D. L.: Monitoring snowpack evolution using interferometric synthetic aperture radar on the North Slope of Alaska, USA, Int. J. Remote Sens., 32, 3985–4003, https://doi.org/10.1080/01431161003801351, 2011. a

Deems, J. S., Painter, T. H., and Finnegan, D. C.: Lidar measurement of snow depth: a review, J. Glaciol., 59, 467–479, https://doi.org/10.3189/2013JoG12J154, 2013. a

Derksen, C., Lemmetyinen, J., King, J., Belair, S., Garnaud, C., Lapointe, M., Crevier, Y., Burbidge, G., and Siqueira, P.: A Dual-Frequency Ku-Band Radar Mission Concept for Seasonal Snow, in: IGARSS 2019 – 2019 IEEE International Geoscience and Remote Sensing Symposium, IEEE, Yokohama, Japan, 28 July–2 August 2019, https://doi.org/10.1109/IGARSS.2019.8898030, pp. 5742–5744, 2019. a

ECCC: Environment and Climate Change Canada - Canadian Climate Normals, https://climate.weather.gc.ca/climate_normals/ (last access: 25 March 2021), 2021a. a

ECCC: Environment and Climate Change Canada – Historical Data – Climate, https://climate.weather.gc.ca/historical_data/search_historic_data_e.html (last access: 9 July 2021), 2021b. a

Engen, G., Guneriussen, T., and Overrein, O.: Delta-K Interferometric SAR Technique for Snow Water Equivalent (SWE) Retrieval, IEEE Geosci. Remote S., 1, 57–61, https://doi.org/10.1109/LGRS.2003.822880, 2004. a, b

Eppler, J.: Code and data for: Snow Water Equivalent Change Mapping from Slope Correlated InSAR Phase Variations. In The Cryosphere Discussions (1.0.1), Zenodo [code and data set], https://doi.org/10.5281/zenodo.6480208, 2022, a, b, c, d

Eppler, J. and Rabus, B. T.: The Effects of Dry Snow on the SAR Impulse Response and Feasibility for Single Channel Snow Water Equivalent Estimation, IEEE T. Geosci. Remote, 60, 1–23, https://doi.org/10.1109/TGRS.2021.3089131, 2021. a

Gruber, S.: Ground subsidence and heave over permafrost: hourly time series reveal interannual, seasonal and shorter-term movement caused by freezing, thawing and water movement, The Cryosphere, 14, 1437–1447, https://doi.org/10.5194/tc-14-1437-2020, 2020. a, b

Grünewald, T., Schirmer, M., Mott, R., and Lehning, M.: Spatial and temporal variability of snow depth and ablation rates in a small mountain catchment, The Cryosphere, 4, 215–225, https://doi.org/10.5194/tc-4-215-2010, 2010. a

Guarnieri, A. M. and Tebaldini, S.: On the Exploitation of Target Statistics for SAR Interferometry Applications, IEEE T. Geosci. Remote, 46, 3436–3443, https://doi.org/10.1109/TGRS.2008.2001756, 2008. a

Guneriussen, T., Hogda, K., Johnsen, H., and Lauknes, I.: InSAR for estimation of changes in snow water equivalent of dry snow, IEEE T. Geosci. Remote, 39, 2101–2108, https://doi.org/10.1109/36.957273, 2001. a, b, c

Hagberg, J. O., Ulander, L. M., and Askne, J.: Repeat-pass SAR interferometry over forested terrain, IEEE T. Geosci. Remote, 33, 331–340, https://doi.org/10.1109/TGRS.1995.8746014, 1995. a

Hanssen, R. F.: Radar interferometry: data interpretation and error analysis, vol. 2, Kluwer Academic Publishers, 2001. a, b

Heginbottom, J.: Canada-permafrost, National Atlas of Canada, publisher: National Atlas Information Service, Natural Resources Canada, 1995. a

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J.: The ERA5 global reanalysis, Q. J. Roy. Meteor. Soc., 146, 1999–2049, https://doi.org/10.1002/qj.3803, 2020. a

Jost, G., Weiler, M., Gluns, D. R., and Alila, Y.: The influence of forest and topography on snow accumulation and melt at the watershed-scale, J. Hydrol., 347, 101–115, https://doi.org/10.1016/j.jhydrol.2007.09.006, 2007. a

Just, D. and Bamler, R.: Phase statistics of interferograms with applications to synthetic aperture radar, Appl. Optics, 33, 4361, https://doi.org/10.1364/AO.33.004361, 1994. a, b

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 Sensing Environ., 215, 242–254, https://doi.org/10.1016/j.rse.2018.05.028, 2018. a

Kokelj, S. V., Palmer, M. J., Lantz, T. C., and Burn, C. R.: Ground Temperatures and Permafrost Warming from Forest to Tundra, Tuktoyaktuk Coastlands and Anderson Plain, NWT, Canada: Permafrost temperatures across tree line northwestern Canada, Permafrost Periglac., 28, 543–551, https://doi.org/10.1002/ppp.1934, 2017. a

Kunzi, K. F., Patil, S., and Rott, H.: Snow-Cover Parameters Retrieved from Nimbus-7 Scanning Multichannel Microwave Radiometer (SMMR) Data, IEEE T. Geosci. Remote, GE-20, 452–467, https://doi.org/10.1109/TGRS.1982.350411, 1982. a

Latifovic, R., Pouliot, D., and Olthof, I.: Circa 2010 Land Cover of Canada: Local Optimization Methodology and Product Development, Remote Sens.-Basel, 9, 1098, https://doi.org/10.3390/rs9111098, 2017. a

Leinss, S., Wiesmann, A., Lemmetyinen, J., and Hajnsek, I.: Snow Water Equivalent of Dry Snow Measured by Differential Interferometry, IEEE J. Sel. Top. Appl., 8, 3773–3790, https://doi.org/10.1109/JSTARS.2015.2432031, 2015. a, b, c, d

Leinss, S., Löwe, H., Proksch, M., Lemmetyinen, J., Wiesmann, A., and Hajnsek, I.: Anisotropy of seasonal snow measured by polarimetric phase differences in radar time series, The Cryosphere, 10, 1771–1797, https://doi.org/10.5194/tc-10-1771-2016, 2016. a

Leinss, S., Löwe, H., Proksch, M., and Kontu, A.: Modeling the evolution of the structural anisotropy of snow, The Cryosphere, 14, 51–75, https://doi.org/10.5194/tc-14-51-2020, 2020. a

Lemmetyinen, J., Derksen, C., Rott, H., Macelloni, G., King, J., Schneebeli, M., Wiesmann, A., Leppänen, L., Kontu, A., and Pulliainen, J.: Retrieval of Effective Correlation Length and Snow Water Equivalent from Radar and Passive Microwave Measurements, Remote Sens.-Basel, 10, 170, https://doi.org/10.3390/rs10020170, 2018. a

Li, L. and Pomeroy, J. W.: Probability of occurrence of blowing snow, J. Geophys. Res.-Atmos., 102, 21955–21964, https://doi.org/10.1029/97JD01522, 1997. a, b

Lin, Y.-n. N., Simons, M., Hetland, E. A., Muse, P., and DiCaprio, C.: A multiscale approach to estimating topographically correlated propagation delays in radar interferograms, Geochem. Geophy. Geosy., 11, Q09002, https://doi.org/10.1029/2010GC003228, 2010. a, b, c, d

Liston, G. E. and Elder, K.: A Distributed Snow-Evolution Modeling System (SnowModel), J. Hydrometeorol., 7, 1259–1276, https://doi.org/10.1175/JHM548.1, 2006. a

Liu, L., Schaefer, K., Zhang, T., and Wahr, J.: Estimating 1992–2000 average active layer thickness on the Alaskan North Slope from remotely sensed surface subsidence, J. Geophys. Res.-Earth, 117, n/a–n/a, https://doi.org/10.1029/2011JF002041, 2012. a

Mackay, J. R. and MacKay, D.: Snow Cover and Ground Temperatures, Garry Island, N.W.T., Arctic, 27, 287–296, https://doi.org/10.14430/arctic2885, 1974. a, b

Mätzler, C.: Applications of the interaction of microwaves with the natural snow cover, Remote Sensing Reviews, 2, 259–387, https://doi.org/10.1080/02757258709532086, 1987. a

Matsuoka, N.: Solifluction rates, processes and landforms: a global review, Earth-Sci. Rev., 55, 107–134, https://doi.org/10.1016/S0012-8252(01)00057-5, 2001. a, b

Morse, P., Burn, C., and Kokelj, S.: Influence of snow on near-surface ground temperatures in upland and alluvial environments of the outer Mackenzie Delta, Northwest Territories This article is one of a series of papers published in this CJES Special Issue on the theme of Fundamental and applied research on permafrost in Canada, Can. J. Earth Sci., 49, 895–913, https://doi.org/10.1139/e2012-012, 2012. a, b

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, https://doi.org/10.5194/tc-14-1579-2020, 2020. a

Nguyen, T.-N., Burn, C. R., King, D. J., and Smith, S. L.: Estimating the extent of near-surface permafrost using remote sensing, Mackenzie Delta, Northwest Territories, Permafrost Periglac., 20, 141–153, https://doi.org/10.1002/ppp.637, 2009. a

Nolan, M., Fatland, D., and Hinzman, L.: Dinsar measurement of soil moisture, IEEE T. Geosci. Remote, 41, 2802–2813, https://doi.org/10.1109/TGRS.2003.817211, 2003. a

Nolan, M., Larsen, C., and Sturm, M.: Mapping snow depth from manned aircraft on landscape scales at centimeter resolution using structure-from-motion photogrammetry, The Cryosphere, 9, 1445–1463, https://doi.org/10.5194/tc-9-1445-2015, 2015. a

Obu, J., Westermann, S., Bartsch, A., Berdnikov, N., Christiansen, H. H., Dashtseren, A., Delaloye, R., Elberling, B., Etzelmüller, B., Kholodov, A., Khomutov, A., Kääb, A., Leibman, M. O., Lewkowicz, A. G., Panda, S. K., Romanovsky, V., Way, R. G., Westergaard-Nielsen, A., Wu, T., Yamkhin, J., and Zou, D.:
Northern Hemisphere permafrost map based on TTOP modelling for 2000–2016 at 1 km^{2} scale, Earth-Sci. Rev., 193, 299–316, https://doi.org/10.1016/j.earscirev.2019.04.023, 2019. a

Oveisgharan, S., Esteban-Fernandez, D., Waliser, D., Friedl, R., Nghiem, S., and Zeng, X.: Evaluating the Preconditions of Two Remote Sensing SWE Retrieval Algorithms over the US, Remote Sens.-Basel, 12, 2021, https://doi.org/10.3390/rs12122021, 2020. a

O'Neill, H. B., Smith, S., and Duchesne, C.: Long-term permafrost degradation and thermokarst subsidence in the Mackenzie Delta area indicated by thaw tube measurements, in: Cold Regions Engineering 2019, 19–21 August 2019, Quebec City, Canada, pp. 643–651, American Society of Civil Engineers Reston, VA, 2019. a

Pachauri, R. K., Mayer, L., and on Climate Change, I. P. (Eds.): Climate change 2014: synthesis report, Intergovernmental Panel on Climate Change, Geneva, Switzerland, 2015. a

Palmer, M.: Vegetation, snow cover, and air and near-surface ground temperature across treeline in the uplands east of the Mackenzie Delta, Northwest Territories, Master of Science, Carleton University, Ottawa, Ontario, https://doi.org/10.22215/etd/2007-08207, 2007. a

Palmer, M., Burn, C., and Kokelj, S.:
Factors influencing permafrost temperatures across tree line in the uplands east of the Mackenzie Delta, 2004–2010 ^{1} This article is one of a series of papers published in this CJES Special Issue on the theme of Fundamental and applied research on permafrost in Canada. 2 Polar Continental Shelf Contribution 03611, Can. J. Earth Sci., 49, 877–894, https://doi.org/10.1139/e2012-002, 2012. a

Rabus, B., Wehn, H., and Nolan, M.: The Importance of Soil Moisture and Soil Structure for InSAR Phase and Backscatter, as Determined by FDTD Modeling, IEEE T. Geosci. Remote, 48, 2421–2429, https://doi.org/10.1109/TGRS.2009.2039353, 2010. a

Rampton, V. N.: Surficial Geology, Tuktoyaktuk Coastlands, District of Mackenzie, Northwest Territories, “A” Series Map 1647A, Scale $\mathrm{1}:\mathrm{500},\mathrm{000}$, Surveys and Mapping Branch, Geological Survey of Canada, 1 sheet, https://doi.org/10.4095/125160, 1987. a

Rott, H., Nagler, T., and Scheiber, R.: Snow mass retrieval by means of SAR interferometry, in: 3rd FRINGE Workshop, 1–5 December 2003, Frascati, Italy, European Space Agency, Earth Observation, pp. 1–6, Citeseer, 2003. a

Rott, H., Cline, D. W., Duguay, C., Essery, R., Etchevers, P., Hajnsek, I., Kern, M., Macelloni, G., Malnes, E., Pulliainen, J., and Yueh, S. H.: CoReH2O, a dual frequency radar mission for snow and ice observations, in: 2012 IEEE International Geoscience and Remote Sensing Symposium, IEEE, Munich, Germany, 22–27 July 2012, https://doi.org/10.1109/IGARSS.2012.6352348, pp. 5550–5553, 2012. a

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, https://doi.org/10.5194/tc-13-3045-2019, 2019. a

Saberi, N., Kelly, R., Flemming, M., and Li, Q.: Review of snow water equivalent retrieval methods using spaceborne passive microwave radiometry, Int. J. Remote Sens., 41, 996–1018, https://doi.org/10.1080/01431161.2019.1654144, 2020. a

Smyth, E. J., Raleigh, M. S., and Small, E. E.: Improving SWE Estimation With Data Assimilation: The Influence of Snow Depth Observation Timing and Uncertainty, Water Resour. Res., 56, https://doi.org/10.1029/2019WR026853, 2020. a

Sturm, M. and Benson, C.: Scales of spatial heterogeneity for perennial and seasonal snow layers, Ann. Glaciol., 38, 253–260, https://doi.org/10.3189/172756404781815112, 2004. a

Sturm, M. and Wagner, A. M.: Using repeated patterns in snow distribution modeling: An Arctic example, Water Resour. Res., 46, https://doi.org/10.1029/2010WR009434, 2010. a

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 Sensing Environ., 115, 3517–3529, https://doi.org/10.1016/j.rse.2011.08.014, 2011. a, b

Walker, D. A.: Hierarchical subdivision of Arctic tundra based on vegetation response to climate, parent material and topography, Glob. Change Biol., 6, 19–34, https://doi.org/10.1046/j.1365-2486.2000.06010.x, 2000. a

Wessel, B., Huber, M., Wohlfart, C., Marschalk, U., Kosmann, D., and Roth, A.: Accuracy assessment of the global TanDEM-X Digital Elevation Model with GPS data, ISPRS J. Photogramm., 139, 171–182, https://doi.org/10.1016/j.isprsjprs.2018.02.017, 2018. a

Zemp, M., Hoelzle, M., and Haeberli, W.: Six decades of glacier mass-balance observations: a review of the worldwide monitoring network, Ann. Glaciol., 50, 101–111, https://doi.org/10.3189/172756409787769591, 2009. a