Snow Water Equivalent Change Mapping from Slope Correlated InSAR Phase Variations

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. Repeat-pass synthetic aperture radar interfer5 ometry (InSAR) as an alternate technique 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 repeatpass InSAR based dry-snow SWE estimation that exploits the sensitivity of the dry-snow refraction-induced InSAR phase to 10 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 spotlightmode dataset near Inuvik, 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 and in areas with high SWE (>150 mm), thus providing 15 complementary coverage to other passiveand active-microwave based SWE estimation methods. Further, our method has the advantage of requiring only a single wavelength band and thus can utilize existing spaceborne synthetic aperture radar systems. In application, a first order analysis of SWE trends within three drainage basins suggests that differences between basin-level accumulations are a function of major landcover types, and that re-vegetation following a forest-tundra fire that occurred over 50 years ago continues to affect the spatial distribution of SWE accumulation in the study area. 20 Copyright statement. 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 2021 1 https://doi.org/10.5194/tc-2021-359 Preprint. Discussion started: 15 December 2021 c © Author(s) 2021. CC BY 4.0 License.


25
Snow melt is one of the greatest sources of water in snow-affected areas, snow accumulation on glaciers is critical to their mass balance and longevity, and snow cover variation is a dominant control on ground thermal regimes in cold regions, 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. Traditionally, SWE variation over space and time is interpolated from sparse point-based measurements, which generally 30 show significant spatial variation due to factors such as topography and vegetation, and temporal bias if the acquisition of each set of points is distributed over a range of dates (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 pointbased 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-35 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, collected intermittently, and with a range of techniques. Snow depth, used to infer SWE when integrated with snow density, has been determined using surface elevation models developed using LiDAR as well as structure-from-motion (optical photogrammetry) systems, either by differential repeat-pass 40 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, widescale 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 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 95 Delta-k method. We then examine 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. Next, we show results from applying our novel SWE retrieval method to the full RADARSAT-2 InSAR dataset. We discuss these results in the context of several in situ SWE transect surveys located within the imaged footprint, as well as a mean SWE temporal history estimated by the ERA5 climate reanalysis model. We then present a first 100 order assessment of the ∆SWE results in the context of landcover type, fire history, and three drainage basins with the study area. Finally, we present our conclusions with a focus on the feasibility of our method for SWE change mapping in the context of existing methods.

Description of 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-105 look complex (SLC) images spanning the period 2012-02-22 to 2021-01-29 at 24-day intervals; nearly continuously with only a few time gaps (Table 1). The study area, situated in Northwest Territories, Canada, at the eastern margin of the Mackenzie Delta, about 120 km south of the Beaufort Sea coast, covers 238 km 2 defined by the SAR imaging footprint and it includes the town of Inuvik. The mean backscatter amplitude image for the dataset is shown in Fig. 1a. Sub-polar or polar shrubland-lichen-moss Sub-polar or polar grasslandlichen-moss Sub-polar or polar barrenlichen-moss Wetland Barren lands Urban Water ALLUVIAL: Floodplains and deltas: silt, clay fine sand, minor gravel, course sand, and organic sediment; sediment in deltas may be more than 20 m thick; includes Aklavik Formation, which constitutes Mackenzie Delta ALLUVIAL: Fans: silt, clay, fine sand, some coarse sand, gravel, generally more than 10 m thick COLLUVIAL: clay, silt, sand, rubble (angular pebble-through boulder sized clasts); generally 0.5 to 3 m thick; materials derived from subaerial weathering; occurs as blankets and veneers on slopes MORAINAL: Hummocky and rolling moraine; generally 4 to 10 m thick MORAINAL : Blanket; variable thickness, but generally 2 to 5 m thick Areas of topography known to have been severely modified by thermokarst; does not include areas of modern thermokarst lakes Nacional para el Conocimiento y Uso de la Biodiversidad, Comisión Nacional Forestal, Instituto Nacional de Estadística y Geografía, and the U.S. Geological Survey, 2015) showing spatial variation of land cover and superimposed with margin of 1968 wildfire (discussed in Sect. 7), (d) surficial geology map of the area from Rampton (1987) © Her Majesty the Queen in Right of Canada, as represented by the Minister of Energy, Mines and Resources Canada, 1987. Positions labelled A-E in panel (a) refer to the snow surveys described in Sect. 6.2.
(DLR). This DEM was derived from source images acquired between 2011-01-13 and 2014-08-08 including snow-season acquisitions so some snow-related bias may affect the estimated elevations. Regarding DEM accuracy, our method described in Sect. 4 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 115 used this DEM for InSAR topographic phase correction and as an input for later snow transport modelling. The DEM shows that 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 m above mean sea level (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 120 to the North American Land Change Monitoring System (NALCMS) 2015 dataset, a 30-metre resolution land classification based on Landsat and RapidEye imagery using the methodology described by Latifovic et al. (2017), there is significant town of Inuvik, segments of the Inuvik-Tuktoyaktuk highway 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 two or more years that is overlain by an active layer that thaws seasonally. The permafrost in the study area is commonly ice rich , and occurs within surficial deposits of several origins 130 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 126 mm water equivalent. For dry-snow conditions, temperatures must remain below freezing; the 135 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 Mackay and MacKay, 1974;Palmer, 2007;Morse et al., 2012). 3 Spatial Variations of Repeat-pass InSAR Dry-snow Phase 140 3.1 Origin of Dry-snow Phase 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 145 contributions will be much less that the primary ground-scattered return. The dry-snow layer has a higher 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 as shown in Fig. 3. Compared to the snow free state, the phase of the SAR signal, Figure 3. Geometry of refracted ray through a dry snow layer over an inclined ground surface compared to the unrefracted (straight) ray trajectory.
150 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 local slope surface normal direction n). Functional dependency on ρ is via the dielectric permittivity , with other parameters in Eq.(1) being the local incidence angle θ and the radar wavelength λ.
In the case of a stratified snowpack, Leinss et al. (2015) have shown that the phase contribution is well approximated by a single layer with mean vertical density, where S is the SWE of the snowpack with depth Z s (measured vertically). Noting that Z s = D s / cos α where α is the local topographic slope angle, Eqs. (1) and (2) can be combined to provide the explicit linear relation between dry-snow phase and SWE, 160 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 density of the removed or added snow layer. As such, ∆S can be either positive 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 one-dimensional flight track 165 (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. The phase contribution due to dry-snow refraction depends on the local incidence angle, which varies along the aperture. The net effect on the focused image phase, 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 170 across-swath variations can be neglected) a horizontally uniform snow layer results in a uniform phase. In contrast, a local spatial region with topographic slope variations will generate dry-snow phase variations corresponding to variations in the local incidence angle and projection of the snow depth onto the slope normal. Figure 4 depicts the topographic modulation of dry-snow 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.

