Rigidity Strengthening of Landslide Materials Measured by Seismic Interferometry

Abstract

Landslides have caused extensive infrastructure damage and caused human fatalities for centuries. Intense precipitation and large earthquakes are considered to be two major landslide triggers, particularly in the case of catastrophic landslides. The most widely accepted mechanistic explanation for landslides is the effective-stress dependent shear strength reduction due to increases in pore water pressure. The Chashan landslide site, selected for the present study, has been intensively studied from geological, geophysical, geodetic, geotechnical, hydrological, and seismological perspectives. Our seismic monitoring of daily relative velocity changes (dv/v) indicated that landslide material decreases coincided with the first half of the rainy period and increased during the latter half of the rainy period. The geodetic surveys before and after the rainy period identified vertical subsidence without horizontal movement. The results from the multidisciplinary investigation enabled us to draw a conceptual model of the landslide recovery process induced by water loading. Where all sliding materials were stable (safety factor > 1.0), unconsolidated landslide colluvium and impermeable sliding surfaces trapped the seepage water to form a water tank, provided that compact forces were acting on the materials below the sliding boundary. The vertical force of compaction facilitates an increase in the cohesion and strength of landslide materials, thereby increasing the landslide materials’ stability. We demonstrated that the recovery process periodically occurs only under the combined conditions of prolonged and intense precipitation and the related stability conditions.

1. Introduction

Landslides represent a major geohazard worldwide but the understanding of landslides remains limited. Landslide occurrence depends on several factors, such as the geometry of the slope (i.e., rocks dipping down a slope), the water content, the pore water pressure, the lithology of the material, and external impulsive forces (e.g., earthquakes and volcanic eruptions). With respect to fluid effects, water plays a crucial role in the mechanism of slope instability, decreasing soil cohesion and increasing pore water pressure on a sliding surface. Deep-seated landslides may be triggered by heavy and intense precipitation. Some studies have conducted a series of stochastic analyses to evaluate the rainfall threshold to be able to provide early landslide warnings [1,2,3]. Physical-based modeling can also highlight how an increase in pore water pressure in sliding material generates liquefaction, which leads to rapid mass movement [4,5,6,7]. Most studies have demonstrated that the mechanism of weakening by fluid is an important trigger in landslides. To further our understanding of the mechanisms behind landslides, geophysical, geotechnical, and geodetic monitoring approaches, including electrical resistivity tomography (ERT) [8], active seismic exploration, borehole drilling [9,10], downhole monitoring of displacement and water level [11,12], and GPS-based surface movement [13,14], have been applied not only to clarify the geometry and geological structures involved but also to ascertain the depth of the multiple sliding interfaces located beneath deep-seated landslides. However, some limitations have led to the poor understanding of landslides: (1) Geophysical profiles, such as active seismic imaging, are usually limited to tens of meters of depth due to the insufficient power of the active source and the complex topography changes in mountain area, (2) only few point-based investigations of surface movement and few borehole drillings have been conducted due to high cost and invasiveness. Therefore, an independent, low-cost, and noninvasive approach is needed for monitoring landslide areas and to provide crucial parameters such as the depth of the sliding interface and temporal changes to properties of the landslide medium—these are especially necessary for landslide forecasting.

Landslide seismology has been applied within minutes after landslide occurrences to comprehensively understand the source parameters, such as event location [15,16], source dynamics [17,18,19], and magnitude [20,21]. However, how landslide seismology can be applied to slope failure prediction is another scientific concern. In the past two decades, seismic noise interferometry has become a key tool for retrieving empirical Green’s functions by using cross-correlating continuous seismic noise recordings from surface station pairs or a single station (or both) to attain crucial information for noninvasively probing the subsurface medium. Applications of this seismic interferometry include monitoring earthquake responses in subsurface materials [22,23], understanding volcanic eruptions [24,25], investigating geothermal activity [26], and (potentially) forecasting landslides [27]. A few studies have presented coda wave interferometry (CWI) as an advanced tool for detecting relatively small velocity changes (daily relative velocity changes (dv/v)) in the medium due to scattered waves corresponding to random coda wave propagations through the medium [28,29]. Monitoring medium changes from earthquake seismograms generated by repeat events is, by contrast, only feasible for much clearer changes in velocity [30]. Mainsant et al. [27] first observed drops in seismic velocity before landslide failure and located the possible sliding interface from the frequency-dependent characteristics of velocity changes, suggesting that the CWI technique may be able to predict landslides.

