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 microbacteria to process more carbon and emit more greenhouse gases. For more background on permafrost carbon cycles, the reader is referred 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, our understanding of permafrost thaw 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, the reader is referred to [7]. For the purpose of analysis in 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) an 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, such as in [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 of around 55 cm, making the RMSE significantly large with respect to the unit of measure. This large error served as motivation for this paper.

ReSALT is based on two main phenomenon. The first is ALT 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_i\right)=N\sqrt{ADDT\left(t_i\right)}$$

(1)

where N is a proportionality factor. This is derived in Appendix A. The second is that soil is fully saturated, meaning all pore space is filled with water. Consequently, when ice thaws into liquid water, the amount of volume reduces, which leads to a subsidence of the ground. This is shown visually in Figure 1 and captured mathematically in (2).

$$\delta \left(t_i\right)=\frac{p_w-p_i}{p_i}{\int }_0^{h\left(t_i\right)}PSdh=f\left(h\left(t_i\right)\right)$$

(2)

where $p_w$ is the density of water, $p_i$ is the density of ice, P is the soil porosity, and S is the soil saturation. The integral captures the volumetric water content. The ratio $\frac{p_w-p_i}{p_i}$ describes the difference in volume height when the ice phase transitions to liquid water per unit volume of ice.2 We can then compute the difference in ground height by multiplying the two quantities, therefore arriving at (2). f is introduced as a convenience to compute subsidence from ALT.

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, and $ADD{T}_{2,i}$ and $ADD{T}_{1,i}$ describe the $ADDT$’s for $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) can be inverted to obtain $h\left(t\right)$. More details on ReSALT’s derivation are shown in Appendix B.

2.2. ReSALT Inconsistency Description

In (4), we see that subsidence difference is modeled as proportional to $\sqrt{ADDT}$. Using (1) and (2), we can see 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

The fundamental problem is that we only measure $ADDT$ and subsidence difference $D={\delta }_{t_i}-{\delta }_{t_j}$, but we know no relationship between the two. We will first demonstrate a simple way to solve this problem for a practical soil moisture model. Then, we will more rigorously analyze the proposed method to understand its properties.

If we ignore R, we can see that the correct inversion LHS (left-hand side) for (3) is ALT difference:

$$h\left(t_i\right)-h\left(t_j\right)=\sqrt{\frac{2k_p}{p_{p}L}}\left(\sqrt{ADDT\left(t_i\right)}-\sqrt{ADDT\left(t_j\right)}\right)$$

(5)

To see if subsidence difference can be related to ALT difference, we analyze (1) as a ratio:

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

(6)

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

(7)

Next, we fix some ratio $Q=\sqrt{\frac{ADDT\left(t_i\right)}{ADDT\left(t_j\right)}}$, create a list of $h_{t_j}$ on some large interval of candidate ALTs, generate corresponding $h_{t_i}$ use (7) to calculate the corresponding subsidences for each $h_{t_j}$ and $h_{t_i}$, and then plot ALT differences and corresponding subsidence differences. See Figure 2 for examples of this plot for two soil porosity models. Each of these models is explained explicitly in [8]. As shown via a “horizontal line test,” a particular subsidence difference can always be directly mapped to an ALT difference with zero ambiguity. Intuitively, this is because (7) enables multiplicative growth of $h_{t_i}$ with respect to $h_{t_j}$. Consequently, we would hope that (2) enables subsidence difference to grow strictly monotonically. Of course, there is no guarantee that this always will be the case. A counterexample is discussed in Section 3.1.

This method is summarized Algorithm 1. We can now rewrite (3) as:

$$\left[\begin{array}{c}{y}_i\\ ⋮\\ 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}A\end{array}\right]$$

(8)

where $y_i$ is the best-matching ALT difference reported for each $D_i$ via Algorithm 1. This now forms a self-consistent modelling of thaw, temperature, and ground subsidence. Hence, we call this new method SCReSALT (Self-Consistent ReSALT). Least-squares inversion gives the best A, which can then be used to find ALT at any point in the thaw season:

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

(9)

Algorithm 1 Self-Consistent ReSALT LHS Retrieval

Require: $Q=\sqrt{\frac{ADD{T}_{t_i}}{ADD{T}_{t_j}}}$ Require: $D=\mathrm{measured}\mathrm{subsidence}\mathrm{difference}$ 1: $N\leftarrow \mathrm{number}\mathrm{of}\mathrm{samples}$ 2: $L\leftarrow \mathrm{lower}\mathrm{end}\mathrm{possible}\mathrm{ALTs}$ 3: $U\leftarrow \mathrm{upper}\mathrm{end}\mathrm{possible}\mathrm{ALTs}$ 4: $step\leftarrow \frac{U-L}{N-1}$             ▹ Calculate step size for linear spacing 5: for$k=1\to N$do 6:     $h_{k,j}\leftarrow L+(k-1)\times step$ 7:     $h_{k,i}\leftarrow h_{k,j}\times Q$ 8:     $y_k\leftarrow h_{k,i}-h_{k,j}$ 9:     $x_k\leftarrow f\left(h_{k,i}\right)-f\left(h_{k,j}\right)$                   ▹ See (2) 10: end for 11: $e\leftarrow$ small neighborhood in which a subsidence difference should be unique 12: if not check_unique$(x,y,e)$ then 13:     raise Error(“x and y have non-unique matching.”) 14: end if 15: Find k such that $x_k$ is closest to D 16: return$y_k$