175
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 sensitivity of the dry-snow phase to a uniform SWE layer can be computed as, where we assume that snow density, ρ, is known. A spatially varying change in dry-SWE during a time-interval spanned by two SAR acquisitions will cause a corresponding spatially variable interferometric phase contribution, modulated by the topographic sensitivity,

185
This effect is depicted in Fig. 5, which shows ξ computed over the RS2-SLA spatial footprint and a late snow-season interferogram which shows spatial phase variations that are correlated with ξ. horizontal (i.e. hilltop) and backslope with a fixed incident radar direction and it shows that the local incidence angle increases from foreslope to backslope. (b) depicts the same hill with the addition of a dry-snow layer and shows the corresponding refracted ray geometry. Refraction increases from foreslope to backslope. Also illustrated is a reduction of the slope-normal projected snow depth, Ds with increasing slope angle, α, while vertical snow depth is constant. 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 both due to aliasing and local incoherence of the InSAR phase. Furthermore, although the 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 dry-SWE phase described in Sect. 3.2 195 allows an absolute estimate without requiring phase unwrapping or knowing the phase offset.
Consider a local estimation window, sufficiently large to contain some topographic slope variation. The phase within this window consists of the superposition of the dry snow phase contribution and other phase components (e.g., due to decorrelation, surface displacement, imperfect topographic correction, soil moisture, or tropospheric delay). For the purposes of estimating dry-snow, we treat the other phase components as a single error phase term, Φ e . Due to the difficulty of phase unwrapping and 200 offset estimation we assume that only the locally de-meaned superposition phase 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)Φ 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 Substituting Eqs. (8) and (9) into Eq.(5) and expanding yields Substituting Eq.(10) into Eq.(6) and noting that ∆ S = ξ = 0 gives, Assuming that the first term, ∆S ξ , is the dominant component, we introduce the following correlation-based estimator for ∆S , This estimator correlates the observable de-meaned phaseΦ with the known topographic variation in dry-SWE phase sensi-220 tivityξ. The estimator relies on the assumption that the horizontal SWE change distribution is uniform within the window and that other interferometric phase components are uncorrelated withξ. The bias of ∆S with respect to ∆S is given by where E(·) denotes the expectation. Therefore, ∆S is biased by horizontal SWE change variations,∆S, that are systematically correlated with eitherξ orξ 2 , and by other interferometric phase components systematically correlated withξ. The significance 225 of the different bias terms is examined in Sect. 5.2 and Sect. 5.3.

