Earthflow reactivation assessment by multichannel analysis of surface waves and electrical resistivity tomography: A case study

2.1 Geomorphological and geological setting

The Tonghua landslide occurred near the Xishan village of Tonghua, Northwest Sichuan Province, China. The typical elevation is approximately 2,300 m, with high mountains and deep valleys in the case-study region. Denudation of vegetation in the region is very high, and the terrain is mostly alpine [24]. A portion of the interim belt extends from the eastern edge of the Qinghai–Tibet plateau to the east Sichuan bowl. Historic earthquakes strongly influence the slope structures of the area; thus, the landslide risk factors are high and may cause safety issues in the future. The most active slope segments in the area are the bulged portions, ranging from 30 to 40°. Ancient superficial deposits on the surface of the slopes include the Devonian system and Quaternary deposits. The basis is generally Devonian crystals containing phyllite and Holocene colluvial deposits [24].

The landslide body is, thus, stratigraphically associated with the collection of Quaternary glacial deposits. The fourth and fifth parties of the Silurian Maoxian establishment control visible divisions in this area (Figure 3), which mainly include metamorphic pillars and diverse deposits of loose sediments. The lithology comprises limestone, sandy slate, phyllite, and metamorphosed stone. The varied deposits of loose sediments are predominantly categorized into three: (i) glacial deposits, (ii) alluvial deposits, and (iii) colluvium. Glacial deposits form a unique residue landform in this region.

Figure 3

(a) Regional geological map of the study area with dominant lithologies, and (b) lithology portion of the research region and adjacent regions acquired from wells [24].

Different geotechnical reports have described the lithology of the various layers. The first layer comprises pebbly clay layers composed of breccia and gravels ranging from 2 to 20 mm. This layer contains 75% clay, making it highly plastic. The second layer, which is gravelly soil, consists of more than 50% gravel ranging from 10 to 50 mm and sand ranging from 0.5 to 2.0 mm. The parent rock is phyllite, and the weathering rate is high. The third layer of the Zagunao River deposits comprises Quaternary alluvial sediments. Large-scale joint fissures with potential flow paths and transportation were first reported in ref. [25], and a dendritic pipe-network seepage system was later identified and categorized in the region [24]. Figure 3(b) shows the borehole lithological log according to previous drill boreholes [24]. The shallow subsurface succession consists of Quaternary diluvium, gravel clay soil, gravel soil, solidly endured phyllite, weathered phyllite, and phyllite bedrock respectively from the top to bottom.

Structurally, the study area is near the Longmenshan fault belt. As shown in Figure 4, the tectonic block is very active and has caused several earthquakes [6]. The China Seismic Network database shows that from 2012, the region within 30 km of the Tonghua landslip has been influenced by more than 28 seismic incidents of surface wave magnitude (Ms) ≥3. Notably, the integrity of the collapsed area was brutally broken by the 2008 Wenchuan Ms-8.0 earthquake as the locale of the landslip is only 66 km from the epicenter of the earthquake. Furthermore, there has been a giant prospective landslide, specifically the Xishan County landslide, in the locality of the Tonghua landslip for many years (Figure 1).

Figure 4

Locations of structural faults and seismic events around the Tonghua landslide modified from ref. [6].

2.2 MASW

MASW is an effective technique for high-resolution examination of surface and subsurface formations using functional resources. It has been used to tackle several tasks in geotechnical and environmental engineering, such as providing a planning basis [26], defining narrow fault regions [27,28,29], sensing cavities [30,31,32,33], seismic microzonation, site reaction experiments, and geotechnical characterization of substances close to the surface [34,35,36], and evaluating landslide strength [37]. Rayleigh wave signs were amassed using a single-ended irritation source layout and a collecting structure with geophones at points of different topographies. Inversion diffusion curves were computed using the frequency–wavenumber (F–K) assay to establish the variation of shear-wave velocity outlines with intensity. The probe strength of MASW is approximately related to V _s/3f, and the shear-wave velocity is susceptible to landslide distortion [35]. Recently, this assay has been used to identify two-dimensional (2D) and/or three-dimensional (3D) landslide structures [38,39].