Taiwan is located at the orogenic belt, and its high elevations and subtropical climate result in a mean annual precipitation of 2500 mm, which combines with the frequent incidence of earthquakes to increase the rapidity with which surface mass wasting occurs [31]. The Central Geological Survey, Ministry of Economic Affairs published a landslide inventory of geologically sensitive areas (https://landslide.geologycloud.tw, last accessed on 15 April 2021). Most landslide sites with high potential for slope failure and records for slope failure have been delineated by satellite and aerial orthoimages and field investigations. Knowing potential landslide sites prompts the question of how to effectively monitor a potential incident area to understand its activity and evaluate the triggers for an eventual slope failure. To address these questions, our study selected the slow-moving Chashan landslide area as the study site and focused on combining the CWI technique with conventional geoengineering-based monitoring. Two short-period seismometers were deployed outside the landslide area to apply the CWI technique, which can provide the relative seismic velocity changes (dv/v) corresponding to changes in material properties. Multidisciplinary geoengineering-based techniques were implemented to comprehensively understand the slow-moving landslide. The study approach consisted of the following elements: (1) Constructing the geological model of the landslide by using borehole drilling, field surveys, downhole and laboratory experiments, and geophysical imaging; (2) determining the landslide movement behavior from observed time-series data of surface-to-borehole displacement, groundwater levels, and seismic ground vibrations; (3) understanding the distribution of inferential water flow through hydrogeological modeling, and (4) identifying the possible failure region by adopting scenario tests for each of the trigger factors. Finally, we produced the end-member model to explain the observations and results of the preceding procedural steps. Thus, we achieved a more comprehensive understanding of possible failure mechanisms of the landslide and made progress toward defining certain thresholds for forecasting landslide occurrences. The thresholds of trigger factors proposed in this study should be validated and dynamically adjusted by using the new observations.

2. Construction of Geological Model

2.1. Background

Chashan, a slow-moving, deep-seated landslide area (Landslide Inventory ID: DS160), was selected as the test site for investigating landslide behavior, especially fluid-related mechanisms. The Chashan landslide, which occurred upstream of the Zhenwen River catchment in southern Taiwan, was triggered in 2009 by Typhoon Morakot. The aspect of the Chashan landslide is southwestern. The bedrock of the Chashan landslide comprises Miocene sedimentary layers of sandstone, shale, and interbedded sandstone and shale; this alternation of sandstone and shale is known as the Changchikeng Formation (Cc) [32]. Adjacent regional geological structures are the Shinmei Anticline and Tatou Thrust Fault, which has a north–south strike to the west of the site (Figure 1a). Aerial photographs depict the historical slope failures and agricultural activities that indicate the distribution of colluvium and deposits. Moreover, topographic features related to the deep-seated landslide can be identified by topography and aerial photography and validated through field survey (Figure S1). The Chashan landslide can be separated into the following three subdivisions and topographic features, which are delineated in Figure 1b: the tension crack zone, erosion gullies, and major/minor scarps.

2.2. Geological Survey, Drill Cores, and Resistivity Profile

Several local folds at the Chashan site (Figure 1b) can be explained on the basis of the attitudes of bedding planes, gleaned through field investigation, and the variation of bedding dipping angles from core interpretations. According to the detailed core interpretations of borehole 01 (BH-01), BH-02, and BH-03, the lithological units were divided into five layers, alternating sandstone and shale (SS/SH), massive sandstone (SS2), shale (SH), sandstone interbedded with shale (SS–SH), and gray sandstone (SS1), as well as displaced rock mass (DRM) and colluvium (Table S1). Field investigation and core interpretation were integrated to produce a geological model, and geological maps and profiles were used to illustrate it in a manner similar to that of Yang et al. (2020) (Figure 1c). Additionally, this study incorporated the electrical resistivity tomography (Profile XX′, Figure 1d) profile into the geological model to provide an inferred material boundary between colluvium and DRM. The Wenner and pole–pole electrical array configuration on a survey line with an electrode spacing of 5 m was used in the data acquisition, and the apparent resistivity data were inverted with two-dimensional inversion software (Geotomo RES2Dinv). The software uses finite element solutions and an iterative scheme to invert for the resistivity structure. A detailed methodology is presented in [33]. The EE′ Profile indicates the shallow portion with low resistivity, which is interpreted as the colluvium. The high resistivity areas include DRM, bedrocks, and low impermeable zones. Additionally, the central bottom of the EE′ Profile displays low resistivity, and its shape fits the anticline of SS1; therefore, the resistivity distribution must be affected by the lithology and structure. Profile XX′ shows the spatial distribution of deposits, strata, and slip surfaces in the slope geological model (Figure 1c). The colluvium is exposed on the ground surface of the Chashan site and covers DRM and bedrock. Scarps S1 to S3 (Figure 1b,c) were identified through topographic analysis and field survey, and their connection to the shear material of cores should allow us to infer the potential slip surface (PSS). Accordingly, three branches combine to form the basal slip surface, and the landslide mass can be divided into DRM1, DRM2, and toe slope failure. The main PSS develops from S1 through the shear zone of BH-02 around the east limb of the anticline (anti dip slope); it then continues along the bedding plane of SS–SH with the west limb of the anticline (dip slope) and passes through the shear zone of BH-03 to the shear-off at the toe of the slope. Moreover, the colluvium covering DRM1 and DRM2 might accumulate groundwater due to the low permeability of DRM and because the shear zones of S2 and S3 form the downslope impermeable boundary. As for the deep-seated Chashan landslide, the bedding plane of the dip slope was favored to slip; however, the joints and shear-off surface prevent the slip occurrence.

3. Monitoring and Measurement Data

3.1. Time-Series Data

A landslide site was equipped with a borehole time domain reflectometer (TDR) [34], which collected samples daily, and a groundwater level (GWL) instrument and a rain gauge, both of which had a sampling interval of 1 hour. An inclinometer, the TDR, and the GWL instrument were installed at BH-01, BH-02, and BH-03, respectively. The inclinometer is one of the most commonly used instruments for landslide observation [35]. The present study used manual measurement to confirm the landslide activity, sliding depth, sliding direction, and sliding rate of the Chashan site. A TDR-based landslide monitoring system was employed according to the similar borehole installation method used by [36], wherein a coaxial cable (Belden RG-8) was installed in the predrilled borehole and backfilled with grout to ensure good contact between the cable and surrounding ground. A TDR pulser (Campbell Scientific TDR100) was used to observe the reflection signals received at the oscilloscope terminal after an electromagnetic pulse was transmitted into the grouted coaxial cable. In the event of coaxial cable deformation due to localized shear displacement, characteristic signals would be observed by the TDR oscilloscope, and the corresponding occurrence depth at the sliding plane could be easily determined with a resolution as precise as the centimeter scale. GWL monitoring was performed by installing a piezoresistive transducer-based water level logger (Model dipperLog NANO) perforated its entire length in the monitoring well for water table measurement. During the observation period, GWL ranged from −36 m to −51 m relative to the borehole surface. The measurement data of the TDR and the inclinometer indicated no apparent sliding surface or lateral displacement. However, in the inclinometer’s data, we did observe an irregular S-shaped buckling from the ground surface to a depth of 26 m (Figure S2). This phenomenon often occurs when ground subsidence causes the inclinometer to be subjected to negative friction. For continuous seismic records, we deployed two short-period (SP) velocity-type seismometers outside the landslide region, which was a much more stable location (Figure 1a). The SP sensor (KINKEI KVS300) measured ground motion with a natural frequency of 2 Hz and an 18-bit digital recorder (EDR-7700). A sampling rate of 100 samples per second was used for seismic data. All instruments were used for the 1-year monitoring period.

3.2. Measurement of Surface Displacement

Real-time kinematic (RTK) measurements were conducted a few times throughout the monsoon season, which was very useful for accurately ascertaining landslide surface movement. The measurement points are displayed in Figure 1b. Three measurements were made in 2018: 10 April, 15 May, and 18 September. During the dry season (from 10 April to 15 May), no movement behavior occurred. By contrast, during the period with extremely heavy rainfall, a significant vertical subsidence of 3–5 cm was observed.

4. Methods

4.1. Hydrogeological Conceptual Model

The surface geometry model was created using a 6 m × 6 m Digital Elevation Model (DEM) of the Chashan site. The profile analyzed from the DEM is indicated by the line AA′ in Figure 1b. The cross-section begins at the ridge, intersects with boreholes BH-01, BH-02, and BH-03, and extends to the toe of the site. The model constructed in GeoStudio is presented in Figure S3. According to the surface geological investigation and interpretation of borehole material, the hydrogeological unit at the Chashan site was divided into the following seven categories: colluvium, DRM, SS/SH, SS2, SH, SS–SH, and SS1. The boundary conditions of the hydrogeological conceptual model are also displayed in Figure S3. The left-side boundary (RA) was set as a no-flux boundary for analyzing infiltration and seepage flow because a crest line had already been established. The right-side boundary (SB) was set as a constant head boundary equal to the water table at the toe of the slope and the adjacent drainage. The lower boundary (AB) was set as a no-flux boundary. The surface of the slope (RS) was then set as a rainfall-infiltration boundary.

Based on laboratory tests, the hydrogeological parameters of each layer are summarized in Table S2. Hydraulic parameters were calibrated by comparing them with monitoring data and the results of a steady-state seepage analysis. According to the steady-state seepage analysis, the constant head value on the right side (SB in Figure S3) was 638 m, and the rainfall-infiltration boundary was 6.85 mm/day (unit flux). Furthermore, the unsaturated soil characteristics were considered in the transient seepage analyses using the Fredlund–Xing equation [37].

4.2. Slope Stability Analysis and Scenario

Generally, the stability of slopes is judged on the basis of a computed safety factor [38,39], an assessment of slope deformation [40], or a partial safety factor [41]. Most computer programs used for slope stability analysis are based on the limiting equilibrium approach for two-dimensional models. The method in the present study was to use GeoStudio software, produced by GEO-SLOPE International Ltd., to perform a slope stability analysis of a two-dimensional model of the landslide. In GeoStudio, the SEEP/W module (finite element analysis module) and the SLOPE/W module (limit equilibrium method module) were used. As much of the modeled slope was unsaturated, the safety factor computed by SLOPE/W was based on the Mohr–Coulomb modified equation that was suggested by Fredlund et al. [42]. Changes in pore water pressure and the subsequent effect on the safety factor of the slope were quantified. Transient analysis results of the SEEP/W module of pore water pressure conditions at various points along the slope were input into the SLOPE/W module, thereby allowing highly irregular saturated/unsaturated conditions or transient pore water pressure conditions to be included in the stability analysis. Furthermore, calculating the change in the safety factor over time was possible.

To effectively quantify the effectiveness of GWL simulations, the mean absolute error (MAE) and the mean relative error (MRE) were evaluated in this study, as shown in Equations (1) and (2), where x_s and x_m are the observed and simulated values, respectively. In the present study, when the MRE was less than 15%, parameter calibration was deemed complete.

$$MAE=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|{\left(x_m-x_s\right)}_i\right|}\times 100\%$$