To extend the formulation to multiple years, we simply need to discount the long-term trend R from the measured subsidence difference $D_i$ before following the above approach. ReSALT jointly inverts for R and E, but a self-consistent model cannot do that. So, we must find a way to estimate R independently of and prior to A. Thankfully, there are many ways around this. Two are discussed here.

• Form interferograms across multiple years at the start of thaw season, ideally when $ADDT=0$ (and perhaps there is some tolerance on how large $ADDT$ can get before it affects the result due to seasonal subsidence). Any observed subsidence can be attributed to R. Compute the best-fit R that solves all $D_i=R(t_i-t_j)$.
• Compute an initial estimate of A by solving (8). A could be solved either per thaw season (which makes $R=0$) or overall. Then, like the previous bullet point, form interferograms across multiple years at the start of thaw season and remove the subsidence contribution due to ALT. This is known because the subsidence contribution due to ALT can be written as: $D_{ALT,i}=f\left(A_{s_{2,i}}\sqrt{ADD{T}_{2,i}}\right)-f\left(A_{s_{i,1}}\sqrt{ADD{T}_{1,i}}\right)$. Note $A_{s_{2,i}}$ denotes A for the thaw season in the 2nd SAR image in the ith interferogram. $A_{s_{1,i}}$ means the same for the 1st SAR image. Compute the best-fit R that solves all $D_i-D_{ALT,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 precisely during the thaw season. One acquisition could be before the thaw season, in which case we can assume $ADDT\approx 0$ 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 )$. The sign change is because $\delta$ is the subsidence underwent when going from $t_j$ to $t_i$. That is expected to be positive if $ADD{T}_{t_i}>ADD{T}_{t_j}$ and negative vice-versa.

3.1. SCReSALT Counterexample Discussion

A more rigorous analysis of SCReSALT is shown in Appendix C. It shows a counterexample with a realistic soil porosity model that does not meet the SCReSALT requirement. Some possible workarounds are described below, but we leave exploring them to future work.

• From multiple interferograms per year and obtain time-series measurements of subsidence difference. Attempt to interpolate the measurements onto an ALT difference vs subsidence difference graph and determine which ALT differences best explain the data.
• Make sure to get a SAR acquisition before the thaw season starts. This allows one to treat the subsidence difference as the true subsidence. Handling this is easy and was described earlier in this paper.
• Incorporating a prior on ALT difference might isolate a particular ALT difference as the most likely candidate.

4. Experiments

To assess the practical improvement of our method, we reproduce and benchmark against the results reported in [9]. [9] applies ReSALT to Utqiagvik (formerly Barrow), Alaska to predict ALT from 2006 to 2010. The paper uses data from the U1 CALM3 plot. It uses ALOS-PALSAR L-band data and generates interferograms using the ISCE2 framework4. The paper defines the ${\chi }^2$ metric:

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

(10)

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.” We have not described uncertainty analysis in this paper, but in brief one could compute uncertainty by assigning an uncertainty to each ReSALT/SCReSALT parameter and either analytically propagating that into an uncertainty in the prediction or more practically using numerical simulation to understand the distribution of predictions as a function of the modeling parameter distributions. Since this paper is primarily concerned with RMSE improvement, and uses ${\chi }^2$ as a way to assess reproduction accuracy, we do not discuss uncertainty explicitly in this analysis. Instead, we take advantage of the fact that [9] states that the ReSALT uncertainty is roughly twice that of the measurement uncertainty. Hence, we use $e_{ReSALT}=e_{SCReSALT}=15.8$ cm. This makes the “great match,” “good match,” and “bad match” categorizations directly comparable. Additionally, [9] uses the same “mixed soil model” from [8]. Further specifics (namely image calibration) on the reproduction are described in Appendix D. With those details explained, results are shown in Table 1.

5. Discussion

As we can see in Table 1, the reproduced ReSALT had much worse results than in the original paper. Private correspondence and a review of the software with the authors of [9] suggested that the reproduction was most likely very similar. Also, it is worth noting that [10] reported similar RSMEs to the reproduced ReSALT (albeit for a broader set of CALM sites). For completeness, we show comparisons between the reproduced ReSALT, the original paper, and SCReSALT. Given that SCReSALT outperforms the original paper on the ${\chi }^2$ metric, we take these results to mean SCReSALT outperforms ReSALT. Of course, this is to be expected. SCReSALT makes no approximations and can perfectly model the “mixed soil model” in a self-consistent way.

6. Conclusion

We present SCReSALT, a modification to ReSALT that enables self-consistency. SCReSALT makes no approximations and reduces to ReSALT under the right soil models. We also rigorously analyze SCReSALT’s requirements and identify a counterexample where the requirements are not met. We then describe possible workarounds but leave exploring that to future work. Experimentally, SCReSALT demonstrates significant improvement in predicting ALTs in Utqiagvik, Alaska. We recommend SCReSALT become a standard tool in permafrost thaw modelling.