The acquisition unit consists of a seismograph controlling 24-spike geophones with a frequency of 4.5 Hz and 2 m spacing, and a 6 kg iron hammer was used as a seismic source to generate wide-bandwidth exterior waves (Table 1) [40]. As shown in Figure 5, four (shear-wave velocity profiles of landslide debris SW-1, SW-2, SW-3, and SW-4) and an adjacent stable section were obtained using the surface wave method [26,34]. With the parameters used in this survey, shear-wave velocities at approximately 30 m depth could be obtained. Although the terrain in this region is rugged due to episodic landslide events, we planted the geophones in almost parallel and straight lines to ensure favorable conditions for the MASW survey. Special consideration was given to the parameters suggested in [41,42,43] to improve the data quality and avoid near-field effects. Data processing was finalized using the “SurfWave® 2020” following three main steps [37,44,45]. First, the dispersion curves were estimated by transforming each Rayleigh-wave shot fold (Figure 6a) to the frequency–phase velocity (f–V) domain for dispersion curve extraction (Figure 6b) [46]. Second, a one-dimensional (1D) S-wave velocity (V _s) profile (Figure 6c) was projected through the inversion process from the dispersion curve using a reiterative reversal algorithm functioning with the least-squares method inversion process to run automatically [40]. Similarly, an initial model, V _s, was generated by a mixture of numerous factors. The original model comprises 10 units, ensuring optimal resolution by iteratively using V _s.

Figure 5

Digital elevation model of the Tonghua landslide overlaid by satellite view and locations of geophysical ERT and MASW data acquisition lines.

Figure 6

(a) Rayleigh-wave shot fold from SW-3, (b) the Rayleigh-wave dispersal image produced from (a), (c) experimental (red dotted line), theoretical (blue dotted line) dispersion curves, and (d) Shear wave velocity (S-velocity) construction using MASW data S-velocity estimate is exact down to the depth of about 26 m.

In the preliminary model, V _s values and depth of access (Z _f) are chosen by assuming that V _s at a depth of Z _f is 1.09 times the measured phase velocity (C _f) at the frequency where wavelength λ _f [34], as expressed in the following equation:

where a is a coefficient that varies only slightly with frequency and is the result of extensive modeling.

Therefore, the maximum access strength imaged using the MASW method was usually monitored by the smallest regularity encountered while generating the initial model. Third, the MASW 1D shear-wave velocity outlines were gathered along with the study lines and then incorporated using the kriging algorithm to create 2D color grid shear-wave velocity profiles [26,31,47,48,49].

According to the National Earthquake Hazards Reduction Program (NEHRP) standards [50,51], V _s of the subsurface must be measured or estimated to a depth of approximately 30 m (V _s30) to study earthquake risks. It can be calculated using the following equation:

where V _s30 is the NEHRP SWV (seismic wave velocity), V _sf is the shear-wave velocity of any layer in m/s, and d _i is the thickness for any layer (0–30 m).

Table 2 lists the details of the landslide and its soil conditions based on the NEHRP standards [50]. The soil is classified into six types (A–F) based on their V _s30.

2.3 ERT

The ERT process produces an image of the allocations of highly saturated and clay-rich areas, showing low resistivity in the landslide, rock weathering, and fracturing [14,52,53,54,55]. Resistivity varies with the mineral elements, groundwater, and electrolytes in the subsurface layers [56]. In this method, a DC electrical current is injected into the subsurface. It runs through conductive media and is received by electrodes planted within the geophysical line [24]. 2D and 3D resistivity tomography methods are used to discover unstable slopes because they provide exact data on the electrical resistivity divisions and physical assets of geological structures [18,57].

For the Tonghua landslide, the electrical resistivity data were acquired using a multielectrode system (56 electrodes) with a SuperStingR8 (AGI-instruments) resistivity meter connected to a multicore cable. In this approach, ERT was used in different array lengths using a dipole–dipole electrode configuration to find variable investigation depths.

Several studies have been conducted using a dipole–dipole array for landslide investigation [13,58,59,60]. The main advantages of a dipole–dipole array include the following: (1) it is suitable for complex and deep-seated landslides; (2) it can determine isometric structures; and (3) it is sensitive to horizontal changes in resistivity [61]. The depth of this investigation was approximately 50 m. Data noise was evaluated using the stacking procedure [62]. Five to ten stacked electrical resistivity measurements were conducted for each point, and the resistivity meter automatically estimated the respective relative standard deviation. Resistivity was calculated using a standard deviation greater than 2%, and noticeable outliers were not considered.

Planning and arrangement of 3D electrode array operations on landslides at the desired depth were complicated and expensive. Although the sliding material was soft, the material distribution was not wide enough to set up a 3D array. A cost-effective alternative to true 3D ERT imaging is the pseudo-3D ERT imaging technique. This involves obtaining parallel and/or orthogonal surface 2D ERT profiles and processing them as though they are true 3D data [63,64].