(1)

$$MRE=\frac{1}{n}{\displaystyle \sum _{i=1}^n\frac{\left|{\left(x_m-x_s\right)}_i\right|}{{\left(x_s\right)}_i}}\times 100\%$$

(2)

Additionally, the present study used RAMMS to understand the runout of potential sliding mass obtained from the previous analysis. RAMMS is a reliable numerical simulation tool that yields the runout distance, flow heights, flow velocities, and impact pressure of hillslope landslides, as well as debris flows, to simulate the movement behavior after slope failure [43]. RAMMS adopts the Voellmy–Salm continuity model as a rheological assumption [44]. Parameters used in the RAMMS model are listed in Table S3.

4.3. Coda Wave Interferometry and Stretching Method

The abundance of natural ambient noise, such as human activities (traffic, farming, hydrological pumps), wind, and ocean waves, means that changes of rigidity in landslide materials can be detected using the seismic ambient noise cross-correlation technique. During data preprocessing, all continuous seismic records for the vertical component were divided into one-day data. After removing the instrument responses and mean and linear trends of each one-day time series, we applied band-pass filtering with frequencies ranging from 2 to 20 Hz. For the retrieval of noise cross-correlation functions (NCFs), we first adopted the phase cross-correlation scheme [45], which can yield stable NCFs without time-domain normalization and spectrum whitening, which are generally used to remove transient time-series items, such as earthquake-induced signals. We then stacked all available daily NCFs over the monitoring period to obtain a reference noise cross-correlation function (RNCF).

