Landslide Thickness Estimated from InSAR-derived 2D Deformation: Application to the Xiongba ancient Landslide, China

1. Introduction

The slow-moving landslide is a typical geological disaster, commonly occurring in the mountainous regions of western China. The failure of slow-moving landslides can cause significant casualties and substantial property losses, while also providing ample material for debris flows, particularly during the rainy season, potentially leading to severe cascading geological disaster events. (Nadim et al., 2006; Petley et al., 2012). Therefore, it is imperative to monitor the surface movement velocity and volume of landslides, with particular emphasis on landslide thickness, as it directly influences the estimation of landslide volume and hazard extent, serving as a critical factor in landslide risk assessment (Guzzetti et al., 2009). However, the detection of landslide thickness poses significant challenges, primarily because it is often not possible to access the underlying base surface of the landslide.

Two primary approaches exist for estimating active thickness are the contact and non-contact methods. Contact methods rely on direct on-site measurements, typically utilizing inclinometer detection (Simeoni et al., 2007; Stark et al., 2008) or ground-penetrating radar techniques (Delbridge et al., 2016). However, contact detection methods are limited to measuring the thickness at specific points rather than across the entire landslide, due to the significant financial and time constraints involved (Mora et al., 2003), making it unsuitable for assessing large-scaleand mega-landslides.

As developed, several non-contact methods are proposed, including the equilibrium cross-section method (Bishop, 1999; Aryal et al., 2015), the dislocation method (Nikolaeva et al., 2014), and the mass conservation method (Booth et al., 2013b; Huang et al., 2017). The equilibrium cross-section method assumes that landslide material behaves as incompressible rigid bodies, neglecting the influence of rheological properties, which leads to a less reliable estimation of landslide thickness (Rowan et al., 1989; Zhang et al., 2014). The dislocation theory method simplifies both the geometric configuration and the physical mechanisms of the landslide. Consequently, its application is limited to the initial stage of landslide development or landslides without significant inelastic deformation (Feigl et al., 1999; Aryal et al., 2015). The mass conservation approach is an emerging approach that integrates surface deformation and rheological properties to estimate the thickness of a moving body. This method has been widely applied in landslide and glacier thickness calculations (Farinotti et al., 2009; Morlighem et al., 2011; Booth et al., 2013; Hu et al., 2018; Handwerger et al., 2021), and is considered a reliable and practical approach for landslide thickness detection (Handwerger et al., 2021; Yang et al., 2021). In this method, the rate of thickness variation is balanced with the divergence of horizontal mass flux across the landslide surface. And the landslide depth can be calculated by integrating along the basal sliding surface with dynamic boundary conditions. However, its implementation requires three-dimensional deformation fields and accurate rheological parameters, which can be difficult to obtain, particularly in remote areas (Booth et al., 2013; Delbridge et al., 2016; Hu et al., 2016; Nguyen et al., 2016).

In this study, an improved mass conservation approach is developed to detect landslide thickness based only on the ascending and descending orbital tracks InSAR observations. The developed method can be used to infer the landslide thickness based on the decomposed 2D deformation along the slope and normal directions. This method is applied to the Xiongba ancient landslide in southwest China. Geological survey and drilling data are introduced to verify the reliability of the proposed method. Additionally, the best-fitting rheological parameters are determined based on the thickness data of sparse points.

2. Study Area and Dataset

2.1. Geologic Background of the Xiongba Ancient Landslide

The Xiongba ancient landslide is located in Shadong Township, Gongjue County, Tibet Autonomous Region, China, on the west bank of the Jinsha River (Figure 1). In this region, the Jinsha River runs predominantly in a north-south direction. The local terrain exhibits significant complexity, with elevations varying between 2600–4000 m. And the area is characterized by a continental plateau monsoon climate, with an average annual rainfall of 480 mm, the majority of which occurs between May and September. In addition, the Xiongba ancient landslide is situated within the Jinsha River seismic zone, which hosts several significant active faults, including the Jinsha River Fault and the Batang Fault, and has a history of frequent seismic activity (Jin et al., 2021).