Hence, seven 2D high-resolution survey profiles, including two longitudinal and five parallel transversal lines spread out on the landslide body, were obtained to investigate the landslide structure (Figure 5). The array length was 132–240 m with 4–5 m electrode spacing to increase the depth coverage, and the investigation depth ranged from 26 to 56 m (Table 1). Finally, the five traversal 2D sections were combined with reasonably accurate images of 3D structures. A statistical analysis of the determined apparent resistivity was performed after data collection to establish the frequency distribution.

The 2D ERT datasets were inverted using ResIPy software [65]. The code for resistivity inversion based on an objective regularized function was integrated with a weighted least-squares evaluation of the misfit (the solution type of Occam) after [66] adjustment for the imposed regularization constraints. To assess the accuracy of the numerical solution, a homogeneous earth structure was used to simulate forward models for different mesh discretization. The 2D ERT data inversion is a nonunique problem that needs constraints on the model structure to minimize possible solutions. An objective function connected to a linear system of equations was minimized to solve the tradeoff between the model constraints and the weighted data misfit [67,68].

Several factors can affect the imaging depth and resolution of the inverse model of electrical methods, including the amplitude of the signal relative to the conditions of ambient noise [69], parameterization type and regularization in the inversion, and geometry of the cables (number of segments, spacing between electrodes, number of electrodes, and type of array) [68].

The structure beneath the depth of investigation (DOI) should not be interpreted geologically. Therefore, estimates of the DOI for DC resistivity surveys are crucial when interpreting models obtained from any inversion. A statistical method developed in ref. [70] was employed to calculate a modified DOI index [71]. Two inversions with homogeneous reference models of different resistivity (0.1 and 10 times the logarithmic mean of the apparent resistivity values) were performed herein. To fit the values of the cumulative distribution of the calculated DOI index, the sum of two normal distributions was used, thus allowing the interpretability index to be calculated for each cell of the model. Another inversion was then performed with the mean of apparent resistivity (a third homogenous reference model). The interpretability index was an “alpha blending” step. To distinguish between poorly and well-constrained cells, the resistivity model was plotted specially.

3.1 Deformation properties and landslide characterization

As shown in Figure 7(a) and (b), two cases of 1D V _s profiles were generated from MASW data (profiles SW3 and SW4) and correlated with existing boreholes using interval velocity functions. The boreholes near the survey lines were drilled by ref. [24], who investigated the Xishan village landslide. The two examples show a good match of the total rapidity from MASW with the stratigraphic logs. The contrast reveals that the top of the bedrock exterior has a shear-wave rapidity of more than 450 m/s. A comparison of MASW-1D profiles helped us to read the 2D V _s sections and debris depth.

Figure 7

Matching 1D velocity profiles obtained from the MASW data to borehole logs situated in the investigated area. (a) 1D velocity profile SW4 was obtained from the stable slope section, and (b) 1D velocity profile SW3 was obtained from landslide debris.

The calculated velocity-depth 2D profiles of the earthflow debris showed a debris depth up to 30 m (Figure 8). After implementing a spatial interpolation approach, we could construct V _s slice maps at different depths [45]. The depth of the bedrock was obtained by comparing the V _s maps with borehole records (Figure 9). Furthermore, the low V _s encountered in the top 5–10 m of the debris indicates soft deposits or high clay content debris. Beneath this layer, the next horizon was a combination of high and low velocities, indicating several earlier episodes of landslides. The second layer exists in all three images (Figure 8(a), (c), and (d)) with high- to low-velocity contrasts. Visual interpretation confirmed the existence of small cavities and high weak zones, causing a collapse under worsening conditions and initiating debris movement. A relatively uniform layer of V _s (600–800 m/s) penetrating at a depth of 24–27 m was interpreted as the termination of debris deposits. It is considered a new slip surface for the existing debris.

Figure 8

V _s 2D sections with the survey profiles. (a) Top of the debris (SW1); (b) at the middle (SW2); (c) at the trailing edge of the landslide debris (SW3); and (d) stable side beside the landslide site (SW4).

Figure 9

(a) 2D S-wave velocity sections projected onto Landsat image of Tongua landslide site (for locations, see Figure 5) and (b) S-wave velocity slicing maps illustrate the distribution of S-velocity at different depths.