To measure the relative seismic velocity change (dv/v) in landslide material, we adopted the stretching method [24] to measure the time shift (dτ) between the daily NCFs and RNCF in coda wave windows for the station pair SP1–SP2. The dv/v was then obtained by dividing dτ by −τ. The τ value represents the center time point of the coda wave window used in the stretch analysis. Generally, a coda wave primarily comprises surface waves [24]. As the NCF is derived from the vertical seismic records, the starting time point of the coda time window can be theoretically determined according to the arrival of the slowest Rayleigh waves. The end of the time window is normally the time point at which the coherence is low overall. However, Obermann et al. [46] demonstrated that the coda signals at relatively later times are dominated by body waves, whereas surface waves contribute to the signals at relatively earlier times. Thus, on the basis of the shear-wave velocity (v_s) structure (dashed black line shown in Figure 2a) obtained by the borehole logging investigation, the arrival time of the slowest Rayleigh wave of the station pair can be predicted (the slowest wave speed shown in Figure 2b). To achieve similar waveforms in the coda wave window between NCF and RNCF, the spectrogram of RNCF was calculated using the S-transformation [15], exhibiting the strong spectral energy of coda wave distribution within a frequency band of 2–10 Hz (Figure 2c). The attenuation effect in seismic wave generation due to the complex topographic variations and dense forest between two stations [47] possibly resulted in a weak spectral energy of higher frequency content (>10 Hz). To consider the results of aforementioned analysis and the coherence of daily NCFs, the time window of the coda wave used in the subsequent stretch analysis ranged from 7 to 10 s, which excited the higher value of power spectral density within the frequency band of 2–10 Hz; thus, 2–10 Hz band-pass filtering was applied to all NCFs and RNCF. Figure S3 displays the filtered (2–10 Hz) daily vertical component NCFs during the monitoring period.

