Application of Resistivity and Seismic Refraction Tomography for Landslide Stability Assessment in Vallcebre, Spanish Pyrenees

Abstract

Geophysical surveys are a noninvasive reliable tool to improve geological models without requiring extensive in situ borehole campaigns. The usage of seismic refraction tomography (SRT), electrical resistivity tomography (ERT) and borehole data for calibrating is very appropriate to define landslide body geometries; however, it is still only used occasionally. We present here the case of a Spanish Pyrenees slow-moving landslide, where ERT, SRT and lithological log data were integrated to obtain a geological three-dimensional model. The high contrasts of P-wave velocity and electrical resistivity values of the upper materials (colluvial debris and clayey siltstone) provided accurate information on the geometry of the materials involved in the landslide body, as well as the sliding surface. Geophysical prospecting allowed us to identify the critical sliding surface over a large area and at a reduced cost and, therefore, gives the geophysical method an advantage over borehole data. The three-dimensional model was used to carry out stability analyses of a landslide in 2D and 3D, which, coherently with previous studies, reveal that the lower part is more unstable than the upper units.

1. Introduction

Landslides are typically complex bodies composed of different geological layers with contrasting and/or gradational physical properties. Landslides are considered the most catastrophic geological hazards, threatening and altering the socioeconomic circumstances of many world countries, and 66 million people are living in high-risk locations [1]. The complex nature of landslides requires a detailed investigation to define the geomorphological conditions of any particular landslide and gather data on stability conditions, soil and rocks properties, type of movement, velocity, etc. Several approaches have been developed to characterize such complex structures, describing landslide intensity and mitigating their effects.

Remote sensing techniques are widely utilized in landslide research, often in combination with laboratory testing. Slope stability modeling, hazard assessment, susceptibility mapping, risk analysis and other methods can also be used to assess landslide threats. [2]. The application of nonintrusive geophysical methods has increased exponentially for the subsurface characterization of landslides, especially to unveil the sliding surface geometry, to evaluate the emergence and growth of fractures, as well as to understand water dynamics and possible landslide reactivation by rainfall [3,4]. When water infiltrates a landslide, it alters the subsurface properties, causing destabilization and triggering slope failure when a critical threshold is exceeded. The changes in properties may be tracked using geophysical methods, and conclusions can be drawn about the destabilizing processes occurring in the subsurface of the landslide system by comparing data acquired over time. The results obtained from geophysical surveys provide cost-effective information on landslide-relevant layers, as well as the changes in water content. This implies that a stability analysis can be performed based on the geophysical investigation supported by geotechnical constraints such as soil unit friction angle and cohesion [5,6].

Geophysical techniques may provide continuous spatial information because the geophysical properties are controlled by porosity, stress state, soil type and degree of saturation [7]. The combination of seismic and electrical resistivity geophysical techniques has proved to be the most successful for studying the structure of landslides and understanding their internal mechanisms [8]. The methods have the ability to offer volumetric and spatial details of the subsurface [9,10,11,12,13]. The electrical resistivity tomography technique (ERT) is a low-cost and fast technique that can provide 2D electrical resistivity cross-sections and 3D resistivity models of the subsurface. Successful identification of subsurface lithology and water content fluctuation is possible thanks to the study and interpretation of their resistivity contrasts. In the last few decades, scholars have employed 2D ERT to unveil landslide geometry and internal structure and localize the sliding surfaces [9,14,15]. Seismic refraction tomography (SRT) is another practical nondestructive technique for determining the inner structure of sediments in most slope areas, particularly bedrock geometry. The first arrivals of seismic signals that travel through a layer with a higher seismic velocity are analyzed and are the base of SRT [5,14,16,17].

If the sliding forces are higher than the retaining forces, a slope might become unstable. As knowledge of geotechnical principles and landslide mechanisms has increased over the last few decades, slope stability analysis has substantially improved. The factor of safety (FS) has been used by many authors to describe landslide stability conditions and is basically the ratio between the resisting forces of a soil mass tending to slide down and its driving forces. There are numerous approaches for the factor of safety estimation, including limit equilibrium methods, which work on quantitatively identifying the slope stability based on momentum and forces analysis [18,19,20].

This paper’s major scope is to improve a three-dimensional geological model of a well-known landslide in the Spanish Pyrenees (Vallcebre landslide). In particular, the objectives are: (i) to apply geophysical methods (involving ERT and SRT) to define the inner features of the landslide, and (ii) to evaluate the landslide’s stability conditions using the new geological model and compare them to previous works.

2. Study Area

The studied site is placed in the Pyrenees range, 110 km west of the Mediterranean Sea. The Vallcebre landslide is located in the Llobregat river basin upper sector, on the occidental flank of the Serra Llacuna. The movement comprises an area of 0.8 km² (about 600 m wide and 1200 m long), with observable ground displacements and cracking. The landslides have been active for at least several centuries, and their exact age is uncertain [21].

The Vallcebre area has a sub-Mediterranean climate, with an average precipitation of 924 mm and 11.8 °C annual temperature. The analysis of the rainfall data has shown an important inter- and intra-annual variability and allows us to validate the seasonal character of the spatial variability of precipitation.

Geologically, the mobilized sediments are composed of a series of claystone, shale and gypsum layers sliding above a thick limestone substrate, which are all from the Upper Cretaceous to the Lower Paleocene period [21,22]. Inverse faults and several associated folds have an effect on the materials. Underneath the landslide, one of its faults has an estimated vertical jump of 10 m. The landslide affects three different slide units: The Lower Unit, the Intermediate Unit and the Upper Unit (Figure 1), and it could be classified as a translational landslide based on the categorization and geomorphological analysis by the authors of [23,24]. All units are defined by a moderate slope surface limited on the downhill edges by a few tens of meters high scarp, with an extension region developing in a graben geological form at the base of each scarp. This is interpreted as the Lower Units moving faster than the Upper Units, which has been verified by the implemented network of monitoring sensors. The landslide’s average slope is approximately 10° [21], and the underlying shearing surface depth reaches 42 m in the Middle Unit and 15 m in the Lower Unit, according to core drill holes (Figure 2), with most of the sliding surface being rather flat [22]. The Lower Unit is the most active section, with the Vallcebre water stream regularly eroding its toe, generating local rotational faults that reduce overall stability. For these reasons, this study will be focused on the Intermediate and the Lower Units (Figure 1).

