SUMMARY

Origins of ancient rockslides in seismic regions can be controversial and must not necessarily be seismic. Certain slope morphologies hint at a possible coseismic development, though further analyses are required to better comprehend their failure history, such as modelling the slope in its pre-failure state and failure development in static and dynamic conditions. To this effect, a geophysical characterization of the landslide body is crucial to estimate the possible failure history of the slope. The Balta rockslide analysed in this paper is located in the seismic region of Vrancea-Buzau, Romanian Carpathian Mountains and presents a deep detachment scarp as well as a massive body of landslide deposits. We applied several geophysical techniques on the landslide body, as well as on the mountain crest above the detachment scarp, in order to characterize the fractured rock material as well as the dimension of failure. Electrical resistivity measurements revealed a possible trend of increasing fragmentation of rockslide material towards the valley bottom, accompanied by increasing soil moisture. Several seismic refraction surveys were performed on the deposits and analysed in form of P-wave refraction tomographies as well as surface waves, allowing to quantify elastic parameters of rock. In addition, a seismic array was installed close to the detachment scarp to analyse the surface wave dispersion properties from seismic ambient noise; the latter was analysed together with a colocated active surface wave analysis survey. Single-station ambient noise measurements completed all over the slope and deposits were used to further reveal impedance contrasts of the fragmented material over in situ rock, representing an important parameter to estimate the depth of the shearing horizon at several locations of the study area. The combined methods allowed the detection of a profound contrast of 70–90 m, supposedly associated with the maximum landslide material thickness. The entirety of geophysical results was used as basis to build up a geomodel of the rockslide, allowing to estimate the geometry and volume of the failed mass, that is, approximately 28.5–33.5 million m3.

1 INTRODUCTION

Slopes in hilly or mountainous regions can present complex morphologies due to erosional, gravitational and seismo-tectonic processes. Gravitational slope failures can appear as large mass movements, such as land- and rockslides or rock falls, generating deep scarps, hummocky surfaces, depletion and accumulation zones. Especially in mountain ranges of seismic regions, numerous deep-seated and voluminous slope collapses can be found (cf. Crosta et al. 2013): for example, Saymareh (Iran; Shoaei 2014), Sary-Chelek, Kara Suu and Beshkiol (Kyrgyzstan; Strom 2010; Havenith et al. 2015), Iskanderkul (Tajikistan; Strom 2010) and Ancient Diexi (China; Fan et al. 2019). Yet, their origin is often not clearly established and must not necessarily be seismic. As a result of weathering and ageing, older collapses can present ‘softened’ morphologies, making it more difficult to understand circumstances of slope development and to trace limits of disturbed rock and bedrock masses. Furthermore, assigning the origin of ancient slope failures to purely climatic or seismic factors is extremely challenging, as already highlighted by Strom & Abdrakhmatov (2018), and requires thorough investigations. The main aspects to be considered are regional precipitation quantities and distributions, local hydrogeological and hydrological regimes as well as the seismic activity, spatial landslide distribution and seismo-tectonic markers of the studied area in relation to ancient and active faults (Crozier et al. 1995; Crozier 1999; Strom & Abdrakhmatov 2018). Also, the shape and location of detachment scarps can be indicative; while shallow mass movements located near the bottom of a slope are more often identified as non-seismic, deep-seated failures close to mountain ridges are more likely to have a seismic origin (Meunier et al. 2008; Strom & Abdrakhmatov 2018).

Geophysical prospection can be used to characterize the internal structure and several material properties of slope collapses. In general, determining the geometry and magnitude of a failed slope (expressed in terms of mass or volume of the release area) allows to estimate initial pre-failure slope conditions as well as to understand possible trigger events through back-analysis. Reviews and case studies of geophysical surveys on landslide bodies all show the advantage of using multiple techniques combined (see Pazzi et al. 2019). Most popular in this field is the use of electrical resistivity tomographies (ERT), spontaneous potential (SP), active and passive seismic surveys, as well as electromagnetic methods such as ground penetrating radar (GPR) and time-domain electromagnetics (TDEM; Godio & Bottino 2001; Meric et al. 2005; Jongmans & Garambois 2007; Jongmans et al. 2009; Thiery et al. 2017; Cody et al.  2020 ; Pazzi et al. 2019).

Electrical resistivity measurements are commonly used to investigate hydrogeological conditions of a slope and possibly to define internal shearing horizons of a landslide body (see, e.g. Lapenna et al. 2003; Danneels et al. 2008; Panek et al. 2008; Perrone et al. 2014; Gance et al. 2016). Commonly combined with electrical tests in landslide reconnaissances are active seismic refraction surveys, such as the seismic refraction tomography (SRT). This technique, described in, for example, Heincke et al. (2006), Grandjean et al. (2007), Göktürkler et al. (2008), Hibert et al. (2012), Uhlemann et al. (2016) and Whiteley et al. (2020), can be used to produce P-wave tomographies or analyse surface wave properties. One method to study the latter is the MASW technique (multichannel-analysis of surface waves), commonly applied to estimate 1-D shear wave velocity profiles (Vs) of the soil. Another method to quantify the shear wave velocity of the subsurface is the analysis of ambient vibration array measurements with multiple seismic stations (Wathelet et al. 2008; Renalier et al. 2010; Hussain et al. 2019; Harba et al. 2019). While the depth and resolution of the aforementioned techniques are constrained by the applicable length of survey lines and the configuration of receivers, single station passive seismic techniques are unbound from such limiting geometries and constitute a powerful tool to retrieve local subsurface properties at great depth. The H/V spectral ratio (HVSR) measurements are increasingly used in the context of slope failure characterization (e.g. Méric et al. 2007; Pazzi et al. 2017; Imposa et al. 2017; Rezaei et al. 2020). The technique facilitates the detection of deep surfaces, especially where two materials are marked by a strong impedance contrast—such as between landslide deposits and underlying bedrock.

In this study, several geophysical prospection methods were applied on an ancient landslide body; the prospective combination with numerical calculations aims at a better comprehension of the slope’s failure history in a possibly seismic context.

2 STUDY AREA

The focus here lies on the bend of the Carpathian Mountains belt in Romania (also called ‘Curvature Carpathians’; see inset figure of Fig. 1), and more precisely, on the seismic region Vrancea-Buzau, where a large number of mass movements are recognized in the relatively smooth relief of the Flysch Carpathians. The majority of the mass movements shown in Fig. 2 consist of rather elongated and most likely climatically induced landslides in form of debris and mud flows. Due to a rather elevated annual precipitation rate of 600–800 mm and thawing periods in late winter and early spring, the Curvature Carpathians show a good correlation between climatic peak events and landslide occurrences (Dragotă 2006; Micu et al. 2013; Ministry of Environment 2018). Only some mass movements in the area can be characterized as deep-seated rockslides. For the latter, we assume that previously existing critical slope conditions in combination with powerful trigger events, such as seismic vibrations, possibly led to repeated or first time failure. In general, the region of Vrancea-Buzau is marked by a high seismicity with intermediate-depth earthquakes; the strongest and best documented Vrancea earthquakes occurred in the years 1802, 1838, 1940 and 1977 with estimated Mw 7.4–7.9 and intensities of I=IX-X MSK (Georgescu 2002; Rogozea et al. 2014; Vacareanu & Ionescu 2016). In Fig. 1, a seismicity map shows the events between the years 1900–2014 during which four events of Mw> 7 occurred in the Vrancea-Buzau region, that is the earthquakes of 1908, 1940, 1977 and 1986 with focal depths of 125, 150, 94 and 131 km, respectively.

Seismicity map of the Vrancea-Buzau region: earthquake events marked in terms of moment magnitude Mw (provided by the Romanian Institutului National de Fizica Pământului) and representation of the focal mechanisms of the 1977 and 1986 events (USGS 2020). The red star marks the location of the Balta rockslide; coordinate system: UTM 35N (in m).
Figure 1.

Seismicity map of the Vrancea-Buzau region: earthquake events marked in terms of moment magnitude Mw (provided by the Romanian Institutului National de Fizica Pământului) and representation of the focal mechanisms of the 1977 and 1986 events (USGS 2020). The red star marks the location of the Balta rockslide; coordinate system: UTM 35N (in m).

Geological map of the study area (modified after IGR 1968) with the outlines of the Balta rockslide and its neighbouring landslides (mapped by Damen et al. 2014); coordinate system: UTM 35N (in m).
Figure 2.

Geological map of the study area (modified after IGR 1968) with the outlines of the Balta rockslide and its neighbouring landslides (mapped by Damen et al. 2014); coordinate system: UTM 35N (in m).

Historic reports of these events concentrate mostly on damages and casualties in densely populated areas. Nonetheless, especially for the 1940 event, geohazards such as groundwater fluctuations, floods, liquefaction, as well as wide and open cracks on slopes, and landslides have been documented (Radulescu 1941; Radu & Spânoche 1977; Georgescu & Pomonis 2012; Mândrescu 1981). As these historic records only cover the period of the last 200 yr in isolated areas, investigations of markers resulting from possible past seismic activity, such as deep-seated slope failures, would contribute to the reconstruction of an extended regional seismic history, and thus support seismic hazard assessment in the Vrancea-Buzau region for longer terms.