Wrapped-phase ('SlopeVar') Estimator
The estimator defined in Eq.(12) requires the demeaned unwrapped-phaseΦ within the spatial window of consideration and we have assumed so far that this can be obtained from the wrapped phase 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 230 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, the phase within the window must either be unwrapped or an alternate wrapped-phase version of the estimator is required.
Here we propose such an alternate estimator (which we denote the 'SlopeVar' estimator) that uses periodogram optimization to iteratively estimate ∆S , using the set of wrapped phases 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∆S = 0, this yields an unbiased estimate, i.e., ∆S w = ∆S , otherwise the ξ ∆ S,ξ∆S and ξ∆ S terms in Eq.(15) will bias ∆S w .
240 Figure 6 shows an example of the SlopeVar estimator being applied to a 1 km square (ground range × azimuth) window for the inset area labelled in Fig. 5. The uncorrected phase is spatially correlated with ξ over the window as seen by comparing Figs. 6a and 6b and by examining the 2D-histogram in Fig. 6d. Figure 6e shows the periodogram with distinct peak at ∆SWE of 28 mm. Figure 6c shows the effect of removing the dry-snow phase 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.

Estimator Implementation
Regarding estimator implementation, 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) multilooked 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 250 computation run-time and was seen to have 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, non-sequential (e.g. 48-, 72-day, . . . ) interfer-255 ograms are known to provide additional information that can be recovered by using a phase-linking estimator (e.g., Guarnieri and Tebaldini, 2008). This could be added to the pre-estimation InSAR processing to reduce estimation variance in the future, but was deemed unnecessary for the current application.
We found the estimator to be sensitive to errors in ξ, especially those resulting from DEM errors and interpolation artifacts.
High spatial-frequency errors in the ξ map tend to bias the estimated ∆SWE magnitude towards zero (discussed in Sect. 5.5).

260
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 blunders in the raw DEM. Estimator performance improved significantly after we (1) manually identified and masked local DEM blunders and interpolated over them; (2) smoothed the DEM using a 2D Gaussian filter with standard deviation of 3 DEM pixels (resulting in approximately 90 m spatial resolution) to reduce high-265 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 the real permittivity of dry snow as (ρ) = 1 + 1.5995ρ + 1.861ρ 3 , according to Leinss et al. (2015) and for simplicity, assumed ρ = 0.3 for all estimates (discussed further in Sect. 5.1).
We then implemented the SlopeVar estimator described in Sect. 4.2. For each candidate value of ∆S w , the periodogram magnitude, e j(φ−∆Swξ) can be computed efficiently using a 2D fast Fourier transform-based rectangular smooth filter 270 over the complex-valued interferogram, allowing for 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 2 grid intervals from the edge of the search range were labeled as invalid. Figure 7a shows the resulting ∆SWE

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 InSAR dry-snow phase sensitivity so that absolute ∆SWE can be estimated without the need for phase unwrapping.

280
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 285 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 phase-sensitivity diversity that exists between separate range-frequency sub-bands yielding dry-snow phase sensitivities of ξ ± ∆ξ/2 where ξ is the sensitivity at the central carrier frequency and ∆ξ = B(1 − b)ξ where B ∈ [0, 1] is the normalized range bandwidth of the SAR and b ∈ [0, 1] is the sub-bandwidth fraction used. Assuming multi-290 looked phase-noise with variance σ 2 Φ , the mean sub-band interferometric phase will have variance of σ 2 Φ /bN and therefore the Delta-K difference phase will have variance σ 2 Φ_DK = 2σ 2 Φ /bN . Delta-K estimates ∆SWE according to ∆S DK = Φ DK /∆ξ and therefore the ∆SWE estimation standard deviation is, For the SlopeVar method, the estimation of ∆SWE corresponds to a linear regression over N samples and therefore, assum-295 ingξ is Gaussian distributed under the estimation window with variance ξ 2 , The relative estimation precision ratio between the two methods is therefore given by, where we have replaced ξ in Eq.(16) with the mean under the estimation window, ξ . We use b = 1/3, which, neglecting 300 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 8 shows the distribution of σ ∆S_DK /σ ∆S_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.
305 Figure 8. Histogram of the predicted ∆SWE estimation precision ratio between the Delta-K and SlopeVar methods (σ∆S_DK/σ∆S_SV) for all 500 m square estimation windows in the RS2-SLA study footprint. The vertical bar denotes the point of equal precision (i.e. σ∆S_DK/σ∆S_SV = 1). This shows that the SlopeVar estimator provides substantially better precision for almost all areas within the study footprint.