The station pair SP1–SP2 has an intrastation distance of 1.2 km that actually samples through the landslide area. However, the lateral resolution along the intrastation path, which can be estimated directly by the radius, R, of first Fresnel zone, is still unclear here. Thus, we further used Equation (3) proposed by Sun and Bancroft [48] to estimate the R-value. Variables associated with Equation (3) are velocity (V), two-way travel-time (t_two) and period (T). The single-way travel-time is picked at the peak amplitude of the envelope of Rayleigh wave coda (coda wave: red rectangles shown in Figure 2c) in causal and anti-causal part. We then calculated the average single-way travel-time, which is about 8.1 s. Based on a distance of 1.2 km, we solve for a velocity (V) of 148 m/s. Spectrogram analysis shows a dominant frequency of Rayleigh wave coda is around 3 Hz, thus, the period (T) is ~0.67 s. The lateral resolution for the station pair SP1-SP2 is approximated as the R-value of 243 m, which is appropriate for studying the dv/v in landslide material of the subdivision B (Figure 1b).

$$R=\frac{V}{2}\times \sqrt{t_{two}\times T}$$

(3)

5. Results

5.1. Model Calibration and Groundwater Level Prediction

The precipitation data from Julian day 195 in 2017 to Julian day 265 in 2018 were used for slope stability analysis. We used the analysis method described in Section 4.1. In the model, rainfall was used as the input data. The GWL at different time steps was obtained by analyzing GeoStudio’s SEEP/W module, and these levels were compared with the actual observation data and adjusted repeatedly until the two exhibited the same trends. The simulation results and observation data are displayed in Figure 3a. Figure 3a displays the differences in the absolute value of GWL, and Figure 3b displays the discrepancies in the amount of GWL change. The simulated and observed GWLs did not sufficiently recover with rainfall during Julian days 195–365 in 2017, and the long-term trend was a GWL decrease. During the rainfall period of Julian days 285–290 in 2017, the GWL rose slightly. Figure 3b shows that the simulation results can indicate the impact of the change in GWL rise.

Equations (1) and (2) were solved to determine that the MAE of the GWL before Julian day 196 in 2018 was 7.7 m and the MRE was 14.1%, which was lower than the acceptable error (15%). Overall, the simulation results had a good reference value before Julian day 196 in 2018. The simulation results’ credibility began to decrease after Julian day 196 in 2018, even when similar rainfall events occurred afterward. Furthermore, the MAE of the GWL after Julian day 196 in 2018 was 8.2 m, and the MRE was 23.7%. GWL did not respond, indicating that some external forces changed the actual formation or physical characteristics. During this time, the initially calibrated model’s credibility was significantly reduced, and retraining was necessary.

5.2. Safety Factors and Failure Scenario

The analytic scenarios used in the slope stability analysis of this study are as follows: (1) the normal condition, under which the GWL was obtained by steady-state seepage analysis, with the effect of seismic force not considered; (2) the 0602 torrential rain event, during which GWL was obtained by transient seepage analysis, with the effect of seismic force not considered; (3) the 0823 torrential rain event, during which the GWL was obtained by transient seepage analysis, with the effect of seismic force not considered, and (4) the earthquake condition, under which the GWL was obtained through steady-state seepage analysis, with the seismic force considered. The peak ground acceleration (PGA) was 280 Gal (cm/s²).

These scenarios considered three inferred potential sliding masses, named B1, B2, and B3 (Figure S5). Additionally, the built-in “Auto Locate” function of GeoStudio was used to automatically search for the sliding masses that were not expected and had the lowest safety factors. The analysis results are listed in Table S4. The analysis results indicated that when the GWL rose, the safety factor of the B3 potential sliding mass at the lower slope was close to 1.0, which indicated a potential for collapse in that area. Therefore, the effects of the failure of the B3 potential sliding mass are worth further discussion. The present study also included a model analysis of the safety factor of the B3 potential sliding mass throughout the observation period, as illustrated in Figure 3b. The GWL of the Chashan site was affected by the 0602 torrential rain event, and the 0823 torrential rain event was associated with a reduction in the safety factor and coincided with a significant increase in the GWL. However, during the period of analysis, the safety factor of the Chashan site was still slightly >1.0.

