Complex dynamics of repeating and river-blocking landslides in Jiangda during 2018

1.   Introduction

Seismic waveform inversion is a convenient way to understand physical mechanisms and kinematic processes of tectonic and non-tectonic earthquakes. Over the past several decades, significant progresses have been made in estimating tectonic earthquake sources. The seismic inversion of non-tectonic earthquakes for both scientific analysis and disaster mitigation has drawn increasing attention. Landslides are one of the most common and destructive disasters that excite seismic signals. During a landslide, friction and other interactions between the sliding mass and ground surface can induce seismic signals, which can be recorded by seismic networks (Kanamori and Given, 1982; Eissler and Kanamori, 1987; Hasegawa and Kanamori, 1987; Kawakatsu, 1989; Fukao, 1995; Allstadt, 2013). By inverting these seismic recordings, the direction and temporal variation of the force exerted on the Earth can be estimated (Allstadt, 2013). From the force history, we can obtain important information of the landslide source, such as the scale, duration, direction, and sliding path, and then estimate some dynamic parameters, e.g., the friction coefficient, based on specific simplified models. So far, seismic inversion is the most effective way to quantitatively describe landslide kinematics and dynamics (Cheng QG et al., 2000). However, landslide inversion has not been performed as widely as earthquake-source inversion. The reason is that landslides are more difficult to be well recorded compared with earthquakes. The energy of landslide waveforms is relatively weaker, making them available at only regional distances. Moreover, many landslides tend to occur in areas with great elevation variation, where regional network stations are sparse, and are generally not sufficient to constrain the landslide sources well.

Southwest China has the most significant terrain differences worldwide; thus, it has frequently experienced severe landslide disasters. Because many streams flow through this area, landslide masses tend to slide into rivers, forming landslide dams and causing large-scale flooding. There have been over 100 landslide dam events in this region, which have caused enormous casualties and economic losses (Chai HJ et al., 1995). The mechanism of these landslides is distinct from that of general landslides. In general, complete river blockages are typically caused by high-speed landslides or collapses (Chai HJ et al., 1998). In a conceptual model of a landslide dam event, the sliding mass moves toward the riverbed after leaving the shear outlet and then stops on the riverbed due to the opposite bank’s obstruction, forming a complete dam body. If the scale is large enough, the frontal material may cross the river, creating a dam body and blocking the river, whereas the rear material stops on the slope due to the dam’s resistance. Once the landslide dam bursts, the rear material may slide to block the river again. However, the above qualitative model for river-blocking landslides has not been confirmed in seismic analysis of actual landslides.

The landslides in October and November 2018 are two examples of typical landslide dam events, which are well recorded by regional broadband networks. The two landslides successively blocked the river at the same location. On October 10, the first landslide blocked the Jinsha River and formed a barrier lake at Jiangda Village, Qamdo (Figure 1). This event was first discovered at approximately 23:00:00 on October 10, 2018 (UTC), and led to over 20,000 people along the river bank being evacuated. Although the barrier lake began to discharge naturally on October 13, the stream was entirely blocked again on November 3, 2018, due to the second landslide in the same location (Figure 1d). The water storage capacity of the barrier lake even reached 300 million cubic meters. Although this water was released through a designated flood relief channel, it still led to the greatest flood ever recorded in Lijiang on November 13, 2018. Before the two landslides, continuous low rate deformation caused small slips in the area (Feng W et al. 2019). Satellite images obtained from Planet Explorer software illustrate geographic changes in this area (Figures 1b-d). Both landslides exhibited a horizontal length of 1.5 km from the top of the slide to the river valley. The horizontal distance and the entire mass-covered area, which ended at the opposite bank, reached 2.1 km and ~1.5 km², respectively. The maximum elevation drop was around 800 m.

Figure  1.  Tomography and satellite images of the study area (a). Red and blue stars, respectively, mark the actual and our located locations of the two landslides. Cyan and green triangles denote stations used to measure the October and November events, respectively, and purple triangles represent stations used for both landslides. The blue line indicates the Jinsha River. Satellite images from Planet Explorer in the landslide area (b) before the October event, (c) after the October event, and (d) after the November event