Since the late 1980s, different monitoring approaches have been tested in the Vallcebre site. For example, surveying and photogrammetry have been used since the monitoring kick-off campaigns, and GPS surveillance started in the middle 1990s. The landslide was equipped with borehole wire extensometers, a network of inclinometers and piezometers over the next few years. Moreover, seven corner reflectors were installed for assessing the effectiveness of the difference interferometric synthetic aperture radar (DIn-SAR) technique in 2006. Ground-based SAR technology was also recently used to monitor the Vallcebre landslide [27].

Most of the first ruptures identified in the slopes near the Vallcebre landslide are the result of the evolution process of the slope itself. However, precipitation is often the most important feature in the reactivation of large landslides, like the one under study, and their acceleration. The Vallcebre landslide is very sensitive to groundwater fluctuation and rainfall [28]. The authors of [22] show that the rate of displacement and the water level are perfectly synchronized. This indicates that subsoil water has a significant impact on the balance of forces that control landslide dynamics. In our case, it is also important to highlight the large volume of water present in the form of waterspouts and streams all along the slope or in the form of water puddles in the counterslope area at the foot of the landslide. These occur around the year and are not limited to the rainy season. The abundant quantity of water—with the presence of gypsum—has generated important karstic dissolutions, giving rise to a significant system of conduits and porosity and causing numerous collapses in the site.

3. Methodology

3.1. Geophysical Methods

Considering that the landslide processes can be extremely complex, especially on slopes with very heterogeneous ground movements, and to better understand these processes, it is required to integrate different complementary approaches. Therefore, two complementary geophysical techniques were adopted to perform the characterization of lateral and vertical lithological changes of the landslide.

Due to the excellent quality of the data obtained and the capacity to generate continuous 2D and 3D subsurface images, the ERT method is one of the most widely used techniques for shallow subsurface surveys to solve environmental and engineering problems [29,30,31,32,33]. This technique allows detailed imaging of the distribution of subsurface electrical resistivity.

The ERT campaign was carried out along the Intermediate and the Lower Units. A Syscal Pro resistivity meter (IRIS Instruments, Orleans, France) with a Wenner–Schlumberger array was used to acquire geoelectrical data. This array has excellent signal strength and moderate properties, allowing it to identify both vertical and horizontal sharps [34,35]. The apparent resistivity measurements were acquired with 48 electrodes separated by 5 m. The electrodes and the ground contact resistance were less than 5 kΩ. The readings were repeated 3 to 6 times to obtain the smallest standard deviation values of the apparent resistivity. In these conditions, a total of 7 profiles were executed in our area of study (Figure 2).

RES2DINV software [36] was used to invert all of the apparent resistivity pseudosections. The authors of [37] describe the integration of the robust model constraint (L1) and smoothness-constrained (L2) optimization norms in two-dimensional resistivity inversion techniques. In practice: the L2 norm should be used if the subsurface resistivity changes are smooth and the L1 norm for sharp boundaries scenarios. In this paper, L1 was used for all data inversion because it was considered more robust with noisy values. Moreover, the L1 norm usually generates models with sharper boundaries, which was more consistent in our site, where large resistivity contrasts are expected between clayey siltstone and limestone units.

Moreover, the depth of investigation (DOI) index was employed in order to discriminate the electrical resistivity values related to geological structures from the ones related to the choice of the initial model. Introduced during the inversion procedure, the DOI index values will be close to 0 in areas that are highly dependent on data. The DOI values increase up to 1, where the gathered resistivity value for the model cell depends greatly on objective function parameters [38].

In this paper, we also provide the results of a four-profile STR survey (Figure 2). Seismic methods have an advantage over electrical methods as there is a direct relationship between porosity and acoustic velocities, though electrical methods are effective in such complicated geological systems for determining the position of the water table.

The seismic refraction method was carried out with a DAQ Link III 24-Channel seismograph (Seismic Source Co, Ponca City, OK, USA) with vertical geophones of 10 Hertz. A total of 24 geophones were installed along a 75-m profile line, with a distance between geophones of 3 m. For each line, nine hammer shots were employed. The 9 shots were distributed equidistantly along the profile, with the first shot at −3 m from the first geophone and the ninth shot at +3 m from the last geophone. To minimize the noise effects and increase the signal, five shots were stacked. The noise was monitored, and shots were taken only in low-noise conditions. A differential GRS1 GPS instrument (Topcon, Itabashi, Japan) was used to establish the geophone and shot coordinates.

To characterize the subsurface, SRT methods generally use a GRID of variable or fixed-sized cells. The travel times between shots and geophone and the ray paths were estimated using forward modeling approaches, such as the finite difference method. Cell P-wave velocities (Vp) were adjusted iteratively until the modeled and observed travel times misfit was minimal and within a tolerable range. The wavepath eikonal traveltime (WET) inversion approach (Morgenstern and Price, 1965) is adopted in this study.

The seismic sections were processed using the Rayfract^® software package (Intelligent Resources Inc. Software, Winnersh, UK). For each geophone and each shot, a manual picking of the first arrival time was performed. An example of P-wave arrivals was recorded, with first arrival times manually picked, and the wave path coverage for the SRT-2 profile is shown in Figure 3. The surface topography was included in the forward and the inversion models.

3.2. Landslide Stability Analysis

The authors of [39,40] started using three-dimensional limit equilibrium methods for slope stability assessment. The use of 3D stability analysis methods is becoming gradually more common due to the appearance of faster computers. The 3D limit equilibrium methods are based on statics, where the equilibrium state applies when the sum of forces applied and moments around a point equals zero. They basically work by dividing the sliding mass into columns with a square cross-section for three-dimensional evaluations. The forces and moment equilibrium equations are resolved for each column, considering the intercolumn forces. Then, the factor of safety is computed.

(1)$\mathrm{Factor}\mathrm{of}\mathrm{Safety}\left(\mathrm{FS}\right)=\frac{\mathrm{Resisting}\mathrm{forces}}{\mathrm{Driving}\mathrm{forces}}$