The encountered rock masses have irregular dispersion with interactions between hard and soft rocks in the stable section. Generally, the fractured mechanical conditions of the geological stratum in this profile are vital for the overburden failure mechanisms. A combination of soft and hard rock layers and their residual bulk characteristics is evident in 2D sections. The overlying layers directly influenced the mechanical strength, failure mode, and residual bulk characteristics of the underlying rock stratum, as deduced from the visual interpretation of the lithological structure. Four geological layers were encountered in the study area based on the 2D sections after comparing with existing borehole logs from nearby areas and interpretation. The top colluvium layer is heterogeneous and contains rock fragments to a depth of 6–7 m. The next layer is 7–8 m thick and comprises altered gravel beds mixed with pebbles, and its velocity structures indicate a discontinuous depositional structure. The third layer is a low-velocity layer inserted between two high-velocity layers and depicts low velocities (150–200 m/s), representing high weathering conditions. This layer is interpreted as weathered phyllite that is 22 m deep. At the bottom, we found rocks with a velocity of 800–1,000 m/s. This is considered a bedrock for the examined slope.

According to the NEHRP standard, V _s30 and soil types are classified into two soil classes (classes D and C) (Table 3). Soil class D was found in the landslide debris stiff soil, whereas class C was found mainly in the stable part of the slope as very dense soil and soft bedrock.

3.2 2D ERT profiles

Several 2D dipole–dipole measurements were conducted on the landslide debris and adjacent stable section. The array covering the stable and debris sections were 132 and 240 m long, respectively. The resistivity inversions achieved optimum estimations in seven repetitions, with less than 4% root-mean-square error. Figure 10a illustrates the 2D tomography ERT-6 profile performed longitudinally above the landslide debris. The debris basement at a depth of 40 m was highly conductive and extended up to 100 m. In the basement, the encountered values were less than 30 Ωm, making the conditions favorable for water saturation. Over the basement, the landslide-induced overlaid sediments contain low-to-high resistive materials of 80–250 Ωm. The debris-flow deposits showed an average thickness of approximately 20–25 m and were thicker in the central part. Intermittent areas with a resistivity of less than 20 Ωm were also found below the freshly accumulated sediment. These anomalous resistivity regions showed the locations of hidden clayey soil beneath the debris sediment. Water infiltration and fissure water runoff in the debris structure were also prominent. The DOI along this section ranged from 20 m at the edges of the section and line ends to 50 m at the profile center (Figure 10b).

Figure 10

(a) ERT-6 profile crossing the main body of the landslide in the center with a spread length of 240 m and (b) DOI sensitivity sections show the adequate investigated depth.

We obtained the 2D ERT profile of the slope adjacent to the debris to understand its failure mechanism. This part of the hill was at the rear under a high debris load but did not move during the debris flow. The ERT-7 profile interpreted this relief as a bedrock outcrop with a resistivity greater than 250 Ωm (Figure 11a). After a 28–30 m-long continuous high-resistivity zone, a 90–95 m-long low-resistivity zone was observed. The low-resistivity region contained intermittent high-resistivity patches, which may be lenses of gravel or weathered phyllite. The stable relief also showed high water content at the toe end and can act as a slip initiator. Electrical resistivity reflects different factors, including the clay content porosity and water content. In general, different types of rocks have a wide range of electrical resistivity. Therefore, the values of electrical resistivity cannot be interpreted based on lithology. Thus, we compared our results with the literature [4,15,16] and nearby exploratory boreholes in the region. Based on the resistivity contrasts, three electrical resistivity ranges were assigned in the stable relief: (a) ρ > 20 Ωm to surface water infiltration and fissure water runoff; (b) ρ = 1–8 Ωm to the water-rich toe of the slope; and (c) ρ > 250 moist gravel and weathered phyllite. In this case, the active landslide material acts as an overburden beneath low-resistivity strata (5–200 Ωm). The active (lower) portion of the relief is classified into two extensive zones. The upper zone comprising displaced, relatively intact blocks (rafts) of bedrock, while the transitions into the lower zone indicate predominantly fractured materials with numerous closely spaced tension cracks and water-rich strata. On this basis, the new susceptible landslide may be classified as a block slide [72] or complex translational rockslide earthflow [73]. The DOI along this section ranges from 18 m at the edges and end lines of the section to 25 m at the profile center (Figure 11b).

Figure 11

(a) Longitudinal ERT-7 profile, with a length of 124 m at the stable landslide section and (b) DOI sensitivity sections show the effective investigated depth.

3.3 Pseudo-3D inverted ERT model