Snow Density Misspecification
Our proposed method relies on knowing ξ, the linear scaling factor between ∆SWE and the corresponding dry-snow phase which, as shown in Eq. (3) , depends on snow density, ρ. Therefore ρ must either be known a priori or an assumed value must be specified. Here we assess the effect of mis-specifying ρ. Since the method estimates ∆SWE from interferometric phase, the 310 required density corresponds to the snowpack change rather than the actual bulk density of the snowpack, such as is the case for fresh snowfall on a dense late-season snowpack. The density of the added layer will be less than the mean density 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. 9a for true and assumed densities over the interval [0, 0.4] and it indicates that assumed densities greater or less than the true density result in a positive or negative bias, respectively. Figure 9a also shows that the bias magnitude is <5% of the true SWE change for the considered density range. In the absence of any ∆SWE, snowpack evolution may result in a bulk density change (e.g. due to settling) which will 320 cause a snow phase contribution. This phase contribution will be spuriously attributed to a SWE change. Consider the case of a snowpack with total SWE of S that undergoes a density change from initial ρ 1 to final ρ 2 . According to Eq.(5) , this will cause dry-snow phase change, Assuming that the temporal mean density of (ρ 1 + ρ 2 ) /2 is somehow known and used for estimation, the resulting ∆SWE 325 error, S e as a fraction of total SWE will be This fractional error is shown in Fig. 9b for initial and final densities over the interval [0, 0.4]. This shows that a density 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 density range, the bias is 330 <5% of the total SWE.
For simplicity we assumed a fixed snow density of ρ = 0.3 for all ∆SWE estimates. This value is likely too large for the early snow-season in Inuvik and, according to Fig. 9a, likely contributes a small positive estimation bias of 1-3% to the results presented in Sect. 6. This bias could be mostly mitigated by assuming a more appropriate density according to the location and time of year (e.g., Sturm and Wagner, 2010).

Violation of ∆SWE Horizontal Homogeneity Assumption
Temporal changes in SWE can be caused by a number of different processes including deposition of new snow, 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 of ∆SWE over a particular time period 340 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. 4.1, the estimation of ∆SWE within a spatial window through the correlation with the topographic variation of dry-SWE phase sensitivityξ, is subject to bias when there are horizontal ∆SWE variations within the window that are correlated with eitherξ itself or withξ 2 . To investigate the significance of this effect, we simulated the seasonal evolution 345 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.

350
The spatially varying vegetation was modeled 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. 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 TanDEM-X DEM used for InSAR topographic For the bare earth case (Figs. 10a, 10b and 10c), a single outlier interval (2012-01-05 to 2012-01-29) 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ξ for that 24-day 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 375 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 380 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-day simulation intervals. The result is a directional cumulative blowing snow index in units of hours (of blowing snow) and is shown in Fig. 11 for all simulated 24-day 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. Figure 11. Directional cumulative blowing snow hours computed using ERA-5 surface parameters applied to the Li and Pomeroy (1997) probability model for all simulated 24-day intervals. The 20120105_20120129 interval appears as an outlier with respect to the cumulative magnitude. The SAR look direction axis is shown for comparison.
For the bare ground case, the temporal mean ∆SWE error map shows mean errors in the +/-2 mm range with a spatial pattern that is correlated with the topography (Fig. 10a). The distribution of all ∆SWE errors (Fig. 10b), 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 change. The histograms of the per-interval spatially averaged errors (Fig. 10c) show that the majority of 24-day intervals have a mean error near zero, suggesting that averaging estimates 390 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. 10d) 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). We compared the relative contributions of the ξ ξ∆ S and ξ 2∆ S contributions to the ∆SWE error which showed that for all 24-day intervals, and both bare earth and NALCMS 2015 cases the RMS of the ξ ξ∆ S component was much greater than that of the ξ 2∆ S component so that with respect to the influence of∆S, the spatial correlation with respect tõ ξ 2 can be neglected.

Correlated Phase Components
As shown in Eq. (12) , additional phase components, other than dry-snow refractive phase, will contribute to ∆SWE estimation error according to their spatial correlation with the dry-snow phase sensitivity: ξ Φ e / ξ 2 . We distinguish between (1) systematic spatial correlation between the dry-snow phase sensitivity and additional phase components that have a topographic dependence and (2) spurious correlations with phase components that are not systematically dependent on topography due 410 to sampling statistics over the finite estimation window. In this section we consider the former by examining several phase components, their spatial characteristics with respect to topographic slope and their expected contribution to ∆SWE estimation error. Uncorrelated phase components are discussed in Sect. 5.4.

Atmospheric delay
The troposphere contributes substantial phase delay during SAR imaging (Hanssen, 2001). This delay can be decomposed into a 415 horizontal mean component and a horizontally variable component which we respectively denote as the 'static' and 'dynamic' components. Temporal differences in the static component between image acquisitions results in a phase component that is modulated by topography. Neglecting any height dependence of the delay over the height range within the estimation window, the phase can be modelled as a simple linear function (Lin et al., 2010),