A factor of safety (FS) larger than 1 theoretically represents a stable slope. For landslides, values between 0.7 and 1.25 are considered a critical slope, and over 1.25 refers to a stable slope [18].

Slide3 and Slide2, modules of the Rocscience commercial software (Rocscience Inc., Toronto, ON, Canada) were utilized to obtain the factor of safety in two- and three-dimensional domains.

The three-dimensional geological model and the initial conditions (water levels) for the stability calculation were obtained from a combination of sources.

• Borehole data: Borehole data were inserted directly in the Slide3 software using the Borehole Manager tool, which allows inserting the coordinates, formations, and water elevation of each borehole;

• Digital elevation model (DEM) from airborne LIDAR data for adding surface topography to the model. DEM has a 2X2 m mesh with an absolute vertical accuracy of 0.15 m [41];

• Electric resistivity tomography profiles: Each material in these profiles has an electric resistivity range that represents it. Data were filtered and divided into four data sets, where each data set represents a material, then finally, the coordinates of the upper limit were saved as XYZ data GRID;

• Seismic refraction profiles: as each seismic refraction velocity range represents a material, velocity data were filtered and divided into four data sets. Each dataset is representing one of the subsurface formations. Finally, the coordinates of the upper limit of each layer were saved as XYZ data GRID.

All the XYZ data derived from the ERT and SRT profiles were integrated to have four data sets, each representing a geological formation of the studied area. The data sets, together with the fracture surface and water level, were inserted into Slide3 software, which has the ability to generate surfaces directly from XYZ data. The resultant surfaces were corrected manually to remove some identified defects before the analysis.

It should be noted that because water is considered a common landslide-triggering force, the highest water level record was used in the analysis, representing the worst-case scenario (lowest FS). Then, all these data, together with the borehole data, were interpreted to build the 3D model and carry out the analysis.

Mohr–Coulomb failure criteria were chosen to calculate the strength properties of the landslide formations, and the data were taken from authors of [21,42] after double-checking that the geometry was properly inserted and interpreted. The introduced materials were assigned to the model, using the section method. The section method allows the user to assign the materials to the model layers using a user-defined model section. The fissured shale layer was considered the weakest layer and the limestone as the bedrock.

According to the authors of [43], Slide3 software implements three-dimensional methods by extending their two-dimensional theories. Slide3 has four available three-dimensional methods: Bishop, Janbu, Spencer and Morgenstern and Price [39,44,45,46].

Bishop’s 3D method is a very accurate approach for circular slip surface stability analysis. The procedure assumes no shear forces between the columns and that the intercolumn forces act horizontally but it considers the moment equilibrium equations.

Janbu’s method adopts that shear forces are equal to zero and neglects the moment equilibrium. The latter calculates the factor of safety (FS) of shaped arbitrarily slip surfaces considering the vertical and horizontal force equilibrium equations of each column, then finally multiplies the calculated factor of safety (f₀) by a correction factor (a) to overcome the fact that the intercolumn shear forces were considered being zero.

(2)$\mathrm{FS}={\mathrm{a}\mathrm{f}}_0$

Spencer’s method considers the force and moment equilibrium equations for each column. The intercolumn forces are considered as well. This method is applicable on circular and arbitrary slip surfaces.

Morgenstern and Price’s method is appropriate for circular and noncircular slip surfaces; it considers force and moment equilibrium equations. The four methods were used to obtain four factors of safety maps which show the variation of FS along the landslide area. Then three sections were created to show the variation of the FS along with the sliding unit. The sections were exported from the Slide3 model Slide2 software for 2D analysis and FS estimation.

In order to perform the geotechnical evaluation and assess the factor of safety, we assigned three main parameters for each formation based on previous studies of the landslide area [21,42]. The parameters are unit weight, cohesion and angle of internal friction (Phi). Table 1 indicates the geotechnical parameters inserted in the Rocscience Slide3 software. Limestone has no cohesion nor angle of friction assigned because it has been considered an infinite-strength material. Limestone represents the bedrock, and the sliding is expected to occur on its upper surface, especially as it is overlaid by a weak bearing capacity thin layer of fissured shales.

The analysis parameters used in Slide3 are presented Table 2. The unit water weight was assumed to be 9.81 kN/m³. To approach convergence, several trials were performed using 50 columns in both X and Y directions, using 50 as a maximum iterations number. The tolerance between initial and final FS was set to be 0.001 and 0.005, and finally, 0.01 was set for the final results. The cuckoo search method [47] was used as a surface search method using a 30 iterations limit.

4. Results and Discussion

4.1. Geophysics

4.1.1. ERT

We selected the inverted cross-sections from the Wenner–Schlumberger (W-S) array for the interpretation of the results. The dipole–dipole (D-D) array inverted cross-sections showed a higher error rate (e.g., in ERT 2 RMS error of the W-S is 2.1; and the RMS error of the D-D is 4.0). A cause could be the considerable length of the profiles, which means the dipoles are separated by a large distance. Moreover, subsurface layers can be better distinguished in W-S cross-sections.

The seven inverted ERT cross-sections acquired in this study using the W-S array showed a root mean square error (RMS) of 1.1 to 2.5%—between measured and calculated apparent resistivity—after five iterations. As we did not have very large resistivity contrasts close to the ground surface, we used the standard model discretization [48]. The DOI index values calculated for all the profiles were less than 0.2. To compare the results of the ERT profiles, all the resistivity sections are displayed using the same color scale (Figure 4). The electric resistivity values range between 15 Ω·m and 750 Ω·m, and corresponding to their electrical resistivity values distribution, we can clearly distinguish two layers in all the cross-sections:

• An upper layer, generally with low resistivity values (<40 Ω·m) and with a thickness varying between 13 m (ERT 7) and more than 40 m (ERT 8).

• A lower layer, with high resistivity values (more than 300 Ω·m), was observed in all the cross-sections.

In some profiles (ERT 1, ERT 4, ERT 5 and ERT 8), the upper conductive layer is very heterogenous with interbedded materials of medium resistivity values (50 to 150 Ω·m) and elongated shapes. In the other profiles (ERT 2, ERT 3 and ERT 6), the upper layers are rather homogeneous.