To explore the influence of draining water out of the landslide debris, five 2D profile measurements were conducted and combined to construct pseudo-3D ERT using the obtained five parallel ERT profiles (Figure 12). The electrode spacing was 4 and 35 m in the X- and Y-axis, respectively. The 2D surveys had deeper exploration depths, lower noise, and simpler array geometry due to the longer array. Specific details about the conductive channels under the debris could be obtained from 3D measurements. Notably, the 2D array data for the inversion were used to construct a model that could summarize the connections of oversaturated drains (Figure 12). As shown in Figure 13(a), the locations with resistivity greater than 250 Ωm represent a lens of stiff deposits with low conductivity, corresponding to the expansion of previous debris movement. The 3D resistivity contour plot shows that the debris consists of three mixed layers: a stiff (dipping) layer with a high range of resistivity; a soft layer (the interim area) with a medium-resistivity range; and a low-resistivity layer, which is a highly saturated stratum. Figure 13(c) shows the 3D volume extraction of resistivity below 30 Ωm and outlines a prominent conductive zone beneath the hard deposits. The clear separation of these two resistivity volumes confirms the water flow paths and results from the perching of rainwater draining through the debris. The unconsolidated layer above the saturated layer exhibits medium resistivity of 80–190 Ωm. A medium-resistivity volume surrounds the phyllite deposits at the interface with compaction. The drilling record confirms that the resistivity of pebbly clay in this region varies from 50 to 200 Ωm [24]. This unconsolidated stratum is a mixture of slit, pebbly clay, and colluvium. The high resistivity is widely distributed within the medium-resistivity region (200–250 Ωm). These deposits cover the slope until the trailing edge and make them vulnerable to the steepness of the slope.

Figure 12

Pseudo-3D inverted ERT model built from combined 3D inversion of parallel ERT-1, ERT-2, ERT-3, ERT-4, and ERT-5 profiles, crossing the debris transversely, almost 25–30 m below the topographic surface. The high distribution of low-resistivity values (blue color) shows the pathways created by infiltrating water.

Figure 13

(a) Represents the inverted results by combining 2D parallel ERT lines to create a 3D slice in the X and Y axes with reasonable RMS error; (b) the 3D inverted image represents the overall status of the debris with core leaking zones of low-resistivity (<20 Ωm) and high-resistivity (>250 Ωm) and represents the volumetric contrasts between extremely hard and saturated debris strata and the volume of debris that is highly saturated with water.

3.4 Landslide development and hydrogeological model

In this study, the spread length of MASW was 46 m, which could limit the shear-wave velocity estimation to a particular depth. The obtained MASW profiles have soil deposits approximately 30 m thick. The landslide was characterized as a zone covering class D bearing, very dense soil, and soft rocks. A visual inspection of the MASW survey revealed well-defined soil and rock characteristics when mapping the shallow 2–30 m deep. At this depth, the shear-wave velocity was attained with fair accuracy in signal processing by preserving the low frequencies (high wavelengths) for deeper penetrations. The geotechnical characteristics of the study area can be used for preliminary design of new buildings. However, laboratory experiments are needed to estimate the friction angles and deformation moduli at certain locations, guided by the conventional relationship between the seepage field and landslide stiffness through the MASW and ERT characteristics. We developed a schematic model of the landslide main body and its hydrogeological features in the study area based on the integrated results of ERT and MASW. The evolution of the landslide was consistent with the hydrogeological features (Figure 14).

Figure 14

A typical situation of the Tonghua Village landslide, the best-fit model, developed after the geophysical investigation, with prominent lithology and water saturation zone.

After the landslide, the adjacent slope was strongly deformed. Potential preferential flow paths were developed and formed by fractured zones, scarps, cracks, and fissures. Thus, much of the recharged groundwater infiltrates and accumulates under the deposited landslide mass. The slope toe is dry and does not drain water as runoff, causing gradual erosion and softening of the debris basement. An increase in the surface runoff and a decrease in the discharge capacity can act as reactivation factors. On the other hand, in the stable slope, the interflow of gravelly soil was not visible in the ERT-7 profile (Figure 11) due to the steepness of the slope. The surface runoff dominates the infiltration replenishment and erodes the gravel top layer in the exposed outcrops of the bedrock slopes. Thus, the absence of overlying colluvium and pebbly clay layers in the surrounding areas requires landslide susceptibility assessment. The sliding beds or moderately weathered bedrock layers on the slopes are greatly affected by the direct infiltration due to the absence of overlying colluvium layers. Hence, the landslide gradually deformed with time and finally collapsed. The trailing edge of the landslide was approximately 200 m wide and followed a steep slope with a depression angle of 25–35°. Currently, the saturated zone comprises a significant portion of the measured volume of the debris. Hence, oversaturation can cause the debris to move again.