DownLoad: Full-Size Img PowerPoint

The October and November landslides in 2018 are two repeating landslides. However, existing studies mainly focused on the October landslide (e.g., Zhang Z et al., 2019; Sheng MH et al., 2020). There are no complete descriptions of the physical mechanisms, sliding processes, and causative relationships of both events. Sheng MH et al. (2020) developed a CAPsf method to estimate the single force of the first landslide. Zhang SL et al. (2020a, 2020b) studied the first landslide initiation mechanism, noting that the long-term creep under the exogenic and endogenic integration should be the initiation, while the active ectogenesis is the primary exogenous factor and the influence of serpentinite and foliation is an endogenic factor. To explore the complexities and possible relationships between the October and November landslides’ physical processes, we investigate their force histories through seismic analysis in this study.

2.   Data and methods

2.1   Data

The seismic signals of the October and November landslides were recorded by regional broadband stations in Sichuan Province and Tibet Autonomous Region (Figure 1a). To balance the azimuthal coverage and data quality, we selected 22 stations and 11 stations for the October and November landslides, respectively (Figure 1a). We showed the data of the nearest usable station JSZJ, 12.6 km from the source, to view the data’s frequency spectrum and their power spectrum density (PSD) more closely (Figure 2). Based on vertical displacements, the October landslide’s maximum amplitude is approximately 50 μm, which is nearly three times that of the November event (17 μm). Compared with the October landslide, more high-frequency burrs are evident in the displacement of the November one. However, both events’ primary frequencies remained in the range of 0.01–0.10 Hz. From the PSD, the October event’s primary energy release is highly concentrated within the first 50 s and 0.01–0.10 Hz. Nevertheless, the November event’s energy release is more scattered in the whole time domain, which may reflect its low signal-to-noise ratio at high frequencies.

Figure  2.  Vertical displacements, frequency spectrums, and power spectrum densities (PSD) of vertical components at station JSZJ for the October (a1–a3) and November (b1–b3) landslides in 2018

DownLoad: Full-Size Img PowerPoint

For waveform inversion, the three-component broadband data were first corrected for their instrument responses, and integrated into the displacement waves. Then, the north and east components were rotated to radial and transverse components, respectively. As the energy of the landslide seismograms was dominant at low frequencies, the waves were filtered by a frequency band of 0.01–0.10 Hz. The frequency band used for filtering usually depends on the signals’ dominant frequency, which is often related to the scale of the landslide and its rate of collapse. A lower cut-off frequency ensures that the frequency band is within the flat range of amplitude responses. An upper cut-off frequency ensures that high-frequency signals, which are difficult to model with a one-dimensional layered crustal structure, are excluded.

2.2   Methods

Unlike the hanging wall and footwall of an earthquake fault, the sliding mass and ground surface are entirely decoupled during a landslide. Because the sliding mass is much smaller than the Earth itself, the latter is assumed to be stationary. Thus, the landslide mechanism is a single force acting on the Earth’s surface (Kanamori and Given, 1982; Eissler and Kanamori, 1987; Hasegawa and Kanamori, 1987; Kawakatsu, 1989; Fukao, 1995; Allstadt, 2013).

According to Allstadt (2013), seismic inversion of a single-force source can be described as follows

where U is ground motion, which can be recovered from seismograms; G is the Green’s function, i.e., the ground motion caused by a unit pulse force, which can be calculated; and F is the force-time function (temporal variation of the force), which is unknown to be resolved. Subscripts n and k represent the components of ground motion and force, respectively. The symbol “*” denotes convolution in the time domain.

In Equation (1), five Green’s functions are required for each station in the landslide inversion. As in Allstadt, (2013), these five Green’s functions can be obtained as follows

where g_rr is the radial component caused by a horizontal force in the radial direction, g_tt is the transverse component produced by a horizontal force in the transverse direction, g_rd is the downward vertical component formed by a horizontal force in the radial direction, g_dd is the downward vertical component caused by a downward vertical force, g_dr is the radial component created by a downward vertical force, and θ is the azimuth from the source to the station measured clockwise from the north.