To better understand deep-seated slope failures of the Vrancea-Buzau region related to dynamic processes, we present the following case study focused on the Balta rockslide (located on the map in Fig. 1, views in Fig. 3), that developed in a Palaeogene schistose sandstone flysch with thick sandstone banks. The south-facing slide is located on a reverse fault and counts numerous, mostly shallower, neighbouring landslides (Fig. 2).

The Balta rockslide shown on Google Earth® images (the direction of failure is marked by the white arrows): (a) outline of the failure with marked subareas; (b) view along the river valley towards the east with the rockslide front spreading towards the north into the valley; (c) photo showing the scarp from the valley (taken in June 2017) and (d) UAV image of the scarp (taken in October 2018) and indication of strata dip.
Figure 3.

The Balta rockslide shown on Google Earth® images (the direction of failure is marked by the white arrows): (a) outline of the failure with marked subareas; (b) view along the river valley towards the east with the rockslide front spreading towards the north into the valley; (c) photo showing the scarp from the valley (taken in June 2017) and (d) UAV image of the scarp (taken in October 2018) and indication of strata dip.

The slope failure can be classified as rock rotational slide (or rock slump, according to Hungr et al. 2014) and presents morphologies that hint at a possible coseismic failure. One of these morphologies is the detachment scarp visible from the valley bottom (see Fig. 3c); remarkable are its geomorphological characteristics, that is in terms of steepness (mean slope angle of 40°, in the upper parts it even exceeds 50°), height (approximately 300 m) and amphitheatre shape, suggesting a powerful (possibly seismic) trigger event. The head scarp furthermore shows a 20–50° inclined antidip slope bedding (see Fig. 3d) and is situated close to the mountain ridge, relatively far from the valley bottom. Note, a more detailed analysis of the structural aspects of the bedding and joint planes of the Balta rockslide, together with a 2-D numerical analysis, is presented in the parallel work of Lemaire et al. (2020). The scarp features, together with the large body of debris that is still in place (in form of crushed landslide material as well as rock boulders), mark the contrast to the more shallow landslides and flows of the region.

The accumulation of the landslide deposits seems to concentrate in the middle part in form of a plateau (see Figs 3a and c), a rather flat area, while the rockslide toe spreads into the valley and supposedly deviated the river Bâsca within (see Fig. 3b; Mreyen et al. 2017). As the opposite slope of Balta is also subject to former landsliding, a plausible historic scenario of simultaneous failure could have implied the formation of an important dam and a temporary barrier lake in the valley. Thanks to the good accessibility and the rather large amount of debris estimated to be on top of in situ rock in the middle part of the slope (plateau), the study site offers good conditions for geophysical explorations.

3 METHODS

For this work, the focus lies on the use of active and passive seismic investigation methods, which are compared with profiles of electrical resistivity. The presented data were collected during three geophysical campaigns (2016–2018) on the Balta rockslide. The analysed geological, geophysical, and seismological outputs were integrated in a 3-D geomodel that was used as tool to interpret the rockslide in terms of morphology and structure, possible sliding horizons (basal and lateral) as well as elastic material properties. A similar approach is presented in Havenith et al. (2018), where the integration of multiple geophysical and geological data in a 3-D geomodel was used to assess the slope stability of an ancient mass movement downstream of a dam construction site.

3.1 Electrical resistivity tomography

The geophysical technique of electrical resistivity tomography (ERT) allows to map the electrical resistivity per unit cell in the subsurface along a profile. Measured changes in electrical resistivities of a subsoil can mark changes in its physical properties of lithological origin. In this context, weathered or fractured rock (often linked to elevated moisture content) will show lower resistivities than intact and non-fractured rock (Jongmans & Garambois 2007; Perrone et al. 2014).

In our case, we used the ERT technique along five sections distributed over the rockslide body in order to detect pronounced conductivity contrasts that could indicate changes in water content and material properties of the debris material. The locations of the ERT surveys are shown in Fig. 4(a). An extended profile (E1, 590 m length) covers the landslide foot to the central part, with a second smaller profile (E2, 395 m length) crossing it approximately perpendicularly. The surveys E3 and E4 are located at the rockslide plateau, while E5 is situated at the outer western part (each with a length of 395 m). The data acquisition was carried out with a Geotom MK8E 1000 high-resolution multi-electrode resistivity system from GEOLOG together with four cables linking 4-m-spaced electrodes (maximum profile length of 395 m); the ERT survey was performed in June 2016. The electrode positions were measured with a differential GPS system (Topcon HiPer V). We applied the Wenner electrode configuration, as it is sensible to vertical resistivity changes (Dahlin & Zhou 2004). The data inversion was performed with the 2-D inversion algorithm of Loke & Barker (1996) implemented in the RES2DINV software. To estimate the true resistivity distribution of the subsurface, a nonlinear smoothness-constrained least-squares inversion was applied to the apparent resistivity data (deGroot Hedlin & Constable 1990). As the 2-D inversion of the latter depends strongly from the applied damping factor, different initial values were tested, that is, 0.05–0.15, before reaching an optimal damping factor of 0.1 by applying the L-curve method (cf. Farquharson & Oldenburg 2004); the latter was increased by a factor of 1.2 with each underlying layer according to Loke (1998).

Location of geophysical surveys (rockslide outlined by hatched line): (a) electrical survey lines E1–E5; (b) seismic survey lines S1–S8 and single station HVSR and (c) onfiguration of the seismic array. Topographical data: TanDEM-X (10.5-m-resolution; German Aerospace Center, DLR); coordinate system: UTM 35N (in m).
Figure 4.

Location of geophysical surveys (rockslide outlined by hatched line): (a) electrical survey lines E1–E5; (b) seismic survey lines S1–S8 and single station HVSR and (c) onfiguration of the seismic array. Topographical data: TanDEM-X (10.5-m-resolution; German Aerospace Center, DLR); coordinate system: UTM 35N (in m).

3.2 Seismic refraction survey

Actively induced seismic energy can be used to determine the elastic properties of the ground by mapping seismic wave propagation in terms of traveltimes.

For the P-wave tomography, first arrival times of the body waves are picked to produce traveltime curves that are used to create an initial velocity model of the survey line followed by a tomographic inversion. Traveltimes are calculated iteratively by ray tracing and are compared with the observed traveltimes until a minimized RMS error is reached. For the data analysis, we used the software package of Seisimager 2D™ developed by Geometrics (more details about the applied data treatment and inversion is given in the supplementary material).

By applying the Multichannel analysis of surface waves (MASW), the Rayleigh-wave energy of the seismic survey is interpreted separately (as surface waves are controlled by dispersion, a property that is not attributed to body waves; Park et al. 1999). The method allows quantifying the 1-D shear wave velocity Vs of the soil (considering surface wave velocity ∼0.9 Vs; Lay & Wallace 1995) without being biased by its water content–contrary to the methods of electrical resistivity and P-wave tomography. The dispersion curve is obtained through a frequency–wavenumber analysis (fk analysis; Rost & Thomas 2002); the inversion of the dispersion curve is performed with the dinver module of Geopsy (see Wathelet 2008). Thousands of velocity models are created by varying the number of layers, their thickness as well as corresponding body wave velocities. Associated theoretical dispersion curves are computed according to Dunkin (1965) and compared with the observed dispersion curve by calculating a misfit value. The generation of models follows a neighbourhood algorithm and the inversion proceeds until the misfit converges to a minimum value (more details concerning the inversion process and corresponding parametrization are included in the supplementary material).

For this study, six seismic profiles were completed on the failed body and two above the rockslide scarp at the hill crest (Fig. 4b). The data acquisition was carried out during two time periods: July 2017 (surveys S1–S5) and October 2018 (surveys S6–S8). The general survey configuration consisted of 24 receivers (4.5 Hz; 5 m spacing, i.e. 155 m profile length for profiles S1–S5) connected to a 16-bit 24 channel seismograph; except for the extended profile S6 of 355 m length, using 48+24 receiver with a 24-bit 48 channel seismograph. Seismic energy was triggered with sledgehammer shots.

3.3 Ambient noise

The analysis of seismic ambient noise is recognized to retrieve the subsurface properties such as shear waves velocities and impedance contrasts (e.g. Wathelet 2008; Del Gaudio et al. 2013; Pazzi et al. 2017). For this purpose, we installed a seismic array of 7 velocimeters GURALP CMG-6TD with lower cut-off frequency of 0.033 Hz on the plateau (in October 2018) and we performed single station measurements in 53 points (during three field campaigns in the years 2016–2018) over the rockslide by using a velocimeter Lennartz 3-D with lower cut-off frequency of 0.2 Hz.

3.3.1 Surface wave dispersion

An array of 7 seismic stations was installed on the plateau of the slide. We used a configuration made of 6 seismic stations on the vertices of two triangles of 60 and 90 m-long medians, respectively, around a central station (see Fig. 4c).

