Permafrost soils in the Arctic store twice the amount of carbon found in the atmosphere (Hugelius et al., 2014;Ping et al., 2008). Over the past decades, warming has led to permafrost thawing (Jorgenson et al., 2006), which may result in the release of stored organic matter into the atmosphere as greenhouse gases and further amplify global warming (Serreze and Barry, 2011;Schaefer et al., 2014;Schuur et al., 2015). In permafrost regions, groundwater flows through the topmost portion of the soil, known as the active layer, that freezes and thaws annually (Woo, 2012;O'Connor et al., 2020). This groundwater flow contains carbon and is important in the export of carbon from land to the ocean and atmosphere (Kling et al., 1991;Stieglitz et al., 2003;Walvoord and Striegl, 2007;Vonk and Gustafsson, 2013;Paytan et al., 2015;Neilson et al., 2018). To understand how thawing permafrost contributes to the global carbon cycle, it is important to understand the hydrologic flow and transport processes in the active layer. Whether the carbon held by the active layer soils will be transformed to carbon dioxide or methane (a more powerful greenhouse gas), or whether it will flow towards rivers and lakes as dissolved carbon in groundwater, depends largely on the wetness or dryness (i.e., how much water is stored) of the active layer (Bond-Lamberty et al., 2016;Taylor et al., 2021).
Most of the arctic permafrost region is hard to access, and in situ observations of water storage and water flow in the active layer are extremely limited. Remote sensing techniques hold promise for local to regional observation of the hydrologic properties and hydrologic states of permafrost. For example, observations from the Gravity Recovery and Climate Experiment (GRACE) mission detect changes in permafrost water mass over a regional scale (Muskett and Romanovsky, 2009), but the spatial resolution is too coarse (∼ 100s of km) to be used in most hydrologic models (Text S1). In comparison, by measuring the phase difference between two paired radar images, Interferometric Synthetic Aperture Radar (InSAR) techniques estimate surface deformation between the two radar acquisition times along the radar Line-Of-Sight (LOS) direction (Rosen et al., 2000;Hanssen, 2001) at the spatial scale (∼ 10s to 100s meters spatial resolution) that overlaps with the scale of hydrologic field measurements and modeling grids. Although spaceborne InSAR has been used for estimating surface deformation associated with solid earth processes since the 1990s (Massonnet et al., 1993;Fialko et al., 2002;Pritchard and Simons, 2002;Shirzaei et al., 2013;Chen et al., 2014), it has only been recently used to estimate surface deformation associated with the seasonal freeze-thaw process of the soil active layer (Liu et al., 2010;Short et al., 2011;Antonova et al., 2018;Strozzi et al., 2018;Rouyet et al., 2019). Because ice density is less than water density (and thus ice volume is greater than water volume), the land surface subsides as the active layer thaws from winter to summer (Liu et al., 2010). Furthermore, InSAR-observed long-term subsidence trend signals over permafrost terrain have been used to study the deepening of the active layer due to wildfires or excessive melt of ground ice (Michaelides et al., 2019;Liu et al., 2014Liu et al., , 2015;;Iwahana et al., 2016;Yanagiya and Furuya, 2020;Abe et al., 2020;Eshqi Molan et al., 2018;Streletskiy et al., 2025).
Existing InSAR permafrost studies tended to associate the magnitude of the InSAR-observed thaw subsidence with the active layer thickness (Liu et al., 2012;Schaefer et al., 2015;Chen et al., 2021). However, the amplitude of the thaw subsidence and frost heave could depend on other factors such as sediment type and local topographic slope (Daout et al., 2017). Chen et al. (2020) found that the amplitude of the seasonal thaw subsidence is proportional to the amount of water stored in the saturated active layer at the end of a thaw season. This is consistent with findings from recent studies that InSAR-derived seasonal subsidence rates reflect spatial soil moisture patterns (Chen et al., 2022(Chen et al., , 2023;;Widhalm et al., 2024). In this paper, we further established a conceptual model that relates InSAR seasonal thaw subsidence observations to soil water storage in the saturated active layer. Our goal is to advance InSAR techniques for the high-resolution mapping of water storage above-permafrost. To demonstrate this, we mapped soil water stored in the saturated active layer using ALOS PALSAR data over a much larger area in the Arctic Foothills than used in Chen et al. (2020). We validated the InSAR results using in-situ soil measurements collected at more than 200 remote sites within ∼ 100 km of the Toolik Field Station as well as optical imagery and land cover maps.
Our results show that InSAR soil water storage estimates derived from two separate satellite frames are consistent with in-situ observations under different vegetation covers. An important new contribution of this work is on uncertainty quantification.
We determine the primary error sources in Toolik ALOS PALSAR Line-Of-Sight (LOS) measurements, and we discuss how errors in InSAR LOS measurements can be linearly related to errors in soil water storage estimates.
In this section, we first describe the conceptual model that relates the soil water storage in the saturated active layer to ground ice melting during summer thaw seasons (Section 2.1). We then explain our InSAR processing strategy for estimating average seasonal thaw subsidence from ALOS PALSAR data (Section 2.2), and discuss key error sources in InSAR measurements 60 (Section 2.3). Finally, we review available field observations and strategies for validating the InSAR results (Section 2.4).
Our study site near Toolik Lake is in continuous permafrost of the upper Kuparuk River basin on the North Slope of Alaska (Figure 1). In 1987, the Toolik Field Station became part of the NSF Long Term Ecological Research program (LTER), which maintains long-term meteorological, ecological, and hydrological observations of the Arctic Foothills (Hobbie and Kling, 65 2014). The availability of the long-term databases of many basic parameters of the permafrost system makes the Toolik area an excellent site for studying how different soils control hydrological dynamics and may change as the climate warms and permafrost thaws.
Based on in situ thaw measurements at Toolik, the active layer starts to thaw in early June, and the maximum seasonal thaw typically occurs in late August (Romanowicz and Kling, 2022). Thawing processes typically slow down around the time of maximum thaw, because (1) thermal diffusivity of ice is larger than that of liquid water; and (2) heat takes much longer to diffuse through a thicker active layer soil column. Due to the density difference between ice and liquid water, the land surface subsides during the thaw season, with the opposite occurring when the active layer refreezes (Short et al., 2011;Painter et al., 2016;Sjoberg et al., 2016;Antonova et al., 2018;Strozzi et al., 2018). The maximum seasonal thaw subsidence (d season ) is proportional to the amount of water stored in the saturated active layer that experiences the ice-to-water phase change in a thaw season (denoted as z water ) following (Liu et al., 2012;Chen et al., 2020):
where ρ w and ρ i are the density of water and ice respectively. Figure 2 illustrates the definition of z water in this work. Here we exclude soil water stored above the water table (tension or unsaturated zone water) in the z water estimation. Because the porosity in the organic soil layers is high (∼ 0.78-0.98), water in the unsaturated zone can expand to fill the empty pore space during freezing without contributing to surface deformation. In this study, we assume the density of water is a constant value of 0.997 g/cm 3 , and the density of ice is a constant value of 0.917 g/cm 3 . Our calculation does not account for variations in subsurface water and ice density due to capillarity associated with surface tension, cation hydration, surface hydration, and interlamellar cation hydration (Zhang and Lu, 2018).
Equation (1) shows that d season is proportional to z water rather than to the Active Layer Thickness (ALT). For example, minimal thaw subsidence signals would be observed over thick but dry active layers (Chen et al., 2020). This means that active layers with higher ice-to-water content are expected to experience larger thaw subsidence, which may have no bearing on ALT.
Furthermore, the active layer (liquid) water storage balance can be defined as:
where ∆S is the change in total soil water storage of the active layer. P , ET , and Q stand for changes in soil water storage due to precipitation, evapotranspiration, and runoff, respectively. A is the amount of soil water change associated with the active layer freeze and thaw process detectable by InSAR. When the active layer thaws during the summer, A > 0; when the temperature drops in autumn and the active layer refreezes, A < 0. In the case that InSAR-observed seasonal thaw subsidence signals are similar over multiple years, the amount of water that experienced the annual freeze-thaw cycle does not change much during this period (the net water drainage P -ET -Q ≈ 0).
We emphasize that many geophysical processes can lead to surface deformation in permafrost terrain detectable by InSAR (Zwieback et al., 2024b). For example, solifluction and other slope creep processes may produce long-term downward deformation trends in regions with large slope angles (Dini et al., 2019). Post-glacial rebound and tectonic motions typically vary at 100-km or larger spatial scales and can be considered as nearly spatially uniform over our study area (Liu et al., 2010;Stephenson et al., 2022). Given that InSAR measures relative deformation with respect to a local reference point, InSAR is only sensitive to spatially varied surface deformation over the study area. Hydrological loading and unloading can produce (2) the water-saturated active layer (thickness of s). The upper surface defines the water table, and the lower surface defines the ice table that separates the thawed active layer from frozen ground (or the permafrost layer at maximum annual thaw depth); and (3) the permafrost (thickness of p), which may or may not be saturated with water as ice. In column (A), zwater represents the amount of water stored in the saturated active layer. While there is tension water in the upper unsaturated zone in (B), zwater is the same in columns (A) and (B). The reason is shown in (C), where the entire soil column has now frozen. The saturated water freezes and the expansion heaves the soil column above and the ground surface (s + 0.09zwater), while the tension water freezes but expands into pores containing soil atmosphere and thus does not contribute to deformation of the ground surface (thickness u does not change).
millimeter-level surface deformation signals (Liu et al., 2010), which is much smaller than centimeter-level freeze-thaw deformation. Furthermore, peat accumulation processes (Jones et al., 2017) may lead to a long-term deformation signal detectable by InSAR, and surface erosion can cause changes in surface scattering properties that decorrelate radar phase measurements (Zebker and Villasenor, 1992). In Section 2.2, we discuss how to extract long-term and seasonal deformation signals from
InSAR observations. The magnitude and characteristics of deformation signals, combined with in-situ observations (Section 2.4), can be used to determine the primary geophysical processes that contribute to the observed deformation patterns.
Interferometric SAR (InSAR) computes the phase difference between two Synthetic Aperture Radar (SAR) images. The resulting interferogram can be used to infer a map of surface deformation between two SAR acquisition times along the radar Line-Of-Sight (LOS) direction (Hanssen, 2001). More specifically, a phase cycle of 2π (in radians) equals to λ/2 of LOS deformation, where λ is the radar wavelength. For L-band ALOS PALSAR data, λ equals 24 cm, and thus a phase cycle of 2π represents 12 cm of LOS deformation occurred between two radar acquisition times.
In a recent proof-of-concept study (Chen et al., 2020), we processed 12 L-band ALOS PALSAR scenes (Table B1) acquired during summer seasons (June to October) between 2006 and 2010 from path 255 frame 1370 over our study region (Figure 1).
Note that we excluded all winter scenes acquired between November and April because the observed phases in winter-winter interferograms are likely related to variations in snow accumulation and snow redistribution, which is not the focus of this study.
We first solved for the long-term LOS deformation trend at a pixel of interest based on a stacking approach (Sandwell and Price, 1998;Lyons and Sandwell, 2003;Rouyet et al., 2019). That is, averaging all interferograms that contain minimal seasonal deformation signals (e.g., a July-to-July pair) and relatively large long-term signals (e.g., span multiple freeze-thaw cycles).
An important finding of this pilot study was that no detectable long-term deformation trend above the InSAR measurement noise level was observed outside the 2007 Anaktuvuk River fire scar (Figure 1) during the study period of 2006 to 2010.
This allowed us to substantially simplify our InSAR processing strategy for reconstructing seasonal freeze-thaw deformation patterns over undisturbed permafrost terrain. We estimated the LOS deformation signatures due to the seasonal active layer freeze-thaw processes between (i) early June and late July, (ii) late July and early September, and (iii) early September and late October by averaging all interferograms that span these periods regardless of how many years those interferograms span.
The averaged LOS deformation between early June and late July was used as an approximation of the maximum seasonal LOS deformation because no ALOS acquisitions were made over the study area around the time of the maximum thaw (late August at Toolik area). This approximation is reasonable given that a few centimeters of late summer thaw (August) of the relatively dry low-porosity mineral layer does not cause surface thaw subsidence detectable by InSAR (O'Connor et al., 2020;Chen et al., 2020). For example, 10 cm of thaw of the saturated organic layer (∼ 90% porosity) would lead to ∼ 0.8 cm of thaw subsidence. In comparison, 10 cm of thaw of the saturated mineral soil (∼ 20% porosity) would lead to < 0.2 cm of thaw subsidence.
We note that InSAR measures the change in distance between the antenna and the ground object, known as the LOS direction.
Assuming the horizontal motion of the land surface is negligible, we converted InSAR seasonal deformation estimates along the LOS direction (d LOS ) to seasonal vertical thaw subsidence estimates d season as:
where e 3 is the vertical component of the radar LOS direction unit vector e = [e 1 , e 2 , e 3 ]. The LOS unit vector e can be computed based on the known satellite position and ground pixel location in the Earth-centered, Earth-fixed (ECEF) coordinate system, and then converted to the local east-north-up (ENU) system (Misra and Enge, 2011). For the ALOS ascending imaging geometry over the Toolik area, e = [0.61, 0.13, -0.78] at the mid-swath, and the variation of e 3 across the entire swath is minimal (less than 3%). This means that ∼ 5 cm thaw subsidence can cause 4 cm positive LOS deformation for the Toolik ALOS PALSAR case.
To confirm that our InSAR processing strategy is suitable for studying the active layer freeze-thaw process over vast areas, here we analyzed an additional 11 L-band ALOS PALSAR scenes (Table B1) acquired during summer seasons (June to October) between 2006 and 2010 from path 255 frame 1380 (Figure 1). We merged interferograms from the same path but two different frames by calibrating the phase differences within the overlapping regions of the two frames. A sample merged interferogram is shown in Figure B1, and the same reference point (68.83 • N, 150.23 • W) as our previous study (Chen et al., 2020) was used to calibrate all interferograms. We chose this reference point because it is in a dry highland area of relatively flat terrain, and the expected seasonal deformation is minimal. Only ∼ 4% of interferograms contain visible phase decorrelation artifacts outside the fire scar, and the overall phase coherence (Figure B2(a)) of the remaining interferograms is comparable to the sample interferogram (Figure B1). We masked out pixels with amplitude dispersion < 0.25 and pixels with phase coherence < 0.2 to exclude water bodies and the area burned by the 2007 Anaktuvuk River fire (Figure B2(b)). A comparable pixel mask can also be generated using the North Slope Science Initiative (NSSI) land cover GIS Data.
Similar to our previous study, we found that the long-term subsidence trend is negligible outside the fire scar (Figure B3(d)).
This allows us to follow the same processing strategy as our previous study to extract seasonal deformation between (i) early
June and late July, (ii) late July and early September, and (iii) early September and late October by averaging all interferograms that span these periods regardless of how many years those interferograms span (Figure B3(a)-(c)). We note that averaging interferograms that contain the common signal of interest (stacking) reduces the impact of random tropospheric turbulent noise by ∼ √ N , where N is the number of independent SAR acquisitions (Sandwell and Price, 1998;Chen et al., 2020). A thaw subsidence pattern similar to the final stacking solution was identified from all individual interferograms that span a common season (e.g., early June to late July). The averaged LOS deformation between early June and late July was used as an approximation of the maximum seasonal LOS deformation because no ALOS acquisitions were made over the study area around the time of the maximum thaw (late August at Toolik area). This approximation is reasonable given that a few centimeters of late summer thaw (August) of the relatively dry low-porosity mineral layer does not cause surface thaw subsidence detectable by InSAR (O'Connor et al., 2020;Chen et al., 2020). For example, 10 cm of thaw of the saturated organic layer (∼ 90% porosity) would lead to ∼ 0.8 cm of thaw subsidence. In comparison, 10 cm of thaw of the saturated mineral soil (∼ 20% porosity) would lead to < 0.2 cm of thaw subsidence.
Based on Equation (1) and Equation (3), we further established a linear relationship between InSAR LOS deformation observations and the amount of water in the saturated active layer that experiences the ice-to-water phase change (z water ) as:
where ρ w and ρ i are the density of water and ice, respectively. e 3 = -0.78 is the vertical component of the ALOS LOS direction
unit vector e = [e 1 , e 2 , e 3 ] as defined in Equation ( 3). This equation shows that InSAR-observed seasonal thaw subsidence is proportional to the active layer water storage z water . For the ALOS Toolik case, 5 cm InSAR LOS deformation measurements (d LOS ) can be related to 70 cm of saturated active layer soil water column (z water ), 1 cm errors in InSAR LOS deformation measurements can lead to 14 cm error in z water estimates. We note that Equation (4) employs the assumption that the horizontal motion of the land surface is negligible. Our study site is a transitional region located between the Coastal Plain and the steep mountains of the Brooks Range, which consists of gently rolling hills and broad exposed ridges that extend along the northern flank of the Brooks Range. Given that the long-term subsidence trend is negligible outside the fire scar and seasonal deformation signatures follow the expected seasonal freeze-thaw patterns (Figure S3), InSAR observations at our study site are primarily related to the volume change associated with water-to-ice phase change rather than slope creep processes. For a 5% slope angle, 1 cm of freeze-thaw deformation perpendicular to the land surface leads to 0.87 mm horizontal deformation and 9.96 mm vertical deformation. For a 10% slope angle, 1 cm of thaw deformation perpendicular to the land surface leads to 1.74 mm horizontal deformation and 9.85 mm vertical deformation. Because the slope angle at most radar pixels is less than 10%, we conclude that the assumption of negligible horizontal motion is reasonable.
To quantify errors in InSAR-based z water estimates, here we evaluate major error sources in InSAR LOS deformation solutions (d LOS ), which can be written as (Zebker and Villasenor, 1992;Zebker et al., 1994Zebker et al., , 1997)):
where λ is the radar wavelength (24 cm for L-band ALOS data), and φ is the average phase of all interferograms that contain the common seasonal deformation signal of interest. The remaining noise terms on the right-hand side are errors due to topography-related artifacts (∆d dem ), phase decorrelation (∆d decor ) and phase unwrapping errors (∆d unwrp ), orbital errors (∆d orb ), atmospheric (∆d atm ) and ionospheric (∆d iono ) artifacts, and other smaller error terms associated with thermal and soil moisture effects (∆d n ).
In the Toolik ALOS InSAR data analysis, we excluded ∼ 4% of interferograms containing visible phase decorrelation and phase unwrapping errors. Long-wavelength phase signatures, varying at spatial scales of tens to hundreds of kilometers and potentially caused by orbital errors and tropospheric or ionospheric noises, were removed as a planar phase ramp from each interferogram (Staniewicz et al., 2020;Wang and Chen, 2022;Zebker et al., 2023). The deramp process does not remove localized freeze-thaw deformation patterns that vary from hilltop ridges to the lowland valleys and riparian zones on the spatial scale of ∼ 100s of meters. Interferograms formed by the SAR scenes acquired on 8 September 2008 (for frame 1380) and SAR scenes acquired on 22 July 2007, 24 July 2008, and 14 September 2010 (for both frame 1370 and 1380) were also excluded because of severely distorted ionospheric artifacts (Gray et al., 2000;Wegmuller et al., 2006;Chen and Zebker, 2012;Fattahi et al., 2017). Because of a cool and dry tundra climate and relatively small elevation variation (∼ 200-300 meters), the stratified tropospheric noise component (Doin et al., 2009) is minimal over the study site. Given that long-wavelength tropospheric and ionospheric noise was removed during the planar ramp removal process, the residual atmospheric noise term (e.g., due to localized temperature or water vapor variations) is mostly random at time scales longer than one day, but is correlated in space and typically increases with distance from the InSAR reference point (Emardson et al., 2003;Staniewicz et al., 2020).
Assuming a 2 cm tropospheric error in each ALOS PALSAR interferogram (Zebker et al., 1997;Emardson et al., 2003),
the turbulent random noise level can be reduced to less than 1 cm after stacking four interferograms formed from four SAR acquisitions. In the remainder of this section, we focus on the dominant error term associated with topography-related artifacts for the ALOS Toolik case.
At a pixel of interest, an error in the Digital Elevation Model (DEM; δ) with respect to the reference pixel can lead to an error (∆d dem ) in the LOS deformation estimates as (Berardino et al., 2002;Werner et al., 2003;Fattahi and Amelung, 2013):
where B perp is the perpendicular component of the InSAR spatial baseline, which can be calculated from known radar imaging geometry. For ALOS interferograms, B perp typically ranges from several hundred to several thousand meters. r is the distance between the radar antenna and the ground pixel, and θ l is the radar look angle. Because the look angle of ALOS PALSAR does not vary much over the ∼ 60 km radar swath, both r ∼ 850 km and θ l ∼ 34 degrees can be approximated as constant values for all ALOS interferograms collected from the same path and frame.
In this study, we removed the topographic phase during interferogram formation using the Arctic DEM (10-meter resolution and resampled to a 30-meter grid) data (Porter et al., 2018), which are widely used in the Arctic community because of its panarctic coverage and high quality (Tozer et al., 2019). Interferograms with comparable quality can be also generated using the Kuparuk River watershed DEM (Chen et al., 2020). While the Kuparuk River watershed DEM has been thoroughly validated and highly accurate (Nolan, 2003b), it does not have complete spatial coverage over the entire study area. It is common to assume that ∆d dem is linearly proportional to B perp . This assumption is valid when δ in Equation ( 6) is introduced by errors in the DEM dataset itself (thus δ is the same for all interferograms). Because thaw subsidence patterns over undisturbed permafrost terrain are expected to be spatially coherent, phase discontinuity in interferograms was visually inspected. If the magnitude of these artifacts is linearly proportional to InSAR perpendicular baseline B perp , they are likely associated with the errors in the Arctic-DEM dataset, given that thaw subsidence signals do not depend on B perp . Furthermore, we discovered that δ can also be introduced by misregistration of the DEM and a SAR image. For example, as shown in Figure 3 steep terrains. In this study, we employed the same image co-registration routine as the standard InSAR processing software such as the InSAR Scientific Computing Environment (ISCE) (Rosen et al., 2012) and GMTSAR (Sandwell et al., 2011). The 2-D cross correlation method for image alignment can achieve sub-pixel accuracy in most cases. However, the alignment can be worse than 1 pixel, because SAR images and DEM data were acquired from sensors with different spatial resolutions and imaging geometries. To better understand these pixel-mismatching artifacts, we approximated the DEM error δ due to pixel misregistration as the difference between the Arctic-DEM and the shifted Arctic-DEM in east/west and north/south directions.
For example, the DEM error δ i,j due to 1 pixel misregistration to the east at pixel (i, j) can be written as:
where h ij is the Arctic-DEM at pixel (i, j). Similarly, we can approximate δ i,j due to 1 pixel misregistration to the south at pixel (i, j) as h i,j -h i,j+1 . We then calculated the expected LOS errors ∆d dem due to δ based on Equation (6) for a given imaging geometry and perpendicular baseline. Results from these numerical experiments were then compared to actual InSAR LOS observations across hill ridges.
Our InSAR thaw subsidence estimates were validated using a relatively large number of field observations collected within ∼ 100 km of Toolik Field Station (Figure 1) in 2018 (August 15 -August 24) and 2019 (July 26 -August 3). Particularly, the amount of water stored in the saturated active layer can be quantified by determining saturated active layer thickness and porosity. Tundra soil in the Toolik area consists of three layers from top to bottom: the acrotelm (peat that contains living plants), the catotelm (peat that contains dead plant materials), and the mineral soil (Figure 4). The thickness of these three soil layers and the depth to the water table were measured at each sampling site. The porosity (ϕ) of each soil core sample was also measured to characterize the water-holding capacity of active layer soils. z water can then be calculated as: where z si and ϕ i are the saturated thickness and the porosity of the i th soil layer. Here, we assume the soil column below the water table is fully saturated. We also note that the mineral soil layer has much lower porosity (thus much less waterholding capability) than organic soil layers. For example, a fully saturated, 10-cm-thick acrotelm layer with a porosity of 0.90 contributes to 9 cm of z water , while a fully saturated, 10-cm-thick mineral soil layer with a porosity of 0.20 only contains 2 cm of z water .
To jointly analyze remote sensing and in situ observations, an exact point-to-point comparison is challenging, if not impossible, because they were collected at very different spatial and temporal scales. A pixel in an InSAR-derived deformation map is ∼ 100-by-100 meter, while field measurements were collected at sites with size ∼ 900 cm 2 (30-by-30 cm plots). To overcome this challenge, we designed a statistical comparison approach. This was done by fitting probability density functions (PDFs)
to the empirical distributions (histograms) of the in-situ soil property measurements, including the thickness and porosity of the acrotelm, the catotelm, and the mineral soil as well as the depth to water table (O'Connor et al., 2020;Chen et al., 2020), and using these distributions to calculate the range of possible thaw subsidence. There are also sources of error in the property measurements, which are (1) errors from reading the measured value, which is typically small (e.g., <0.5 cm for thaw depth measurements from probing), and (2) in situ measurements varying due to the sub-meter-scale heterogeneity of arctic soils.
To reduce estimation bias, we targeted specific vegetation cover types and soil layers needing larger sample sizes over time to improve statistical robustness. The PDF fitting results did not change much after a second year of sampling, indicating that the sample size in this study is sufficiently large to capture the statistical characteristics of soil properties. We drew random samples from the PDFs of soil properties, and calculated the statistical distribution of z water following Equation (8).
Finally, we validated InSAR-observed z water using field-based predictions of z water at different vegetation types. For the purpose of studying active layer soil properties, we grouped various subclassifications of vegetation types over the study area (Walker and Walker, 1996;Stow et al., 2004;Walker et al., 2017) into four primary land cover types: "sedge", "tussock", "woody shrub", and "sparse vegetation" (O'Connor et al., 2020). The sedge land cover typically occurs in wet to saturated sites (e.g., riparian zones) and may occasionally mix with shrub mounds on slightly elevated ground. The tussock land cover is distributed broadly from ridges to riparian zones. The term "woody shrub land cover" refers to areas dominated by woody-stem plants, which include both woody shrubs along the water tracks and heath vegetation on ridges. Because the soil in water tracks typically consists of well-drained acrotelm with underlying gravel and boulders, we did not collect soil samples in the water tracks. As a result, this study focuses on soil measurements collected over three land cover types: sedge, tussock, and woody shrub on ridge-tops (referred to as "heath"). Photographs of the land cover types are shown in Figure B4. We also identified the land cover type of each InSAR pixel using the North Slope Science Initiative (NSSI) Land Cover Map (Payne et al., 2016), which does not distinguish between woody shrubs within water tracks and heath on hill ridges.
3.1 InSAR-estimated soil water storage in the saturated active layer InSAR-observed average seasonal thaw subsidence estimates (2006)(2007)(2008)(2009)(2010) between early June and late July from two independently processed ALOS PALSAR frames are consistent with no visual discontinuity or artifacts (Figure 5(a)). This confirms that our InSAR analysis is robust for reconstructing thaw subsidence over permafrost terrain. Ninety-five percent of the observed thaw subsidence ranged from 0 cm to 5.4 cm, which correlates with the topography as well as the watershed and river network morphology (Figure 5(b)). The drier ridge-top areas usually show less than 2 cm thaw subsidence, while the wetter valleys and riparian zones show up to 6 cm subsidence. Thaw subsidence of ∼ 4 cm is observed near the transition zone as the steeper hilly terrain (south) transitions to flatter plains (north). Based on Equation (1), 1 cm thaw subsidence (∼ 0.78 cm LOS deformation) is caused by an ∼ 11 cm z water . Ninety-five percent of z water estimates range from 0 to 62 cm in the Toolik area, with up to 75 cm z water observed in the wettest riparian zone after removing less than 3% of outliers (Figure 5(c)). Here, pixels are marked as outliers if they are larger than the upper adjacent value, which is, by definition, the largest observation that is less than or equal to the threshold located at the 1.5 Inter-Quartile Range (IQR) above the upper quartile (Q3). The spatial variation of z water is consistent with groundwater flows and accumulation from the higher ridges to the flatter riparian zones and valleys (e.g., Figure 6). Large z water values are often observed in wet local low regions.
Because the amount of soil water influences the type of vegetation that can grow, the spatial pattern of InSAR-observed z water correlates with land cover types (Figure 5(d)). We found that land cover type indicates the characteristics of soil stratigraphy (Figure B5), and each soil layer possesses different characteristics (e.g., porosity and thickness) that influence the waterholding capability of the active layer. For example, water-loving sedges tend to grow on wet soils with a thick porous catotelm layer and a shallow water table, while heath vegetation is often found on dry hill ridges with a thin catotelm layer and a deep water table. To further illustrate the spatial correlation between the amount of soil water and land cover types, Figure 7 (a)-(c)
shows a zoomed-in area near the Toolik Field Station from frame 1370 (with the location outlined by the purple dashed line in Figure 5(c)), where all three major land cover types are present. We found that (1) sedges are often distributed over regions 1). Figure 7 (d)-(f) shows another zoomed-in area from frame 1380 (with the location outlined in a red dashed line in Figure 5(c)), where the terrain transitions from rolling hills to coastal plains. This region is wetter than the Toolik Lake area (Figure 7 (g) and Table 1). Here, tussock is the dominant land-cover type, and water-loving shrubs and sedges are distributed along the water tracks (visible in the optical image).
To validate InSAR z water estimates, the expected distribution of z water was also calculated from field measurements col-315 lected near Toolik (Figure 1) following Equation ( 8). We found that z water estimated from field and satellite observations is statistically consistent (Figure 8(a)-(d)). The median z water values derived from field data, ALOS PALSAR path 255 frame 1370 data, and ALOS PALSAR path 255 frame 1380 data are 30.8 cm, 28.4 cm, and 28.1 cm, respectively. The standard deviation values of z water derived from field data, ALOS PALSAR path 255 frame 1370 data, and ALOS PALSAR path 255 frame 1380 data are 21.2 cm, 18.9 cm, and 17.3 cm, respectively. Both InSAR and field observations are also consistent over three major land cover types (Figure 8 (e)).
Based on in situ data, z water has a median of 19.5 cm, 24.1 cm, and 34.9 cm for heath, tussock, and sedge land covers. This is consistent with InSAR observations over three land cover types: 20.0 cm for woody shrubs, 28.8 cm for tussocks, and 35.9 cm for sedges (Table 1). InSAR observations for each land-cover type generally show a larger variation of z water compared to field observations. This is likely because InSAR pixels were classified using the NSSI land cover map, which is less accurate than field-based land-cover classification at each sampling site. We also note that the land cover map used for classifying InSAR pixels does not distinguish woody shrubs located near the water tracks and those on dry ridge tops, while during field data collection, we only sampled dry heath land covers on the ridges. This is another reason that InSAR woody shrub observations show a larger variation compared to the other two land cover types.
Due to the remote nature of the study area, the number of available field observations is limited. We acknowledge that most field sites are located within the coverage of ALOS PALSAR path 255 frame 1370. Nevertheless, both frames exhibit similar land cover type combinations, suggesting similar climatic and geological settings (Figure 5 (d)). While some of our field sites are located outside the radar footprint (Figure 1), field observations at these sites follow similar statistical distributions as those sites located within the radar footprint. Furthermore, InSAR-observed average seasonal thaw subsidence estimates (2006)(2007)(2008)(2009)(2010) from two independently processed ALOS PALSAR frames are consistent with no visual discontinuity or artifacts at the frame boundary (Figure 5 (a)). This indicates the InSAR processing strategy produced consistent thaw subsidence estimates. We observed surface subsidence due to the thawing of the active layer from early June to late July and a net surface uplift between late July and late October resulting from the refreezing of the soil (Figure B3(a)-(c)). Minimal long-term subsidence trends were observed outside the fire scar (Figure B3(d)). These observations confirm that the observed InSAR seasonal deformation signals at our study site are primarily related to the volume change associated with water-to-ice phase change. We translated
InSAR measurements into soil water storage in the saturated active layer following Equation (4), which does not require additional information on soil properties. We used in-situ soil measurements as an independent validation for InSAR results in regions wherever it is possible, and our goal is to develop a remote sensing technique that can fill the observational gaps in remote Arctic areas with no in-situ observations.
We also acknowledge that field observations were collected in 2018 and 2019, while ALOS PALSAR InSAR data were used to estimate the average seasonal thaw subsidence between 2006 and 2010. In-situ thaw depth measurements show that the August 11 thaw depth (∼ 40 cm) at the Toolik long-term monitoring site has increased very slightly since 1990. At the Imnavait site, the August 11 thaw depth increase between 2006 and 2010 is ∼ 5 cm. At both sites, we did not observe any long-term subsidence trend above the InSAR noise level (Chen et al., 2020). This is because a 5 cm thaw of the low porosity (thus less water-holding capacity) mineral soils was unlikely to cause any soil water content increase that is detectable by InSAR. Therefore, our study focuses on the comparison between InSAR average seasonal thaw subsidence estimates (2006)(2007)(2008)(2009)(2010) and recent field observations over undisturbed permafrost terrain (relatively stable with long-term changes undetectable by InSAR). Due to the limited ALOS PALSAR temporal sampling rate, the investigation of inter-annual variability of InSAR thaw subsidence patterns is outside the scope of this work. Future work can focus on studying how the signal magnitude of seasonal thaw subsidence changes over multiple years using Sentinel-1 data collected with 6-12 day revisit cycles (Zwieback and Meyer, 2021;Zwieback et al., 2024a).
Because DEM error is the dominant error source in the ALOS Toolik InSAR data, here we discuss errors in thaw subsidence estimates associated with errors in the DEM data. When there is an error in the DEM data, a similar signature may be observed in the InSAR surface deformation observations. For example, Figure 9 shown in Fattahi and Amelung (2013), DEM errors in the LOS measurement are linearly proportional to the perpendicular baseline for a fixed error in the Arctic-DEM (Equation ( 6)).
To better illustrate that our observations are consistent with existing InSAR DEM error studies, we analyzed all 51 interferograms from path 255 frame 1380 and identified a linear relationship between the InSAR perpendicular baseline and the thaw subsidence errors at P1-P2 across the discontinuity line (marked in Figure 9 (c) and (d)). The observed linear slope (Figure 9 (f)) suggests a 1.16-meter error in the Arctic-DEM, which is consistent with the ∼ 1-2 m DEM discontinuity observed in the actual Arctic-DEM data (Figure 9 (e)). Existing InSAR studies tend to assume non-negligible DEM artifacts are typically observed in areas with steep terrain (Li et al., 2014;Staniewicz et al., 2020;Zhou et al., 2020). However, our results show that more than 1.5 cm LOS errors associated with inaccurate DEM can be observed over relatively flat areas in ALOS PALSAR interferograms with B perp values > 4000 m. This error is on the same order of magnitude as the observed centimeter-level thaw subsidence signals; thus, it is not negligible. These artifacts caused by errors in the Arctic-DEM can be mitigated by ( 1)
excluding interferograms with larger B perp or (2) estimating and removing a phase component that is proportional to B perp from all interferograms (Berardino et al., 2002). In our case, the discontinuity line is no longer observable after applying the stacking technique (as shown in Figure 5(c)), thus having minimal impact on the final z water estimates.
An important finding of this study is that pixel misregistration between the DEM and a SAR image can lead to DEM-related errors in InSAR LOS measurements. As an example, Figure 10(a) shows an interferogram formed by SAR images acquired on 8 June 2008 and 30 July 2010, and Figure 10(b) shows an interferogram formed by SAR images acquired on 3 June 2006 and 27 July 2009. Because the long-term subsidence trend is negligible, similar early June to late July thaw subsidence patterns are present in these two interferograms. To quantify InSAR LOS errors in areas with larger percent slopes (> 7.5%), we zoomed in to a hilly area near Imnavait Creek with a slope between 8.4% and 11.0%, and calculated the phase difference (a phase difference of 2π is equivalent to 12 cm LOS distance difference for L-band ALOS PALSAR data) between points P E on 390 the east-facing slope and P W on the west-facing slopes across a hill ridge. Because water flows away from ridges with very small catchment areas, we expect to observe minimal freeze-thaw deformation on either side of the dry hill ridge. However, the phase difference between P E and P W is 1.23 rad (an equivalent LOS deformation error of 2.3 cm) for the interferogram that spans 8 June 2008 and 30 July 2010 (Figure 10(d)), and 0.92 rad (an equivalent LOS deformation error of 1.7 cm) for the interferogram that spans 3 June 2006 and 27 July 2009 (Figure 10(e)). Although the perpendicular baselines of these two ALOS interferograms across ridges (Figure 11). These phase artifacts are most noticeable in interferograms formed using one of the three SAR images (acquired on 8 June 2008, 24 October 2008, and 27 October 2009), which likely suffer more severe pixel misregistration errors than other SAR scenes. By contrast, an error in the DEM dataset itself can lead to LOS errors that are linearly proportional to the perpendicular baseline (B perp ) in all interferograms (Section 3.2), while long-term deformation 400 trend signals are proportional to temporal baselines (e.g., related to various slope processes as discussed in (Dini et al., 2019)).
To confirm that the observed InSAR phase errors between east-facing and west-facing slopes are indeed associated with DEM-SAR misregistration, the Kuparuk River watershed DEM data (Nolan, 2003a) were shifted to the east by 1 pixel (∼ 12 m). The difference between the original and shifted DEM was used as an approximation of the DEM error (δ) caused by 1-pixel misregistration to the east as described in Equation ( 7). In this case, a positive DEM error on the east-facing slope with respect to the west-facing slope was observed (Figure 12 (a)). Similarly, a negative DEM error was observed on the east-facing slope with respect to the west-facing slope when the original DEM was shifted by 1 pixel to the west (Figure 12 The land-surface slope is another factor that could affect the magnitude of DEM errors δ at different pixels. We classified radar pixels into four groups based on their percent slope. For each group, phase errors associated with pixel-mismatching (1 radian phase error is equivalent to 1.9 cm LOS error) were calculated for the case that the location of the ridge is off by 1 pixel (∼ 12 m) to the east and a perpendicular baseline of 5104 m. This scenario can be considered as the error upper bound because
(1) the perpendicular baselines of L-band ALOS Toolik data are typically less than 5104 meters, and (2) the amount of pixel misregistration in the standard InSAR processing software packages is on the order of sub-pixels. We found that topographic artifacts associated with DEM-SAR pixel misregistration are most noticeable in areas with a slope larger than 10%, and the majority of the surface area with a low slope (0-5%) show negligible phase errors due to DEM-SAR pixel misregistration (Figure 13). To further illustrate this, Figure 14 shows the amount of the LOS errors (in cm) in all interferograms at a steep area and a flat area. Up to ∼ 6 cm LOS errors associated with pixel misregistration were observed in the steep area, while < 1 cm LOS errors were observed in the flat area. 8.4% to 11.0% (blue dots) and around P f with percent slopes 2.0% to 1.8% (orange dots). The location of these pixels is shown in Figure 10.
In summary, we found that (1) the DEM error δ increases as the amount of pixel misregistration increases for a given pixel location (Figure 12); (2) the DEM error δ increases with local slopes at different pixel locations for the same amount of pixel misregistration (Figure 13 and Figure 14); and (3) the relationship between the LOS error due to δ and the perpendicular baseline is non-linear. We emphasize that both the amount of pixel misregistration and the slope influence the DEM error δ. For example, δ equals 0 if there is no pixel misregistration. At a given pixel location, δ increases as the amount of pixel misregistration increases. For a fixed amount of pixel misregistration, δ increases as the slope increases at different pixel locations. This means that the perpendicular baseline is not the only factor that controls the observed DEM artifacts in InSAR LOS measurements ∆d dem . It is difficult to fully correct the pixel misregistration because SAR images and DEM data were acquired from sensors with different spatial resolutions and imaging geometries. For example, the generation of the ArcticDEM using multiple imagery data acquired at different times can introduce distortions, which makes it challenging to precisely quantify the propagation of this effect in the misregistration. Additionally, pixel misalignment could also be influenced by atmospheric distortions in optical imagery. Given that these pixel misregistration artifacts are mostly observed in a small subset of pixels with relatively large slope angles, we did not develop a misalignment correction algorithm in this study. Nonetheless, our approach provides a method to estimate spatial characteristics and upper bound of InSAR phase errors due to DEM-SAR pixel misregistration in individual interferograms.
InSAR-estimated seasonal surface thaw subsidence can be related to the amount of water stored in the saturated soil active layer above permafrost, which can be used to constrain hydrologic models and water mass budgets. In the Toolik area, 95% of z water estimates range from 0 to 62 cm, and the spatial distribution of z water correlates with elevation and vegetation cover types.
The amount of error in InSAR-estimated z water is linearly proportional to the error in InSAR LOS deformation measurements.
Although most InSAR measurement noises have been mitigated during the processing procedure, errors in the Arctic-DEM data and DEM-SAR misregistration can lead to visible InSAR LOS measurement errors. In the ALOS Toolik case, a 1-2 meter error in the Arctic-DEM data can lead to a LOS error larger than 1.5 cm when the perpendicular baseline is larger than 4000 m.
Errors associated with the DEM-SAR misregistration are determined by the amount of pixel misregistration, the local slope, and InSAR perpendicular baselines. In the Toolik area, these pixel-mismatching artifacts are mainly observed in regions with a steeper slope (> 10%) in interferograms formed using a subset of SAR scenes with noticeable misregistration issues. Most pixels in our study area have percent slopes smaller than 5%, and the LOS measurement error is generally smaller than 1 cm (equivalent to z water errors smaller than 14 cm.). As the landscape near Toolik Lake on the North Slope of Alaska transitions from hilly terrain to the south to flat plains to the north, DEM-SAR misregistration no longer produces visible phase artifacts in InSAR LOS observations. Our study shows that InSAR is an effective and powerful technique for accurately monitoring the status of and changes in hydrological characteristics in active-layer soils above continuous permafrost. InSAR estimates of soil water depth are statistically consistent with in situ observations, and the advantages of InSAR estimates include broader spatial coverage, higher spatial resolution, and the ability to map spatial patterns. Author contributions. Y.W. was responsible for writing the original draft, while the review and editing process involved contributions from
interferograms are similar (∼ 1500 meters), the observed errors are different. These artifacts were observed in many Toolik