Directly solving Equation (1) will lead to a solution in which the three force components may have different time histories, which is reasonable because the motion path of the sliding mass is not necessarily a straight one. In this study, the Green’s functions in Equation (2) were computed using the Qseis (Wang RJ, 1999) based on the CRUST1.0 model (Laske et al., 2013). Similar to data processing, the calculated velocity of the Green’s functions were also integrated into displacements, and then bandpass filtered using 0.01–0.10 Hz.

3.   Seismic inversion results

3.1   Determination of location and origin time

Although we can precisely locate the landslide through satellite images in some cases, we still need to locate the landslide with seismic recordings. This is because in some cases, landslides may occur far away from the crowd, or satellite data cannot be transmitted in real-time. We must perform quick inversions to estimate the scale, sliding path, duration, etc. For landslide waves, the P and S wave arrival times are too unclear to be identified. In turn, we used a two-dimensional grid search to optimize the location, in which waveform inversion was performed to solve Equation (1) for each latitude and longitude. The data fittings of each inversion were calculated, and the location was obtained from the maximum fitting. For the October event, the location was found at (98.69°E, 31.02°N). The origin time was identified as 14:05:33 (UTC), on October 10, 2018, from the distance-time curve following the procedure described in Allstadt (2013) and Ekström and Stark (2013). Based on the October event’s location, the November landslide’s origin time was identified as 09:21:25 (UTC), November 3, 2018. Our location result is about 7 km away from that revealed from satellite images (98.71°E, 31.08°N).

3.2   Force-time functions

With the optimized location and origin time, the force-time functions of the two landslides were estimated (Figure 3). The correlation coefficients of synthetic data were restored from the estimated force-time functions, and the observed signals for both landslides were all larger than 0.65 (Figures 4 and 5). We divide the sliding process into three stages or sub-events using the complex force-time functions: S1, S2, and S3. S1 and S2 represent the acceleration and deceleration stages in which the vertical force is positive and negative, respectively. After S2, we distinguish S3 events when the forces are weak and the sliding direction is unreasonable. These three-component force-time functions reveal differences between the two landslides in terms of both scale and duration (Figure 3). The peak force value for the October event was approximately three times that of the November event. Moreover, the events’ effective durations were 57 s and 35 s, respectively, when considering only S1 and S2. For both landslides, the deceleration sub-event (S2) displayed greater amplitude but shorter durations than the acceleration sub-event (S1). Particularly for the October event, the S2 sub-event was approximately three times greater than the S1 sub-event in scale, but it lasted only 24 s, which was 9 s less than the S1 sub-event. These relatively strong and short decelerations may imply that the sliding masses did not terminate naturally by friction. Instead, they were blocked or slowed down within a short time period. We conjecture that the sliding mass may have collided with the riverbed, leading to a much stronger deceleration. In contrast, the November event’s sliding mass may have been blocked by the landslide dam, which can explain the weaker S2/S1 ratio compared to that of the October event.

Figure  3.  Three-component force-time functions (FTF), temporal variation of slope angle, and horizontal sliding direction for October (a1–a3) and November (b1–b3) landslides

DownLoad: Full-Size Img PowerPoint

Figure  4.  Observed and synthetic waves of the October landslide. Correlation coefficients are marked on the left, and waves were filtered using 0.01–0.10 Hz. T, R, Z indicate tangential, radical and vertical directions, respectively

DownLoad: Full-Size Img PowerPoint

Figure  5.  Same as Figure 4, except for the November landslide

DownLoad: Full-Size Img PowerPoint

By integrating the force-time functions and taking the arc tangent value of the ratio between the vertical and horizontal momenta, we obtained a series of angles. These angles are only meaningful for the acceleration stage (S1) because deceleration possibly did not occur on the slope. The October landslide exhibited a nearly uniform slope angle of approximately 30° during sub-event S1. In contrast, the November event’s slope angle gradually decreased from a large one (~60°) to approximately 30°. This result suggests that in the November event, mass slid along a curved listric surface during the acceleration stage.

