A Modification to the Remotely Sensed Active Layer Thickness (ReSALT) Algorithm to have Self-Consistency

1. Introduction

Climate forecasting studies estimate that permafrost thaw could contribute anywhere between 20 to 500 Gt of CO₂-eq by 2100 ([1,2,3,4]). These estimates varied in the phenomenon they modeled, how that phenomenon was modeled, and assumptions of warming trends, therefore explaining the wide range of predictions. Regardless of these differences, a crucial variable in all forecasts is the active layer thickness (ALT), which describes the topmost layer of permafrost subject to annual freeze and thaw. Warmer temperatures lead to deeper active layers. Organic carbon is distributed throughout the soil column, meaning deeper active layers permit microbes to process more carbon and emit more greenhouse gases. Hence, there is a strong interest in developing modeling and forecasting algorithms for ALT. For more background on permafrost carbon cycles, refer to [5].

2. Background

2.1. ReSALT Description

The remotely sensed active layer thickness (ReSALT) algorithm was introduced to model ALT with broad spatial coverage while preserving heterogeneous predictive power [6]. Without remote sensing, permafrost thaw monitoring is constrained to sparse in-situ observations. Given the vast extent of permafrost, it would be incredibly useful if remote sensing could accurately monitor ALT. ReSALT leverages a technique known as interferometric synthetic aperture radar (InSAR) to create its predictions. For background information on InSAR, refer to [7]. Within the context of this paper, an InSAR measurement should be thought of as having three properties: i) the date of the reference SAR acquisition, ii) the date of the secondary SAR acquisition, and iii) a geocoded image containing ground subsidence1 per pixel. Specifically, if the ground position is measured as $\delta$ in an unchanging, earth-centered reference frame and the positive direction points towards the Earth’s center, then the image will consist of $D={\delta }_{reference}-{\delta }_{secondary}$ per pixel.

ReSALT has been used in several studies since its conception ([6,8,9,10]). While it has demonstrated promising results, there is room for improvement. For example, the Permafrost Dynamics Observatory project reported a root mean square error (RMSE) of 15.18 cm in ReSALT predicted vs in-situ ALTs [10]. The in-situ data had an average ALT between 50 and 60 cm, making the RMSE somewhat large in magnitude. This error served as motivation for this paper.

ReSALT is based on two main phenomenon. The first is that thaw depth h is proportional to $\sqrt{ADDT}$, where $ADDT$ stands for Accumulated Degree Days of Thaw. This is the well-known Stefan equation and is shown in (1).

$$h\left(t\right)=N\sqrt{ADDT\left(t\right)}$$

(1)

where N is a proportionality factor. This is derived in Appendix A. The second phenomenon is that when ice thaws into liquid water, the amount of volume reduces. This leads to ground subsidence as shown visually in Figure 1 and mathematically in (2).

$$\delta \left(t\right)={\displaystyle \frac{{\rho }_w-{\rho }_i}{{\rho }_i}}{\int }_0^{h\left(t\right)}PSdh=f\left(h\left(t\right)\right)$$

(2)

where ${\rho }_w$ is the density of water, ${\rho }_i$ is the density of ice, P is the soil porosity, and S is the soil saturation. The integral captures the volumetric liquid water content. The ratio ${\displaystyle \frac{{\rho }_w-{\rho }_i}{{\rho }_i}}$ captures the (positive) difference in volume height when liquid water phase transitions to solid ice.2 The ground subsidence can be computed by multiplying the two quantities, therefore arriving at (2). f is introduced as a convenience to compute subsidence from thaw depth.

Any particular InSAR image contains subsidence difference $D=\delta \left(t_i\right)-\delta \left(t_j\right)$. Temperature is typically known at a coarse resolution. Hence, ReSALT strives to find a way to relate these two quantities and ultimately retrieve $h\left(t\right)$. It does so by posing the following least-squares inversion problem:

$$\left[\begin{array}{c}{D}_1\\ ⋮\\ D_N\end{array}\right]=\left[\begin{array}{cc}{t}_{2,1}-t_{1,1}& \sqrt{ADD{T}_{2,1}}-\sqrt{ADD{T}_{1,1}}\\ ⋮& ⋮\\ t_{2,N}-t_{2,N}& \sqrt{ADD{T}_{2,N}}-\sqrt{ADD{T}_{1,N}}\end{array}\right]\left[\begin{array}{c}R\\ E\end{array}\right]$$