To ensure optimal azimuthal coverage, the medians of the two triangles are rotated relatively to the other by an angle of 60°. A simultaneous record duration of at least 75 min was ensured for the 7 seismic stations. The signals are analysed by using a frequency-wavenumber technique (e.g. Rost & Thomas 2002) applied on windows of 50 s with an overlap of 50 per cent. The signal windows are first selected through a STA/LTA approach (Earle & Shearer 1994) in order to eliminate contributions due to transient events (STA = 1 s; LTA = 30 s; threshold = 3.5). We derive the kinematic properties of waves crossing the array from the inter-station phase differences derived from the cross-spectral matrix. By exploiting a high-resolution technique (Capon 1969), the slowness vector is estimated at each frequency step from the cross-spectral matrix components by maximizing a beam function. The dispersion curve is computed by converting the slowness vector into apparent phase velocities. We determine the domain of validity of the dispersion curve based on the array geometry by computing the array response (Rost & Thomas 2002) and seismic noise energy content. The observation of the dispersion curve provides direct information on the resolution and maximum investigation depth of the velocity profile; empirical laws suggest that these latter correspond to approximately one third of the minimum and maximum measured wavelength, respectively (Lay & Wallace 1995; Park et al. 1999).

As for the dispersion curves obtained from MASW, the inversion of the Rayleigh dispersion properties derived from array measurements is performed by using the dinver module of Geopsy (see Wathelet 2008).

3.3.2 HVSR

The horizontal-to-vertical-spectral-ratio method (HVSR or H/V; Nakamura 1989) permits to identify impedance contrasts at depth by the estimation of the fundamental frequency of resonance of the investigated site. The technique exploits the three-component recordings at each measurement point by computing the spectral ratios between the horizontal and the vertical components. The fundamental resonance frequency is associated with the H/V peak measured at low frequency.

The signal is firstly pre-processed with a band-pass filter between 0.2 and 20 Hz. The H/V curves are computed on 50-s-long windows sliding the entire recording time (averaged recording time is about 20 min). Similarly to the surface wave dispersion analysis, we previously applied an STA/LTA filter in order to avoid the contribution of transient signals in the analysis. We applied this method to the seismic measurements recorded at seismic stations reported in Fig. 4(b) including the seismic stations installed as an array as shown in Fig. 4(c). In the simple geological configuration where an impedance contrast is encountered at depth, like a loose surface layer overlying a bedrock, the resonance frequency fr is related to the thickness of the surface layer thickness h through the following expression (Ibs-von Seht & Wohlenberg 1999):
(1)
where Vs is the S-wave velocity in the surface layer.

The results of our geophysical exploration, with a special focus on HVSR, will be integrated and analysed with the help of the 3-D geomodelling software Leapfrog Geo by Seequent©. The software allows compiling surface as well as subsurface data, for example geophysical tomographies, geological sections or drillings, and to interpret these in terms of lithological units as well as geological and tectonic structures in a 3-D space. As such, 1-D and 2-D subsurface data can be linked using the in the software integrated radial basis function algorithm (RBF; cf. Powell 1992; Horowitz et al. 1996). Furthermore, the geomodel allows for the creation of subsurface units - here, used to estimate the volume of rockslide mass.

4 RESULTS

4.1 ERT

The five electrical surveys carried out on the slide deposits show resistivity values between 10 and 350 Ωm. A maximum depth of penetration of 60–70 m could be reached in their central part. The surveys E3 and E4 (see Fig. 5a) situated on the plateau, that is, close to the detachment scarp, indicate resistivity values up to 350 Ωm, while the relatively higher values were detected near the surface. For the surveys E5, E1 and E2 (see Figs 5b and c, respectively) on the lower part of deposits, that is, at the rockslide foot, values generally do not exceed 150–200 Ωm. Compared to the upper part of the rockslide body, a clear trend towards lower resistivity values can be noted at the foot. The intersecting surveys E3–E4 and E1–E2 show a good coherence in the measured electrical resistivities at their respective crossing points.

Electrical resistivity measurements on the rockslide body: (a) intersecting surveys E3 and E4 at the plateau, view from SE; (b) survey E5 at the eastern part of rockslide mass and (c) intersecting surveys E1 and E2 at the landslide foot, view from NW.
Figure 5.

Electrical resistivity measurements on the rockslide body: (a) intersecting surveys E3 and E4 at the plateau, view from SE; (b) survey E5 at the eastern part of rockslide mass and (c) intersecting surveys E1 and E2 at the landslide foot, view from NW.

4.2 Active and passive seismics

In the following, results of the active (P-wave refraction) and passive (HVSR and array) seismic surveys are presented.

4.2.1 P-wave velocity Vp

In Fig. 6, we show the results of the active seismic refraction survey as P-wave (Vp) tomographies created with the tool package Plotrefa of Seisimager 2D™, Geometrics. The colour scale was fixed in the range of Vp < 500 to 3200 m s–1as a result of the minimum and maximum Vp values detected during ray tracing of the tomographic inversion (with exception of the survey S6 where higher Vp values were detected close to the scarp, see explications below).

Seismic refraction surveys showing profiles of P-wave velocity [m s–1]: (a) S1 at the central foot; (b) S8 at the hill crest; (c) S2 at the western foot; (d) S3 at the toe; (e) S6 on the plateau; (f) and (g) S4 and S5 at the eastern limit of deposits, respectively. The dashed lines indicate the lower limit of ray paths outlined during tomographic inversion of the respective survey. Note, survey S7 is not represented due to insufficient precision of DGPS collected elevation data (profile located in a dense forest area).
Figure 6.

Seismic refraction surveys showing profiles of P-wave velocity [m s–1]: (a) S1 at the central foot; (b) S8 at the hill crest; (c) S2 at the western foot; (d) S3 at the toe; (e) S6 on the plateau; (f) and (g) S4 and S5 at the eastern limit of deposits, respectively. The dashed lines indicate the lower limit of ray paths outlined during tomographic inversion of the respective survey. Note, survey S7 is not represented due to insufficient precision of DGPS collected elevation data (profile located in a dense forest area).

In general, at the surface of the rockslide body, very low velocities <1000 m s–1 were detected in the upper 20–30 m. The survey S1 performed near the central part of the landslide foot (Fig. 6a) shows a relatively low velocity gradient, reaching Vp of 2000–2500 m s–1 in a depth of 50 m, while the survey S3 at the rockslide toe (Fig. 6d) demonstrates the same velocity range in shallower depths of 25–30 m. Similar to the latter, survey S4 (Fig. 6f) shows very low velocities at the surface and a relatively constant velocity gradient with depth. The surveys S2 and S5 (Figs 6c and g, respectively) show high velocities starting from 45 m depth up to Vp > 3000 m s–1 close to the designated western and eastern limits of the main rockslide body, respectively.

At the plateau, the extended survey S6 (Fig. 6e) demonstrates relatively low velocities in the northern half of the profile, that is, the downhill part of the plateau area (Vp < 1500 m s–1 over the first 25 m below surface). However, in the southern half of the profile, that is, the flat area that is close to the main scarp, velocities increase strongly, attaining Vp > 3000 m s–1 in 35–40 m depth (compared to 1800–2300 m s–1 at the same depth in the central and northern part of the survey line) and a maximum Vp of 4195 m s–1 in 50 m depth. The maximum depth of the S6 profile is reached in the central part with 71.4 m below surface (at 230 m profile length; corresponding Vp = 2814 m s–1).

The surveys at the hill crest, south of the rockslide detachment scarp, demonstrate very low Vp at the direct surface, gradually increasing with depth, for example, S8 (Fig. 6b) presents low Vp < 1100 m s–1 in the upper 10 m below surface and a strong velocity increase over 15 m to finally attain Vp > 2300 m s–1 in 25–30 m depth.

4.2.2 S-wave velocity Vs

The Vs profiles obtained from MASW are shown in Fig. 7 and are defined approximately up to 50 m below the surface. On the lower part of the slide (S2, S3; Fig. 7a) low velocities of approximately 180 m s–1 are measured at the surface. These velocities follow a common increase up to the depth of 25 m where S3 is subject to a high velocity contrast (Vs = 1300 m s–1) while S2 reaches 650 m s–1 at 36 m. The S1 profile situated on the foot presents higher velocities of about 320 m s–1 at the surface increasing to 710 m s–1 at a depth of 40 m. The S6 survey at the plateau presents velocities increasing from 410 to 510 m s–1 in its upper 40 m. On the eastern profiles S4 and S5, we observe an increase in velocity from 300 and 180 m s–1 at the surface respectively to a velocity of 660 m s–1 at a depth of 50 m (Fig. 7b). At the crest, higher velocities are measured at depth for the profiles S7 and S8: velocities starting from 250 to 300 m s–1 at the surface reaching approximately 840 m s–1 at a depth of 35 m as shown in Fig. 7(c).

Best velocity models resulting from the inversion of the MASW dispersion curves: (a) surveys at the toe (S3), foot (S1 and S2) and plateau (S6); (b) surveys on the edges (S4 and S5) and (c) surveys at the hill crest (S7 and S8).
Figure 7.

Best velocity models resulting from the inversion of the MASW dispersion curves: (a) surveys at the toe (S3), foot (S1 and S2) and plateau (S6); (b) surveys on the edges (S4 and S5) and (c) surveys at the hill crest (S7 and S8).

In Table 1, we list Poisson’s ratio ν for two respective depths per seismic survey. Note, the P-wave velocity values Vp were thereby extracted from the 2-D refraction profiles centre, while S-wave velocity values Vs originate from the 1-D MASW analysis. The Poisson’s ratio is defined by the following relation:
(2)
The Poisson’s ratio constitutes an important parameter in elastic rock material analysis and is often used in landslide monitoring, slope stability assessment and modelling (e.g. Uhlemann et al. 2016). It allows us to classify rock types and granular soils, as well as to estimate the water content of natural ground, with ν ranging from 0 to 0.5 and ν = 0.5 referring to saturated cohesive soils (Gercek 2007). Here, values of ν range from 0.29 to 0.48, that is, rather elevated values referring to an important water content in a soft material, as we encounter in the fractured landslide material of Balta.
Table 1.