By double-integrating the force-time functions in the first two stages (S1 and S2), we obtained the sliding paths of the two landslides (Figure 3). The October event moved toward the east-northeast, which is consistent with the satellite image observations (Figure 1). In contrast, the November event moved to the north-northeast, which differs slightly from the October event and from the satellite image observations. The October landslide’s sliding direction was nearly constant, but the acceleration and deceleration of the November landslide suggested different sliding directions. This supports the finding that the November event’s sliding surface was more complicated than the October event’s. In addition, the acceleration and deceleration path lengths were comparable for the October event. However, the November event exhibited a longer path for deceleration than for acceleration. This matches the change in slope angles discussed earlier.

4.   Discussion

4.1   Location

Our location result is about 7 km away from the location mapped in satellite images. The location difference may be caused by the fact that we locate the landslide through low-frequency waveform inversions. The velocity of S wave in the target area is about 3.5 km/s, and the minimum period we used is 10 s, making the corresponding wavelength of the S wave be 35 km. According to the Rayleigh criterion, the maximum location difference we can distinguish is 35/4≈9 km. Our location result is about 7 km away from the satellite location, which is less than the maximum location error acceptable. On the other hand, the location difference of 7 km corresponds to a time difference of 2 s for the S wave. This is much smaller than the minimum period we used in the waveform inversion (10–100 s); thus, it would not biasly afftect the dynamic models. Besides, we have made time shift for observed waves in the actual inversion process by searching the maximum correlation coefficients between the observed and synthetic waves of the preliminary inversion. This practice further minimizes the effect of the time difference caused by the location error; that is to say, even if we cannot get the precise location information through ways such as field investigation and satellite imaging, we can still obtain reliable seismic dynamic models of landslides.

4.2   Sub-events

The occurrence of additional sub-events after the deceleration period (S2) should not be ignored, especially for the November event, where a second acceleration period was observed to the east (positive vertical and negative eastern forces) after the deceleration. We also noticed that the vertical components after decelerations are also significant in other studies (Zhang Z et al., 2019), which were thought to be the gradual accumulating process of landslide mass. To demonstrate the contributions of the sub-events to the seismogram fitting, we compare the observed seismograms with the corresponding synthetic waves of the two landslides and their sub-events (Figure 6). To quantitatively verify the reliability and rationality of these subsequent sub-events, we calculated the L2 norms of the sub-event synthetic seismograms and the inversion residuals (Figure 7).

Figure  6.  Comparison between the observed (black) and synthetic (red) seismograms at five stations that were used for the inversions of both October and November events. Synthetic waves of sub-events S1 (orange), S2 (blue), and S3 (green) are also plotted in each sub-graph. The value in the upper left corner of each subgraph represents the maximum amplitude of the observed waveform. The waves were filtered with 0.01–0.10 Hz

DownLoad: Full-Size Img PowerPoint

Figure  7.  L2 norms of the sub-event synthetic seismograms and inversion misfits for the October (a) and November (b) events. The stations are ranked from left to right by their distances from the landslide source

DownLoad: Full-Size Img PowerPoint

Comparison between the observed (black) and synthetic (red) seismograms at five stations that were used for the inversions of both October and November events. Synthetic waves of sub-events S1 (orange), S2 (blue), and S3 (green) are also plotted in each sub-graph. The waves were filtered with 0.01–0.10 Hz.

As illustrated in Figure 6, the synthetic seismograms of deceleration sub-events (S2) are significantly stronger than the inversion residuals. In contrast, the synthetic seismograms of the acceleration sub-events (S1) are comparable with the residuals. The subsequent sub-events (S3) are the weakest, and the L2 norm of their synthetic seismograms is approximately one order of magnitude smaller than that of the residuals. This result indicates that the S3 sub-events contributed little to the data fitting, i.e., they were poorly constrained in the inversion. Therefore, what we can distinguish are just the accelerations and decelerations. In this case, the effective durations of the two landslides were approximately 57 s and 35 s, respectively.

