Excess Ground Ice Profiles in Continuous Permafrost Mapped From InSAR Subsidence

Excess ground ice formation and melt drive surface heave and settlement, and are critical components of the water balance in Arctic soils. Despite the importance of excess ice for the geomorphology, hydrology and biogeochemistry of permafrost landscapes, we lack fine‐scale estimates of excess ice profiles. Here, we introduce a Bayesian inversion method based on remotely sensed subsidence. It retrieves near‐surface excess ice profiles by probing the ice content at increasing depths as the thaw front deepens over summer. Ice profiles estimated from Sentinel‐1 interferometric synthetic aperture radar (InSAR) subsidence observations at 80 m resolution were spatially associated with the surficial geology in two Alaskan regions. In most geological units, the estimated profiles were ice poor in the central and, to a lesser extent, the upper active layer. In a warm summer, units with ice‐rich permafrost had elevated inferred ice contents at the base of the active layer and the (previous years') upper permafrost. The posterior uncertainty and accuracy varied with depth. In simulations, they were best (≲0.1) in the central active layer, deteriorating (≳0.2) toward the surface and permafrost. At two sites in the Brooks Foothills, Alaska, the estimates compared favorably to coring‐derived profiles down to 35 cm, while the increase in excess ice below the long‐term active layer thickness of 40 cm was only reproduced in a warm year. Pan‐Arctic InSAR observations enable novel observational constraints on the susceptibility of permafrost landscapes to terrain instability and on the controls, drivers and consequences of ground ice formation and loss.

. The susceptibility to near-term instability is strongly controlled by the excess ice characteristics of the upper permafrost (Jorgenson et al., 2015).In addition to damaging infrastructure, permafrost-thaw-induced deformation can radically alter the movement and storage of water, ecosystem composition and the cycling of carbon and nutrients (Jorgenson et al., 2010;C. Ping et al., 2015).
The impacts of seasonal excess ground ice in the active layer are more subtle but no less pervasive (Mackay, 1984;Romanovsky et al., 2008).Seasonal heave and subsidence of up to ∼15 cm, driven by the formation and destruction of excess ice, can be detrimental to infrastructure and vegetation (Palmer & Williams, 2003; C. L. Ping et al., 2008;Walker et al., 2008).During freeze up, the movement of water toward growing ice lenses near the surface and the active layer-permafrost boundary (assuming no talik) promotes formation of a desiccated zone in between (C.L. Ping et al., 2008).These processes are implicated in solute and organic matter redistribution; changes in hydraulic and mechanical properties; steep gradients in biogeochemical conditions; and elevated perennial excess ice contents in the upper permafrost (Mackay, 1983;Shur et al., 2005; C. L. Ping et al., 2008).
The abundance of excess ground ice varies greatly in continuous permafrost landscapes (Jorgenson et al., 2010;Morse et al., 2009).The large spatiotemporal variability of seasonal excess ice is poorly constrained and difficult to model, owing to complex interactions between factors such as soil texture, water availability, and meteorological forcing (Gruber, 2020;Nicolsky et al., 2008;Rempel, 2010).Perennial ground ice in permafrost below the active layer varies vertically and horizontally.Ice wedges, a common form of massive ice formed by repeated cracking and water infiltration, illustrate the small-scale (∼10 0 -10 1 m) horizontal variability, with excess ice contents varying by an order of magnitude across polygonal networks.On landscape scales (∼10 2 -10 4 m), the large spatial variability of ice in the upper permafrost is the result of an intricate balance of ice formation and melt over decades to millennia, influenced by such factors as drainage conditions, topography, and disturbance (French & Shur, 2010;Shur et al., 2005).
There is a need for spatially resolved estimates of near-surface excess ground ice profiles on regional scales, but existing methods have limited spatial coverage.Coring is the most reliable method but expensive and spatially restricted (Kanevskiy et al., 2013;Morse et al., 2009).Local-scale geophysical observations such as electrical resistivity techniques also provide valuable constraints (De Pascale et al., 2008;Ross et al., 2007;Yoshikawa et al., 2006).Alternatively, surface displacement observations of the cumulative seasonal heave or subsidence are proxies for the excess ice content integrated over the given year's active layer (Gruber, 2020;Harris et al., 2011;Little et al., 2003;Romanovsky et al., 2008).Profile information can be obtained from temporally resolved observations (Gruber, 2020;Harris et al., 2011;Overduin & Kane, 2006;Smith, 1985), albeit hitherto mostly qualitatively.The recent availability of pan-Arctic subsidence time series from satellite interferometric synthetic aperture radar (InSAR) enables spatially resolved data products, but we lack a method and quantitative assessments.

Contributions
Here, we introduce a method to estimate near-surface excess ice profiles from InSAR subsidence time series from the thaw period (Figure 1).The central tenet is that the instantaneous subsidence rate is proportional to the excess ice content at the depth of the thaw front (Harris et al., 2011;Overduin & Kane, 2006).The idea, then, is to probe the ice content at increasing depths as the thaw front deepens over summer (Zwieback & Meyer, 2021).
We pursue the following objectives: 1. Devise a Bayesian excess ice profile estimation method 2. Conduct a synthetic performance assessment 3. Estimate and assess ice profiles from InSAR We meet the first objective by inverting a forward model that simulates subsidence for large ensembles spanning a range of plausible soil (and excess ice) stratigraphies (Figure 2).The Bayesian inversion combines the simulation results with the InSAR subsidence observations to yield posterior distributions of excess ice profiles.It accounts Figure 1.The idea is to relate the progression of subsidence over the thaw period to the excess ice profile.For instance, at the time of low subsidence rate indicated by the brown line on the left, the thaw front (brown line on the right) is progressing through materials with very little excess ice.Later, the subsidence rate increases when the thaw front progresses down to ice-rich strata.
10.1029/2023WR035331 3 of 19 for the direct and the thaw-front-mediated dependence of subsidence on the excess ice profile.
For the second objective, we apply the inversion to simulated InSAR observations.Comparison to the input excess ice profiles quantifies the location, width and fidelity of the posterior distribution as a function of depth and observational accuracy.We additionally perform a resolution analysis to determine how accurately ice-enriched layers at variable depths can be reconstructed.
For the third objective, we map excess ice profiles over two study regions from Sentinel-1 InSAR observations.We assess the results by comparison to excess ice profiles determined independently from cores at two sites within the first region (Brooks Range Foothills).We also quantify the year-to-year variability, which in the future may constrain hydrological and biogeochemical processes.In Northwestern Alaska, our second study region, we estimate excess ice contents at the top of (the previous years') upper permafrost in a very warm year when the thaw front penetrated materials that had previously been perennially frozen.This second analysis highlights the method's potential for mapping the susceptibility to terrain instability triggered by permafrost thaw.

Forward Model
Our forward model couples a thaw subsidence module with a lumped energy balance module.

Thaw Subsidence
We consider surface subsidence prompted by the melting of subsurface excess ice, with consolidation ensuing from meltwater drainage (Morgenstern & Nixon, 1971).
In our one-dimensional model, the instantaneous rate of surface subsidence ds(t)/dt equals the thickness of melted excess ice per unit time (Figure 3).Excess is assumed to melt at the depth of the thaw front, y f (t), so that where e(y) is the volumetric excess ice content profile.The fraction of excess ice volume within a volume of frozen soil, e is dimensionless but may be reported in m 3 m −3 .
The cumulative subsidence from the beginning of the thaw period t = 0 to time t, s(t), equals the total amount of excess ice above the position of the thaw front at time t The time t = 0 is taken to be the start of the thaw period y f (0) = 0.
The y coordinates in Equation 2 refer to the depth below the surface before the onset of thawing.Owing to the subsidence from thaw consolidation, the depth of the thaw front with respect to the surface at time t, ν(y f (t)), will be less than y f (t).

Assumptions
Equations 1 and 2 rest on simplifying assumptions.The model assumes that subsidence arises exclusively from the expulsion of excess meltwater (Morgenstern & Nixon, 1971).This mechanism corresponds to classical primary consolidation, according to which compaction results from the dissipation of excess pore pressure, typically induced by an external load (Biot, 1941).Our model implicitly postulates a soil medium of incompressible particles and fluids, whose void space decreases from the pre-thaw to the post-consolidation stage by the excess ice content.It does not account for the dependence on the effective stress (history) or the minor volume changes of quasi-saturated ground arising from melt-induced porewater volume expansion (Dagesse, 2010).It further neglects deviations from poroelastic compaction such as creep (Fowler & Noon, 1999;Liu & Borja, 2022).
Excess porewater pressures induced by meltwater release are assumed to dissipate faster than the subweekly time scales of interest, τ.Within Terzaghi's one-dimensional consolidation theory, we require Δν 2 /c v ≪ τ, where Δν is the distance from the pressurized stratum to unsaturated materials or the surface (Morgenstern & Nixon, 1971).
The consolidation coefficient c v scales with the hydraulic conductivity (Biot, 1941), rendering fine-grained materials most susceptible to delayed consolidation.However, we caution that thaw consolidation under natural conditions at ∼100 m scales is poorly understood owing to such phenomena as cryoturbated organic pockets (Figure 8c) and macropores.
The one-dimensional model does not capture lateral variability in ground ice contents, thaw depths, and subsidence (Jorgenson et al., 2015).Neither can it capture complex three-dimensional stress fields (e.g., cavity formation) (Osterkamp et al., 2000).

Modeling the Thaw Front
Linking the observed subsidence as a function of time to the sought excess ice as a function of depth necessitates knowledge of thaw front progression.

Lumped Energy Balance
We describe the thaw front progression using a lumped model for the energy balance of the ground above the thaw front.We assume that the thaw front is at the bulk melting temperature T m = 273.15K, the temperature at which the excess ice melts.The freezing point depression (Romanovsky & Osterkamp, 1997) of the interstitial ice is not considered.
We consider four terms in the energy balance (Figure 3).The surface heat flux Q s drives the thawing and warming of the ground.The major energy sink is the melting of ice at the thaw front y f .Denoting the (extensive) phase-change enthalpy of ground ice by H i , its time derivative is proportional to the rate of thaw front deepening and the (intensive) volumetric enthalpy of phase change L, which increases with the ice content (Lunardini, 1991;Romanovsky & Osterkamp, 1997).We also include minor terms in the energy balance, namely the heat flux into the frozen ground, Q f , and the change in enthalpy of the thawed ground, H t .The energy balance then reads A simplified energy balance will be used for demonstration purposes and to rapidly obtain initial estimates.As in the Stefan solution (Lunardini, 1991), the energy required to warm the thawed ground and the heat flux into the frozen substrate are neglected.Thus, The accuracy is expected to increase as the initial temperature of the frozen substrate approaches the melting point T m , and as the Stefan number  St = c (s − m)∕ approaches zero.This is the case when phase change is the main heat sink, as then the volumetric heat capacity    times the deviation of the surface temperature T s from T m is negligible compared to L.

Energy Balance Terms
The surface heat flux is modeled using a quasi-steady approximation of Fourier's law of conduction through (Lunardini, 1991) where k (y f ) is the effective thermal conductivity of the thawed ground, and the surface temperature T s (t) is modeled on daily time scales using an n-factor approach (Cable et al., 2016;Klene et al., 2001).The denominator is the depth of the thaw front relative to the subsided surface.The assumed quasi-steady temperature profile will be appropriate when the diffusion length in the thawed materials exceeds y f for time scales in excess of one day.
The full model accounts for temperature changes in the thawed ground by tracking  t( f ) = ∫  f 0 ct() ( )d , where   t is the thawed heat capacity of the consolidated soil per unit volume of frozen ground (Zwieback et al., 2019).A steady-state temperature profile T (t, y) is assumed.
The heat flux into the frozen ground, Q f , is estimated using a heat-balance integral approach (Goodman, 1958;Lunardini, 1991).It postulates a quadratic temperature profile between the thaw front y f and the temperature penetration front y p .At y p (t), the temperature equals the initial pre-thaw temperature T 0 < T m with vanishing depth derivative for all t.The heat flux Q f is given by Fourier's law just below y f , assuming a frozen thermal conductivity k f .A detailed description can be found in Section 1.

Numerical Solution
Our fast algorithm retrieves y f (t) and thus s(t) at daily time steps using a perturbation approach.We first solve the simplified energy balance based on Equation 4b to obtain a first-guess y f (t).We then refine y f (t) iteratively by estimating the minor energy sinks from the previous y f (t) and subtracting them in Equation 4a.
The first-guess y f (t) is obtained from integrating the separable ordinary differential equation that results from using Equation 5 in Equation 4b, viz.
Root finding yields y f (t) for each day.
The iterative refinement evaluates the minor energy balance terms in Equation 4a using the previous y f (t) estimate.It then adjusts y f (t) by subtracting the minor terms in Equation 7. A single refinement step was used for our study area owing to the large water content and not too cold permafrost.

Inversion
We first characterize the dependence of the subsidence time series s(t) on the excess ice profile e(y), before introducing our Bayesian inversion for estimating e(y) from subseasonal InSAR subsidence observations.

Sensitivity From Linearization
The subsidence predicted from Equation 2 depends on the excess ice profile e(y) directly and also through y f (t).The Gateaux variation, reflects a balance between (a) the direct influence ∫δe(y) and (b) a thaw-front mediated term e (y f )δy f .The direct influence (a) manifests as larger subsidence for increased excess ice above the (fixed) y f .However, a concomitant δy f < 0 can reduce the total subsidence through (b).We assume that δe(y) induces δL(y) = λe(y).The numerical value λ, a volumetric enthalpy of phase change, will depend on what constituent δe "replaces." The sensitivity of s(t) to e(y) for the simplified energy balance of Equation 4b and assuming k(y) = k can then be described by where the integral kernels K e (y) and   f () capture the direct and thaw-front-mediated sensitivity, respectively.The dependence of the kernels on t is omitted.The direct influence kernel K e (y) = 1 states that the total change in excess ice above y f (t) induces an equal change in subsidence, for constant y f .Using scales with e (y f ).It consists of two contributions.The first term in the parentheses is nonpositive and proportional to λ: the energy expended to melt additional excess ice, δe > 0, retards the thaw front, thus decreasing s(t).The second term is nonnegative: additional excess ice induces more subsidence, thus keeping icy materials closer to the (subsided) surface, which in turn enhances subsidence (integration by parts makes the physical link more readily apparent).

Simulated Sensitivity
Simulations illustrate the dependence of the predicted subsidence on excess ice profiles.The predictions were obtained using the model from Section 2.3, rather than the simplified energy balance, for a hypothetical uniform excess ground ice profile of e = 0.05.The two sensitivity scenarios had a 10 cm layer centered at 20 and 45 cm, respectively, enriched by δe = 0.1.
The predicted subsidence accelerated when the thaw front hit the enriched layers (Figure 4a).The subsidence rate ds/dt increased relative to the baseline, resulting in greater subsequent cumulative subsidence s(t).After melt of the enriched layer, ds/dt was similar to the baseline, yielding a quasi-parallel s(t) trajectory.The legacy of the enriched layer is visible in the thaw front position y f (Figures 4b and 4c), which initially lagged that of the baseline due to the increased energy sink but (partially) recovered.

Bayesian Inversion of InSAR Observations
We propose a fast approach to estimate e(y) from InSAR, exploiting the dependence of subsidence on e(y).It takes InSAR subsidence observations  d and produces a Bayesian posterior distribution for e(y).To represent the uncertainty of parameters such as the thermal conductivity, an ensemble of stratigraphies, including e(y) is assembled.The forward model applied to each ensemble member yields a subsidence time series s(t).By weighting the ensemble members according to how compatible their s(t) is with the uncertain InSAR observations, an ensemble-based posterior distribution of e(y) is obtained.The method is fast because the same forward model runs can be applied across many similar pixels.
The posterior distribution captures the uncertainty in e(y) in dependence of the error covariance of  d and the parameteric uncertainty in stratigraphies.Through use of the forward model, it naturally accounts for the direct and the thaw-front-mediated dependence of subsidence on excess ice.Conversely (Overduin & Kane, 2006), estimated ice content profiles from in-situ subsidence observations by prescribing the thaw front progression using in-situ data, obviating the need to account for the coupling between e(y) and y f (t).

Bayesian Inference
We evaluated the posterior distributions of e(y) (and also y f (t)) conditional on the InSAR subsidence observations using importance sampling.More precisely, we obtained an approximate sample-based estimate of the posterior density (Vehtari et al., 2019).The procedure requires a prior ensemble of N stratigraphies β n , where β n encompasses all the parameters in the model.We can think of a prior β n as a plausible configuration of the excess ice profile e(y) but also of parameters like k(y) or c(y).These latter depend on variables such as saturation and organic layer thickness; here, we propose simple parameterizations and prior distributions introduced in Section 3.2.2.The estimation also requires temperature forcing data.
The staring point for the importance sampling is the prior ensemble (β 1 , …, β N ), from which the predicted subsidence time series s n (t) for each β n was obtained using the forward model.The method then evaluated for each candidate ensemble member how consistent its s n (t) was with the time-lapse InSAR observation  d .The consistency was quantified using the likelihood of β n ,   (   ; d) , accounting for the characteristics of the observations (to be introduced in Section 3.2.3).Each ensemble member was assigned an unnormalized weight , thus emphasizing those whose subsidence matched the observations (Hesterberg, 1995).Formally, this corresponds to importance sampling using the prior distribution as proposal, similar to the original particle filter (Gordon et al., 1993).Importance sampling can suffer from large variance induced by a handful of samples being assigned exceptionally large weights.To reduce the variance, we used Pareto smoothing (Vehtari et al., 2019), yielding adjusted final weights (w 1 , …, w N ).The posterior density of any variable θ-for instance e at depths y 1 and y 2 -could then be approximated by Rather than dealing with the Dirac delta δ directly, we computed expectations such as means or probabilities.

Prior Ensemble of Stratigraphies
For our study areas, we modeled stratigraphies on a regular grid to a depth of y max = 1.5 m with a grid spacing Δy = 2 mm.
The prior profiles of excess ice before thaw onset were represented by smoothing splines with 12 knots.The knot spacing increased quadratically from 5 cm at the surface, attaining 15 cm at 0.5 m, a typical active layer depth in our region.The logit of the excess ice contents at the knots, l e , was drawn from an autoregressive normal distribution with mean    = −3.0 , standard deviation    = 3.0 , and auto-correlation    = 0.7 .The knot spacing and auto-correlation were selected to achieve profiles whose smoothness reflects expectations on the InSAR scale of ∼80 m and stabilizes the estimation using time-lapse observations.The logits l e were interpolated between the knots and then converted to e profiles.The parameter values yield average e of 0.17 with a standard deviation of 0.25, thus not imposing strong prior assumptions.
The thermal properties were modeled based on the stratigraphy.A detailed description can be found in in Section 2. To summarize, the volume not occupied by excess ice at each depth was allocated to the constituents mineral, organic, non-excess ice and air.The organic layer thickness was modeled with a uniform distribution, as were the saturation of the organic and mineral layers.From the volumetric fractions of all constituents, the thawed thermal conductivity was modeled following (Cosenza et al., 2003).Finally, the n-factor was drawn from a scaled beta distribution.The distributions of these quantities were assumed to be mutually independent.

Stochastic Model for InSAR Observations
Let the InSAR subsidence observations be denoted by  d , a P-dimensional vector derived from P + 1 time-lapse InSAR scenes acquired over a thaw season.The element  d encodes the displacement at t i , the acquisition time of scene i + 1, relative to the time t 0 , the first scene of the season.The displacement is measured along the line of sight to the satellite rather than along the vertical.Assuming the actual displacement to be along the vertical, the observation relates to the predicted s(t) as where θ is the incidence angle and ɛ i the observation error for scene i.
We modeled the observation error using a normal distribution with autocorrelated errors.The InSAR displacement errors had zero mean, as we did not consider systematic errors from soil moisture or vegetation changes (Zwieback & Hajnsek, 2016).Their covariance matrix Σ = Σ s + Σ a was partitioned into speckle and atmospheric errors, respectively (Hanssen, 2001).Speckle induces quasi-random fluctuations in the InSAR phase (Zwieback & Meyer, 2021).Their magnitude increases the lower the interferometric coherence, a measure of the temporal stability.The errors propagate to the displacement estimates, often inducing substantial temporal autocorrelation.Here, we estimated the covariance matrix Σ s using the full Fisher-Information (Zwieback & Meyer, 2022) evaluated using the regularized InSAR magnitude matrix (Zwieback, 2022).Second, atmospheric errors affect the relationship between the observed phases and inferred displacement.They are spatially correlated.We modeled the atmospheric error covariance matrix Σ a based on temporally uncorrelated errors, one for each scene (Hanssen, 2001).Thus, the log-likelihood of parameters β giving rise to a predicted d(β) was described by where |⋅| denotes the determinant and c is a constant of no consequence in normalized importance sampling.

Rationale
To explore the estimation in a controlled setting, we conducted a synthetic twin experiment and a resolution analysis.
The twin experiment provided constraints on the real-world estimation performance.In contrast to real-world conditions, the twin experiment was based on the premises that the forward model provided a perfect representation of nature and that the uncertainty in the stratigraphy parameters was accurately captured.For instance, even though the organic layer thickness and thus k(y) was not assumed to be known in the inversion, its variability as captured by the ensemble matched the variability of the synthetic truth.To isolate errors associated with the stratigraphy and InSAR observations, we assumed the temperature forcing was known.
A complementary resolution analysis was conducted to determine the depth resolution achievable in the presence of parametric and observational uncertainty.The goal was to identify how the resolution varies with depth and illustrate how parametric uncertainty in the stratigraphy relates to a loss of resolution.

Synthetic Twin Experiment
In a first step, we simulated the synthetic truth, that is, the subsidence time series s(t), for prescribed stratigraphies drawn from the prior distribution.Second, we generated synthetic noisy InSAR observations  d .Third, we applied the Bayesian inversion to these observations to estimate e(y) and compared it to the prescribed profile.

Synthetic Truth
The synthetic truth was obtained by running the forward model from Section 2.3 for N sim = 500 stratigraphies β n drawn independently from the prior distribution (Table S1 in Supporting Information S1).No additional synthetic noise was applied to the predicted subsidence time series.

Synthetic Observations
Synthetic InSAR observations  d were generated from Equation 13for an incidence angle θ of 30° and regular 12-day sampling.For each of the N sim stratigraphies, we generated N rep = 100 sets of replicates with independent and identically distributed errors ɛ.
We considered three scenarios for the observation error ɛ.First, a baseline scenario based on the actually observed covariance Σ at our Ice Cut (IC) site (see Section 5.1).Second, a low-accuracy scenario where Σ was scaled by 16.Third, a high-accuracy scenario where Σ was 1/16th of the baseline.These scenarios are intended to capture low coherence and high coherence, respectively.

Assessment
We compared the posterior distribution of the discretized e(y) to the prescribed e(y) for a given y.We retrieved for each synthetic observation  d using the same atmospheric forcing as in the forward model.
The error of the posterior mean of e(y) was quantified through the mean absolute deviation (MAD) with respect to the truth.
The posterior uncertainty was characterized through its width and calibration.The width was summarized by the credible interval I 10−90 bounded by the 10th and 90th percentile of the posterior  p(()) .The decrease in size, expressed through the half width HW(I 10−90 ), with respect to the corresponding prior interval is an internal measure of the insight gained from the observations.The interval's independence of the truth implies that a small half-width may be of little value if the posterior is not calibrated (Gneiting et al., 2007).The calibration of the posterior distribution was evaluated using the coverage C(I 10−90 ) of the 10%-90% credible interval.Theory predicts a coverage of 80%, as the simulated truth was drawn from the prior distribution (Dawid, 1982).Deviations from 80% are diagnostic of issues with the approximate Bayesian inference.

Resolution Analysis
The resolution analysis considered the same scenarios with ice-enriched layers as the sensitivity analysis in Section 3.1.While the abrupt breaks in excess ice profiles are difficult to reconcile with the smooth prior, they illustrate how accurately the depth and thickness of ice-enriched layers can be estimated.Forward modeling, generation of synthetic InSAR observations and the ensemble-based inversion were identical to the synthetic twin experiment.
Analysis of the ensemble members with high posterior weight sheds light on the sources of uncertainty.By comparing the ice profiles with the depth of thaw across ensemble members, we can relate the depth resolution to the parametric uncertainty of the soil stratigraphy.Consider two ensemble members with contrasting thaw depth trajectories that are equally compatible with the uncertain synthetic observations according to Equation 14.For instance, one may have shallower thaw due to a thicker organic layer.We expect its ice-enriched zone to be shallower as well.We visualized the depth variability by plotting the e profile of those members most compatible with the synthetic observations, along with their thaw depth at the first and last InSAR observation.

Synthetic Twin Examples
Synthetic inversion results for two contrasting stratigraphies are shown in Figure 5.One is ice-rich only near the surface (a), while the second features an abrupt increase in excess ice at 40 cm (b).The associated subsidence trajectories are reconstructed with fidelity by the Bayesian inversion of the InSAR subsidence.
The accuracy of the inverted excess ice profiles varied with depth.In the central active layer (15-30 cm), the posterior mean followed the prescribed with deviations of ≲0.05.The deviations were greater near the surface, most strongly so for "unprobed" layers that thawed prior to the first InSAR scene (Figures 5b and d).The accuracy also decreased toward the base of the active layer, near ∼40 cm.
The posterior uncertainty, as quantified by the 10%-90% credible interval in Figure 5, was largest near the surface or near the base of the active layer, especially when these strata were ice rich (panel b).

Synthetic Twin Assessment
Across all stratigraphies in the synthetic twin experiment, Figure 6a shows the accuracy of the posterior mean from N = 10 5 ensemble members was best in the central active layer.Between 15 and 35 cm, the MAD was around 0.05.It increased in either direction.For depths exceeding ∼60 cm, it approached the MAD of the prior, as these layers remained frozen throughout.Throughout the profile, the accuracy was virtually identical when only 10 4 ensemble members are used (Figure S1 in Supporting Information S1).
The accuracy of the InSAR observations had a moderate impact on the estimation performance as measured by the MAD.The MAD improved by up to 0.02 in the high-accuracy scenario, for which the variance of the observation errors was reduced by a factor of four.In the low-accuracy scenario (quadrupled variance), the MAD increased by up to ∼0.05.For both scenarios, the changes with respect to baseline were greatest in the central active layer.
The estimated uncertainty was likewise smallest in the central active layer (Figure 6c), with an average width of the 10%-90% credible interval of around 0.06.It approached the prior one of ∼0.2 at the surface and below 60 cm.The sensitivity to the observational accuracy was comparable to that of the MAD.
The 10%-90% credible intervals are shown to be well calibrated in Figure 6c.The uncertainty intervals were slightly too wide in the active layer, with a coverage of 0-5 percentage points below the 80% target.The deviations were greater for the high-accuracy scenario and for N = 10 4 ensemble members (Figure S1 in Supporting Information S1).

Resolution Analysis
Figure 4d shows the inversion results for the baseline profile.The posterior mean matched the prescribed truth between 10 and 50 cm depth.It deviated from the prescribed profile in the upper ∼8 cm, which largely thawed prior to the first observation on May 21 for the ensemble members (dark markers).
The shallow ice-enriched layer was identified in the posterior mean (Figure 4e), but at reduced resolution.The high-weight ensemble members featured an enriched layer of ∼10 cm thickness, but the peak location varied by ∼15 cm.Averaging over the ensemble members degraded the resolution, similar to the example in Figure 5a.
The deep ice-enriched layer was reconstructed with substantial uncertainty (Figure 4f), with elevated posterior means ≳10 cm above the layer (cf. Figure 5b).The highest-weight ensemble members captured the increase in excess ice content, albeit with substantial depth dispersion.Those ensemble members with shallow (deep) ice-enriched layers also had shallow (deep) thaw at the last SAR observation (light marker).The depth dispersion induced a loss of resolution when averaging across ensemble members with disparate stratigraphies.The ice-poor materials beneath the enriched layers were not resolved, as they thawed after the last observation in most ensemble members.

InSAR Inversion
Our Sentinel-1 InSAR analyses focused on two regions in the continuous permafrost zone of Alaska (Figure 7a).

Brooks Foothills
Our first study region is located along the Dalton Highway, a rolling tundra landscape between the Brooks Range and the Beaufort Sea (Figure 7c).Mean annual air temperatures of ∼−7°C sustain continuous permafrost.Sentinel-1 observations from path 131, frame 363 were made every 12 days in the summers of 2019 and 2022.We focused on two sites, Happy Valley (HV) and IC.Each contains a 100 m × 100 m plot at which we obtained ground ice and thaw depth data (Section 5.2.2).Their vegetation cover includes sedges, most notably the tussock-forming Eriophorum vaginatum, and dwarf shrubs (Figure 7d).The soils at both are silty loams with some larger clasts.The gently sloping HV site is moderately well drained due to its crest-proximal location, and it does not contain recognizable ice wedge troughs.The IC site is on a flat plateau; while no surficial signs of ice wedges are apparent within the plot, disturbed locations nearby feature well-developed troughs.

Kivalina, Northwestern Alaska
Our second study region is near the town of Kivalina in northwestern Alaska.The tundra landscape is underlain by warm continuous permafrost.We focused on the Tatchim Isua location, a candidate site for relocating the flood-prone town (Tryck Nyman Hayes, 2006).Geotechnical analyses conducted in 2005 revealed that a narrow rocky bench was ice poor whereas the upper permafrost of the surrounding hillslope was rich in segregated and massive ground ice, limiting the location's suitability for relocation (Shannon & Wilson, Inc, 2006).
We studied Sentinel-1 data from the extremely warm summer of 2019, when thawing degree days exceeded the long-term average by more than 50%.A strong association between elevated late-season subsidence and independent estimates of the ice content at the top of the permafrost was previously documented by Zwieback and Meyer (2021) across the wider landscape.Here, we illustrate the profile inversion from Sentinel-1 (path: 15, frame 367) for the Tatchim Isua location, with the goal of generating a quantitative ground ice map.The absence of in-situ excess ice measurements precludes a quantitative assessment.

Forward and Inverse Modeling
The forward modeling yielded an ensemble of N = 10,000 predicted displacements d at the Sentinel-1 observation times.The daily temperature forcing for the Brooks Foothills was taken from NOAA station USS0048U01S, located 30 km north of HV.For Kivalina, we used the MERRA-2 re-analysis product (Gelaro et al., 2017).
The InSAR displacements  d were derived in four steps.First, we multilooked the single-look complex observations to obtain the interferometric covariance matrix at a ground-equivalent resolution of 80 m.Second, from each pixel's covariance matrix we estimated the phase history using the Hadamard regularization from Zwieback (2022).Third, the phase history was unwrapped sequentially with SNAPHU (Chen & Zebker, 2001) and converted to an equivalent displacement.Fourth, we locally referenced the displacement to a location (e.g., rocky outrcrop) assumed to be stable.From the second step, we also obtained an estimate of the covariance due to speckle, Σ speckle .An estimate of the total covariance Σ was obtained by addition of atmospheric errors.Our study areas' spatial extent of only ∼2 km implies small atmospheric errors.We assumed a displacement-equivalent standard deviation of 4 mm, resulting from the difference in the optical path lengths through the atmosphere.
The inversion used the importance sampling from Section 3.2.Each InSAR pixel was treated independently.Only one ensemble was used for all pixels, with stratigraphies generated using parameters from Table S1 in Supporting Information S1.While the parameter values are plausible for the two Brooks Foothills sites with coring data, their applicability to other locations such as recent floodplains or the rocky bench is questionable.For instance, underestimation of the thaw depth at sites with thinner organic layers will compress the inverted profiles.

Comparison With Cores
The quantitative assessment was focused on the two locations in the Brooks Foothills.We compared the InSAR estimates with the site-averaged excess ice profiles derived from cores.
We collected cores of the frozen active layer and the upper permafrost from May 02 to 05, 2022.At the HV and IC sites, we took 11 cores each within the 100 by 100 m transect area.Using a 7.5 cm diameter SIPRE corer, we drilled down to depths of 90-150 cm but only show results for the top 60 cm.Two example cores are shown Figure 8.In core HV-C1 (subfigure a), visible ground ice is found primarily at the top of the mineral soil (12-20 cm; dark) and below 45 cm, where it appears transparent, gray or dark.The previous year's active layer thickness as measured from the surface position in spring is unknown but was presumably ∼45 cm.Conversely, the central active layer (25-40 cm) is characterized by a massive soil structure without visible excess ice.The core IC-D4 in Figure 8c also features an ice-poor central active layer and ice-rich strata below 40 cm, which are, however, interrupted by organic soil.
We determined each core's excess ice profile in the laboratory.The profiles were obtained in 5-10 cm increments by sawing the frozen core into segments (Figure 8a).The segment length varied because the cores already fragmented while drilling.For each segment, we first determined the frozen volume V f by submersion in water.The sample was subsequently thawed in a beaker, homogenized and allowed to settle (Figures 8b-8d and 8h).Following Morse et al. (2009), we determined the volume of supernatant water V s (by weighing) and estimated the excess ice content as  = (ws)∕(if ). (15) The estimates tend to be conservative because air inclusions in ground ice (Morse et al., 2009) and uptake of excess ice meltwater by unsaturated mineral soil reduce V s .Wedge ice-encountered in five IC cores-was assumed to consist entirely of excess ice.
Site-level e profiles were determined on a 1 cm grid from the arithmetic mean across cores.An 80% confidence interval was derived from percentile bootstrapping using 1,000 samples, assuming independence between cores.
According to the site-level profiles at HV and IC (Figures 8e and 8f, respectively), average excess ice below the surface organics are low ≲0.10, decreasing to zero below 20 cm and gradually increase below 35 cm to values ≳0.25 at 50 cm.The observed variability was substantial at all depths with excess ice (10-20 cm; below 35 cm), in contrast to the central active layer.
We measured the 2022 late-season thaw depth at 150 transect points each at IC (August 15) and HV (August 17) by manual probing (Figure 7d).The thaw depths were measured from the subsided surface, thus corresponding to υy f .The average values were 41 and 40 cm, respectively.Thaw was shallow in 2022 compared to preceding years such as the warm summer of 2019 (Figure 8g), according to mid-August thaw depth measurements from the Circumpolar Active Layer Monitoring (CALM) 1 km grid (121 points) at HV (Nyland et al., 2021).
In the Kivalina study region, we qualitatively compared the excess ice estimates at 55-65 cm depth, an approximate active layer depth for vegetated locations (Shannon & Wilson, Inc, 2006), to a binary classification of the upper permafrost (ice-rich/ice-poor) based on 8 cores collected in 2005 (Shannon & Wilson, Inc, 2006;Zwieback & Meyer, 2021).

Comparison to Cores, Brooks Foothills
The 2022 coring excess profiles largely fell within the 10%-90% credible interval of the InSAR estimates, but the posterior mean deviated notably from the coring profile, especially in the lower active layer (Figures 9a-9d).Near the surface, the large uncertainty and concomitantly elevated mean mirrored the simulations.The ice-enriched layers near 15 cm, corresponding to subsidence of ∼5 mm, were not resolved by the inversion (Figures 9a-9d).
The agreement was better in the central (15-35 cm) active layer.At both sites, the coring estimates captured the low excess ice contents e < 0.05 in the central active layer (15-35 cm), with narrow (<0.1) uncertainty intervals.This corresponds to ∼0 cm InSAR subsidence from mid-June to late July (Figures 9c-9f).
The inversion underestimated the increase in e at the depth of the long-term active layer.The inferred probed depth, the estimated y f on September 14 (final Sentinel-1 acquisition), was 48 and 44 cm at IC and HV, respectively.The posterior mean e at 45 cm of ∼0.1 was smaller than the core-derived values of ∼0.2 at both sites.Furthermore, the core results show a steep increase of average excess ice content at and below this depth.Testing whether the low ice contents inferred from limited late-season subsidence of ∼1 cm were due to the thaw front barely penetrating the ice-rich materials would require field observations from the end of the thawing season.Instead, we were only able to assess the inferred probed depths in mid-August 2022.
At IC, the inversion-based inferred probed depth of y f = 41 cm was close to frost-probing-derived υ(y f ) of 42 cm.At HV, y f = 37 cm is to be compared to υ(y f ) of 40 and 38 cm from our transects and the CALM grid, respectively.
For 2019, the InSAR-derived e below 40 cm matched the coring data from 2022 reasonably well, capturing the increased ice contents at depth.The discrepancy in years limited quantitative comparisons of e in the active layer, but the shallow thaw in the intervening years (Figure 8g) suggests comparisons below ∼45 cm remained meaningful.Thaw was deep in 2019, with the mean mid-August υ(y f ) at the CALM grid of 46 cm exceeding the 2022 value by 8 cm.The inferred probed depth y f of 46 cm was comparable to the CALM observations on August 12, increasing to 51 cm at the day of the last Sentinel-1 acquisition.At both sites, the InSAR posterior mean profiles captured but smoothed the steep increase in ice contents below 40 cm depth.The identification of ice-enriched layers at depth reflected the substantial late-season subsidence of ∼4 cm observed at both sites in this warmer summer.

Ice Cut, Brooks Foothills
The estimated e profiles in 2022 varied little across the landscape, with consistently low (≲0.1)posterior means (Figure 10).Exceptions include the floodplain with pockets with elevated e ≳0.15 near the surface and also at the base of the active layer.These areas do not align neatly with apparent deposit age and vegetation cover.The interannual variability in estimated e profiles was most apparent where the upper permafrost is interpreted to be ice rich.In the warmer summer of 2019, large posterior means ≳0.3 at 50 cm prevailed across older and more ice-rich geological units, corresponding to ∼3 cm late-season subsidence (Figure S2 in Supporting Information S1).Conversely, low estimates at 50 cm were largely confined to (in)active floodplains.Elevated e estimates in the central active layer were associated with larger ice contents at depth, consistent with the limited depth resolution found in the simulations.At 10 cm, the inferred excess ice contents were even lower than in 2022.

Happy Valley, Brooks Foothills
In the cool summer of 2022, low InSAR-derived excess ice contents in the active layer predominate.Minor, and patchy, exceptions were found on the floodplain (near the surface and at ∼0.5 m) and in the uplands (at 0.5 m).
The warmer summer of 2019 yielded elevated InSAR-derived excess ice at depth across units where the upper permafrost is thought to be ice rich.Across transect T1, these include abandoned floodplain deposits and soil-covered hillslopes with visible ice wedge polygons.Conversely, low excess ice contents were inferred for inactive floodplains and rocky outcrops, corresponding to ∼0 cm late-season subsidence (Figure S3 in Supporting Information S1).

Tatchim Isua, Northwestern Alaska
Vast parts of the landscape exhibited large estimated excess ground ice at depths of ∼0.6 m.While the posterior means of the excess ice contents e(y) were near 0.1 at 10 cm depth and smaller at 25 cm (Figures 12a and 12b), values exceeding 0.4 were prevalent at 60 cm (Figure 12c) on the hillslopes and in the drained lake basin (lower left).The most notable exceptions with e ∼ 0 were rocky outcrops, including the bench, and the floodplain in the upper left corner, corresponding to ∼0 cm late-season subsidence (Figure S4 in Supporting Information S1).
Intermediate values of ∼0.2 were most conspicuous along the margins of the drained lake basin.
At the Tatchim Isua rocky bench (transect T1, Figure 12e), the InSAR-derived excess ice contents at ∼60 cm increased sharply up and downslope.The contrast in ice content agrees with the 2005 cores (Figure 12c), whose sparsity precluded an assessment of the sharpness of the transition.Upslope of the bench (550-750 m distance in T1), where field and high-resolution imagery observations did not detect conspicuous manifestations of ice-rich permafrost (Shannon & Wilson, Inc, 2006), the InSAR-derived e profiles indicate ice contents comparable to below the bench and further up the hillslope, where ice wedge troughs are recognizable.
The estimated ice contents of the upper and central active layer are generally lower.In particular, posterior means in the central active layer at 25-35 cm were close to zero in many areas (e.g., T1).Exceptions with elevated e(y) ≳ 0.2 above 30 cm were clustered in locations disposed to fine-grained soils and abundant water supply.For instance, the mature inactive floodplain at the start of transect T2 (Figure 12f) had large reconstructed ice contents near ∼20 cm but not at depth.

Performance Assessment
The inversion method provides spatially resolved estimates of near-surface excess ice profiles from satellite InSAR subsidence observations.It exploits the close physical connection between the ice content at the thaw front and the subsidence triggered by the melt of this ice.Its applicability to broad spatial scales is an advantage over coring or surface geophysical methods.The method largely yielded plausible excess ice profiles across landscape units.For instance, during the warm summer of 2019 it captured the desiccated central active layer underlain by ice-rich materials at our coring sites (Figure 9), while nearby ice-poor units such as rocky outcrops or (in)active floodplains had low estimated excess ice contents (Figure 11).
Our quantitative assessment revealed consistent estimates near the surface and greater deviations at depth.In the upper 40 cm, the inferred profiles largely overlapped with the coring-derived profiles from 2022.One shortcoming is the limited ability to capture the moderate ice contents in the upper mineral soil, but the uncertainty intervals of the InSAR and coring-derived excess ice profiles overlap.The uncertainty intervals were wide (∼0.5) near the surface and narrow (≲0.1) at 30 cm, similar to the simulations (see Figure 5).Below 40 cm, the inversion did not reproduce the increase in ice content at the base of the active layer and upper permafrost in the cool summer of 2022, owing to limited late-season subsidence from InSAR.The inconsistency between coring and InSAR-derived ice contents below 40 cm may partially be due to unresolved variability.For instance, the ice-rich layer in core IC-D4 (Figure 8c) may be closer to the surface than on average partially because thaw tends to be shallow, potentially confounding the comparison.Conversely, errors in y f are less plausible a reason because the internally estimated average depth of probing appeared reasonable (Figure 9).To conclusively tackle this challenge, co-located thaw depth and coring observations are required, but this is challenging because coring is destructive.The method identified ice-rich materials in the (previous years') upper permafrost in Northwestern Alaska and the Brooks Foothills in the warm year of 2019.At the Tatchim Issua bench in Northwestern Alaska, estimated e ≈ 0.4 at 60 cm corresponded to ice-rich permafrost and e ≈ 0.0 to ice-poor materials.While the comparison to (old) cores was only qualitative, the inferred landscape-scale variability was plausible.Similarly, in the two Brooks Foothills locations, large excess ice contents at 50 cm were inferred in 2019 across units known to be ice-rich (e.g., visible ice wedge troughs), and low excess ice contents at ice-poor locations such as recent floodplains.

Limitations and Future Improvements
Inferring excess ice contents from subsidence is intrinsically limited to the depth of thaw.The materials below y f at the last observation remain "unprobed."To infer how ice-rich the (long-term) upper permafrost is, a year with deep thaw is needed.It may also happen that even in a warm year, deeper ice-rich materials remain "unprobed" because they are too deep (for instance, due to previous disturbance).Conversely, when the thaw front sufficiently penetrates ice-rich materials at depth, simulated and InSAR inversions can detect these materials, albeit with reduced depth resolution (Figure 12).The reduced depth resolution is associated with the dispersion of predicted thaw depths (Figure 4), suggesting performance gains from tighter prior constraints on soil stratigraphies.
The subresolution variability in stratigraphies within the ∼80 m InSAR pixel is a major challenge for the method.
As the forward model is one-dimensional, it does not represent the ∼0.2 variability in e (Figure 8), the ∼25% dispersion in y f (Nyland et al., 2021), and the association between y f , ground ice and such factors as organic layer thickness.The ensemble inversion yields an effective e profile.Potential remedies would be ensembles of three-dimensional model runs or an ensemble of subresolution tiles, but at substantial computational cost (Abolt et al., 2018;Martin et al., 2021).
Extending the method to regional scales requires improvements to the prior stratigraphies and the InSAR error modeling.Improved stratigraphies should account for spatially variable soil characteristics, especially the organic layer thickness and n-factor.For instance, our default ensemble assumed too thick an organic layer on the rocky bench in Figure 12e, thus presumably underestimating the thaw depth; the impacts on the estimated ice profiles were minor because the bench did not subside notably.Regional analyses further require data-driven compensation of the atmospheric errors, which can reach ∼2 cm at 50-100 km (Hanssen, 2001).Finally, uncertainty in the temperature forcing is a challenge, in particular in high-relief terrain.

Conclusion
Based on synthetic and Sentinel-1 subsidence observations, we showed how satellite InSAR can map near-surface profiles of excess ice in favorable conditions.The InSAR estimates' variability across the landscape aligned with the surficial geology, with soils unconducive to excess ice formation generally exhibiting low ice contents.The estimated profiles predominantly had low ice contents in the central active layer and, to a lesser extent, the upper active layer.Elevated ice contents near the surface were largely restricted to floodplains.Inversion of InSAR observations from the very warm summer of 2019 yielded large ice contents ≳0.3 at the base of the active layer and the previous years' upper permafrost across geological units known to host ice-rich permafrost.
The accuracy and posterior uncertainty varied with depth.In simulations, they were best (≲0.1) in the central active layer, deteriorating (≳0.2) toward the surface and permafrost.In the Brooks Foothills, the estimates compared favorably to coring-derived excess ice profiles down to 35 cm, while the increase in excess ice contents near the base of the (average) active layer were reproduced in a warm but not a cool summer.In Northwestern Alaska, the extreme summer of 2019 with deep thaw facilitated detection of ice-rich materials in what previously had been the upper permafrost, as the retrievals matched qualitative coring observations.
The pan-Arctic availability of suitable InSAR data enables novel insights into water fluxes and ground ice storage changes.The rates, controls and drivers of ground ice loss and formation are poorly constrained but of critical importance to the stability, hydrology and biogeochemical functioning of permafrost landscapes.Empirical constraints from our InSAR method and complementary observations portend improved monitoring and predictive capabilities across the rapidly changing Arctic.

Figure 4 .
Figure 4. (a-c) Sensitivity of modeled subsidence to excess ice profiles.Compared to the baseline with uniform e = 0.05, the shallow and deep scenarios featured ice-enriched layers indicated by the shaded bands in panel (b).For the shallow layer shown in purple, subsidence in panel (a) increased in late May when the thaw front hit the enriched layer below 15 cm shown in panel (b), followed by a sustained reduction in thaw depth relative to baseline, plotted in panel (c).(d-f)Resolution analysis with synthetic observations.The prescribed excess ice profile is shaded, the posterior mean shown with a thick line, and the highest-weight ensemble members with thin lines, with markers to indicate y f at the time of the first and last observation.

Figure 6 .
Figure 6.Synthetic performance assessment for N = 10 5 ensemble members and three settings of the observational accuracy.(a) Mean absolute deviation of the posterior mean vs. the prescribed excess ice content; (b) average width of the posterior 10%-90% credible interval; (c) coverage of the 10%-90%, with the 80% target shown in gray.

Figure 5 .
Figure 5. Excess ice profile estimates for two synthetic stratigraphies (columns).(a-b) True excess ice profile in dark blue, posterior mean in brown, with the shaded region indicating the 10%-90% credible interval; estimated thaw front depth y f at time of last interferometric synthetic aperture radar (InSAR) scene is shown in gray.(c-d) Dark line shows the true subsidence simulated with the forward model with prescribed stratigraphy, the brown line the posterior mean obtained by Bayesian inversion based on synthetic InSAR observations converted to subsidence (brown circles).

Figure 7 .
Figure 7. (a) Map of Alaska; Landsat true-color composite of (b) the Kivalina, Northwestern Alaska and (c) the Brooks Foothills study region, respectively; (d) thaw depth measurements at Happy Valley, characterized by tussock and dwarf shrub tundra.

Figure 8 .
Figure 8. analysis results from Happy Valley (HV) and Ice Cut (IC), Brooks Foothills.(a, c) Photograph of the cleaned cores stitched together from pictures taken in the lab.Depth increases toward the right; the segments and their excess ice content e are indicated by the purple boxes (photograph c was made before the segments were cut).(b, d, h) Beakers containing the thawed samples, with supernatant water in panels (b, h) appearing yellow.(e, f) Site-level excess ice profiles (heavy line: mean, shaded area: 10-90% confidence interval, thin lines: individual cores).(g) Mid-August thaw depth at the HV circumpolar active layer monitoring (CALM) grid approximates the active layer thickness.

Figure 9 .
Figure 9.Comparison of interferometric synthetic aperture radar (InSAR) and coring profiles for Ice Cut (top row) and Happy Valley (bottom row) in the Brooks Foothills.(a, d) Coring profile from 2022 in blue (line: mean, shaded: 80% confidence interval) versus 2022 InSAR profile in brown (line: posterior mean, shaded: 10%-90% credible interval); the transparency of the InSAR profiles is proportional to the probability of ground remaining frozen at the time of the last InSAR acquisition, and the inferred probed depth y f (50% probability of remaining frozen) is marked with a horizontal line.(b, e) InSAR profiles from warmer summer 2019 vs. 2022 cores.(c, f) InSAR time series converted to vertical subsidence (marker: point estimate, vertical bar: ±standard deviation).

Figure 10 .
Figure 10.Ice Cut, Brooks Foothills: InSAR-derived excess ice profiles for 2022 (cool) and 2019 (warm).(a-f) Posterior means averaged over 10 cm layers of increasing depth, with the deepest (45-55 cm) being at the base and below a typical active layer thickness in organics-covered locations.(g) Sentinel-2 near-infrared-red-green image with elevation contours and annotated coring location (circle), reference point (cross) and transect.(h-i) Posterior mean profiles along T1 (in direction of arrow), with inferred probed depth (estimated y f at last interferometric synthetic aperture radar observation) shown in white.

Figure 12 .
Figure 12.Tatchim Isua, Northwestern Alska: InSAR-derived excess ice profiles from the exceptionally warm summer of 2019.In panel (c), binary classification of upper permafrost ice content from cores taken near the rocky bench, visible in the NIR-R-G Sentinel-2 composite (d).Otherwise same as Figure 10.