Poisson’s ratio ν at given depths for the active seismic surveys S1–S8 (* indicates values originating from 2-D Vp profiles with strong lateral variations in P-wave velocities and thus resulting in uncertain ν estimations).

SurveyDepth (m)Vp (m s–1)Vs (m s–1)ν (–)
S11012554850.41
4019967090.43
S2*108502780.44
3719536580.44
S31117743450.48
33249313590.29
S498773190.42
4022896690.45
S5*912503240.46
3022974980.47
S6*911914100.43
3024935100.48
S7912824650.42
3822338770.41
S8911824010.43
3225248920.42
SurveyDepth (m)Vp (m s–1)Vs (m s–1)ν (–)
S11012554850.41
4019967090.43
S2*108502780.44
3719536580.44
S31117743450.48
33249313590.29
S498773190.42
4022896690.45
S5*912503240.46
3022974980.47
S6*911914100.43
3024935100.48
S7912824650.42
3822338770.41
S8911824010.43
3225248920.42
Table 1.

Poisson’s ratio ν at given depths for the active seismic surveys S1–S8 (* indicates values originating from 2-D Vp profiles with strong lateral variations in P-wave velocities and thus resulting in uncertain ν estimations).

SurveyDepth (m)Vp (m s–1)Vs (m s–1)ν (–)
S11012554850.41
4019967090.43
S2*108502780.44
3719536580.44
S31117743450.48
33249313590.29
S498773190.42
4022896690.45
S5*912503240.46
3022974980.47
S6*911914100.43
3024935100.48
S7912824650.42
3822338770.41
S8911824010.43
3225248920.42
SurveyDepth (m)Vp (m s–1)Vs (m s–1)ν (–)
S11012554850.41
4019967090.43
S2*108502780.44
3719536580.44
S31117743450.48
33249313590.29
S498773190.42
4022896690.45
S5*912503240.46
3022974980.47
S6*911914100.43
3024935100.48
S7912824650.42
3822338770.41
S8911824010.43
3225248920.42

We remind that the MASW Vs surveys, in contrast to Vp, represent an average of the shear-wave velocities over the entire profile length as a function of depth, and thus 2-D effects might lead to inaccurate estimates of the Poisson’s ratios. An example for the latter are the surveys S5 and S6, presenting strong lateral Vp variations and the resulting elevated ν of 0.47. For the surveys S1, S3, S4 and S8, however, velocities increase homogeneously with depth along the profiles and estimations of Poisson’s ratios are more reliable.

The Rayleigh dispersion curve obtained through the analysis of ambient noise vertical signals recorded at the seismic array on the plateau is shown in the supplementary material, Fig. S1. The array response limits the validity of the curve to the frequency range of 2.7–6 Hz. As the survey S6 and the seismic array are colocated (see Fig. 4b), we gathered both dispersion results. The dispersion curves obtained from the MASW of survey S6 are defined over a 4.5–12.7 Hz frequency band for the fundamental mode of propagation, and a 7–14 Hz frequency band for the first higher mode of propagation (see Fig. S2). The fundamental modes of the dispersion curves found by array analysis and MASW are in good agreement as the phase velocities are consistent at 4.5 Hz (see Fig. 8); we thus decided to join the MASW survey S6 (4.5–12.7 Hz and 7–14 Hz) and the array results (2.7–4.5 Hz) in order to define a combined fundamental mode and a first higher mode that can be used as target in the inversion process (see Fig. 8).

Dispersion curves from seismic ambient noise analysis (from 2.7 to 4.5 Hz) and MASW (from 4.5 to 14 Hz), fundamental mode and 1st higher mode of propagation are picked as target for the inversion process.
Figure 8.

Dispersion curves from seismic ambient noise analysis (from 2.7 to 4.5 Hz) and MASW (from 4.5 to 14 Hz), fundamental mode and 1st higher mode of propagation are picked as target for the inversion process.

The minimum and maximum wavelengths measured on the resulting dispersion curve are λmin = 28 m and λmax = 310 m, respectively. We therefore estimate the resolution of the velocity profile at about 9.3 m and its maximum investigation depth at 103 m.

At the array location, the HVSR measurements indicate an averaged resonance frequency of 1.35 Hz. We integrated this information with the combined dispersion curve and applied a joint inversion process. The resulting velocity profile is shown in Fig. 9(b), while Fig. 7(a) only precises the upper 50 m of the best model obtained. Over the first ten meters below surface, a Vs of 410 m s–1 was measured; the velocities increase gradually to 510 m s–1 up to a depth of about 62 m. Below this depth, we observe a velocity contrast: in the half-space, the shear wave velocity reaches approximately 850 m s–1. Due to decreasing resolution with depth, we encounter larger uncertainties associated with the depth of the interface: models with misfit lower than 0.05 present interfaces with depths spanning the range between 62 and 92 m (see Fig. 9b). A sensitivity analysis was applied following Luo et al.(2011) and is reported in Fig. S3. This latter highlights that the array measurements principally control the depths below 60 m, while the MASW measurements define results from surface to 60 m depth.

Results of the inversion of Rayleigh wave dispersion curves: (a) dispersion curve target of the inversion shown by black dots, while the coloured curves are the dispersion curves associated with the generated models for the inversion. (b) Generated velocity models associated with the dispersion curves in (a). For both, (a) and (b), the colour scale represents the misfit.
Figure 9.

Results of the inversion of Rayleigh wave dispersion curves: (a) dispersion curve target of the inversion shown by black dots, while the coloured curves are the dispersion curves associated with the generated models for the inversion. (b) Generated velocity models associated with the dispersion curves in (a). For both, (a) and (b), the colour scale represents the misfit.

4.2.3 HVSR

The HVSR measurements are distributed over different sections on the rockslide: (i) the failed body, that is, on the plateau, representing the main body, the foot and the toe; (ii) in the proximity of the failed body and on the hill crest, that is, above the main scarp. Four examples of H/V curves taken from the toe, the foot, the plateau and the crest are shown in Fig. 10.

Examples of H/V curves and corresponding spectral ratios measured at four different locations of the landslide: the toe, the foot, the plateau and the crest.
Figure 10.

Examples of H/V curves and corresponding spectral ratios measured at four different locations of the landslide: the toe, the foot, the plateau and the crest.

An overview of the recorded resonance frequencies f0 over the area is given in Fig. 11. On the plateau (i; approximately 610 m a.s.l.), the resonance frequency shows a peak at 1.3 Hz. On the lower part of the rockslide body (approximately 450 m a.s.l.), in the centre of the foot, we measured f0 values ranging from 1.2 to 2 Hz. Higher values are observed on the side parts of the failed body mass with 2–5 Hz resonance frequencies; values >2.5 Hz were generally recorded at the borders and front of the rockfall deposits. These relatively elevated frequencies could either indicate shallower depths of detected impedance contrasts between two materials or could be due to 2-D–3-D effects (e.g. close to the edges of the main body). On the hill crest (ii; approximately 950 m a.s.l.), measurements show resonance frequencies between 1 and 1.3 Hz.

Resonance frequency f0 (Hz) results of the single station H/V survey on the Balta rockslide. The colours of the measurement points designate the range of measured resonance frequency, while the shapes refer to the surface wave velocity Vs used per area (averaged over the whole thickness of the near-surface layer—which is increasing from the toe of the rockslide to the plateau area).
Figure 11.

Resonance frequency f0 (Hz) results of the single station H/V survey on the Balta rockslide. The colours of the measurement points designate the range of measured resonance frequency, while the shapes refer to the surface wave velocity Vs used per area (averaged over the whole thickness of the near-surface layer—which is increasing from the toe of the rockslide to the plateau area).

5 INTEGRATION AND DISCUSSION OF RESULTS

Though the origin and the age of the Balta rockslide are unknown, its morphological and structural characteristics indicate a significant (probably unique) trigger event. From a geomorphic-geological point of view, the failure presents a deep-seated shape (cf. Lemaire et al. 2020) together with a rather large body of landslide deposits. The main slope exposed to the north could indicate a partial climatic origin of the landslide (as it is supposedly affected by a higher moisture content). Considering the relatively low altitude of the study site, that is, 600–900 m a.s.l., and the consistent vegetation cover of the area, we primarily consider the slope as initially stable.

5.1 Characterization of the rockslide body based on geophysical results

The detected range of resistivity values from 15 to 350 Ωm indicates a subsurface with relatively low resistivity contrasts. However, we could identify a general trend showing higher resistivity values at the plateau (surveys E3 and E4) than at the lower part of the slide deposits (surveys E1, E2 and E5). This could indicate an increase in the moisture content towards the valley bottom and, especially, close to the river at the rockslide front (toe). The SE end of E5 survey, even though located at the border of the plateau area, does not present elevated resistivity values as E3 and E4, probably also due to elevated soil moisture (as a small rivulet source was observed in that area during our field campaigns—possibly as natural drainage of the slope). The high conductivity of the soil due to water-filled pores could explain the measured values (Loke 2004, suggests a resistivity range of approximately 10–100 Ωm for fresh groundwater). In relation to that, we assume a variation in the degree of fragmentation over the landslide body; the material near the toe of the rockslide seems to be more crushed than the debris forming the plateau; though also singular very large sandstone blocks (>7 m in diameter) could be found in the lower part of the deposits.