The present study was conducted according to previous analysis results under specific rainfall and earthquake conditions for the B3 potential sliding mass. Furthermore, we assessed the rainfall- and earthquake-induced landslide threshold at the Chashan site. The rainfall conditions were as follows: a rainfall intensity of 30 mm/h, rainfall duration of 72 h, uniformity in the rainfall pattern, and a total rainfall of 2160 mm. Figure 4a displays the variation curve of the safety factor of a B3 potential sliding mass with accumulated rainfall. The analysis results indicated that the B3 potential sliding mass’s safety factor significantly decreased as the accumulated rainfall increased. When the accumulated rainfall reached 880 mm, the safety factor was reduced to 1.0, which indicated that the B3 potential sliding mass could collapse under such conditions. When the accumulated rainfall exceeded 1000 mm, the GWL reached the ground surface; thus, the safety factor could no longer decrease. As for the set seismic conditions, the PGA was between 0 and 560 Gal. Figure 4b presents the relationship between the PGA and the safety factor. The results suggested that the safety factor of the slope decreased as PGA increased. When the PGA reached 530 Gal, the safety factor dropped to 1.0, indicating that the B3 potential sliding mass could collapse under such conditions.

Above results indicate that the B3 potential sliding mass may collapse when the accumulated rainfall reaches 880 mm or the PGA reaches 530 Gal. Additionally, slope stability analysis revealed that the landslide volume was 53,142 m³. In the present study, this volume was used as the initial volume in the RAMMS analysis, and the remaining input parameters are listed in Table S3. RAMMS analysis results revealed that the maximum accumulation depth was approximately 3.08 m (Figure 4c) and the average accumulation depth was approximately 0.7 m. The deepest accumulation area was located on the stream bed (approximately 480 m long) in the middle of the creek. The affected area after the landslide was approximately 4.4 hectares (calculated based on the accumulation depth being >10 cm). As the accumulation was still 1.8 km away from the downstream Chashan tribe, the destruction of the B3 potential sliding mass was estimated to have had no direct impact on downstream settlement preservation objects.

5.3. Relative Seismic Velocity Changes

We applied the seismic interferometry technique and the stretching scheme to continuous seismic recordings that had daily temporal resolutions. Our results revealed that the relative seismic velocity variations (dv/v) from a station pair with a trajectory crossing the landslide area can be measured within ±1.0% variance; this indicates the background level of dv/v measurements within ±0.2%, except for the period affected by external environmental forces, such as earthquakes, water penetration, and/or property changes to the basal sliding interface (Figure 5a). Overall, we observed a significant decrease in coherence during the rainfall events. Therefore, we propose that the loss of coherence is related to the changes in the subsurface’s medium [26,46] (e.g., the scattering properties of the medium) (Figure S7). However, changes in the local noise sources can also partly contribute this de-coherence. Bontemps et al. [49] highlighted the temporally combined effect between earthquakes and precipitations on slow-moving landslide movement. For example, they proposed the cumulated forcing of the frequent small earthquakes combined with a high water content in soil materials that prevents the landslide from healing its rigidity. In the present study, the environmental force-related dv/v changes can be characterized into three periods: several rainfall episodes without earthquake activity (period 1, P[1]), large effect of earthquake forcing (period 2, P[2]), and prolonged and intense precipitations (period 3, P[3]) (Figure 3a) without earthquake effect, which would be very helpful to investigate the effect of individual forcing in dv/v observations. In P1, an obvious pattern is the typical cycle of reduction and recovery (RAR) in dv/v in response to material changes related to water infiltration and in-filled water running out (Figure 3a). Comparing daily dv/v observations with changes in the GWL responses to precipitation penetrating the landslide material indicated a velocity drop could be caused by rises in GWL under wet weather conditions. By contrast, a falling GWL results in an increasing dv/v at a specific recovery rate (arrows in Figure 3a). For example, after an intense precipitation event accompanied by the predicted water elevation increase (P1 in Figure 3a), the velocity apparently recovered by 1.2% in 30 days. A series of RAR (dashed arrows during P1) was strongly sensitive to rainfall episodes that had <100 mm cumulative precipitation. We suspect that the recovery rate of dv/v was mainly controlled by the drainage properties of aquifers and the intensity of a preceding rainfall event.