The ratios between the L2 norms of sub-event (S1 and S2) synthetic seismograms differ for the two landslides. As displayed in Figure 7, the L2 norm of the synthetic seismograms of S2 is approximately three times that of the October S1 sub-event. However, for the November event, the corresponding value is slightly greater than twice of the S1 sub-event’s. This consequence is consistent with the result of the force-time functions, and confirms that the two landslides exhibited different terminations.

4.3   Landslide schematics

A more explicit description of the landslide sliding processes is displayed in Figure 8. From an aerial perspective, the October landslide’s sliding area was substantially more extensive than that of the November landslide’s; however, the two horizontal sliding distances were comparable (Figure 1). This indicates that only materials from a small part of the slope slid during the November event, leading to the smaller scale of the November landslide. The October event’s horizontal sliding direction was approximately N75°E, which was consistent with the extension direction of the sliding area (right panels in Figure 1). For the November event, however, the average sliding direction was approximately N45°E. The deceleration sub-event (S2) was directed further north than the acceleration sub-event (S1) (Figure 3). The main reason for this may be that the November landslide sliding body was blocked by a barrier in the southeast corner of the sliding area (Figure 8), and was thus forced to move northward. The barrier marked in red in Figure 8 has been confirmed by field investigations (Deng JH et al., 2019), and can be identified faintly in the Planet Explorer images (Figure 1). The slope angles of the two landslides were also different (Figure 3). The uniform slope angle of the October event implies that the sliding body moved along a straight slope with an angle of approximately 30°, which agrees with the result of field observations (Deng JH et al., 2019). For the November event, however, the slope angle varied with time and decreased from approximately 60° to 30°, which may indicate that the sliding surface was steep in the upper area but flatter in the lower parts; i.e., it had a curved sliding surface (see the cross section in Figure 8). Photographs of the area provide strong support for this feature. For example, a photo taken by Xu DL (2018) reveals a constant sliding angle for the October landslide (Figure 8e). In contrast, a photo taken by Lu X (2018) shows a much sharper angle at the top of the November event (Figure 8f). Previous field investigations also confirmed that the November landslide occurred in the back region of the first landslide source area with a steep angle of slightly over 65° (Xu Q et al., 2018; Deng JH et al., 2019; Feng W et al., 2019).

Figure  8.  Schematics of the October and November 2018 landslides from an aerial perspective (a, b) and along a cross-section (c, d). Turquoise and pale brown areas indicate the sliding areas of the October and November events, respectively. Brown arrows highlight directions from inverted force history. Blue area denotes the Jinsha River. Red area indicates the barrier that may have altered the sliding direction of the November event. (e) October and (f) November photographs are reproduced from Xu DL (2018) and Lu X (2018), respectively

DownLoad: Full-Size Img PowerPoint

Allstadt (2013) proposed that the vertical force-time function amplitude is typically much lower and more susceptible to noise, which reduces the reliability of the vertical component, while landslides’ slipping angles are small. However, the high consistency between the inverted angles and post-field investigation illustrates that the vertical force-time function is reliable when the landslide slope angle is sufficiently large. Combining the inversion results of the two landslides, we can obtain a complete description of the processes involved. On October 10, approximately 75% of the mass slid along a straight slope, collided with the riverbed, and formed a dam blocking the Jinsha River. The river was later naturally released, but the dam was not wholly destroyed. On November 3, the remaining 25% of the mass slid along a curved surface, accumulated on the dam, and blocked the river again. The November event’s listric curved sliding surface may have been caused by the October landslide, which scoured the ground surface along its entire path, leading to higher surface slope angles.

4.4   Friction coefficient

Low basal friction for massive landslides with a long runout has long been a heated debate subject (e.g., Hsu, 1975; Kilburn and Sorensen, 1998). The apparent frictional coefficient was simply derived using the sliding durations and previously calculated slope angles. Using the simplest model, where the entire mass slides from the top to the bottom, we have