420
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ξ. Note that the horizontal variation of ∆K is gradual and therefore ∆K can be assumed constant within the estimation window. Therefore, the spatial variation of Φ sa within the ∆SWE window (assumed to be ≤ 1 km) can be attributed solely to topographic height variation.
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 km (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 as well as ξ prior to ∆SWE estimation, similar to the method employed by Lin et al. (2010) for estimating ∆K from interferograms.
430 Figures 12a and 12d show the expected ∆SWE estimation error for the case of ∆K = 10 -5 (i.e. 10 mm km −1 ). For this case the RMS error is 1.1 mm SWE. However, the distribution is centered around zero so that the mean error is only 0.04 mm SWE.
Therefore, spatial averaging of the estimates over a sufficient area will tend to mitigate the error due to static atmospheric delay. (c) and (f) show the error for reference solifluction case (Vsf = 37 mm y −1 · tanα, corresponding to surface velocity of 10 mm y -1 for a slope of α = 15 • , evaluated for a 24 day interval).
Regarding the temporal distribution of ∆K, given that it is a pairwise differential quantity, it has zero expectation on an annualized time-scale. However, there is some seasonal dependence on tropospheric delay (Cong et al., 2012) and therefore 435 ∆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 kilometers (Lin et al., 2010;Bekaert et al., 2015), or from climate reanalysis data (Cong et al., 2012). Therefore, interference by static atmospheric phase is not considered a limiting factor for the proposed 440 method, since the effect can be largely mitigated.

Surface displacement
Surface displacements may result from freezing that causes upward heave and thawing that causes downward subsidence, either near the surface in soil layers that freeze and thaw seasonally or deeper in the case of permafrost aggradation or degradation in periglacial regions. Surface displacements associated with freezing or thawing are predominantly in the direction of the local 445 surface normal and therefore spatial variation of the InSAR phase signal induced by such displacement 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-450 winter related to freeze-back to the top of permafrost, which is followed by static conditions until the spring when air temperatures rise above 0°C. Then the active layer thaws continually throughout the summer with gradual subsidence until the mean daily surface temperature remain below 0°C again. 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 such as ground ice loss (subsidence) due 455 to climatic warming trends or surface disturbance that cause permafrost degradation, or ground ice gain (heave) due to climatic cooling trends or changing surface conditions that cause permafrost aggradation. Annual displacement amplitudes for seasonally frozen terrain typically occur in the range of 20 to 120 mm (Gruber, 2020), and are normally greater in fine-grained, frost-susceptible soils (e.g., morainal deposits with a high silt content), and lesser in coarse-grained non frost-susceptible soils (e.g., glaciofluvial deposits with a high sand content).

460
To assess the significance of heave as an interfering factor we consider a reference case where d h = 10 mm over a 24-day interval. Figures 12b and 12e show the expected ∆SWE estimation error for this reference case. The error is always positively valued and corresponds to an approximately +9 mm ∆SWE bias over the 24-day estimation interval which is very substantial.
Hence, there is potential for displacement processes, if present, to significantly interfere with the proposed estimation method.
One mitigation approach is to limit use of the estimator spatially and temporally to correspond to known areas and periods 465 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.

470
We consider slow moving, gravity-driven surface soil flux, i.e., solifluction, as a potential interfering phase 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 environmental conditions that include vegetation type, soil moisture content, and soil grain-size distribution (Mat-475 suoka, 2001), and it is common in periglacial regions. The flow direction is principally in the local downslope direction and therefore the resulting InSAR phase is approximated by where V sf is the surface velocity magnitude due to solifluction, ∆t is the InSAR pair temporal baseline, and L and D are the unit magnitude SAR look and downslope direction vectors. Assuming homogenous soil conditions within the estimation 480 window, spatial variation of V sf can be modelled as a function of topographic slope angle only. A commonly applied model is given by Matsuoka (2001), where β sf is a scaling factor, which we assume is constant over the window. In order to assess the significance of solifluction as an interfering factor we consider a reference case where V sf = 10 mm y -1 for a slope of α = 15 • . This corresponds to β sf = 37 485 mm y -1 . Figures 12c and 12f show the expected ∆SWE estimation error for this reference case. The error is always positively valued and corresponds to a 2.2 mm SWE bias for a 24-day estimation interval.
It is noteworthy that both the heave (upward normal-to-surface) and solifluction (downward tangential-to-surface) result in positively valued biases. In the case of symmetric surface displacement (heave followed by equal magnitude subsidence along surface normal -so called 'retrograde' motion) the net ∆SWE bias over the cyclical period (including snow-free thaw-season 490 interferograms) is zero. However in the more general case of heave followed by a combination of retrograde motion and solifluction, e.g. as demonstrated by Harris et al. (2008), the result is net positive ∆SWE bias.

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 495 slope aspect on solar heating. As such, there may be some correlation with the dry-snow phase sensitivity 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 the ∆SWE estimation by considering that the soil moisture phase magnitude, Φ sm is limited to some maximum 500 value, Φ sm_max according to sensor wavelength and soil type (De Zan et al., 2014). Considering a case where ξ ∈ [ξ min , ξ max ], the spatial distribution of Φ sm producing maximum possible ∆SWE estimation bias is which corresponds to perfect correlation between Φ sm and ξ. According to Eq.(12) this yields a ∆SWE bias limit ofΦ sm_max /(ξ max − ξ min ). Considering our study case, ξ ∈ [0.22 rad mm −1 , 0.28 rad mm −1 ] (see Fig. 5a) and assuming Φ sm_max = 90°, which is 505 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 SWE which is relatively small especially when considering that any actual spatial correlation is likely much less than this ideal upper limit.

Uncorrelated Phase Components
Phase components not systematically correlated withξ will contribute to the ∆SWE estimation error due to finite-sampling over 510 the estimation window. These include decorrelation phase noise and the non-ξ correlated parts of the components described previously in Sect. 5.3. 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.
Decorrelation occurs in repeat-pass InSAR because of spatial differences in the imaging geometry, temporal changes in the 515 scattering environment including snow wetting/melt and system effects including thermal noise and SAR processing errors (Zebker and Villasenor, 1992). Furthermore, we expect some dependence of decorrelation on the true ∆SWE magnitude due to variations in the refractive dry-snow phase at both the sub-resolution scale and at the sub-estimation window scale.
The contributions include those due to horizontal snow depth variations caused either by varying ground surface heights as described by Rott et al. (2003) and especially in the early snow-season when short scale undulations 'fill-in', and spatial depth 520 variations due to snow surface structure at scales larger than the radar wavelength.
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 ∆SWE estimation variance is made difficult by the fact that the 525 distribution ofξ 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. This requires knowledge of the space and time varying uncorrelated phase error distribution. We approximate this as the distribution of 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, 530 γ. We estimate γ as the ensemble coherence of the residual phase after removing the dry-snow phase components using the ∆SWE estimate, As such, this approach can be used to estimate a unique ∆SWE estimation variance map for each source interferogram.

540
Errors in the DEM used for analysis contribute to the ∆SWE estimation error due to both residual topographic phase (as a component of total phase error, Φ e ) and error in the DEM-derived phase-sensitivity,ξ, used by the estimator.
The topographic phase error 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 B ⊥ E ξ ∆h 545 for a single interferogram and with E (B ⊥ ) E ξ ∆h 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.
Furthermore, the effect can also be neglected if one assumes that that ∆h is uncorrelated withξ, which may or may not be the case depending on the DEM used.
Next, we consider the effect of errors inξ caused by DEM errors. The estimator defined in Eq.(12) can be better expressed 550 using the DEM derived dry-snow phase sensitivity,ξ D =ξ +ξ ∆h , instead of trueξ whereξ ∆h is the error inξ D . Furthermore, for simplicity we assume that homogenous dry-snow change is the dominant phase contribution such thatΦ ∼ = ∆S ξ . Therefore, If we assume thatξ ∆h is uncorrelated withξ, which corresponds to a multiplicative bias towards lower amplitudes that is more pronounced in areas with low ξ diversity, i.e. low ξ 2 . We have not attempted to quantify ξ 2 ∆h for our TanDEM-X DEM derived ξ map but we made efforts to reduce it by smoothing the DEM (see Sect. 4.3), noting that this approach has its limits since over-smoothing will tend to increase ξ 2 ∆h .

∆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 ∆SWE estimation maps using an estimation window of 500 m × 500 m (ground range × azimuth). This corresponds to 46 snow-season interferograms over 2 partial and 8 complete snow-seasons and 49 non-snow 565 interferograms over 9 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. 5.4 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 13 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).