(3)

where $D_i$ is the measured subsidence difference in the ith interferogram, $t_{2,i}$ and $t_{1,i}$ describe the times of the two SAR acquisitions that were used to create the ith interferogram, $ADDT$ is the accumulated degree days of thaw, and $ADD{T}_{2,i}$ and $ADD{T}_{1,i}$ describe the $ADDT$’s at $t_{2,i}$ and $t_{1,i}$. Note that all $D_i$ refer to the same geographic location, but may differ in pixel location from image to image. R and E are the unknown variables and are solved for. R captures a long-term subsidence trend. More information can be found in [6,11]. E captures the relationship between $\sqrt{ADDT}$ and subsidence difference, which can then be used to identify subsidence at any arbitrary time t via:

$$\delta \left(t\right)=E\sqrt{ADDT\left(t\right)}$$

(4)

Then, (2) is inverted to obtain $h\left(t\right)$. More details on ReSALT’s derivation are shown in Appendix B.

2.2. ReSALT Inconsistency Description

In (3) and (4), subsidence is modeled as proportional to $\sqrt{ADDT}$. Using (1) and (2), it becomes evident that this only happens when P is constant with depth. Hence, any application of ReSALT that uses a non-constant porosity model is self-inconsistent. Examples of self-inconsistent applications of ReSALT are [8] and [9], which use the “mixed soil model” shown later in this paper.

3. Method

Observe that if R is ignored, the desired left-hand side (LHS) for (3) is thaw depth difference:

$$h\left(t_i\right)-h\left(t_j\right)=N\left(\sqrt{ADDT\left(t_i\right)}-\sqrt{ADDT\left(t_j\right)}\right)$$

(5)

However, only subsidence difference and $ADDT$ are measured. To see if subsidence difference can be related to thaw depth difference, (1) is analyzed as a ratio:

$${\displaystyle \frac{h\left(t_i\right)}{h\left(t_j\right)}}=\sqrt{{\displaystyle \frac{ADDT\left(t_i\right)}{ADDT\left(t_j\right)}}}$$

(6)

$$h\left(t_i\right)=h\left(t_j\right)\sqrt{{\displaystyle \frac{ADDT\left(t_i\right)}{ADDT\left(t_j\right)}}}$$

(7)

For convenience, $K:=\sqrt{{\displaystyle \frac{ADDT\left(t_i\right)}{ADDT\left(t_j\right)}}}$ is introduced. Next, a list of candidate $h_{t_j}$ on some large interval of possible thaw depths is generated, and corresponding $h_{t_i}$ is computed using (7) for a fixed K. The subsidence for each candidate $h_{t_j}$ and $h_{t_i}$ can be computed using (2). A plot of candidate thaw depth differences and corresponding subsidence differences can then be constructed. See Figure 2 for examples of this plot for two soil porosity models. Each of these models is explained in [8]. As shown in the bottom two figures, a particular subsidence difference can unambiguously be mapped to a thaw depth difference. Intuitively, this is because (7) enables multiplicative growth of $h_{t_i}$ with respect to $h_{t_j}$. Consequently, the hope is that (2) enables subsidence difference to grow strictly monotonically. Put another way, for each interferogram, K (which is known) acts as a constraint parameter that helps search the space of possible thaw depth pairs $h\left(t_j\right)$ and $h\left(t_i\right)$. Corresponding subsidences can be computed and a graph can depict whether subsidence difference can be related to thaw depth difference. Of course, there is no guarantee this will work for all soil porosity models. A counterexample is discussed in Section 3.1. Also, note that in Figure 2 the bottom images demonstrate the effect of K on the method’s robustness. For the constant porosity model, there is no effect because the thaw depth difference is linearly related to the subsidence difference. For the “mixed soil model,” the slope is higher in the bottom-left image, making it easier to identify a particular thaw depth difference from a given subsidence difference. This is because larger K make the thaw depth pair have smaller values for the same thaw depth difference, which then exhibit larger subsidence gradients, as the top-right plot shows.

This method is summarized in Algorithm 1. (3) can be recast as:

$$\left[\begin{array}{c}{y}_1\\ ⋮\\ y_N\end{array}\right]=\left[\begin{array}{c}\sqrt{ADD{T}_{2,1}}-\sqrt{ADD{T}_{1,1}}\\ ⋮\\ \sqrt{ADD{T}_{2,N}}-\sqrt{ADD{T}_{1,N}}\end{array}\right]\left[\begin{array}{c}N\end{array}\right]$$

(8)

Algorithm 1:Self-Consistent ReSALT LHS Retrieval

Require: $K=\sqrt{{\displaystyle \frac{ADD{T}_{t_i}}{ADD{T}_{t_j}}}}$  Require: $D=\mathrm{measured}\mathrm{subsidence}\mathrm{difference}$  Require: $N=\mathrm{number}\mathrm{of}\mathrm{samples}$  Require: $L=\mathrm{lower}\mathrm{end}\mathrm{of}\mathrm{possible}\mathrm{thaw}\mathrm{depths}$  Require: $U=\mathrm{upper}\mathrm{end}\mathrm{of}\mathrm{possible}\mathrm{thaw}\mathrm{depths}$  1: $step\leftarrow {\displaystyle \frac{U-L}{N-1}}$▹ Step size for linear spacing  2: for$k=1\to N$do  3:     $h_{k,j}\leftarrow L+(k-1)\times step$  4:     $h_{k,i}\leftarrow h_{k,j}\times K$  5:     $y_k\leftarrow h_{k,i}-h_{k,j}$  6:     $x_k\leftarrow f\left(h_{k,i}\right)-f\left(h_{k,j}\right)$▹f is defined in (2)  7: end for  8: if not check_monotonically_changing$(x,y,K)$ then  9:     raise Error(“Fails SCReSALT requirement.”)  10: end if  10: This can be done with binary search, which enables fast solving for hundreds of pixels. Using a root-finding approach like Newton’s method would be slower.  11: Find k such that $x_k$ is closest to D  12: return$y_k$

where $y_i$ is the best-matching thaw depth difference reported for each $D_i$ via Algorithm 1. This now forms a self-consistent modeling of thaw, temperature, and ground subsidence under variable soil porosity models. Hence, this new method is called SCReSALT (Self-Consistent ReSALT). Least-squares inversion gives the best N. (1) can then be used to find the thaw depth at any point in the thaw season.

To extend the formulation to multiple years, the long-term trend R must be discounted from $D_i$ before following the above approach. ReSALT jointly inverts for R and E, but a self-consistent model cannot do that. So, R must be estimated independently of and prior to N. Solving this is out-of-scope for this paper, but here are some ideas.

(i)

Form interferograms across multiple years at the start of the thaw season, ideally when $ADDT=0$. Any observed subsidence can be attributed to R. Compute the best-fit R that solves all $D_i=R(t_i-t_j)$.

(ii)

Compute an initial estimate of N by solving (8). N could be solved either per thaw season (which makes $R=0$) or overall. Then, form interferograms across multiple years at the start of the thaw season and remove the subsidence contribution due to any initial thaw. This is known because the subsidence contribution due to initial thaw can be written as: $D_{h,i}=f\left(N_{s_{2,i}}\sqrt{ADD{T}_{2,i}}\right)-f\left(N_{s_{i,1}}\sqrt{ADD{T}_{1,i}}\right)$. Note $N_{s_{2,i}}$ denotes N for the thaw season in the 2nd SAR image in the ith interferogram. $N_{s_{1,i}}$ means the same for the 1st SAR image. Compute the best-fit R that solves all $D_i-D_{h,i}=R(t_i-t_j)$.

Once R is known, the multi-year problem can be solved using (8) but instead of giving Algorithm 1 the observed subsidence difference $D_i$, give it $D_i-R(t_i-t_j)$.

For completeness, it is worth noting that both ReSALT and SCReSALT do not need both SAR acquisitions to be during the thaw season. One acquisition could be before the thaw season, in which case $ADDT\approx 0$ can be assumed for that acquisition. Hence, any measured subsidence difference can be modeled as the true subsidence at that point. ReSALT can directly use (3). SCReSALT would need a slight modification because it operates using ratios. Algorithm 1 could be modified such that if $ADD{T}_{t_j}<ϵ$, simply report $y=f^{-1}\left(\delta \right)$. Similarly, if $ADD{T}_{t_i}<ϵ$, report $y=f^{-1}(-\delta )$.