During the dry season, as indicated in P2 (Figure 3a), the possible effects of rainfall penetration on the observed dv/v could be neglected as the significant earthquakes occurred in P2. Indeed, a clear scatter pattern in dv/v measurements is evident. Studies have suggested that seismic velocity variations may be caused by strong ground-shaking from earthquakes [22,49,50]. In this study, we obtained the Central Weather Bureau’s (CWB) earthquake catalog to understand whether significant ground motion induced by regional earthquakes could have affected the landslide material. During the monitoring period, seismicity within 150 km of the landslide site was limited to the local magnitude (M_L), which ranged from 4.1 to 6.2, and a focal depth between 5.0 km and 35.6 km (Figure S6). A total of 8 earthquakes with the seismic intensity higher than II (intensities III and IV on the CWB intensity scale, corresponding to the PGA ranges of 8–25 and 25–80 Gal, respectively) were reported by the CWB earthquake center. In our case, we suspect that the frequent strong ground-shaking caused dv/v measurements to tend to scatter. This coseismic dv/v pattern can probably be explained by the earthquake-forced redistribution of the cracks and fractures embedded in the landslide body. We propose that this nonlinear structural weakening of the colluvium layer and the rock materials can explain the observed GWL fluctuations that were caused by the hydraulic property changes, resulting in a change in ground water conditions (e.g., change in preferential water infiltration paths) (P2 in Figure 3a) [49,50]. In Taiwan, the annual monsoon/typhoon season is generally accompanied by prolonged and intense rainfall between June and September (Julian days 152–273, 2018; P3). In our case, a cycle of RAR in dv/v exhibited a longer cycle duration of approximately 80 days during P3 that coincided with frequent rainfall episodes. We further noticed that the dv/v increased with intense precipitation, but the response of the observed GWL was conspicuously uncorrelated with rainfall data. A large discrepancy in GWL change rates between modeled and observed data was also observed, which implied the occurrence of temporal changes in rainwater infiltration properties (Figure 3a).

6. Discussion and Conclusions

During P3, we did not expect to observe a clear increase in dv/v amid heavy precipitation at this landslide site. To investigate the mechanisms that fundamentally cause dv/v increases during high precipitation rate and intense rainfall, we first located the possible sources of dv/v around the landslide area, which corresponded to the changes to the bulk elastic properties of landslide material (Figure 5a), such as the GWL [24], ground-shaking [22,49,50], and subsurface pressure variations [25]. Finally, a water-load-based compacting model was proposed to support our observations and results.

6.1. Factors Influencing Daily Relative Velocity Changes

Recent studies have highlighted the continuous decreases in dv/v along the sliding interface for a few days prior to landslide failure [27]. Obermann et al. [46] noted that frequency-dependent dv/v measurement could be a powerful tool for locating the depth of velocity anomaly. Based on the available PS-logging shear-wave velocity (v_s), the assumption of coda waves dominated by Rayleigh waves, and the Computer Programs in Seismology package [51], a layer of sliding material of BH-03 (a possible shear zone at a depth ranging from 51 to 54 m) with a velocity drop (ΔV) was tested to compute the relative surface-wave phase velocity differences between PS-logging v_s and perturbed velocity models (Figure 2a). In a synthetic test involving 3 m thickness (H) and 67% ΔV (wherein initial velocity = 1.5 × perturbed velocity), a peak of phase velocity drop of 2.7% was observed at a frequency of 5 Hz (Figure 2b). We calculated the Rayleigh wave sensitivity kernels as a function of depth for the fundamental mode, which resulted in sensitivity depths ranging from a few meters to 60 m at 2–10 Hz frequencies; in our calculations, we sampled both the response range of the GWL and the sliding material. However, no shear sliding/deformation was observed at borehole displacement sensors, such as the TDR and the inclinometer, during the monitoring period; thus, we concluded that the 2–10 Hz dv/v measurement was mainly caused by the GWL changes and bounded within ±1.0% variance. However, only GWL-related dv/v fluctuations could not explain the abnormal cycle of RAR observed during a series of intense and prolonged precipitation events (P3 in Figure 3a), therefore, a different mechanism must be involved.

6.2. Evidence to Support the Water-Load-Based Compacting Model