The water content of the slope supposedly also affects the results of the seismic refraction survey, with water saturated layers showing higher P-wave velocities than the vadose zone (due to the general relation of soil moisture and Vp; Gregory 1976). To verify the impact of moisture content, the Vp results should be related to the MASW results from the same data analysed in terms of Vs. The Poisson’s ratios (Table 1) indicate that the material is highly saturated, as it is also shown by the ERT surveys at the rockslide foot (e.g. ν of 0.48 at 11 m depth of the S3 profile). However, these values should be relativized in terms of pure rock properties due to the assumed high water content of the slope. As given by the ERT results, we expect the foot part of the debris body to be even more affected by an elevated moisture content than the upper part close to the detachment scarp.

In general, we notice very low P-wave velocities at the surface of each seismic refraction survey, on the rockslide body as well as south of the detachment scarp on the hill crest. These values supposedly indicate a general weathering of the shallow soil and rock material, at depths of 5–15 m. For surveys showing relatively low velocities, we assume the presence of a larger amount of fractured and even crushed rocks at depth, than for surveys where larger values of Vp and Vs could be measured. The surveys S1 and S2 at the foot of the slide demonstrate a lower increase of P-wave velocity with depth than the survey S3 at the toe. This supposedly indicates a rather large accumulated debris volume in the central part that decreases towards the toe of the slope. As already shown by the Vp results, the corresponding MASW Vs curve clearly indicates a strong velocity increase at medium depth (20–30 m) for the S3, in contrast to the S1 and S2 surveys where no such contrast could be found (cf. Fig. 7a). This agreement between elevated Vp and Vs, together with the elevated H/V frequencies of 2.5– 7.5 Hz in its vicinity, supposedly confirms that only a shallow layer of soft material is in place at the S3 survey location, and consolidated in situ rock was detected at shallow depth. The latter can be supported by the decrease of Poisson’s ratio from 0.48 to 0.29, 11–33 m below surface, respectively, possibly referring to a change in soil composition and related to that, soil moisture. For the lateral limits of the rockslide, the western section of survey S2 shows P-wave velocities increasing to higher values over shallower depth compared to its eastern section. The latter possibly indicates the thickness of debris flattening towards the western end of the profile. The same can be observed for the survey S5 in the eastern part of the survey with supposed layers of fractured material thinning towards the eastern limit of the rockslide body.

At the plateau, the survey S6 reveals the detection of relatively elevated Vp values close to the scarp. These high velocities (Vp > 4000 m s–1 in 50 m depth), reached in the southern part of the survey, suggest a rather consolidated subsoil. However, the combined (seismic array and MASW) Vs study at this location highlight a constant increase of elastic properties with depth, until a marked velocity contrast is reached at a depth of about 70–75 m depth that can be associated with an interface (cf. Figs 7a and 9b). The northern part of the survey on the other hand is comparable to the survey S1 in the middle part of the slope with similar Vp and Vs detected in similar depths (cf. Figs 6a, e and 7a). The measured values of the hill crest survey lines S7 and S8 demonstrate for both, P and S wave, constantly increasing velocity gradients with depth, reaching higher velocities in depth than for the surveys located on the rockslide debris (cf. Figs 6b and 7c). For this site, we suppose that the in situ rock is largely weakened at the surface due to climatic effects, possibly enhanced by strong amplification effects near the crest.

Note, for the surveys S2, S4 and S5, the dispersion curves (see Fig. S3) are identified down to a frequency of 3.5 Hz, that is, below the natural geophone frequency of 4.5 Hz. This implies that between 3.5 and 4.5 Hz, the signals can be affected by amplitude damping and signal distortion due to the instrumental response, leading possibly to erroneous estimations of the surface waves phase velocities. However, the fk method applied on the data is based on the computation of phase delays between pairs of signals, that are relative measurements. The possible phase delay introduced by the instrument is therefore eliminated by calculating the differences between the signals under analysis, as the 48 geophones are identical and share the same distortion effects. For our purpose, we considered the distortion on the signals below the natural frequency of instruments as negligible and extended the Vs analysis down to 3.5 Hz for these three surveys.

On the plateau, we combined the results obtained from active seismic and passive seismic measurements. Havenith et al. (2007) previously showed that the combined dispersion curve helps to take advantage of both techniques. By using the combined frequency band, we exploit the high resolution given by the phase velocities at high frequencies (from MASW) and the larger investigation depth given by the phase velocities at low frequencies (from array). At 4.5 Hz, a good compliance in phase velocities is reached. Between 4.5 and 6 Hz, the phase velocities present higher uncertainties at the array; these latter might be due to the contributions of body waves present in the seismic noise, inferring higher phase velocity measurements. Also, as noted on the velocity–frequency diagram in Fig. S1, aliasing effects affect the array measurements above 5 Hz. For the given observed uncertainties between 4.5 and 6 Hz at the array, we therefore adopted the MASW measurements on the 4.5–12.7 Hz frequency band for the fundamental mode and limited the array results between 2.7 and 4.5 Hz.

Parametrization in the inversion process of surface waves dispersion properties is an important aspect to take into account, being a non-linear problem and subject to non-uniqueness solutions. The number of model parameters should be kept as low as possible to avoid ill-posed problems on one side, and sufficient enough to represent the complexity of the measured dispersion curves (Wathelet et al. 2004). The sensitivity analysis applied on the model at the Balta plateau shows that S-wave velocity is a dominant inversion parameter in respect to the P-wave velocity (see Fig. S4); this result is in agreement with Xia et al. (1999) and Foti et al. (2018) for near-surface studies. In this paper, we adopted a simple parametrization by considering an increasing number of horizontal layers varying in thickness and S-wave velocities, and kept the minimum number of layers with misfit of the best solutions converging to a minimum value.

Note, the presented geophysical data was collected during three field campaigns in the years 2016–2018. As the ERT data were collected during one single summer campaign, the moisture content of the soil is supposed to be steady throughout all surveys. For the SRT data, the surveys S1–S5 were performed in summer, whereas the surveys S6–S8 were performed in early autumn in the consecutive year. It should be kept in mind that the moisture content is thus not identical for both survey parts; however, the Poisson’s ratios (Table 1) do not indicate any marked differences.

5.2 Depth of the main shearing horizon

Based on the passive seismic tests combined with the MASW from the active seismic surveys, seismic impedance contrasts were derived for each HVSR measurement point. In the area of the main body of failure, such a contrast can be used as indicator for the vertical limit of debris material over in situ rock. To estimate the thickness of debris per measurement point, the main body was subdivided into three zones of averaged S-wave velocity as, depending on the location on the slope and the associated subsurface material, the detected Vs values vary (cf. Fig 7). The HVSR measurement point locations were linked to these representative S-wave velocities of the area (Fig. 11).

For the rockslide toe and foot, we defined a Vs of 330 m s–1 (averaged from the MASW of S2 and S3), while for the middle slope, that is, upper part of the foot area, a higher Vs of 430 m s–1 was averaged (from the MASW of S1, S4 and S5). Here, impedance contrasts are observed in the depth range of approximately 40–70 m by using f0 values in the range of 1.2–2 Hz. Close to the plateau area, depth ranges of material contrasts detected with the HVSR method can go down to 90 m. The mean Vs value designated to the plateau area is 510 m s–1 (averaged from the MASW of S6 and the seismic array). Combining the latter with the measured central resonance frequency value of 1.3 Hz, impedance contrasts at a maximum depth of 98 m are estimated. Considering both, the results of the dispersion curves analysis and the recorded HVSR resonance frequencies, we assume the maximal depth of the shearing horizon in approximately 70–90 m below the surface.

5.3 Volume estimation of the failure

For the further analyses, we used a total of 58 HVSR measurements on the main body that present clear resonance frequency peaks. Together with their associated S-wave velocity, that is, 330, 430 or 510 m s–1, the depth to the detected impedance contrast was calculated (see eq. 1). Here, this depth is inferred from the thickness value that was calculated for the fractured debris material, subtracted from the surface elevation value at the measurement point. In Fig. 12, we show an exemplary cross section of the study area with the modelled 1-D HVSR logs in its vicinity (located at maximum 30 m lateral distance from the profile line). As described before, the distance from the surface down to the detected impedance contrast, here expressed as logs, diminishes from the scarp to the toe. By assembling a large number of 1-D measurement points, these thicknesses can be spatially linked and interpolated to a subsurface layer.

Cross section A–B (location marked by red line in the view on the left). Estimated thicknesses of failed mass shown as logs from the HVSR measurements H11, H17, H26, H09, H20 and H30 along the section; the green line marks the interpolated basal shearing plane. The double arrows mark the defined Vs zones, respectively (cf. Fig. 11).
Figure 12.

Cross section A–B (location marked by red line in the view on the left). Estimated thicknesses of failed mass shown as logs from the HVSR measurements H11, H17, H26, H09, H20 and H30 along the section; the green line marks the interpolated basal shearing plane. The double arrows mark the defined Vs zones, respectively (cf. Fig. 11).