Lastly, in Appendix A, the Stefan equation is shown to also depend on a constant porosity model. Appendix D introduces SCReSALT-NS which relaxes this constraint. `NS’ stands for `non-Stefan’.

3.1. SCReSALT Counterexample Discussion

In Appendix C, a counterexample that does not meet the SCReSALT requirement is derived. Some possible solutions are described below.

(i)

If a subsidence difference vs thaw depth difference graph can uniquely identify a subset of thaw depth differences, one could use the known thaw depth differences and the expectation that thaw depth only increases with time to select the right thaw depth.

(ii)

Make sure to get a SAR acquisition before the thaw season starts. The measured subsidence difference is now true subsidence. Handling this was described earlier.

Note that any soil porosity model can be tested to see if it meets the SCReSALT requirement using Algorithm 1.

4. Experiments

To assess the practical improvement of SCReSALT, it is benchmarked against the results reported in [9]. [9] applies ReSALT to Utqiagvik (formerly Barrow), Alaska and uses the “mixed soil model” to predict ALT from 2006 to 2010. Metrics are extracted from the paper and the data product they release, which can be found at https://daac.ornl.gov/ABOVE/guides/ReSALT_InSAR_Barrow.html. [9] uses data from the U1 CALM3 plot. It uses L-band ALOS-PALSAR data and generates interferograms using the ISCE2 framework4. The paper defines the ${\chi }^2$ metric:

$${\chi }^2={\displaystyle \frac{{(AL{T}_{pred}-AL{T}_{gt})}^2}{e^2}}$$

(9)

where e is the uncertainty in the measurement. In [9], $e=7.9$ cm. If ${\chi }^2<1.0$, it is called a “great match.” If the $AL{T}_{gt}$ is within the uncertainty of $AL{T}_{pred}$, it is called a “good match." Otherwise, it is called a “bad match.” Uncertainty analysis is not described in this paper, but in brief, it could be computed by assigning an uncertainty distribution to each ReSALT/SCReSALT parameter and either analytically propagating that into a prediction uncertainty or more practically using numerical simulation. This work uses the fact that [9] stated that the ReSALT uncertainty is roughly twice that of the measurement uncertainty. Hence, $e_{ReSALT}=e_{SCReSALT}=15.8$ cm. Other experimental details are described in Appendix E. Results are shown in Table 1.

5. Results

The metrics reported in [9] closely match the results computed from the data product. The “good match” category is a little worse because $e_{ReSALT}$ is not precisely known. That then influences the “bad match” category. Otherwise, the alignment is high, meaning there is high confidence in the derived Pearson R, mean absolute error (MAE), and RMSE values. The author also attempted to reproduce ReSALT but generally obtained worse numbers. It is difficult to get precise replication without source code, which is not available for [9]. Regardless, SCReSALT is shown to outperform ReSALT on virtually all metrics and handles outliers better given its reduced RMSE and ${\chi }^2$. If compared to the reproduced ReSALT, SCReSALT’s improvement is even stronger. Furthermore, SCReSALT-NS demonstrates improvement over SCReSALT, suggesting that incorporating non-Stefan thaw dynamics can improve modeling accuracy. Also, the author was puzzled by the ReSALT data product’s low Pearson R. It is not clear why, but given the high measurement uncertainty, it is possible to have many “great matches” despite no correlation between predicted and ground-truth ALTs.

6. Conclusion

This work presents SCReSALT and SCReSALT-NS, two algorithms that modify ReSALT to enable self-consistent modeling of thaw, temperature, and ground subsidence in permafrost. It describes the conditions where SCReSALT can be used and shows its applicability to previously developed soil models. It is also shows that SCReSALT doesn’t work on all soil models. Possible solutions are discussed, but exploring them is left to future work. Experimentally, SCReSALT demonstrates significant improvement in estimating ALTs in Utqiagvik, Alaska. The author recommends that SCReSALT become a standard tool in permafrost thaw modeling.

Conflicts of Interest: The author declares no conflicts of interest.