The Xiongba ancient landslide presents a chair-shaped structure, and can be divided into three distinct regions according to the topographic feature (Figure 1b). The rear wall (Region I) of the landslide is characterized by steep terrain, with elevations ranging from 3359 to 3800 meters. The slope gradient generally varies between 40° and 45°, but occasionally exceeds 60–70°. Previous remote sensing interpretation and field geological surveys indicate that the rear wall region is generally stable (Yan et al., 2021). The central platform (Region II) features relatively flat terrain, primarily distributed at elevations between 3326–3534 m, with slopes ranging from 10° to 15°, and generally exhibits a stable disposition. The frontal deformed area (Region III) is the main portion of the ancient landslide, with elevations varying between 2693–3216 m, and is currently in a creeping state. The landslide extends ~2.13 km in length and ~2.68 km in width, covering an area of ~5.33 km²;, and the volume is estimated to be between 2.6×10⁸–6.0×10⁸ m³ (Li et al., 2021). Furthermore, the topographic features and deformation distribution indicate that the Xiongba ancient landslide consists of two secondary landslides of H1 and H2, as shown in Figure 1b. (Jin et al., 2021).

The bedrock of the Xiongba ancient landslide is mainly composed of mica-quartz schists of Permian-Lower Triassic Gangtoite Formation (PT_1g). The sediment of the Xiongba landslide is mainly composed of Quaternary broken stone soil (Q₄^del), which contains approximately 50% crushed stone with particle sizes of 2~10 cm, 20% block stone with block diameter of 20~40 cm, and silty clay (Jin et al., 2021). And riverbed sediment of the Jinsha River is mainly composed of Quaternary alluvial-proluvial sand and gravel (Q₄^al+pl). Previous study suggested that the Xiongba landslide may fail along weak surfaces under the influence of localized heavy rainfall, potentially triggering a chain of disaster events, including river blockage and dam failure, which could seriously threaten the safety of surrounding villages (Li et al., 2021).

2.2. Ground Deformation of the Xiongba Landslide Extracted by InSAR Time-Series Analysis

In this study, 59 ascending track and 68 descending track radar images (Coverage and detailed parameters can be found in Figure 1a and Table S1) captured by the Sentinel-1 satellite were collected to map the deformation fields of the Xiongba landslide. The Precise Orbit Ephemerides (POD) data provided by European Space Agency (ESA) was employed to refine satellite orbit vectors. The SRTM Digital Elevation Model (DEM) data with 30 m resolution was applied to remove the topography phase component from the interferograms. And the minimum cost flow algorithm was adopted to unwrap the interferogram (Costantini et al., 1999; Pepe et al., 2011). The Generic Atmosphere Correction Online Service for InSAR (GACOS) data were used to correct for atmospheric delays.

The SBAS-InSAR technique was applied to extract the surface InSAR velocity fields (Figure 2) of the Xiongba landslide based on the collected ascending and descending SAR images, in which a combination of high-pass and low-pass filtering operations were performed to further mitigate the atmospheric delay errors (Berardino et al., 2002). The temporal and spatial baselines were set to 120 days and 300 meters, respectively, resulting in 445 and 587 interferometric pairs for the ascending and descending InSAR tracks. Considering that the GNSS data is lacking, the absolute mean of the InSAR observations in the relatively stable zone, far away from the landslide, is calculated to quantitatively assess the uncertainty of the observed InSAR deformation. And it is found that the uncertainties of 5.1 mm/yr and 5.7 mm/yr are respectively for the ascending and descending InSAR tracks. The millimeter-level accuracy mentioned above indicates that the extracted InSAR observations are sufficiently reliable for further revealing the surface motion of the Xiongba-H2 landslide.

Figure 2a,b show notable differences between the ascending and descending InSAR velocity fields of the Xiongba-H2 landslide in terms of deformation magnitude and spatial distribution, which could be attributed by the difference in imaging geometry between the ascending and descending InSAR tracks. The velocity of the ascending InSAR varies from -172 mm/yr to 8 mm/yr, and the descending InSAR deformation rate varies from -10 mm/yr to 78 mm/yr. However, some similarity can still be found in the ascending and descending InSAR velocity fields. First, two regions with different deformation patterns can be addressed in both ascending and descending velocity field, including a frontal zone with significant deformation rates marked by the Region 2 and a relatively stable zone (Region 1). Then, there is no rapid acceleration shown in the displacement time series for T1 from ascending track and T2 from descending track, suggesting that the Xiongba-H2 landslide is likely in a relatively stable status with no immediate risk of failure. However, it should be noted that the landslide risk assessment relies on long-term deformation observations, along with a comprehensive analysis of regional rainfall, topography, geomorphology, and the rock-soil characteristics of the landslide.