The latter was done within a 3-D geomodelling software (Leapfrog Geo, Seequent) where the 1-D H/V basal impedance contrasts were interpolated to an isosurface, that is, the estimated basal shearing plane of the rockslide. While the Balta failure is estimated to cover a superficial area of 1.47 × 106 m2, a volume of 33.5 × 106 m3 of debris mass in place could be calculated by using the created basal shearing plane together with the surface elevation model of the actual topography. In Fig. 13, we show a general view and two cross sections of the created geomodel; the displayed HVSR logs define the estimated depth and shape of the rockslide plane, as well as the calculated volume of debris (represented in green). For the seismic refraction surveys at the landslide foot and toe, plotted within the geomodel (see Fig. 13), it can be observed that the estimated basis of debris is located at depths where P-wave velocities of 2000–2500 m s–1 are reached. On the plateau on the other hand, the elevated Vp measured at the southernmost part of S6, near the source area, probably indicate that the sliding surface was reached (located at a depth of less than 50 m). In addition, at this location, the survey possibly indicates that part of the in situ rocks were less fragmented, also due to the minor sliding distances, and thus allowed the build-up of the flattened plateau area. However, towards the north (beyond the centre of the S6 line), the sliding surface is located beyond the SRT penetration depth (>60 m).

3-D geomodel of the Balta rockslide with volume of rockslide debris (coloured in green), integrated SRT (coloured sections), and HVSR (black logs): (a) full view; (b) E–W cross section and (c) N–S cross section.
Figure 13.

3-D geomodel of the Balta rockslide with volume of rockslide debris (coloured in green), integrated SRT (coloured sections), and HVSR (black logs): (a) full view; (b) E–W cross section and (c) N–S cross section.

Considering the phenomenon of bulking, that is, the expansion that occurs during slope collapse, the here estimated volume of 33.5 × 106 m3 is most probably overestimated (e.g. with a bulking factor of 15 per cent, the originally detached rock volume would equal to approximately 28.5 × 106 m3 as computed similar to recommendations by Hungr & Evans 2004; Jaboyedoff et al. 2019). Another uncertainty in terms of debris volume quantification is a probable mass loss due to river erosion or dam breach, provided that the failure actually blocked the river valley.

The geophysical results, as well as the geometrical estimations of the sliding surface depth and the debris volume, present important information for further analyses, for example, a numerical back-analysis. The latter is able to examine various parameter, such as initial slope morphology, structural aspects, geomechanical parameter, and water content of the slope (for a purely static analysis), combined with seismic scenarios (for a dynamic analysis) and helps to understand the possible factors that could have led to its actual state. A follow-up paper is currently in preparation, where we will back-analyse the massive rock slope failure through several scenario computations, involving also seismic inputs.

6 CONCLUSIONS

An extended geophysical survey of the deep-seated failure of Balta allowed us to characterize the elastic properties of its accumulated debris mass, as well as of the hill crest area above its detachment scarp.

The active seismic profiling completed over the landslide body revealed mostly homogeneously increasing Vp values at larger depths; averaged shear wave velocities could be determined through application of the 1-D MASW technique. The analysis of Poisson’s ratio suggests, consistent with the electrical resistivity results, the presence of particularly fractured and wet rock material, filled by clayey material, that originated from the flysch sandstones and schists.

The seismic surveys at the hill crest also revealed weak rock material, probably due to climatic and possibly seismic fatigue, but in contrast to the surveys on the landslide deposits, the low quality of rock appears only at shallow depths.

The flat area in direct vicinity of the detachment scarp, the plateau, was used as location for a passive seismic array configuration. By combining the dispersion curve resulting from the latter together with the MASW dispersion curve of the colocated active seismic profile, a pronounced contrast could be found at a depth of about 70–75 m below the surface. By combining passive 1-D HVSR measurements with the processed S-wave velocity profiles, the transition of fragmented rock to in situ rock could be quantified in terms of depths by detecting impedance contrasts between the two materials. With this method, we could define 1-D structures of the soft landslide material deposits; we assume the northern part of the plateau to be the area where maximum thicknesses of rockslide deposits are reached. The combination of seismic and microseismic methods indicate a maximum landslide thickness of 70–90 m.

With the help of a geomodel, the 1-D information on the surface layer was interpolated to create the body of the deposits, which allowed us to compute their full volume (estimated at 28.5–33.5 million m3). Such a geomodel is fundamental for further analyses of the rockslide; by quantifying the failed volume we are able to remodel a possible initial pre-failure state of the slope that can be used as input for a numerical back analysis.

SUPPORTING INFORMATION

Figure S1. Results of seismic ambient noise analysis performed with the seismological array: diagram of phase velocity as a function of frequency. The colors represent the occurrence of phase velocity measured at a given frequency. The surface wave dispersion curve is clearly identified between 2 and 4.5 Hz. The white dashed lines represent the domain of validity of the dispersion curve, computed on the basis of the array geometry.

Figure S2. Results of MASW applied to the S6 profile: diagram of phase velocity as a function of frequency. The colors represent the occurrence of a phase velocity to be measured at a given frequency. A fundamental mode and a first higher mode are identified, integrated with the dispersion curve obtained from array measurements and the site resonance frequency to be included in the inversion process.

Figure S3. Dispersion curves obtained through MASW analysis and fk analysis for the fundamental mode of propagation. These curves serve as target in the inversion process to retrieve the Vs structure along the profiles.

Figure S4. (a) Sensitivity analysis of S-wave velocities on Rayleigh phase wave velocities. (b) Sensitivity analysis of P-wave velocities on Rayleigh phase wave velocities.

Figure S5. Sensitivity analysis applied for each MASW results.

Table S1. Parameter space defined for the inversion of the Rayleigh wave dispersion properties obtained from MASW, the [] indicates the additional layer in case of 5 (instead of 4) layers overlying an half-space.

Table S2. Parameter space defined for the inversion of the Rayleigh wave dispersion properties obtained from array measurements and MASW applied on S6.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

We would like to thank Philippe Cerfontaine (University of Liege, Belgium) who highly contributed to this work on the field and with the related data processing. Furthermore, we thank the FRIA-FNRS for the research fellowship granted to Anne-Sophie Mreyen. Finally, we acknowledge the reviewer’s helpful comments that improved the quality of this work. This research was partly funded by the WBI (with a bilateral project ′Evaluation des risques long-termes liés aux mouvements de masse déclenchés par les séismes dans la région de Vrancea, Roumanie′), as well as by the F.R.S.–FNRS Belgium in the frame of the Belgian-Swisscollaboration project ‘4D seismic response and slope failure’. The data underlying this article will be shared on reasonable request to the corresponding author.

REFERENCES

Capon
 
J.
,
1969
.
High-resolution frequency-wavenumber spectrum analysis
,
Proc. IEEE
,
57
(
8
),
1408
1418
..

Cody
 
E.
,
Draebing
 
D.
,
McColl
 
S.
,
Cook
 
S.
,
Brideau
 
M.-A.
,
2020
.
Geomorphology and geological controls of an active paraglacial rockslide in the New Zealand Southern Alps
,
Landslides
,
17
,
755
776
.

Crosta
 
G.
,
Frattini
 
P.
,
Agliardi
 
F.
,
2013
.
Deep seated gravitational slope deformations in the European Alps
,
Tectonophysics
,
605
,
13
33
..

Crozier
 
M.
,
Deimel
 
M.
,
Simon
 
J.
,
1995
.
Investigation of earthquake triggering for deep-seated landslides, Taranaki, New Zealand
,
Quater. Int.
,
25
,
65
73
..

Crozier
 
M.J.
,
1999
.
Prediction of rainfall-triggered landslides: a test of the antecedent water status model
,
Earth Surf. Process. Landforms
,
24
(
9
),
825
833
..

Dahlin
 
T.
,
Zhou
 
B.
,
2004
.
A numerical comparison of 2D resistivity imaging with 10 electrode arrays
,
Geophys. Prospect.
,
52
(
5
),
379
398
..

Damen
 
M.
,
Micu
 
M.
,
Zumpano
 
V.
,
Van Westen
 
C.
,
Sijmons
 
K.
,
Balteanu
 
D.
,
2014
.
Landslide mapping and interpretation: implications for landslide susceptibility analysis in discontinuous data environment
, in
Proceedings of the International Conference Analysis and Management of Changing Risks for Natural Hazards
,
18–19 November 2014, Padua, Italy
, pp.
177
186
.

Danneels
 
G.
,
Bourdeau
 
C.
,
Torgoev
 
I.
,
Havenith
 
H.-B.
,
2008
.
Geophysical investigation and dynamic modelling of unstable slopes: case-study of Kainama (Kyrgyzstan)
,
J. geophys. Int.
,
175
(
1
),
17
34
..

deGroot Hedlin
 
C.
,
Constable
 
S.
,
1990
.
Occam’s inversion to generate smooth, two-dimensional models from magnetotelluric data
,
Geophysics
,
55
(
12
),
1613
1624
..

Del Gaudio
 
V.
,
Wasowski
 
J.
,
Muscillo
 
S.
,
2013
.
New developments in ambient noise analysis to characterise the seismic response of landslide-prone slopes
,
Nat. Hazards Earth Syst. Sci.
,
13
(
8
),
2075
.