We proposed a hypothetical model based on the water-load compacting force, and it is displayed in Figure 5b. During Julian days 150–198 in 2018 (Stage 1 during P3; Figure 3a and Figure 5b), the GWL exhibited a strong correlation with the precipitation rate but had a delayed reaction to the rainwater inputs, coinciding with a gradual velocity reduction over 50 days. At Stage 1, waterflow from rainfall episodes penetrates through the landslide material, and the GWL at the top of the main shear zone increases (sliding surface is in the water table), resulting in the effective normal stress (σ_n) decreasing. After Julian day 198, Stage 2 of P3 (Figure 3a), an obvious out-phase correlation between observed GWL and precipitation rate was observed, which implied that the hydraulic conductivity had probably already changed. The response of measured GWL was less sensitive to rainfall, resulting in large discrepancies in the GWL rate compared with modeled results. We suspect that water trapped in the landslide’s colluvium-formed tank compacts the material beneath the sliding interface which could have reduced the capability of rainwater penetration and led to increases in dv/v accompanied by groundwater running out, cohesion, and σ_n increasing. In the aforementioned case, the dynamic closures of cracks and fractures would enhance the rigidity strength of landslide material. A hypothetical model for a water load forcibly inducing dv/v increases is presented in Figure 5b. In summary, the formation of a water tank during the rainfall/typhoon seasons requires satisfaction of the following four landslide material preconditions: (1) The impermeable interface (shear zone) acts as an aquitard; (2) preexisting landslide colluvium located (aquifer) at the top of aquitard is required to trap the rainwater; (3) the safety factor for the storm period is >1.0, and (4) materials underneath the aquitard boundary are compressible. The first two preconditions can be clarified by the geological model. During the monitoring period, the triggering factors did not reach the thresholds of landslide failure that we proposed in our stability analysis. The uniaxial compression strength test of BH-02 in laboratory testing indicated that the axial failure strain of the SS–SH rock sample was 1.0%. The results indicated that although the rock material was brittle, the compressibility of the rock material was still sufficiently high. The landslide site in the present study satisfied the aforementioned four preconditions; thus, the increases in the rigidity of the landslide body being caused by a water-related compacting force was strongly supported by not only the dv/v results but also the RTK (Figure 3b) and inclinometer measurements (Figure S2). RTK measurements revealed a vertical subsidence with a range of 3–5 cm without horizontal movements in response to extreme rainfall events. As for practical applications, loading the top of a slope may significantly influence stability; thus, water loading should not be at the top of a slope. Notably, during a late rainfall period with a heavy precipitation rate (Stage 3 in P3, 290 mm/day displayed in Figure 3a), the rainwater penetrated into the vicinity of a high-permeability region and rapidly flowed in and out, resulting in the impulse-wave pattern of observed GWL. This rapid return of GWL could possibly be related to deep groundwater from a specific flow path, and the combined effect of GWL increases and entrained subsurface water loading would lead to the absence of a concurrent dv/v drop (Stage 3 in Figure 5b). Truly discriminating the details of the flow path is difficult due to scant independent data. After precipitation peaked, the dv/v exhibited a gradually increasing trend (Stage 4 in Figure 5b).

6.3. Implications

To conclude, we applied seismic interferometry to monitor the temporal velocity fluctuations in over a year’s worth of continuous ambient noise recordings, which responded to seasonal precipitation, typhoon activity, and seismicity in Taiwan. At the landslide site, the material rigidity increases being related to the surface water loading explained the approximately 0.5% increase in dv/v that concurrently occurred with a series of extreme rainfall events (Stage 2 of P3, Figure 5b). A temporal change in hydraulic conductivity, caused by vertical compacting force, was proposed; we also proposed that the compact process induced by water forcing could yield a better understanding of the fluid-related landslide mechanism and cause the rigidity of landslide material to tend to strengthen. We can also conclude that the use of temporal changes in hydraulic parameters in numerical modeling is urgently required for studying landslide sites in a world with relatively intense and prolonged rainfall events.

The possible strengthening of landslide material can occur periodically during intense rainfall events, causing the landslide to tend toward stability, which is one reason why the thresholds of trigger factors for early landslide warnings should be dynamically adjusted. In our case, an estimated velocity reduction of 2.7% was related to basal sliding over a give thickness (H = 3 m) and could potentially be used as an additional and independent threshold for early warning of landslide failure. However, a certain threshold of seismic velocity reduction for our landslide area is needed to further investigate. For practical applications, the combination of a landslide seismology technique and hydrogeological modeling could prove feasible in understanding the temporal changes in hydrological parameters, which would help in recalculating the thresholds of landslide trigger factors to better predict landslides. In the present study, a cumulative rainfall of 880 mm and a PGA of 530 Gal were considered to be the lower bound of triggering thresholds for landslide forecasting, since aforementioned trigger factors should be recalculated after observed material strengthening.