The electric resistivity variations along the ERT profiles allow interpretation of the landslide main layers. The electric resistivity threshold values of each layer are known based on the correlation with borehole lithological data available in the area of study (Figure 2). Where the borehole log is missing, we extrapolated the stratigraphy boundaries over ERT results using the same threshold values. The upper low resistivity layer is correlated with the clayey siltstone registered in all the boreholes, whereas the deep-seated high resistivity layer is attributed to the limestone response. It is worth mentioning, firstly, that the layer of intermediate resistivity values registered in some profiles is associated with colluvial debris logged in boreholes S6–S7 and S10–S14. Secondly, the fissured shales crossed by several boreholes are not clearly identified in the ERT profiles. This is probably due to the electrical response of this level having values among the conductive shallow layer and the resistive deeper layer, so it is reflected as a transition area that cannot be detected with the ERT method.

In the ERT 8 cross-section—mainly carried out on the Intermediate Unit—it was possible to distinguish values between 100 and 150 Ω·m at the beginning of the cross-section, which corresponds to the colluvial debris sediments. The thickness of this layer is variable, with maximum values at the beginning of the profile and decreasing as we approach the cliff (almost at the end of the profile), where it disappears, transitioning to the clayey siltstones. In the deepest part of this profile, and unlike the rest of the profiles, the high resistivity values are only registered from the beginning to 100 m, which indicates that at least the contact zone with the limestone is deeper than 45 m in this part of the cross-section. From 100 m onwards, high values begin to be recorded, even at depths of less than 15 m in the escarpment zone (180 to 200 m from the beginning).

The ERT 5 profile also identifies, and more clearly, the scarp zone in the central part of the cross-section (120 m from the beginning). In this case, the dip in depth is very evident. The top of the limestone basement goes from less than 15 m deep before the escarpment to more than 30 m after the escarpment, thus indicating the boundary between the Intermediate and the Lower Units. It is also observed that the upper layer of the Intermediate Unit presents 100 to 150 Ω·m resistivity values corresponding to the colluvial debris, and its thickness decreases progressively as we move towards the escarpment.

ERT 2, ERT 3 and ERT 7 show similar resistivity distribution typical of the Lower Unit. In these cross-sections, we identify low resistivity values in the upper part and high resistivity values in the deeper zone, which correspond to clayey siltstone and limestone, respectively.

4.1.2. SRT

The datasets of the four P-wave SRT profiles were processed to characterize the lateral variation of the materials, estimate the thickness of layers—especially the clayey siltstone and the limestone contact—and evaluate the limestone’s top depth.

The final models were obtained after 20 WET iterations and a good fitting (RMS mean error < 5%) with the raw data was achieved (Figure 5). The ray coverage—which describes the ray paths model that passes below the surface—showed high values of ray density g in a large area of the profiles (see for instance Figure 3) and indicates good survey resolution [49]. The investigation depth for all the inverted datasets ranged from 17 m to 35 m using the ray coverage at the last iterations (Figure 3) and was hence suitable to characterize the clayey siltstone and the limestone layers contact, except SRT 3, where it was not possible to reach the depth of this contact.

The interpreted Vp models for cross-sections SRT 1 to SRT 4 are presented in Figure 5. A deep study of all cross/sections exposes that the P-wave tomography unveils the occurrence of a shallow medium velocity layer Vp (<2500 m/s), which punctually shows low Vp values (<1200 m/s) at a distance of 30, 59 and 68 m from the beginning in the case of the profile SRT 1 and a distance of 55 m from the beginning in the case of cross-section SRT 3. The thickness of this layer is largely variable, with a maximum depth for profile SRT 1 (about 23 m) and a minimum depth for profile SRT 4 (about 8 m). The identified geometry of the overburden—clay and colluvial debris—compares favorably with the 2500 m/s Vp contour depth, as we can observe in boreholes in the surrounding area of the cross-sections SRT 1 and SRT 2 (Figure 2). A second layer is defined by higher velocities (2500 m/s < Vp < 3000 m/s), and it is placed below this first layer. The thickness of this layer is also variable with an average of 6 m (SRT 1 and SRT 4) up to more than 10 m (SRT 2). We interpret this intermediate layer as the clayey siltstone according to the same borehole characteristics. It is important to highlight that this layer covers almost the entirety of the profile SRT 3.

At greater depth, P-wave velocity ranges from 3000 m/s to more than 4000 m/s in the seismic cross-sections SRT 1, SRT 2 and SRT 4. The isocontour of this high-velocity layer displays strong variability, with a depth lower than 20 m in the case of profile SRT 2 and SRT 4 and more than 25 m deep in the case of profile SRT 1. Considering the closest boreholes in the area, this layer was correlated with limestone recorded in the deepest part of these lithological logs. The fact that no high velocities (greater than 2500 m/s) are recorded in profile SRT 3 indicates that the limestone layer could not be reached. This agrees with the results of the nearest borehole (S7), where the limestone layer was recorded at deeper than 35 m and exceeds the maximum depth of SRT 3.

We also compare qualitatively the results obtained from refraction seismic and electrical resistivity campaigns, adding the closest lithological log data in ERT 2-SRT 2 cross-sections, and showing a good fit among them (Figure 6). The values of both seismic P-wave velocity and electrical resistivity that we assigned to limestones and siltstones are analogous to the values reported by the authors of [50,51] in similar settings.

4.2. Landslide Stability Analysis

Figure 7 shows the factor of safety map applying the four three-dimensional limit equilibrium methods available in Slide3 software. The map shows the variation of the FS along the studied area. The four methods showed quite similar results. The factor of safety ranges between 0.7 and 2, which clearly represents the instability of the area. The graben area (scarp) between the Intermediate Unit and Lower Unit shows the lowest values of FS (less than 1), while the Lower and Intermediate Unit showed higher FS values. There are some slight differences between the four methods (about 10%). Bishop, Morgenstern and Price (GLE), and Janbu showed another low FS anomaly in the upper limit of the Intermediate Unit, which is overlaid by the Upper Unit of the landslide area. Spencer showed another low FS anomaly on the right side of the model and a small anomaly in the Lower Unit. All of the methods showed a sliding trend towards the Northwest: 281.5°, 281.5°, 310.1° and 281.6°, for Bishop, Janbu, Spencer and GLE, respectively, which is the true sliding trend of the Vallcebre landslide.