Dragotă
 
C.S.
,
2006
,
Precipitaţiile excedentare în România
,
Academiei Române, Bucureşti
, pp.
53
57
.

Dunkin
 
J.W.
,
1965
.
Computation of modal solutions in layered, elastic media at high frequencies
,
Bull. seism. Soc. Am.
,
55
(
2
),
335
358
.

Earle
 
P.S.
,
Shearer
 
P.M.
,
1994
.
Characterization of global seismograms using an automatic-picking algorithm
,
Bull. seism. Soc. Am.
,
84
(
2
),
366
376
.

Fan
 
X.
,
Yunus
 
A.P.
,
Jansen
 
J.D.
,
Dai
 
L.
,
Strom
 
A.
,
Xu
 
Q.
,
2019
.
Comment on ‘Gigantic rockslides induced by fluvial incision in the Diexi area along the eastern margin of the Tibetan Plateau’ by Zhao et al. (2019) geomorphology 338, 27–42
,
Geomorphology
,
in press, doi:10.1016/j.geomorph.2019.106963
.

Farquharson
 
C.G.
,
Oldenburg
 
D.W.
,
2004
.
A comparison of automatic techniques for estimating the regularization parameter in non-linear inverse problems
,
J. geophys. Int.
,
156
(
3
),
411
425
..

Foti
 
S.
 et al. ,
2018
.
Guidelines for the good practice of surface wave analysis: a product of the interPACIFIC project
,
Bull. Earthq. Eng.
,
16
(
6
),
2367
2420
..

Gance
 
J.
,
Malet
 
J.-P.
,
Supper
 
R.
,
Sailhac
 
P.
,
Ottowitz
 
D.
,
Jochum
 
B.
,
2016
.
Permanent electrical resistivity measurements for monitoring water circulation in clayey landslides
,
J. appl. Geophys.
,
126
,
98
115
..

Georgescu
 
E.
,
2002
.
The partial collapse of Coltzea Tower during the Vrancea earthquake of 14/26 October 1802: the historical warning of long-period ground motions site effects in Bucharest
, in
International Conference Earthquake Loss Estimation and Risk Reduction
,
October 24–26, 2002, Bucharest, Romania
.

Georgescu
 
E.-S.
,
Pomonis
 
A.
,
2012
.
Building damage vs. territorial casualty patterns during the Vrancea (Romania) earthquakes of 1940 and 1977
, in
15th World Conference on Earthquake Engineering
,
24–28 September, Lisbon, Portugal
, pp.
12504
12514
.

Gercek
 
H.
,
2007
.
Poisson’s ratio values for rocks
,
Int. J. Rock Mech. Min. Sci.
,
44
(
1
),
1
13
..

Godio
 
A.
,
Bottino
 
G.
,
2001
.
Electrical and electromagnetic investigation for landslide characterisation
,
Phys. Chem. Earth, Part C.
,
26
(
9
),
705
710
..

Göktürkler
 
G.
,
Balkaya
,
Ç.
,
Erhan
 
Z.
,
2008
.
Geophysical investigation of a landslide: the Altındağ landslide site, İzmir (western Turkey)
,
J. appl. Geophys.
,
65
(
2
),
84
96
..

Grandjean
 
G.
,
Malet
 
J.-P.
,
Bitri
 
A.
,
Méric
 
O.
,
2007
.
Geophysical data fusion by fuzzy logic for imaging the mechanical behaviour of mudslides
,
Bulletin de la société géologique de France
,
178
(
2
),
127
136
..

Gregory
 
A.
,
1976
.
Fluid saturation effects on dynamic elastic properties of sedimentary rocks
,
Geophysics
,
41
(
5
),
895
921
..

Harba
 
P.
,
Pilecki
 
Z.
,
Krawiec
 
K.
,
2019
.
Comparison of MASW and seismic interferometry with use of ambient noise for estimation of S-wave velocity field in landslide subsurface
,
Acta Geophysica
,
67
,
1875
1883
.

Havenith
 
H.-B.
,
Fäh
 
D.
,
Polom
 
U.
,
Roullé
 
A.
,
2007
.
S-wave velocity measurements applied to the seismic microzonation of Basel, Upper Rhine Graben
,
J. geophys. Int.
,
170
(
1
),
346
358
..

Havenith
 
H.-B.
,
Strom
 
A.
,
Torgoev
 
I.
,
Torgoev
 
A.
,
Lamair
 
L.
,
Ischuk
 
A.
,
Abdrakhmatov
 
K.
,
2015
.
Tien Shan geohazards database: earthquakes and landslides
,
Geomorphology
,
249
,
16
31
..

Havenith
 
H.-B.
,
Torgoev
 
I.
,
Ischuk
 
A.
,
2018
.
Integrated geophysical-geological 3D model of the right-bank slope downstream from the Rogun dam construction site, Tajikistan
,
Int. J. Geophys.
,
2018
.

Heincke
 
B.
,
Maurer
 
H.
,
Green
 
A.G.
,
Willenberg
 
H.
,
Spillmann
 
T.
,
Burlini
 
L.
,
2006
.
Characterizing an unstable mountain slope using shallow 2D and 3D seismic tomography—seismic survey of an unstable mountain
,
Geophysics
,
71
(
6
),
B241
B256
..

Hibert
 
C.
,
Grandjean
 
G.
,
Bitri
 
A.
,
Travelletti
 
J.
,
Malet
 
J.-P.
,
2012
.
Characterizing landslides through geophysical data fusion: example of the La Valette landslide (France)
,
Eng. Geol.
,
128
,
23
29
..

Horowitz
 
F.G.
,
Hornby
 
P.
,
Bone
 
D.
,
Craig
 
M.
,
1996
.
Fast multidimensional interpolations
, in
26th Proceedings of the Applications of Computers and Operations Research in the Minerals Industry
,
Soc. for Mining, Metallurgy and Exploration Inc.
, pp.
53
56
.

Hungr
 
O.
,
Evans
 
S.
,
2004
.
Entrainment of debris in rock avalanches: an analysis of a long run-out mechanism
,
Bull. geol. Soc. Am.
,
116
(
9–10
),
1240
1252
..

Hungr
 
O.
,
Leroueil
 
S.
,
Picarelli
 
L.
,
2014
.
The Varnes classification of landslide types, an update
,
Landslides
,
11
(
2
),
167
194
..

Hussain
 
Y.
,
Cardenas-Soto
 
M.
,
Uagoda
 
R.
,
Martino
 
S.
,
Sanchez
 
N.P.
,
Moreira
 
C.A.
,
Martinez-Carvajal
 
H.
,
2019
.
Shear wave velocity estimation by a joint inversion of HVSR and fk curves under diffuse field assumption: a case study of Sobradinho landslide
,
Anuário do Instituto de Geociências
,
42
(
1
),
742
750
.

Ibs-von Seht
 
M.
,
Wohlenberg
 
J.
,
1999
.
Microtremor measurements used to map thickness of soft sediments
,
Bull. seism. Soc. Am.
,
89
(
1
),
250
259
.

IGR
,
1968
.
Harta geologica 1:200000, Covasna sheet
,
Institutul Geologic al Romaniei
.

Imposa
 
S.
,
Grassi
 
S.
,
Fazio
 
F.
,
Rannisi
 
G.
,
Cino
 
P.
,
2017
.
Geophysical surveys to study a landslide body (north-eastern Sicily)
,
Nat. Hazards
,
86
(
2
),
327
343
..

Jaboyedoff
 
M.
,
Chigira
 
M.
,
Arai
 
N.
,
Derron
 
M.-H.
,
Rudaz
 
B.
,
Tsou
 
C.-Y.
,
2019
.
Testing a failure surface prediction and deposit reconstruction method for a landslide cluster that occurred during Typhoon Talas (Japan)
,
Earth Surf. Dyn.
,
7
(
2
),
439
458
..

Jongmans
 
D.
,
Garambois
 
S.
,
2007
.
Geophysical investigation of landslides: a review
,
Bulletin de la Société géologique de France
,
178
(
2
),
101
112
..

Jongmans
 
D.
,
Bièvre
 
G.
,
Renalier
 
F.
,
Schwartz
 
S.
,
Beaurez
 
N.
,
Orengo
 
Y.
,
2009
.
Geophysical investigation of a large landslide in glaciolacustrine clays in the Trièves area (French Alps)
,
Eng. Geol.
,
109
(
1–2
),
45
56
..

Lapenna
 
V.
,
Lorenzo
 
P.
,
Perrone
 
A.
,
Piscitelli
 
S.
,
Sdao
 
F.
,
Rizzo
 
E.
,
2003
.
High-resolution geoelectrical tomographies in the study of Giarrossa landslide (southern Italy)
,
Bull. Eng. Geol. Environ.
,
62
(
3
),
259
268
..

Lay
 
T.
,
Wallace
 
T.C.
,
1995
.
Modern Global Seismology
,
Academic Press
.

Lemaire
 
E.
,
Mreyen
 
A.-S.
,
Dufresne
 
A.
,
Havenith
 
H.-B.
,
2020
.
Analysis of the influence of structural geology on the massive seismic slope failure potential supported by numerical modelling
,
Geosciences
,
10
(
8
),
323
.

Loke
 