3. Methods

3.1. Decomposition of the Derived InSAR LOS Deformation Rates

Figure S2 shows the landslide motion coordinate system, in which OK, OT and OI respectively indicate the slope, perpendicular to the slope and the normal directions of the landslide. The mathematical relationship between the slope coordinate system and geographic coordinate system could be depicted as

$$\left[\begin{array}{l}{D}_E\\ D_N\\ D_U\end{array}\right]=\left[\begin{array}{ccc}\cos \alpha \sin \beta & -\cos \beta & \sin \alpha \sin \beta \hfill \\ \cos \alpha \cos \beta & \sin \beta & \hfill \sin \alpha \cos \beta \\ -\sin \alpha & 0& \hfill \cos \alpha \end{array}\right]\left[\begin{array}{l}{D}_K\\ D_T\\ D_I\end{array}\right]$$

(1)

Where, $D_E$, $D_N$ and $D_U$ represent deformation along the East, North, and Vertical directions in the geographic coordinate system. $D_K$, $D_T$, and $D_I$ indicates the deformation along OK, OT and OI directions, respectively. $\alpha$ and $\beta$ are the slope and aspect angles of the landslide.

According to the SAR imaging geometry, the LOS deformation could be expressed by the $D_E$, $D_N$ and $D_U$ deformations as

$$D_{LOS}=\left[ -\sin \theta \cos \phi \sin \theta \sin \phi \cos \theta \right]{\left[D_E D_N D_U\right]}^T$$

(2)

Where, $D_{LOS}$ denotes the deformation along LOS direction. $\theta$ and $\phi$ are respectively the incidence angle and azimuth angle. Combining the Equations (1) and (2), the InSAR LOS deformation can be expressed by the 3D deformation under the slope coordinate system, shown as Equation (3).

$$\begin{array}{l}{\mathrm{D}}_{\mathrm{LOS}}=\left[\mathrm{a} \mathrm{b} \mathrm{c}\right]\left[\begin{array}{c}{\mathrm{D}}_{\mathrm{K}}\\ {\mathrm{D}}_{\mathrm{T}}\\ {\mathrm{D}}_{\mathrm{I}}\end{array}\right]\\ a =–\cos \theta \sin \alpha +\sin \theta \sin \phi \cos \alpha \cos \beta – \sin \theta \cos \phi \cos \alpha \sin \beta \\ b =\sin \theta \sin \phi \sin \beta +\sin \theta \cos \phi \cos \beta \\ c =\cos \theta \cos \alpha +\sin \theta \sin \phi \sin \alpha \cos \beta – \sin \theta \cos \phi \sin \alpha \sin \beta \end{array}$$

(3)

Where, $a$, $b$, and $c$ are projection coefficients in the three directions.

The 3D deformation in the landslide motion coordinate system can be solved when three or more orbital tracks of InSAR observations are available. However, it is generally difficult to simultaneously extract the ground deformation from three independent InSAR tracks. Considering that the slow creeping of the landslide is controlled by the gravity of the shallow rock-soil body, the displacements along the slope and normal directions are expected to be the main deformation components. Therefore, the deformation perpendicular to the slope direction (along the OT axis in Figure S2) could be ignored due to the negligible magnitude. And above-mentioned assumption has been widely used to estimate the surface deformation of the landslide (Raucoules et al., 2013; Meng et al., 2019; Meng et al., 2021; Dai et al., 2022). Fortunately, the Sentinel-1 satellites are capable of providing SAR images from both ascending and descending orbital tracks for most regions globally with a high temporal resolution (12-day revisit time). It suggests that the landslide surface deformation rates could be probably obtained using the Sentinel-1 ascending and descending acquisitions. Subsequently, the 2D deformation, comprising the slope and normal components in the slope coordinate system, can be estimated based on observations from two InSAR tracks, further enhancing the applicability of InSAR technology for analyzing landslide motion.