The FS maps also showed fair agreement with [52], who studied the Lower Unit by using the infinite slope theory implemented in a geographic information system (GIS) and considered the landslide as infinite slices, assigning geotechnical parameters (unit weight, angle of friction and cohesion) for each as a rigid block, adding the fracture surface, water level variation, then finally determining the factor of safety to produce the FS map of the Lower Unit. The algorithm used has neglected the interslice forces, the 3D effect, and the fact that the landslide area consists of four formations with different geotechnical parameters. In the current research, the method used is quite different, as stated previously; this explains the differences between FS values in both articles. However, the values of FS were close in general as both methods showed convergence to unity, but there are differences in the location of low FS anomalies.

Figure 8 shows three geological cross-sections created with Slide3 software. The picture was edited to highlight the variation of FS along the sliding units (above the limestone formation). The sections show that the middle area has a lower factor of safety (fewer than 1), which represents the area between the Intermediate Unit and the Lower Unit. Section II shows the lowest values of FS. Afterward, these sections were exported to Slide2 software for two-dimensional analysis (Figure 9).

The four methods used for the two-dimensional analysis showed similar results, with a difference of about 5–10%. The reason relies on the fact that the two-dimensional model does not consider the variation of geological models (sections) along the slope area, while the three-dimensional analysis does [53]. In this research, the analysis performed using Slide2 (2D analysis) software showed about 10–15% differences from the analysis performed using Slide3 software (3D analysis).

Modeling exclusively using 2D may lead to either overestimation or underestimation of the factor of safety. The variation of FS along the sliding units in the 2D analysis showed slightly higher values (about 10%) than the 3D analysis for Section I and Section II, while it showed lower values than the 3D analysis for section III (about 15%). This slight difference between 2D and 3D analysis is similar to studies carried out on other sites, such as the authors of [54], who had a difference of about 15%. The authors of [55] studied five different cases to compare 2D and 3D analysis and showed an average difference of 26%, 0.5%, 12%, 2%, and 442.25%, respectively, between 2D and 3D analysis using the same commercial software Slide2 and Slide3.

The difference in 2D and 3D results is consistent with results obtained by other authors [53], as the outcomes of the 2D evaluation are always conservative, while the 3D analysis tends to enhance the safety factor [56].

Performing slope stability analysis using three-dimensional analysis would approximate the real condition better, as slope failures are three-dimensional in most cases [54,55]. However, three-dimensionality has its drawbacks. It is still a hot new topic that may require further improvements, and it needs a lot of work to build the geometry, especially when it comes to complex large-scale structures, as in our case. It also requires a substantial amount of computational power, especially for complex large-scale projects [55], which explains the several collapses of the computer used for analysis in this project.

Vallcebre landslide is of a transitional nature and continuous movement; thus, limit equilibrium analysis can give an insight into the nonbalanced forces magnitude and evaluate the spatial distribution of the instability. A full assessment of the Vallcebre landslide case would require dynamic models, including the time domain [21].

The values of FS shown in Figure 6 are quite similar to those obtained by the authors of [24], who developed a dynamic analysis approach highlighting the effect of mass variation on the stability of landslide and then applied it to the Vallcebre landslide using velocity, displacement and groundwater data obtained in borehole S2 for two years (21 November 1996–29 October 1998). The FS values from the four 3D models in this research show a value of 0.7–0.8 around the S2 borehole. The FS value calculated by the authors of [24] is about 0.67 for the highest groundwater level recorded on 1 January 1998. The slight difference in FS values is due to the different stability analysis methods used in each study.

5. Conclusions

The fundamental purpose of slope stability analyses is to determine the location and geometry of the principal slip surface. In this respect, geophysical techniques combining electrical resistivity and seismic refraction tomography proved to be highly successful, since it was feasible to distinguish the critical sliding surface over a large area of the study zone and at a reduced cost. This fact represents an advantage of the geophysical method over borehole data and complemented remote sensing analysis (photogrammetry and GPS surveillance techniques) used in the Vallcebre site since the early 1980s.

According to the results of the electrical resistivity tomography, we were able to differentiate between the upper layer of clayey siltstone, with minimal resistivity values, and the deeper layer of limestones, where the resistivity values are high. Between these two layers, we identified the transition zone with relatively lower resistivity values than the deep layer and which we related to the fissured shale layer. It is in this layer that the landslide zone was identified.

The seismic data indicate an increase in velocities with depth, as we move from the clayey siltstone at the surface to the limestone layer at depth. These velocity changes coincide with the interlayer boundaries that we recorded with the ERT.

The Depth of Investigation (DOI) parameters that we used to conduct the analyses can be seen in the Supplementary document.

The integration of all the data granted us to better identify the geometry of the landslide necessary for the computation of the factor of safety in the area. In this sense, the factor of safety was calculated using four different methods, but the results are quite similar. The analysis of the potential sliding surface showed an area with the highest risk of occurring a landslide in the escarpment zone between the Lower and the Intermediate Units. The results of the 2D modeling show a slight difference from the 3D modeling. Finally, we consider this case study as an interesting example of the combination of geophysical and geotechnical data to evaluate landslide stability.

Author Contributions

Conceptualization, M.H., A.S., M.E. and C.A.; methodology, M.H., L.R. and. G.P.D.; formal analysis, M.H., A.C. and A.S.; data curation, M.H., M.A., A.U. and R.L.; writing—original draft preparation, M.H., M.A. and A.C.; writing—review and editing, C.A., R.L., A.U., L.R., G.P.D., M.E., M.A., A.C., A.S. and M.H. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. (a) Location sketch and geological map of the studied zone. Modified from [25]; (b) Geomorphological sketch (modified from [21]) and digital elevation model of the Vallcebre landslide obtained using LIDAR data from Catalan Geological Survey [26].