where a represents sliding acceleration, µ shows apparent friction coefficient, g is the acceleration due to gravity (9.8 m/s²), θ is slope angle, d represents the sliding distance, and t is duration. The two landslides’ sliding distances are approximately 1.7 km (horizontal and vertical distances of approximately 1.5 km and 0.8 km, respectively). Considering that the durations (t) for the October and November events were 57 s and 35 s, respectively, and their average slope angle (${{\theta }}$) was 30°, we estimate their apparent friction coefficient to be approximately 0.45 and 0.25, respectively. These estimates are comparable with those obtained by Brodsky et al. (2003). The difference in frictional coefficients of the two landslides agrees with their successive occurrences. On October 10, the dominant portion of the mass (75%) slid first, with 25% remaining on the slope. However, the October landslide scoured its entire path, leaving an unstable free surface for the subsequent November landslide (Xu DL, 2018). Rain infiltration may have reduced the frictional coefficient to approximately 0.25, causing an unstable zone on the uppermost slope and finally triggering the second event.

4.5   Limitations of the seismic inversion method

Although we obtained the three-component force-time functions for the 2018 Jiangda landslides through the seismic inversion, this method has certain limitations. First, only the force-time function is obtained, which is the mass plus acceleration over time; thus, the calculated mass and sliding distance have a trade-off. Additionally, the centroid sliding distances are not easily determined. If we use the landslide’s horizontal length from satellite maps as the horizontal sliding distance, which is larger than the actual value, the landslide mass is underestimated. However, because only the driving part of the sliding mass induces seismic waves (Allstadt 2013; Ekström and Stark 2013; Zhang et al. 2019), the inverted mass only represents a small portion of the total mass.

We considered the temporal complexity (force-time) for the seismic inversion in this study; therefore, the dimension was neglected, and the landslide source is assumed to be a point. This practice may lead to oversimplification of the landslide process, i.e., the mass on the slope may exhibit different behaviors at different locations, whereby some parts may accelerate while others may remain stable or even decelerate. However, in a point-source inversion, the inversion recovers only the total sliding mass, which can underestimate the landslide scale. In this case, inversion of the spatial-temporal process of a landslide is required, which is similar to a finite-fault inversion of an earthquake source. This detailed inversion demands a higher frequency in both observations and Green’s functions. Furthermore, because the landslide seismograms are dominated by surface waves, topographic effects should be considered when modeling using high-frequency signals. Also, the sliding velocity of a landslide is typically less than 100 m/s. This is significantly slower than seismic velocities, resulting in extremely weak Doppler effects in the observations. Therefore, it is difficult to use far-field waveform data to constrain the sliding direction and dimension of a landslide. Instead, near-field data should be considered. A densely distributed seismic network around potential landslide areas is required to capture this type of data.

5.   Conclusions

In this study, we estimated the force histories of two co-located Jiangda landslides that occurred on October 10 and November 3, 2018 by inverting regional broadband seismograms. The two landslides exhibited significant differences in scale, sliding direction, and sliding surface. Specifically, the deceleration processes of the two river-blocking landslides were stronger than their accelerations, which may be due to the collisions between the sliding mass and the Earth. Because the October event’s deceleration/acceleration ratio was stronger than the November’s, we infer that the October event was terminated by the collision between the sliding mass and riverbed, and the November event was blocked by the dam formed in the October event. The October landslide brought materials along its sliding path into the river valley and steepened the surface of the mountain slope in the area. Certain unstable elements then slid along a curved slope, leading to the November landslide. A barrier in the southwest corner of the sliding surface may have caused the November event to slide further north.

Acknowledgments

The authors acknowledge Planet Earth for permission to use their satellite data and China Seismological Bureau for permission to use their seismic data. We also thank Prof. Xinghui Huang, Dr. Jiaqi Li, and Mr. Shenjian Zhang for their help in calculating Green’s functions. This work was supported by the National Natural Science Foundation of China (Nos. 42074058, 41822401, and 42021003) and the National Key Research and Development Program of China (No. 2018YFC1503705).