M.
,
1998
.
RES2DINV version 3.3: rapid 2D resistivity and IP inversion using the least-squares method
,
Computer Disk and Manual, Penang, Malaysia
.

Loke
 
M.
,
2004
.
Tutorial: 2-D and 3-D Electrical Imaging Surveys
,
Available at: www.geoelectrical.com
.

Loke
 
M.
,
Barker
 
R.
,
1996
.
Practical techniques for 3D resistivity surveys and data inversion 1
,
Geophys. Prospect.
,
44
(
3
),
499
523
..

Luo
 
Y.
,
Xia
 
J.
,
Xu
 
Y.
,
Zeng
 
C.
,
2011
.
Analysis of group-velocity dispersion of high-frequency Rayleigh waves for near-surface applications
,
J. Appl. Geophy.
,
74
(
2–3
),
157
165
..

Mândrescu
 
N.
,
1981
.
The Romanian earthquake of march 4, 1977; aspects of soil behaviour
,
Revue roumaine de géologie, géophysique et géographie. Géophysique
,
25
,
35
56
.

Meric
 
O.
,
Garambois
 
S.
,
Jongmans
 
D.
,
Wathelet
 
M.
,
Chatelain
 
J.-L.
,
Vengeon
 
J.
,
2005
.
Application of geophysical methods for the investigation of the large gravitational mass movement of Séchilienne, France
,
Can. Geotech. J.
,
42
(
4
),
1105
1115
..

Méric
 
O.
,
Garambois
 
S.
,
Malet
 
J.-P.
,
Cadet
 
H.
,
Gueguen
 
P.
,
Jongmans
 
D.
,
2007
.
Seismic noise-based methods for soft-rock landslide characterization
,
Bulletin de la Societe Geologique de France
,
178
(
2
),
137
148
..

Meunier
 
P.
,
Hovius
 
N.
,
Haines
 
J.A.
,
2008
.
Topographic site effects and the location of earthquake induced landslides
,
Earth planet. Sci. Lett.
,
275
(
3–4
),
221
232
..

Micu
 
M.
,
Bălteanu
 
D.
,
Micu
 
D.
,
Zarea
 
R.
,
Raluca
 
R.
,
2013
.
Landslides in the Romanian Curvature Carpathians in 2010
, in
Geomorphological Impacts of Extreme Weather
, pp.
251
264
.,
Springer
.

Ministry of Environment B.
,
2018
.
Raport anual privind Starea Mediului în Romania pe anul 2018
,
Ministerul Mediului Bucuresti
,
1
678
., .

Mreyen
 
A.-S.
,
Micu
 
M.
,
Onaca
 
A.
,
Cerfontaine
 
P.
,
Havenith
 
H.-B.
,
2017
.
Integrated geological-geophysical models of unstable slopes in seismic areas
, in
Workshop on World Landslide Forum
, pp.
269
279
.,
Springer
.

Nakamura
 
Y.
,
1989
.
A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface
,
Railway Tech. Res. Inst., Quart. Rep.
,
30
(
1
).

Panek
 
T.
,
Hradecky
 
J.
,
Šilhán
 
K.
,
2008
.
Application of electrical resistivity tomography (ERT) in the study of various types of slope deformations in anisotropic bedrock: case studies from the Flysch Carpathians
,
Studia Geomorphologica Carpatho-Balcanica
,
42
,
57
73
.

Park
 
C.B.
,
Miller
 
R.D.
,
Xia
 
J.
,
1999
.
Multichannel analysis of surface waves
,
Geophysics
,
64
(
3
),
800
808
..

Pazzi
 
V.
,
Tanteri
 
L.
,
Bicocchi
 
G.
,
D’Ambrosio
 
M.
,
Caselli
 
A.
,
Fanti
 
R.
,
2017
.
H/V measurements as an effective tool for the reliable detection of landslide slip surfaces: case studies of Castagnola (La Spezia, Italy) and Roccalbegna (Grosseto, Italy)
,
Phys. Chem. Earth, Parts A/B/C
,
98
,
136
153
..

Pazzi
 
V.
,
Morelli
 
S.
,
Fanti
 
R.
,
2019
.
A review of the advantages and limitations of geophysical investigations in landslide studies
,
Int. J. Geophys.
,
2019
, doi:10.1155/2019/2983087

Perrone
 
A.
,
Lapenna
 
V.
,
Piscitelli
 
S.
,
2014
.
Electrical resistivity tomography technique for landslide investigation: a review
,
Earth-Sci. Rev.
,
135
,
65
82
..

Powell
 
M.J.
,
1992
.
The theory of radial basis function approximation in 1990
,
Adv. Numer. Anal.
,
II
,
105
210
.

Radu
 
C.
,
Spânoche
 
E.
,
1977
.
On geological phenomena associated with the 10 november 1940 earthquake
,
Revenue Romanian Geology Geophysics et Geographic.–Geophysique
,
21
,
159
165
.

Radulescu
 
N.A.
,
1941
.
Considérations géographiques sur le tremblement de terre du 10 november 1940
,
Comptes Rendus de Séances de LAcadémie des Sciences de Roumanie
,
5
(
3
),
243
269
.

Renalier
 
F.
,
Jongmans
 
D.
,
Campillo
 
M.
,
Bard
 
P.-Y.
,
2010
.
Shear wave velocity imaging of the Avignonet landslide (France) using ambient noise cross correlation
,
J. geophys. Res.
,
115
(
F3
), doi:10.1029/2009JF001538.

Rezaei
 
S.
,
Shooshpasha
 
I.
,
Rezaei
 
H.
,
2020
.
Evaluation of landslides using ambient noise measurements (case study: Nargeschal landslide)
,
Int. J. Geotech. Eng.
,
14
(
4
),
409
419
..

Rogozea
 
M.
,
Marmureanu
 
G.
,
Radulian
 
M.
,
Toma
 
D.
,
2014
.
Reevaluation of the macroseismic effects of the 23 January 1838 Vrancea earthquake
,
Rom. Rep. Phys.
,
66
(
2
),
520
538
.

Rost
 
S.
,
Thomas
 
C.
,
2002
.
Array seismology: methods and applications
,
Rev. Geophys.
,
40
(
3
),
2
1
..

Shoaei
 
Z.
,
2014
.
Mechanism of the giant Seimareh Landslide, Iran, and the longevity of its landslide dams
,
Environ. Earth Sci.
,
72
(
7
),
2411
2422
..

Strom
 
A.
,
2010
.
Landslide dams in Central Asia region
,
J. Jpn. Landslide Soc.
,
47
(
6
),
309
324
..

Strom
 
A.
,
Abdrakhmatov
 
K.
,
2018
.
Rockslides and Rock Avalanches of Central Asia: Distribution, Morphology, and Internal Structure
,
Elsevier
.

Thiery
 
Y.
,
Reninger
 
P.-A.
,
Lacquement
 
F.
,
Raingeard
 
A.
,
Lombard
 
M.
,
Nachbaur
 
A.
,
2017
.
Analysis of slope sensitivity to landslides by a transdisciplinary approach in the framework of future development: the case of la Trinité in Martinique (French West Indies)
,
Geosciences
,
7
(
4
),
135
.

Uhlemann
 
S.
,
Hagedorn
 
S.
,
Dashwood
 
B.
,
Maurer
 
H.
,
Gunn
 
D.
,
Dijkstra
 
T.
,
Chambers
 
J.
,
2016
.
Landslide characterization using P-and S-wave seismic refraction tomography—the importance of elastic moduli
,
J. appl. Geophys.
,
134
,
64
76
..

USGS
,
2020
.
U.S. Geological Survey - Earthquake Hazards Event Page. Available at: https://www.usgs.gov/natural-hazards/earthquake-hazards/earthquakes (accessed 2 September 2020)
.

Vacareanu
 
R.
,
Ionescu
 
C.
,
2016
,
The 1940 Vrancea Earthquake. Issues, insights and lessons learnt
, in
Proceedings of the Symposium Commemorating 75 Years from November 10, 1940 Vrancea Earthquake
,
Springer
.

Wathelet
 
M.
,
2008
.
An improved neighborhood algorithm: parameter conditions and dynamic scaling
,
Geophys. Res. Lett.
,
35
(
9
), doi:10.1029/2008GL033256.

Wathelet
 
M.
,
Jongmans
 
D.
,
Ohrnberger
 
M.
,
2004
.
Surface-wave inversion using a direct search algorithm and its application to ambient vibration measurements
,
Near Surf. Geophys.
,
2
(
4
),
211
221
..

Wathelet
 
M.
,
Jongmans
 
D.
,
Ohrnberger
 
M.
,
Bonnefoy-Claudet
 
S.
,
2008
.
Array performances for ambient vibrations on a shallow structure and consequences over Vs inversion
,
J. Seismol.
,
12
(
1
),
1
19
..

Whiteley
 
J.
 et al. ,
2020
.
Landslide monitoring using seismic refraction tomography—the importance of incorporating topographic variations
,
Eng. Geol.
,
268
, doi:10.1016/j.enggeo.2020.105525.

Xia
 
J.
,
Miller
 
R.D.
,
Park
 
C.B.
,
1999
.
Estimation of near-surface shear-wave velocity by inversion of Rayleigh waves
,
Geophysics
,
64
(
3
),
691
700
..

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Supplementary data