Figure 2. Location of the geophysical surveys and simplified lithological logs of and the geological cross-section A’—A boreholes (Vallcebre area). Modified from [21].

Figure 3. (a) Refraction seismic shots recorded during the SRT-2 acquisition with the corresponding Picking analysis; (b) Wavepath coverage for the SRT-2.

Figure 4. ERT inverted model resistivity cross-sections and their interpretation based on the stratigraphic boundaries from borehole lithological logs and electrical resistivity threshold values.

Figure 5. Inverted seismic refraction tomography pseudo-sections with the corresponding interpretation.

Figure 6. Overlaying of SRT, ERT and lithological logs results.

Figure 7. Factor of safety variation map using different three-dimensional limit equilibrium methods: Bishop method, Morgenstern and Price method (GLE), Janbu method and Spencer method.

Figure 8. Factor of safety variation along the sliding materials of three different cross-sections.

Figure 9. The factor of safety values of Sections I (a), II (b) and III (c) using two-dimensional limit equilibrium methods (Slide2 software). The analysis was performed using Bishop, Janbu, Spencer and GLE; (d) location of Sections I, II and III.

Table 1

Geotechnical properties of materials used in Rocscience Slide3 software [16,32].

Unit Name	Unit Weight (kN/m3)	Cohesion (kPa)	Phi (°)
Limestone *1	20	- *1	- *1
Clayey Siltstone	18	0	14.7
Fissured Shales	19	0	11.8
Debris of Gravel	20	0	14.7

*¹ Limestone does not have cohesion and friction value because it has been considered an infinite-strength material [21,42].

Table 2

Analysis parameters used in the Rocsience Slide3 software.

Group	Parameter	Value
	Unit weight of water	9.81 kN/m3
Convergence	Number of columns in X or Y	50
	Maximum Iterations	50
	Tolerance	0.01
	Intercolumn force function	Half Sine
Surface options	Search method	Cuckoo Search
	Maximum Iterations	30
	Iteration tolerance	0.0001

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs14246333/s1.

References

1. Kjekstad, O.; Highland, L. Economic and Social Impacts of Landslides. Landslides—Disaster Risk Reduction; Sassa, K.; Canuti, P. Springer: Berlin/Heidelberg, Germany, 2009; pp. 573-587. ISBN 978-3-540-69970-5