3.2. Estimation of Landslide Thickness Based on 2D Deformation

For a landslide with constant density, the principle of mass conservation dictates that

$${\displaystyle \frac{\partial h}{\partial t}}+\nabla \cdot (h \overline{v})=0$$

(5)

Where, $h$ denotes landslide thickness, $t$ is time. $\overline{v}=(v_x,v_y)$ indicates the depth-averaged horizontal velocity vector, written as $\overline{v}=f\cdot v_{surf}$. $v_{surf}$ refers to the surface horizontal velocity, $f$ represents the rheological parameter ranging from 0 to 1 (Delbridge et al., 2016).

It can be assumed that the changing of landslide sliding base is insignificant during the deformation observed period. Then, the changing rate of landslide thickness is equal to the landslide deformation velocity along normal direction $v_z$, thus the Equation (5) could be re-written as:

$$v_z=-f \nabla \cdot (h{v}_{surf})$$

(6)

Furthermore, as mentioned above, the motion perpendicular to the slope-dipping direction can be ignored, as the slow-moving landslide is predominantly controlled by the deformation components along the slope and normal directions. Employing the finite difference method, Equation (6) could be further re-written as

$$v_z(i,j)=-{\displaystyle \frac{f{v}_x(i,j)}{2\Delta x}}\cdot h_{i+1,j}+{\displaystyle \frac{f{v}_x(i,j)}{2\Delta x}}\cdot h_{i-1,j}-f({\displaystyle \frac{v_x(i+1,j)-v_x(i-1,j)}{2\Delta x}})\cdot h_{i,j}$$

(7)

Where, $\Delta x$ indicates the spatial sampling interval in slope direction, $v_x$ denotes the deformation velocity along the slope direction. It suggests that the inversion of landslide thickness can be achieved based on the estimated 2D deformation from the combined ascending and descending InSAR observations.

3.3. Determination of Optimal Rheological Parameter

It can be seen from Equation (7) that the rheological parameter is crucial for estimating landslide thickness. This parameter can be assumed to be spatially uniform, with its value depending on the landslide rheology, ranging from 0 to 1 (Delbridge et al., 2016). To determine the best-fitting rheological parameter of the Xiongba-H2 landslide, a comparison was made between the landslide thickness measured through geophysical surveys and drilling data, and the thickness estimated through geodetic data inversion. First, we collected the geophysical exploration and drilling data of the Xiongba-H2 landslide. Then, the landslide thickness along the geological profiles I-I’ was extracted from the geophysical surveys and drilling data (Figure 1c). The optimal rheological parameter was determined from a testing of the experimental values of the rheological parameter varying from 0 to1 with a step of 0.01. For each given f, the landslide thickness along the profile I-I’ was estimated using Equation (7). The model misfit between the estimated and measured landslide thickness was derived. Figure 3d shows the model misfit versus f, and the best-fitting rheological parameter of 0.75 is selected when the model misfit reaches its minimum.

4. Results and Discussions

4.1. Two-Dimensional Deformation Fields of the Xiongba-H2 Landslide and Field Investigation

As described in Section 3.1, the two-dimensional deformation fields can be derived using the observed InSAR LOS deformation rates of the Xiongba-H2 landslide (Figure 2a,b), with the results shown in Figure 2c,d. It can be easily seen that significant deformation is mainly concentrated in the Region 2. The deformation rates along the slope direction range from 10.7 mm/yr to 195.1 mm/yr, while those in the normal direction range from -47.9 mm/yr to 60.5 mm/yr. Focusing on Figure 2c, significant deformation with a magnitude greater than 70 mm/yr is mainly concentrated in the area marked by ‘A’, located in the middle of the Region 2. However, the deformation in normal direction is mainly concentrated at the front edge (marked by ‘B’, Figure 2d) of the Xiongba-H2 landslide. The observed positive normal deformation within the ‘B’ area suggests a significant surface uplift motion. Furthermore, considering that the surface deformation rates along the slope direction are relatively low in the ‘B’ area, thus the surface uplift at the front edge of Xiongba-H2 landslide could probably be resulted from the pushing of the rock-soil mass in the ‘A’ area.