Comparison of ∆SWE Estimates with In Situ Measurements
In order to verify the ∆SWE estimates we compiled a number of mid-to-late winter snow transect surveys that are summarized in and their positions A-E are labelled in Fig. 1a. These surveys are comprised by: 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 GNWT-ENR at two-week intervals at a site with mixed trees and shrubs that is adjacent to the Inuvik Satellite Station Facility.
To compare with these absolute SWE measurements, we generated snow season cumulative SWE maps by integrating the sequential SlopeVar estimated ∆SWE maps through time; starting with the first dry-snow interferogram identified for each snow season. As discussed in Sect. 4.3, some ∆SWE estimates are labelled as invalid due to lack of a distinct periodogram 580 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 SlopeVar estimated 585 total SWE maps.
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. The in situ measurements were each made on a partic-590 ular date, whereas the InSAR cumulative SWE maps are produced at 24-day 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. correspond to areas of low topographic variation as indicated by the rightmost column of Table 2, and, as shown in Fig. 14,

600
(locations 'A', '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.

Comparison of ∆SWE Estimates with Predicted Temporal Change
The estimator results were examined temporally by comparing the spatial mean of each estimated ∆SWE map, after projection to ground range geometry, with the corresponding SWE change predicted by the European Centre for Medium-Range Weather

605
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 to 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-day intervals. The additional five intervals   each span 48-days due to missed SAR acquisitions. For the temporal comparison, each of these five intervals were split into 610 two 24-day intervals with half the estimated SWE change arbitrarily assigned to each interval. Figure 15a compares the continuous ERA5 24-day interval ∆SWE history time-series with the SlopeVar estimator results. Figure 15b shows the same data as a scatter-plot comparison. Both plots are labelled according to 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 615 dry-snow period (Oct-Mar), 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 Oct-Dec and Jan-Mar estimates because these correspond with the dry-snow period when the estimator is expected to provide useful results. Figure 15b shows that these do follow the 1:1 correspondence line. However, the Oct-Dec points appear to have a significant positive bias whereas the Jan-Mar points do not. This is quantified in Table 3, which shows that, in aggregate, the Oct-Mar time period has a correlation of 0.57 620 with respect to ERA5, which is moderately strong. There is an 8.6 mm per-interval ∆SWE bias for Oct-Dec but only 1.7 mm for Jan-Mar. This is consistent with the expected impact of early freeze-season heave, which, as shown in Sect. 5.3 is expected to yield a positive bias until freeze back of the active-layer completes. Considering the result from Fig. 12e where 10 mm of heave results in a mean ∆SWE bias of 9 mm, the Oct-Dec bias can be explained by (8.6 mm/0.9)(92 days/24 days) = 37 mm of mean upward displacement amplitude, which is a plausible value for the area during active-layer freeze back (Gruber, 2020).

625
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 unexpected since the snow-free period is expected to correspond to subsidence and therefore a negative bias with magnitude similar to but less than Oct-Dec 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. The RMSE of the 630 no-snow set is 4.2 mm which serves as a measure of the estimation uncertainty due to factors other than SWE inhomogeneity.
The Apr-Jun intervals are included for completeness even though these are not expected to agree because this corresponds to the snow melt period. This is especially apparent for the two intervals with large negative ERA5 ∆SWE values corresponding to periods of intense snow melt. 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.

Spatial Distribution of the Estimated ∆SWE
Spatial variations in the estimated ∆SWE 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 estimated ∆SWE maps since we expect the bias to differ between seasons. Figure 16   640 shows the temporal mean ∆SWE maps for the early snow season (Oct-Dec), late snow season (Jan-Mar) and all snow-free intervals. These show a moderate degree of negative correlation between the Oct-Dec and snow-free maps, especially in the upland area east of the delta-margin. The Pearson correlation coefficients between the three maps are: {R Oct-Dec/Jan-Mar = 0.28; R Oct-Dec/Snow-free = -0.52; R Jan-Mar/Snow-free = -0.17}. These are consistent with a dominant bias contribution from surface normal heave and subsidence contributing to the Oct-Dec and snow-free subsets but with opposite polarity. The weaker correlations 645 between these and the Jan-Mar 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.
Since the snow-free maps have 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)  and all snow-free periods (temporal subsets described in Table 3).
the top ξ diversity quartile of spatial samples (corresponding to threshold of ξ 2 1/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.

Correcting for the Dry-Snow Phase Component
The main focus of this article is estimation of ∆SWE. However, for repeat-pass InSAR applications other than ∆SWE estima-655 tion (e.g. surface deformation estimation), the dry-snow phase can be considered an interfering factor. In this case, knowing the spatial ∆SWE distribution for a given interferogram allows for removal of the snow phase component 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 the dry-snow phase corresponding to this initial estimate; and then (3) attributing the remaining InSAR phase to additional (high-spatial frequency) ∆SWE which can then be estimated per res-660 olution cell according to Eq.(5) . We performed a preliminary investigation of this by scaling the estimated ∆SWE maps by the dry-snow phase sensitivity map to generate spatial maps of estimated dry-snow phase 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 estimated ∆SWE maps. This approach therefore requires further investigation because, in order to maintain spatial phase smoothness required for subsequent phase-unwrapping, the ∆SWE maps likely 665 require imposition of smoothness constraints.

Spatial Resolution
Regarding choice of estimator window size, larger windows improve precision by mitigating the effect of spatially random noise. However, as discussed in Sect. 5, 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. 5.4, the precision of the SlopeVar estimator depends on both the coherence level and the spatial diversity of the SWE sensitivity, summarize by ξ 2 1/2 . These factors are spatially variable and therefore in order to generalize 675 the achievable spatial resolution to other locations one needs to consider both the degree of topographic slope variation and conditions affecting temporal 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.

Landcover and Drainage Basin Analysis
A detailed interpretation of the SWE results is not within the scope of this paper, but to demonstrate the advance in estimating SWE over large areas we present a first order assessment of the ∆SWE results in the context of landcover type, fire history, and three drainage basins. To simplify the assessment, we created four major landcover classes by merging the two NALCMS needleleaf forest classes into 'needleleaf forest', the broadleaf deciduous forest and mixed forest classes into 'decid-685 uous/mixed', the two shrubland classes into 'shrubland', and the four grassland or barren classes into 'grassland/barren' (Table   4). We disregard 'wetland', 'urban' and 'water'.  (Fig. 17a), and the cumulative mean SWE history per major landcover class (Fig. 17b)  Using the SAGA Terrain Analysis tool (Conrad et al., 2015), we demarcated three drainage basins located within our scene that have similar areas ranging between 22 km 2 and about 35 km 2 (Fig. 17c). Landcover class fractions per basin are shown in Table 4. To compare landcover class areas between drainage basins, we averaged the SWE estimates over the non-water areas (~89 to 95% total basin area) within each basin. Cumulative mean SWE curves (Fig. 17d) show that at the end of the winter,

700
Basin 3 typically has distinctly less SWE than Basins 1 and 2. Of the three basins, Basin 3 has the highest proportion of oldgrowth needleleaf forest, whereas Basin 1 and Basin 2 both have over double the fraction of deciduous/mixed forest. Notably, the dense deciduous/mixed forest east of Inuvik developed following a major forest-tundra fire in 1968 (extent outlined in Fig.   17c) that burned the widely spaced needleleaf forest (Mackay, 1995), resulting in dense growth of tall alder (Alnus crispa) and willows (Salix spp.) and regrowth of partly-burned paper birch (Betula papyrifera). The 1968 fire completely burned the 705 forested portion of Basin 1, the forest in Basin 2 was partially burned, but the forest in Basin 3 was beyond the fire limit (Fig.   1c). The SWE data suggest that over half a century since the initial event, the forest-tundra fire near Inuvik is an important influence on variation in hydrological dynamics in the area, as it is on variation in ground thermal conditions (Palmer et al., 2012), via ecological succession and its influence on snow cover. The mechanism is not investigated here, but may be related to net radiation and snow ablation in early spring due to the large influx of solar radiation and differences in albedo between 710 forest and non-forest areas (Rouse, 1984), as at this time of the year the tall shrubs are covered almost completely by snow and the needleleaf trees stand above the snow cover (Palmer et al., 2012).
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 dry-snow phase sensitivity map over a suitably 715 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.
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 by simulating an evolving snow-pack with the Snow-Model snow transport model and found that surface morphology caused sufficient ∆SWE inhomogeneity to generate modest 720 (few mm of SWE) estimation errors whereas spatial variations in vegetation height (modelled using NALCMS 2015) induced much larger errors (tens of mm). This leads us to conclude 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 by correlation with the dry-snow phase sensitivity map. Similarly, we also investigated solifluction and soil moisture phase, two bias sources expected to only affect snow-free interferograms used as zero ∆SWE 725 controls for method validation. We found that static atmospheric delay and soil moisture are likely not significant sources of bias and that both surface deformation and solifluction, if present, may lead to substantial estimation bias, generating ∆SWE errors at scales comparable to the deformation magnitude. We also considered the effect of DEM error and found that DEM errors at scales smaller than the estimation window size lead to damping multiplicative ∆SWE bias and therefore DEM filtering, pre-processing, or both should be performed to minimize this effect. 730 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 735 correlation and 7.9 mm RMSE (per 24-day period) over the Oct-Mar period. Seasonal analysis revealed a substantial early winter (Oct-Dec) bias (8.5 mm per 24-day period) that is largely absent in late winter (Jan-Mar) 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 on both terrain factors affecting the coherence 740 and on the degree of slope variation, both of which are location dependent.
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 745 sparsely map SWE change in otherwise flat areas. In other areas the constraints are similar in that they both require dry-snow