2. Zhao, C.; Lu, Z. Remote Sensing of Landslides—A Review. Remote Sens.; 2018; 10, 279. [DOI: https://dx.doi.org/10.3390/rs10020279]

3. Ivanov, V.; Arosio, D.; Tresoldi, G.; Hojat, A.; Zanzi, L.; Papini, M.; Longoni, L. Investigation on the Role of Water for the Stability of Shallow Landslides—Insights from Experimental Tests. Water; 2020; 12, 1203. [DOI: https://dx.doi.org/10.3390/w12041203]

4. Pazzi, V.; Morelli, S.; Fanti, R. A Review of the Advantages and Limitations of Geophysical Investigations in Landslide Studies. Int. J. Geophys.; 2019; 2019, 2983087. [DOI: https://dx.doi.org/10.1155/2019/2983087]

5. Caris, J.P.T.; Van Asch, T.W.J. Geophysical, geotechnical and hydrological investigations of a small landslide in the French Alps. Eng. Geol.; 1991; 31, pp. 249-276. [DOI: https://dx.doi.org/10.1016/0013-7952(1)90011-9]

6. Whiteley, J.S.; Chambers, J.E.; Uhlemann, S.; Wilkinson, P.B.; Kendall, J.M. Geophysical Monitoring of Moisture-Induced Landslides: A Review. Rev. Geophys.; 2019; 57, pp. 106-145. [DOI: https://dx.doi.org/10.1029/2018RG000603]

7. Di Maio, R.; De Paola, C.; Forte, G.; Piegari, E.; Pirone, M.; Santo, A.; Urciuoli, G. An integrated geological, geotechnical and geophysical approach to identify predisposing factors for flowslide occurrence. Eng. Geol.; 2020; 267, 105473. [DOI: https://dx.doi.org/10.1016/j.enggeo.2019.105473]

8. Havenith, H.-B.; Torgoev, I.; Ischuk, A. Integrated Geophysical-Geological 3D Model of the Right-Bank Slope Downstream from the Rogun Dam Construction Site, Tajikistan. Int. J. Geophys.; 2018; 2018, 1641789. [DOI: https://dx.doi.org/10.1155/2018/1641789]

9. Chambers, J.E.; Wilkinson, P.B.; Kuras, O.; Ford, J.R.; Gunn, D.A.; Meldrum, P.I.; Pennington, C.V.L.; Weller, A.L.; Hobbs, P.R.N.; Ogilvy, R.D. Three-dimensional geophysical anatomy of an active landslide in Lias Group mudrocks, Cleveland Basin, UK. Geomorphology; 2011; 125, pp. 472-484. [DOI: https://dx.doi.org/10.1016/j.geomorph.2010.09.017]

10. Di Maio, R.; Piegari, E. Water storage mapping of pyroclastic covers through electrical resistivity measurements. J. Appl. Geophys.; 2011; 75, pp. 196-202. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2011.07.009]

11. Grandjean, G.; Gourry, J.C.; Sanchez, O.; Bitri, A.; Garambois, S. Structural study of the Ballandaz landslide (French Alps) using geophysical imagery. J. Appl. Geophys.; 2011; 75, pp. 531-542. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2011.07.008]

12. Himi, M.; Casado, I.; Sendrós, A.; Lovera, R.; Rivero, L.; Casas, A. Assessing preferential seepage and monitoring mortar injection through an earthen dam settled over a gypsiferous substrate using combined geophysical methods. Eng. Geol.; 2018; 246, pp. 212-221. [DOI: https://dx.doi.org/10.1016/j.enggeo.2018.10.002]

13. Jongmans, D.; Garambois, S. Geophysical investigation of landslides: A review. Bull. Société Géologique Fr.; 2007; 178, pp. 101-112. [DOI: https://dx.doi.org/10.2113/gssgfbull.178.2.101]

14. Hack, R. Geophysics for slope stability. Surv. Geophys.; 2000; 21, pp. 423-448. [DOI: https://dx.doi.org/10.1023/A:1006797126800]

15. Perrone, A.; Lapenna, V.; Piscitelli, S. Electrical resistivity tomography technique for landslide investigation: A review. Earth-Sci. Rev.; 2014; 135, pp. 65-82. [DOI: https://dx.doi.org/10.1016/j.earscirev.2014.04.002]

16. Hibert, C.; Grandjean, G.; Bitri, A.; Travelletti, J.; Malet, J.-P. Characterizing landslides through geophysical data fusion: Example of the La Valette landslide (France). Eng. Geol.; 2012; 128, pp. 23-29. [DOI: https://dx.doi.org/10.1016/j.enggeo.2011.05.001]

17. Jongmans, D.; Bièvre, G.; Renalier, F.; Schwartz, S.; Beaurez, N.; Orengo, Y. Geophysical investigation of a large landslide in glaciolacustrine clays in the Trièves area (French Alps). Eng. Geol.; 2009; 109, pp. 45-56. [DOI: https://dx.doi.org/10.1016/j.enggeo.2008.10.005]

18. Cahyaningsih, C.; Riau, U.I.; Asteriani, F.; Crensonni, P.F.; Choanji, T. Safety factor characterization of landslide in Riau-West of Sumatra highway. Int. J. Geomate; 2019; 17, pp. 323-330. [DOI: https://dx.doi.org/10.21660/2019.63.69945]

19. Chowdhury, R.; Flentje, P. Geotechnical analysis of slopes and landslides: Achievements and challenges Geologically Active. 11th IAEG Congress of the International Association of Engineering Geology and the Environment; Taylor and Francis: Oxford, UK, 2010; Volume 31.

20. Johansson, J. Impact of Water-Level Variations on Slope Stability. Licentiate Thesis; Luleå University of Technology: Luleå, Sweden, 2014.

21. Corominas, J.; Moya, J.; Ledesma, A.; Lloret, A.; Gili, J.A. Prediction of ground displacements and velocities from groundwater level changes at the Vallcebre landslide (Eastern Pyrenees, Spain). Landslides; 2005; 2, pp. 83-96. [DOI: https://dx.doi.org/10.1007/s10346-005-0049-1]

22. Gili, J.A.; Moya, J.; Corominas, J.; Crosetto, M.; Monserrat, O. Past, Present and Future Monitoring at the Vallcebre Landslide (Eastern Pyrenees, Spain). Appl. Sci.; 2021; 11, 571. [DOI: https://dx.doi.org/10.3390/app11020571]

23. Cruden, D.M.; Varnes, D.J. Landslides-Investigation and Mitigation, Special Report No. 247; National Academy Press: Washington, WA, USA, 1996; ISBN 9780309061513

24. Ferrari, A.; Ledesma, A.; González, D.A.; Corominas, J. Effects of the foot evolution on the behaviour of slow-moving landslides. Eng. Geol.; 2011; 117, pp. 217-228. [DOI: https://dx.doi.org/10.1016/j.enggeo.2010.11.001]

25. ICGC. Mapa Geològic de Catalunya 1:50,000. Available online: www.icgc.cat (accessed on 3 October 2021).

26. ICGC. Topography Map 1:5000 and DEM 2X2. Institut Cartogràfic i Geològic de Catalunya. Available online: www.icgc.cat (accessed on 2 December 2020).

27. Crosetto, M.; Gili, J.A.; Monserrat, O.; Cuevas-González, M.; Corominas, J.; Serral, D. Interferometric SAR monitoring of the Vallcebre landslide (Spain) using corner reflectors. Nat. Hazards Earth Syst. Sci.; 2013; 13, pp. 923-933. [DOI: https://dx.doi.org/10.5194/nhess-13-923-2013]

28. Moya, J.; Corominas, J.; Gili, J.A.; Ledesma, A.; Lloret, A. Twenty years of continuous monitoring of the Vallcebre landslide (Eastern Pyrenees): Lessons learned. Proceedings of the JTC1 Workshop on Advances in Landslide Understanding; Barcelona, Spain, 24–26 May 2017; pp. 101-104.

29. Cho, I.-K.; Yeom, J.-Y. Crossline resistivity tomography for the delineation of anomalous seepage pathways in an embankment dam. Geophysics; 2007; 72, pp. G31-G38. [DOI: https://dx.doi.org/10.1190/1.2435200]

30. Loke, M.H.; Dahlin, T. A comparison of the Gauss–Newton and quasi-Newton methods in resistivity imaging inversion. J. Appl. Geophys.; 2002; 49, pp. 149-162. [DOI: https://dx.doi.org/10.1016/S0926-9851(01)00106-9]

31. Griffiths, D.H.; Barker, R.D. Two-dimensional resistivity imaging and modelling in areas of complex geology. J. Appl. Geophys.; 1993; 29, pp. 211-226. [DOI: https://dx.doi.org/10.1016/0926-9851(93)90005-J]

32. Loke, M.H.; Barker, R.D. Practical techniques for 3D resistivity surveys and data inversion. Geophys. Prospect.; 1996; 44, pp. 499-523. [DOI: https://dx.doi.org/10.1111/j.1365-2478.1996.tb00162.x]

33. Sjödahl, P.; Dahlin, T.; Johansson, S.; Loke, M.H. Resistivity monitoring for leakage and internal erosion detection at Hällby embankment dam. J. Appl. Geophys.; 2008; 65, pp. 155-164. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2008.07.003]

34. Bedrosian, P.A.; Burton, B.L.; Powers, M.H.; Minsley, B.J.; Phillips, J.D.; Hunter, L.E. Geophysical investigations of geology and structure at the Martis Creek Dam, Truckee, California. J. Appl. Geophys.; 2012; 77, pp. 7-20. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2011.11.002]

35. Martorana, R.; Fiandaca, G.; Casas Ponsati, A.; Cosentino, P.L. Comparative tests on different multi-electrode arrays using models in near-surface geophysics. J. Geophys. Eng.; 2009; 6, pp. 1-20. [DOI: https://dx.doi.org/10.1088/1742-2132/6/1/001]

36. Loke, M.H.; Barker, R.D. Rapid least-squares inversion of apparent resistivity pseudosections by a quasi-Newton method. Geophys. Prospect.; 1996; 44, pp. 131-152. [DOI: https://dx.doi.org/10.1111/j.1365-2478.1996.tb00142.x]

37. Loke, M.H.; Acworth, I.; Dahlin, T. A comparison of smooth and blocky inversion methods in 2D electrical imaging surveys. Explor. Geophys.; 2003; 34, pp. 182-187. [DOI: https://dx.doi.org/10.1071/EG03182]

38. Oldenburg, D.W.; Li, Y. Estimating depth of investigation in DC resistivity and IP surveys. Geophysics; 1999; 64, pp. 403-416. [DOI: https://dx.doi.org/10.1190/1.1444545]

39. Morgenstern, N.R.; Price, V.E. The Analysis of the Stability of General Slip Surfaces. Géotechnique; 1965; 15, pp. 79-93. [DOI: https://dx.doi.org/10.1680/geot.1965.15.1.79]

40. Anagnosti, P. Three dimensional stability of fill dams. Proceedings of the 7th International Conference on Soil Mechanics and Foundation Engineering; Mexico City, Mexico, 4–8 May 1969; International Society for Soil Mechanics and Geotechnical Engineering: London, UK, 1969; pp. 275-280.

41. ICGC. Model d’Elevacions del Terreny de Catalunya 2x2metres (MET-2); version 1.0 Institut Cartogràfic i Geològic de Catalunya: Barcelona, Spain, 2016.

42. Corominas, J.; Moya, J.; Lloret, A.; Gili, J.A.; Angeli, M.G.; Pasuto, A.; Silvano, S. Measurement of landslide displacements using a wire extensometer. Eng. Geol.; 2000; 55, pp. 149-166. [DOI: https://dx.doi.org/10.1016/S0013-7952(99)00086-1]

43. Cheng, Y.M.; Yip, C.J. Three-Dimensional Asymmetrical Slope Stability Analysis Extension of Bishop’s, Janbu’s, and Morgenstern–Price’s Techniques. J. Geotech. Geoenvironmental Eng.; 2007; 133, pp. 1544-1555. [DOI: https://dx.doi.org/10.1061/(ASCE)1090-0241(2007)133:12(1544)]

44. Bishop, A.W. The use of the Slip Circle in the Stability Analysis of Slopes. Géotechnique; 1955; 5, pp. 7-17. [DOI: https://dx.doi.org/10.1680/geot.1955.5.1.7]

45. Janbu, N.; Bjerrum, L.; Kjaernsli, B. Soil Mechanics Applied to some Engineering Problems; Norwegian Geotechnical Institute: Oslo, Norway, 1956.

46. Spencer, E. Thrust line criterion in embankment stability analysis. Géotechnique; 1973; 23, pp. 85-100. [DOI: https://dx.doi.org/10.1680/geot.1973.23.1.85]

47. Wu, A. Locating General Failure Surfaces in Slope Analysis via Cuckoo Search. Available online: https://static.rocscience.cloud/assets/verification-and-theory/Slide2/Cuckoo-Search.pdf (accessed on 3 November 2022).

48. Loke, M.H. RES2DInv; Rapid 2D Resistivity and IP Inversion Using the Least Squares Method Geotomo Software: Tanjung Penaga, Malaysia, 2004.

49. Nurul Nida, H.; Asta Yogantara, M.; Otolomo, M.A.P.; Sanny, T.A. Application of Refraction Seismic Tomography Method Using Wavepath Eikonal Traveltime Inversion for Modelling the Subsurface. Proceedings of the IOP Conference Series: Earth and Environmental Science; Bali, Indonesia, 7–10 August 2018; Sule, R. Institute of Physics Publishing (IOP): Bali, Indonesia, 2019; Volume 318, 012016.

50. Caylak, C.; Kocaslan, A.; Gorgulu, K.; Buyuksarac, A.; Arpaz, E. Importance of ground properties in the relationship of ground vibration–structural hazard and land application. J. Appl. Geophys.; 2014; 104, pp. 6-16. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2014.02.003]

51. Yilmaz, S.; Okyar, M.; Gurbuz, M. Electrical resistivity imaging for investigation of seepage paths in the Yukari Gökdere Dam, Isparta, Turkey. Bull. Eng. Geol. Environ.; 2021; 80, pp. 115-125. [DOI: https://dx.doi.org/10.1007/s10064-020-01952-2]

52. Santacana, N. Análisis de la Suceptibilidad del Terreno a la Formación de Deslizamientos Superficiales y Grandes Deslizamientos Mediante el uso de Sistemas de Información Geográfica. Aplicación a la Cuenca alta del Río Llobregat. Ph.D. Thesis; Universitat Politècnica de Catalunya: Barcelona, Spain, 2001.

53. Arief, M. Studi pengaruh geometri Lereng Pada analisis Kemantapan Lereng 2D Dan 3D Dengan Metode Kesetimbangan Batas. Indones. Min. Prof. J.; 2020; 2, pp. 51-56. [DOI: https://dx.doi.org/10.36986/impj.v2i1.27]

54. Bar, N.; McQuillan, A. 3D Limit Equilibrium Slope Stability Analysis for Anisotropic and Faulted Rock Masses in Australian Coal and Iron Ore Mines. Proceedings of the ISRM International Symposium-10th Asian Rock Mechanics Symposium; Singapore, 29 October–3 November 2018.

55. Firincioglu, B.S.; Ercanoglu, M. Insights and perspectives into the limit equilibrium method from 2D and 3D analyses. Eng. Geol.; 2021; 281, 105968. [DOI: https://dx.doi.org/10.1016/j.enggeo.2020.105968]

56. Xie, M.; Esaki, T.; Zhou, G.; Mitani, Y. Geographic Information Systems-Based Three-Dimensional Critical Slope Stability Analysis and Landslide Hazard Assessment. J. Geotech. Geoenvironmental Eng.; 2003; 129, pp. 1109-1118. [DOI: https://dx.doi.org/10.1061/(ASCE)1090-0241(2003)129:12(1109)]