In addition, the field investigation indicates that a transverse tensile crack is found in the upper ‘A’ area, trending NE 2°, with a length of approximately 350 m and widths of 10-50 cm, as shown in Figure S3f. And the Unmanned Aerial Vehicle (UAV) optical images show that multiple longitudinal cracks (marked by the red solid lines in Figure S3) are present on the slope surface, with lengths ranging from 100 to 500 m. Since displacements in both the slope and normal directions are not significant in the C area (Figure 2), these surface fractures could result from the forward movement of the rock-soil mass in the A area.

4.2. Estimated Thickness of the Xiongba-H2 Landslide

The distributions of the landslide thickness were estimated using different rheological parameters of 0.5, 0.75 and 1.0 for comparison, as shown in Figure 3a–c. It can be found that the three results obtained from different rheological parameters show a similar spatial distribution pattern. However, the estimated landslide thickness gradually decreases as the value of the rheological parameter increases. Meanwhile, three regions (I, II, and III) with relatively greater thickness were identified from the results. Notably, these regions are characterized by abnormally elevated local topography (Figure S3), which is likely due to the accumulation of rock-soil mass from past landslide events. In addition, Region III is thicker than other parts of the landslide’s front edge along the Jinsha River, indicating an ample sediment source for the potential generation of landslide hazard. Furthermore, the front edge of the Xiongba-H2 landslide is characterized by steep topography, sparse vegetation cover, and abundant water resources, all of which suggest that Region III has a high potential for future failure.

Figure 3e shows the comparison between the estimated landslide thickness by geodetic inversion and the measured thickness by geophysical surveys and drilling data, along the profile P-P’. It could be found from Figure 3e that the estimated thickness distribution along the P-P’ profile differs for each rheological parameter, however the overall trend in thickness change remains consistent across the results. The estimated thickness reaches its maximum at the gently sloping platform located at the front edge of the Xiongba-H2 landslide, approximately 950 m from the starting point P. It can be easily seen that the estimated thickness profile with the rheological parameter of 0.75 has a better consistency with the measured thickness. And Figure 3d also indicates that the optimal rheological parameter is 0.75, where the misfit reaches its minimum. The average absolute difference is calculated to be less than 10 m, which is negligible compared to the average thickness of ~90 m along the profile P-P’. Such insignificant difference suggests the reliability of the proposed method of estimating the landslide thickness based on InSAR-derived 2D deformation. Above-mentioned difference between the estimated and measured thickness should result from a joint contribution of uncertainties in InSAR deformation caused by atmospheric delay and phase unwrapping errors, as well as the neglect of deformation along the direction perpendicular to the slope. The results estimated with the best-fitting rheological parameter of 0.75 (Figure 3b) show that the derived thickness of the Xiongba-H2 landslide is ranging from 10 to 100 m, corresponding to the active volume of 6.17×10⁷ m³.

5. Conclusions

In this study, a practical method of detecting the landslide thickness based on the 2D deformation derived from two-track InSAR observations is developed and applied to estimate the thickness of the Xiongba-H2 landslide. The proposed method originates from the principle of mass conservation. Considering the challenges of simultaneously acquiring InSAR observations from three or more orbital tracks, our proposed method is more feasible, as landslide thickness estimation relies solely on two-track InSAR observations. The comparison between the landslide thickness estimated by inverting the InSAR-derived 2D deformation and the measured thickness by geological and drilling methods verifies the reliability of our proposed method.

The result shows that two regions with different deformation patterns, including a frontal zone with significant deformation and a relatively stable zone, are identified in the InSAR LOS deformation rate fields. The decomposed 2D deformation indicates that the surface uplift motion at the front edge of Xiongba-H2 landslide could probably be resulted from the pushing of the rock-soil mass in the area with the most significant deformation along the slope direction. The inversion result suggests that the thickness of the Xiongba-H2 landslide is ranging from 10 to 100 m, corresponding to the active volume of 6.17×10⁷ m³. Addition, it suggests that a thicker region at the front edge along the Jinsha River, characterized by steep topography, sparse vegetation cover, and abundant water resources, poses the potential of further failure. The presented approach, which is low-cost and easy to implement, is expected to enhances the availability of InSAR technology for landslide motion analysis and hazard assessment.