Introduction

On 12 May 2008, the Wenchuan earthquake of moment magnitude (Mw) 7.9 devastated the mountainous area of Longmen Shan in south-western China, which separates the Tibetan Plateau from the Sichuan Basin, causing over 87,000 victims. The earthquake triggered about 200,000 landslides (Xu et al. 2014), which were responsible for about one third of the earthquake fatalities (Huang and Fan 2013). Among the collateral effects of this event, one that has drawn special attention for its damaging potential is the occurrence of river damming by large landslides, which can induce severe flooding downstream in case of breaching (Fan et al. 2018 and references therein). Therefore, the study of landslide dam structure and mechanical properties is of the utmost importance, since it can allow understanding the evolution of their stability conditions (cf. Wang et al. 2014, 2018a, b).

The present work investigates the landslide dam deposit in a gully named Yang Jia Gou (Beichuan County), originated from a rock avalanche triggered by the Wenchuan earthquake. Although the landslide dam has been breached within weeks following the 2008 earthquake, the remaining deposits are marginally stable, and their erosion during monsoon seasons has led to a number of destructive debris flows (Li et al. 2021; Wasowski et al. 2021).

Due to the considerable thickness of the deposits, difficulties were expected in investigating its structure with active seismic techniques, whose implementation would require the use of a very strong energizing source of explosive type (unlikely to be authorized near marginally stable slopes) and the deployment of long linear arrays (hindered by the presence of a rough topography). Therefore, we tested an application of single-station passive seismic techniques, supported by ERT (electrical resistivity tomography), along a profile following the dam deposits exposed above the left (NE) gully wall. In particular, we used passive techniques, relying on single-station three-component recordings of seismic ambient noise, to derive curves of the ratios between the amplitudes of horizontal and vertical components of ground vibrations as function of frequency. These curves present pronounced peaks at site resonance frequencies, i.e., frequencies for which the ground response to seismic shaking presents a maximum, which can be interpreted in terms of thickness and rigidity of soft deposits overlying a stiffer substratum.

In the present study, ambient noise recordings were analysed using the standard Nakamura’s technique (Nakamura 1989) and the innovative technique HVIP (Del Gaudio 2017). However, the interpretation of the outcome of these kinds of analyses is not straightforward. There are limitations related to the sensitivity of data inversion to a number of unknown variables, which implies a wide interpretative ambiguity (e.g. Castellaro and Mulargia 2009). It is therefore necessary to combine such passive seismic surveys with other types of investigations providing some additional constraints, for instance on the geometry of the deposited layers. In particular, we used ERT surveys in order to add independent constraints useful for the modelling of the uppermost layers of the landslide dam deposit. Moreover, ERT surveys provided higher resolution information on lateral variations of mechanical properties within the first 10 m of depth of the rock avalanche deposit.

Geological setting

Regional setting and landsliding generated by the 2008 Wenchuan earthquake

The area of interest is located in the Beichuan County, in the central part of Longmen Shan, a northeast trending mountain belt in the Sichuan Province, China (Fig. 1). The literature reports clear evidence of the recent (late Pleistocene and Holocene) tectonic activity of Longmen Shan (e.g. Densmore et al. 2007). There are four major range-parallel fault zones, which represent seismogenic structures capable to generate large magnitude earthquakes (Densmore et al. 2007; Ran et al. 2013).

Fig. 1
figure 1

Location and geological features of the study area in the Sichuan province, China (upper-right) (modified after Wasowski et al. 2021); background image (2019) from Google EarthTM. Note Yangjia gully (YJG) and Tangjiawan landslides (both rock avalanches) triggered by 2008 Wenchaun earthquake, and 2008 surface rupture on nearby Yingxiu-Beichuan (Y-B) fault (rupture location from Dai et al. 2011); lower-right inset shows geological formations (based on 1:200,000 geological map of the China Geological Survey CGS 2001), and source area of YJG rock avalanche (white circle). Legend: 1 Maoxian Fm. Silurian; 2 Loureping and Shamao Fm. - Silurian; 3 Longmaxi Fm. Lower Silurian; Baota Fm.—Middle Ordovician; 5 Qingping Fm.—Cambrian; 6 Upper Sinian; 7 fault, dashed when inferred; 8 headscarp of 2008 rock avalanche

Full size image

The May 12, 2008 Mw7.9 Wenchuan earthquake was the most recent example of the seismic activity of the region. The event involved two major faults, which generated over 300 km long northeast trending surface rupture and peak ground acceleration (PGA) locally approaching 1 g (Li et al. 2008; Wen et al. 2010). In situ investigations of the surface rupture indicated both reverse and strike-slip components in the Beichuan area, which is distant about 135 km away from the epicentre (Liu-Zheng et al. 2009).

The 2008 Wenchuan earthquake induced numerous co-seismic hazards, including tens of large (volume > 1 million m3) rock avalanches (Qi et al., 2010) and in total about 197,000 landslides (Xu et al. 2014). Moreover, the event produced several hundreds of landslide dams (Xu et al. 2009; Fan et al. 2012). More information on the geological hazards caused by the 2008 earthquake and additional relevant literature can be found in an overview article by Fan et al. (2018).

The Beichuan County was among the hardest-hit areas by the Wenchuan earthquake. There, the Yingxiu-Beichuan fault surface rupture produced PGA values that exceeded 0.4 g (Li et al. 2008). In addition, few strong aftershocks (Mw > 6) occurred in the Beichuan area. Tang et al. (2011) studied regional-scale landsliding generated by the 2008 earthquake in the north-eastern part (414 km2) of the Beichuan County crossed by the Yingxiu-Beichuan fault. They documented high co-seismic landslide frequency amounting to 5.4/km2.

Local setting

The area studied is about 20 km to the northeast of the town of Beichuan demolished by the 2008 earthquake (Fig. 1). The 2008 earthquake surface rupture (Fig. 1), mapped on the west valley side of the Duba River (Dai et al. 2011), is only 1.5 km away. The area is drained by the Yangjia gully torrent, a local tributary of the Duba River, which flows in a large valley that follows the Yingxiu-Beichuan fault zone.

The landscape around the Yangjia gully is that of moderately high mountains and deeply incised river valleys with the elevations ranging from about 700 m to over 1300 a.s.l. Slopes are typically steep and locally very steep (> 45°).

Local geology bears the imprint of the northwest-dipping Longmen Shan thrust zone and the nearby active Yingxiu-Beichuan fault (Fig. 1). The faults and associated major fractures are locally superimposed on the stratigraphic discontinuities (former bedding) or boundaries between different lithological units, which are affected by low-grade metamorphism. There seem to be two predominant orientations of the faults reflecting the deviation of the Yingxiu-Beichuan fault toward northerly direction from the general NE-SW trend of the Longmen Shan thrust belt. The main and secondary lithologies are, respectively, greyish-beige to black slates and light-grey meta-limestones belonging to Lower Paleozoic age formations (Fig. 1). The rocks are typically intensely fractured and sheared. This and the presence of the high relief make the slopes susceptible to seismically triggered failures.

The Yangjia gully rock avalanches

The impact of the 2008 earthquake on the local landscape has been discussed in detail by Wasowski et al. (2021), who documented the occurrence of 68 landslides in an area of 3 km2 surrounding the Yangjia gully. The Yangjia gully rock avalanche was the largest among the 2008 co-seismic slope failures occurred in the study area. The avalanche deposit dammed the Yangjia gully stream and the source of the avalanche was at the top of the north-facing slope of the Weijia Mountain, at an elevation of 1250 m a.s.l., while the slope base (Yangjia gully base level) was at ~835 m a.s.l. (Fig. 1). The maximum width and length of the rock avalanche deposit were, respectively, 400 and 880 m, and the overall area amounted to about 170,000 m2 (Wasowski et al. 2021). In situ observations of the exposures along the Yangjia gully indicated that the average thickness of the landslide dam deposit could be around 20–30 m. However, this estimate is uncertain, because the base of the 2008 rock avalanche is only locally exposed at its distal eastern-most portion (Fig. 2), and the reconstruction of the pre-failure landscape is hampered by the lack of access to pre-earthquake topographic maps. Considering the tectonic activity of the area and the occurrence of large magnitude earthquakes, and hence high rate of geomorphic processes, we suspect that the morphology of the base of the 2008 rock avalanche deposit could locally be complex. Moreover, the 2008 rock avalanche overlaps at least in part the deposit of an older (perhaps pre-historic) rock avalanche (Wasowski et al. 2021). The contact of the base of the older rock avalanche with slates that form the local bedrock is only exposed in one isolated outcrop in a small stream channel, near its confluence to the Yangjia gully (Figs. 1 and 2). The slate bedrock moderately dips to the west and is not visible along the Yangjia gully walls upstream of the confluence. However, the slate bedrock crops out at the north-western distal part of the 2008 rock avalanche deposit (see Fig. 3). In this area, the bedrock dips steeply (>45°) to the south-southwest.

Fig. 2
figure 2

(upper) Downstream, eastern portion of the landslide dam near the confluence of a small stream channel to the Yangjia gully. Slate bedrock dips to the west and disappears under the old avalanche deposit (no. 1), which is overlaid by the 2008 co-seismic rock avalanche deposit (no. 2). Note several piles and parts of foundations of two check dams destroyed by the June 2018 flood-debris flow; (lower) Landslide dam deposits along the left (NE), up to 60 m high Yangjia gully wall; photo shows the thickest exposed portion of the dam coinciding, to the left, with the downstream end (B) of the ERT profile, and to the right with ambient noise station YIG8 (cf., Fig. 3)

Full size image
Fig. 3
figure 3

2018 UAV image showing the upstream, north-western part of the landslide dam deposit exposed along the left (NE) Yangjia gully wall with the positions of ambient noise measurement stations (identified by yellow triangles), the reference station (identified by a blue point) and the ERT profile A-B (red line). Note the landslide dam lake (upper left) and slate bedrock steeply dipping to S-SW

Full size image

The material properties of the rock avalanche dam reflect the characteristics of the source rocks, i.e. intensely sheared and fractured slates. On its rapid decent downslope, the rock avalanche mass disintegrated into coarse and fine debris and formed a poorly consolidated and loose landslide dam deposit. Most recently, Li et al. (2021) reported the particle size distributions of the Yangjia gully landslide dam deposit, which indicate the predominance of gravel size angular clasts and sand; the silt size fraction amounts to approximately 5%. Hence, this deposit can be described as cohesionless material. The landslide dam deposits cropping out along the gully walls show marginal stability and high susceptibility to erosion (Figs. 2 and 3).

Data acquisition and analysis methodology

After a preliminary test carried out in October 2017, noise measurements were repeated in three campaigns (May and October 2018, October 2019), using a set of 3 tromographs Tromino, produced by MoHo s.r.l. and, in particular, two models ENGY PLUS and one model 3G, equipped with velocity meter sensors operating in a frequency range of 0.1–300 Hz at its maximum sensitivity (10-9 m/s) and with a full scale of 1.2 mm/s (for more details see https://moho.world). These tromographs are all-in-one instruments including a 24-bit acquisition system with a selectable sampling frequency. The latter was set to 128 Hz in our measurement campaigns, and data were acquired keeping the 3G sensor, having a longer autonomy, in continuous recording at a reference station, named YJG0. The other tromographs were used for 30–40-min recording sessions at 8 different stations, named YJG1-8, located along a 325-m-long profile on the left bank of the river (see Fig. 3). The arrangement of a reference station made it possible to distinguish between signal differences related to station-specific site conditions from the variations that may reflect temporal changes in the properties of the noise sources during the recordings.

Acquired data were first analysed using the standard Nakamura’s technique, also known as HVNR (Nakamura 1989), which calculates, on a series of temporal intervals of few tens of seconds, the average spectral ratios H/V between the horizontal and vertical component of noise recording. In particular, we adopted a subdivision of recordings into time windows of 20 s, processed following the guidelines recommended by the SESAME project (SESAME 2004), which implies the removal of time windows with transient signals and the smoothing of spectra, according to the method proposed by Konno and Ohmachi (1998). Mean HVNR values were calculated setting as numerator the amplitude of horizontal components rotated at 10° intervals, in order to analyse the directional variations of spectral ratios.

The HVIP technique (Del Gaudio 2017) aims at identifying, within the recordings of ambient noise, wave packets in which Rayleigh waves are dominant. They are recognised by means of an instantaneous polarization analysis based on the analytic transformation:

$$\begin{aligned}\overrightarrow{{u}_{c}}(t)=\overrightarrow{u(t)}+j \kern 0.03000 em \kern 0.03000 em\overrightarrow{\widehat{u}(t)}=\overrightarrow{A(t)}{e} \kern 0.03000 em \kern 0.03000 em \kern 0.03000 em^{j\Phi (t)}\end{aligned}$$
(1)

where \(\overrightarrow{\widehat{u}(t)}\) is the Hilbert transform of the 3-component noise recording \(\overrightarrow{u(t)}\) and j is the imaginary unit.

Morozov and Smithson (1996) demonstrated that ground motion can be described at each instant as consisting of portions of elliptical trajectory whose major and minor semi-axes, \(\overrightarrow{a}(t)\) and \(\overrightarrow{b}(t)\), have direction and length that can be calculated at each instant from the analytic transformation (1) through the formulae:

$$\overrightarrow{a}(t)=\mathit{Re}\left[{e}^{-j{\Phi }_{o}}\cdot {\overrightarrow{u}}_{c}(t)\right]$$
(2)
$$\overrightarrow{b}(t)=\mathit{Re}\left[{e}^{-j\left({\Phi }_{o}+\frac{\pi }{2}\right)}\cdot {\overrightarrow{u}}_{c}(t)\right]$$
(3)

where

$${\Phi }_{o}=\frac{1}{2}\mathit{arg}\left[\frac{1}{2}{\sum }_{k}{\left({u}_{k}+j{\widehat{u}}_{k}\right)}^{2}\right]$$
(4)

The HVIP technique exploits the results of this instantaneous polarization analysis to identify coherent Rayleigh wave packets, for which the plane of the elliptical trajectory is close to vertical and the major semi-axis \(\overrightarrow{a}(t)\) is close to horizontal (for H/V > 1) or to the vertical (for H/V < 1). Operatively, a three-component noise recording (Fig. 4a–c) is first passed through a narrow-band filter centred on different frequencies (e.g. Fig. 4d–f), then is subjected to the analytic transformation to identify the time intervals during which the instantaneous elliptical trajectory is close to vertical, and the major semi-axis \(\overrightarrow{a}(t)\) is close to horizontal (or to vertical), as effect of the presence of Rayleigh wave packets prevailing on other wave types. It is then possible to restrict the calculation of the ratio between the instantaneous amplitude of horizontal and vertical component of ground motion just to such intervals, reducing the influence of other wave types on the measurement. The application of this technique has shown that the condition of Rayleigh wave dominance actually occurs for very short intervals (see Fig. 4g, h), typically in the order of 1–2% of the recording (Del Gaudio et al. 2017). Nevertheless, the estimate of the H/V ratio can be averaged on a large number (in the order of thousands) of instantaneous values. The use of an analysis providing H/V estimates for each recording sample allows a higher resolution in the separation of Rayleigh from other waves. This distinguishes the HVIP method from other techniques (e.g. Jurkevics 1988; Hobiger et al. 2009), which require processing a number of consecutive data samples for Rayleigh wave identification.

Fig. 4
figure 4

Rayleigh wave identification through HVIP analysis: a, b, c component of ambient noise recording; d, e, f result of narrow-band filtering centred at 2.95 Hz; g instant by instant estimates of deflection dipp of the plane of ground motion elliptical trajectory from the vertical; h inclination dipa of the ellipse major semi-axes. Rayleigh wave samples (marked by grey stripes) are recognized when both dipp and dipa are simultaneously close to 0°

Full size image

The resulting average of H/V ratios, calculated for different frequencies, can be interpreted as representative of the variations of Rayleigh wave ellipticity as function of frequency. The ellipticity curves are characterized by pronounced peaks at site resonance frequencies caused by the impedance contrast between surface soft layers and a stiffer substratum. These curves can provide information on properties of subsoil material, i.e. (i) layer thickness and velocity from the resonance frequencies, (ii) impedance contrast between surface layer and bedrock from the H/V peak amplitudes and (iii) deposit mechanical anisotropy from directional variation of the H/V peak amplitude. In simple cases, i.e. in presence of a single soft layer overlying a stiffer substratum, the information on (i) can be obtained from the fundamental resonance frequency through the so-called quarter wave-length law, expressed as (Haskell 1960)

$${\mathrm{f}}_{0}=\frac{{\mathrm{V}}_{\mathrm{s}}}{4\mathrm{H}}$$
(5)

where f0 is the fundamental resonance frequency identified from the H/V peak, VS is the velocity and H the thickness of the surface layer. In presence of more layers, causing multiple resonance peaks, the entire ellipticity curve can be used to constrain a subsoil velocity model. In the inversion of the H/V curve, the peak amplitudes provide additional constraints to the velocity model, since they are sensitive to the S-wave velocity contrast between different layers and to the Poisson’s ratio of the same layers.

In comparison to the HVNR method, the HVIP technique proved to provide more stable results in the determination of H/V peak amplitudes and more details on site resonance frequencies (Del Gaudio 2017; Del Gaudio et al. 2018, 2019, 2021). The resulting H/V curve can be more reliably interpreted in terms of Rayleigh wave ellipticity curve in comparison to the curve derived from H/V spectral ratios by the Nakamura’s method. Indeed, in the HVNR curve, one cannot distinguish the contribution of different types of waves to the horizontal and vertical spectrum, which can change at the same site in different conditions of ambient noise generation.

However, the inversion of Rayleigh wave ellipticity curves in terms of sub-soil velocity model is not univocal, since one can obtain the same curve by different combination of S-wave and P-wave velocities. Therefore, to solve data interpretation ambiguity, independent constraints are needed. In the present study, ERT surveys were therefore carried out as support to the interpretation of ambient noise data.

ERT technique is a well-established methodology commonly applied in landslide investigations (e.g. McCann and Forster 1990; Gallipoli et al. 2000; Hack 2000; Perrone et al. 2014; Bellanova et al. 2016). It can provide useful information on the geometrical characteristics of landslide deposits and on surrounding potentially instable slope areas.

Resistivity measurements are typically carried out by deploying an array of steel electrodes connected to a geo-resistivity meter, a multichannel apparatus which controls the injection into the ground, through couples of electrodes, of a direct current delivered by an about 100 Watt power source, and which measures the potential drop at another couple of electrodes by a voltmeter. The multichannel system can select combinations of arrays of current and measurement electrodes according to different configuration (e.g. Wenner, Schlumberger, dipole-dipole), differing for the relative position of the two electrode pairs and having different performance in terms of investigation depth and lateral resolution. For each combination, an apparent resistivity value ρa is calculated from the equation

$${\rho }_{\mathrm{a}}=\mathrm{K}\frac{\Delta \mathrm{V}}{\mathrm{I}}$$
(6)

where I is the intensity of the injected current, ΔV is the potential drop and K a geometric coefficient depending on electrode distances.

The quantity (6) obtained for each electrode combination is considered representative of the subsoil resistivity at a point located below a middle position of the electrode array, at a depth depending on the electrode maximum distance. Several measurements are carried out for each point, averaging the results after removing the extreme (maximum and minimum) values.

In this study, to overcome the intrinsic limitations of each electrode configuration due to the trade-off between investigation depth and horizontal resolution, both Alpha Wenner and dipole–dipole configurations were used (Perrone et al. 2014); in both cases, an array of 120 electrodes was used. During two different survey campaigns, carried out in 2018 and 2019, we tested comparatively the results obtained by arrays with different spacing (2–5 m). Figure 3 shows the position of the ERT profile A-B.

To obtain a subsurface image of the electrical resistivity, the apparent electrical resistivity data have to be inverted into true electrical resistivity values by means of specific inversion software. For this purpose, we used the software package RES2DINV (Loke 2001) to obtain 2D electrical resistivity images of the subsurface. The inversion routine is based on the smoothness-constrained least-squares inversion method, implemented by using a quasi-Newton optimisation technique (Sasaki 1992; Loke and Barker 1996). The optimisation method adjusts the 2D electrical resistivity model trying to iteratively reduce the difference between the calculated and measured apparent resistivity values. The root-mean-squared (RMS) error provides a measurement of this difference.

Results

Ambient noise recording analyses

Noise data at the reference station YJG0 were acquired for several hours in all the three measurement campaigns (3 h 48 m, 3 h 14 m, 5 h 36 m, respectively); thus, it was possible to verify the stability of the results of noise analysis carried out with the two techniques (HVNR, HVIP). Figure 5 shows diagrams summarizing the H/V spectral ratios and Rayleigh wave ellipticity as function of frequency and azimuth, averaged over each measurement campaign and over the entire set of measurements. While both HVNR and HVIP analysis provided consistent evidence of a significant resonance effect on a frequency band of 3–4 Hz, possibly related to the combination of two resonance peaks at two close frequencies, the two methods showed differences in the peak amplitude estimations and in the stability of the results of these estimations over different measurements.

Fig. 5
figure 5

Polar diagrams of H/V spectral ratios obtained by the HVNR method (top) and of the Rayleigh wave ellipticity obtained by the HVIP technique (bottom), from the analysis of noise recorded at the reference station YJG0 in three measurement campaigns (May and October 2018, October 2019). H/V ratios are represented by a colour scale as a function of radially plotted frequencies, along different directions. The polar diagram of H/V values averaged over the entire set of recordings are also reported, together with plots of the curves of the HVNR or HVIP values along the directions of major peaks resulting from the analysis of data of each measurement campaign. Curves corresponding to differently oriented major peaks derived from the average of all the data from three measurement campaigns are also shown. Legends specify date and azimuth to which the curves are referred

Full size image

The major HVNR peak has an amplitude varying from a minimum of 5 to a maximum of 12, whereas the peaks resulting from the HVIP analysis are comprised between 5 and 6. However, the mean amplitude is quite similar (5.96 for HVNR, 5.76 for HVIP), which confirms our previous observations (Del Gaudio et al. 2018, 2019, 2021) about the lower reliability of the HVNR amplitude estimates derived from data acquired in a single campaign with respect to the HVIP technique.

The greater variability of the HVNR results does not seem to depend on the effect of random short-term changes of noise wavefield properties. Indeed, the consistency of results obtained in different periods does not improve by averaging H/V ratios over long time intervals, which in our study were longer than those recommended by the SESAME guidelines (about half an hour). Thus, it is likely these variations are caused by changes of environmental conditions influencing the contribution of different type of waves (P, S, Rayleigh, Love) to the noise wavefield. This makes more instable the results derived from techniques (e.g. HVNR) that do not separate the contributions of different waves.

As additional difference between the two techniques, the HVIP analysis provided a clear evidence of a significant secondary resonance effect at higher frequency (amplitude ~ 4 at 10 Hz: see the grey curve in the bottom right diagram of Fig. 5), which would be considered negligible based on the HVNR results alone (maximum amplitude of only 2.5 at 10.8 Hz). This confirms the better ability of HVIP to reveal secondary resonance effects.

Considering the above observations, the data interpretation that follows refers to the results of the HVIP analysis. Figure 6 shows the results obtained at eight stations, from YJG1 to YJG8, as average of the measurements acquired in all the campaigns, together with the curves of ellipticity along the azimuth of major peaks observed in each campaign. All these results show a common pattern with a major peak at a relatively lower frequency F0 (between 2 and 4 Hz) and a secondary peak at a higher frequency F1 (from about 5 to 12 Hz).

Fig. 6
figure 6

Polar diagrams of HVIP values averaged over the results from the three campaigns of noise recordings at stations YJG1 through YJG8 and diagrams of HVIP curves along the directions of major peaks. Curves observed in each measurement campaign (colour curves) and on the average of all the measurements (black thick curve) are reported. Diagram legends specify date and azimuth to which curves are referred

Full size image

The main characteristics of the major and secondary peaks are summarized in Table 1. Along the measurement array, frequency F0 first shows a gradual decrease from 4.05 at YJG0 to a minimum of 2.20 Hz at YJG5, followed by an increase up to 3.55 Hz at YJG7 (Fig. 7a). A similar trend is observed for the amplitude A0 of the peak of ellipticity found at this frequency, even though with some local oscillations (Fig. 7b).

Table 1 Main peaks of Rayleigh wave ellipticity derived from the HVIP analysis of ambient noise recorded at the measurement stations located on the rock avalanche deposit in Yangjia gully. Peak parameters include the frequency in Hertz for lower (F0) and higher (F1) resonance frequencies, the amplitude of ratio H/V, the azimuth of the direction of maximum and a directivity index (Idir) given by Eq. (7)
Full size table
Fig. 7
figure 7

Variation, along the noise measurement array, of a the lower resonance frequency (F0) and b the amplitude A0 of the H/V peak at frequency F0. For comparison, the topographic profile of the investigated section is also shown

Full size image

The frequency F0 should reflect site fundamental resonance caused by a relatively deep impedance contrast. According to Eq. (5), the variation of resonance frequency is inversely correlated with the thickness of surface layer and directly correlated with its S-wave velocity (Vs). On the other hand, the peak amplitude variations are expected to be correlated with the velocity contrast between surface layer and substratum. Thus, the simultaneous decreasing trend of both F0 and A0 observed from YJG0 to YJG5 cannot be explained as effect of lateral decrease of Vs in the surface layer alone, since it would imply also an increase (instead of decrease) of A0. Therefore, an increase in the thickness of the rock avalanche deposit provides the most likely explanation. This is also consistent with the increase in the surface elevation in the landslide accumulation zone, as one can notice from a comparison with the topographic profile (see Fig. 7a). In the stations from YJG5 to YJG8, the increase of F0 and A0 could be due to a decrease of the deposit thickness and/or lateral variations of velocity both in the surface layer and in the substratum.

Frequency F1 also shows apparently gradual variations (Fig. 8a), with two relative maxima (~12 and 9 Hz at YJG1 and YJG6, respectively) and two minima (6.0 and 5.5 Hz at YJG5 and YJG7, respectively), whereas the corresponding ellipticity peak amplitude changes in a somewhat irregular way (Fig. 8b). The frequency F1 is likely related to an impedance contrast between a thin surface layer of material overlying a deeper and stiffer part of the rock avalanche deposit. Thus, it is likely heterogeneous materials with lateral variations of thickness and/or velocity characterise this layer.

Fig. 8
figure 8

Variation, along the noise measurement array, of a the higher resonance frequency (F1) and b the amplitude A1 of the H/V peak at frequency F1. For comparison, the topographic profile of the investigated section is also shown

Full size image

The directional properties of site resonance can be evaluated from the polar diagrams of Figs. 4 and 5 and the directivity index Idir reported in Table 1. This index is calculated as

$${I}_{dir}=\frac{{H/V}_{max}\left({F}_{x}\right)}{{H/V}_{min}\left({F}_{x}\right)}$$
(7)

where H/Vmax(Fx) and H/Vmin(Fx) are the maximum and minimum ellipticity, respectively, observed at the resonance frequency Fx (with x = 0 or 1) among different directions. Del Gaudio et al. (2008) used a similarly defined parameter in the analysis of HVNR results and found it to be indicative of a significant directivity in site response when such a ratio exceeds 1.5. According to this criterion, the resonance F0 presents a significant directivity at five stations, from YJG0 to YJG5. Apart from YJG0, the other stations also revealed a consistent orientation of maximum direction (with azimuth in the range 55°–85°). The reasons for lack of significant directivity of F0 resonance at stations YJG6–YJG8 are unclear. This section, however, coincides with the thickest exposed portion of the 2008 landslide dam and the collision zone where the 2008 co-seismic avalanche collided against the remnant of a pre-existing landslide dam formed by an ancient rock avalanche (Wasowski et al. 2021). Here, the inferred overlap surface of the two landslide dam deposits dips at moderate angles towards south, whereas elsewhere (e.g. at stations YJG1-YJG5), the contact is sub-horizontal.

The resonance F1 shows a major variability of directivity properties, with an almost equal number of sites with and without evidence of site response directivity and with a major variability of directions of H/V maxima. Considering that directivity can reflect anisotropy of mechanical properties, this variability is indicative of heterogeneity of such properties in the most surficial layer of the rock avalanche deposit.

Overall, the ambient noise data are consistent with the presence of superimposed layers, each characterized also by lateral variations of geometrical and/or mechanical properties. However, without additional independent information (e.g., from ERT surveys), it is not possible to identify the exact nature of these variations (thickness, stiffness controlling S-wave velocity or both).

ERT data

Figure 9a, b show the results of the ERT surveys carried out in 2018 and 2019 using the Alpha Wenner array configuration. The subsoil resistivity models were obtained by inverting the apparent resistivity, taking into account topographic correction, through the software RES2DINV using the standard least squares smoothness model bound. The maximum investigation depth reached in the central portion of the profile is about 40 m.

Fig. 9
figure 9

2D resistivity profiles resulting from the inversion of ERT data acquired in two surveys (a 2018, b 2019) carried out along the Yang Jia Gou river profile A-B (see Fig. 3). Black vertical lines represent the thickness of the uppermost surficial layer inferred from ERT cross section at the sites of ambient noise recordings. The black dashed line represents the streambed profile whose direction is sub-parallel to the ERT layout

Full size image

The results of the two surveys show a similar resistivity variation pattern along the investigated profile. The uppermost portion of the profiles contains more resistive (greater than 50 ohm · m) and heterogeneous material. This is consistent with a rock avalanche origin. Rock avalanche deposits often show the presence of coarser rocky material at their top (Dufresne et al. 2016). The deeper portion of both profiles shows lower resistivity values (less than 20 ohm · m), which could be related to water-bearing finer sediments or to the presence of a higher percentage of clayey material. Indeed, resistivity of sub-soil material depends on a series of factors such as porosity, degree of water saturation and concentration of clayey minerals, and the decrease of resistivity can reflect the passage to more porous deposits with higher water content or a change in material composition.

It is apparent that the upper boundary of the conductive layer does not correspond to the streambed level (see Fig. 9). Therefore, it is likely that the low resistivity can be due to the infiltration of water migrated laterally from the outer (NE) to inner parts of the gully (cf. Fig. 3). Moreover, resistivity data alone do not allow us to detect the bedrock depth, but suggest a possible presence of a surface separating deposit’s layers with different characteristics. As the investigations have always been carried out in periods following the rainy season, it is apparent that the presence of water in the subsoil does not help identifying the possible separation surface between the different layers of deposit.

Additional resistivity data were acquired adopting the dipole-dipole electrode configuration, which can provide a better lateral resolution in the characterization of properties of surface layers. This is obtained at the expense of a decrease in the investigation depth to only 20 m (Fig. 10). The results show, within the first 5–10 m, strong lateral variations of resistivity with the nuclei of high resistivity material within a matrix of less resistive deposit. This is consistent with the presence, within the uppermost part of the 2008 rock avalanche deposit, of chaotically distributed blocks of rocky material immersed in a “matrix” of finer fragmented material.

Fig. 10
figure 10

2D dipole-dipole resistivity profile resulting from the inversion of ERT data acquired along the Yang Jia Gou river profile A-B (see Fig. 3)

Full size image

This uppermost layer overlies a relatively more continuous and more conductive layer. There is an apparent increase of resistivity in the deepest part of the section. However, such feature could be an artefact due to the poorly constrained definition of subsoil properties near the limit of the investigation depth, as no resistivity increase at depth resulted from the Wenner profile inversion (see Fig. 9).

Numerical modelling of ellipticity curves

Further insights on the physical-geometrical characteristics of the rock avalanche deposit were obtained by carrying out the inversion of Rayleigh wave ellipticity curves derived from the HVIP analysis. In order to search models consistent with the experimental data, for each station, a mean ellipticity curve was used as a target of inversion, exploiting data acquired during all the measurement campaigns.

Sub-soil velocity models consistent with the ellipticity curves were then obtained using the dinver module of the Geopsy open source software (Wathelet et al. 2020). It adopts a Neighborhood Algorithm approach to explore efficiently the parameter space, subdivided according to Voronoi cells, in search of misfit minimization around a number of best solutions found at each iteration of the inversion process. The code finally provides a set of 1D models compatible with observation within a certain level of measurement uncertainty (for more details, see Wathelet 2008). The use of a 1D modelling of the sub-soil velocity is an approach commonly followed to interpret the results of single-station methods of ambient noise analysis (cf. Castellaro and Mulargia 2009). Indeed, the ratio between the amplitude of horizontal and vertical components of ground motion mainly reflects the terrain seismo-stratigraphic properties within a limited area around the measurement site.

A problem for the use of 1D modelling arises, at some sites, from the evidence of directional variations of H/V curves, related to possible anisotropy of the mechanical properties of the rock avalanche deposit. All H/V curves show two major peak values (see Table 1), one centred on a relatively lower frequency (2–4 Hz) and the other on higher frequencies (5–12 Hz), which supports the hypothesis of stratification inside the deposit. The amplitudes of peaks can show a considerable azimuthal variability and generally, at the same site, the two peaks reach a maximum along different azimuths. Therefore, in order to obtain velocity models representing average properties of the ground below each measurements site, a single target curve was defined for each station by averaging the HVIP values of the three campaigns along the different azimuths, regardless of the site response directivity properties. The central frequency and amplitude of the two peaks that can be recognized in such average curves are summarized in Table 2.

Table 2 Characteristics of the resonance peaks obtained at each noise measurement station from curves of HVIP values averaged over different azimuths and the three measurement campaigns. F0 and F1 values identify the lower and the higher resonance frequency at each station, respectively
Full size table

This kind of data adjustment introduces some variations in the amplitude and frequency trend along the measurement profile. In particular, A0 tends to increase along the profile (Fig. 11b), instead of showing a minimum in the central part of the profile (at YJG5: see Fig 7b). This depends on the presence of a more pronounced directivity affecting the low frequency response in the first half of the profile, so that the averaging of amplitudes over different directions causes a considerable decrease of the resulting A0 value, in comparison to the more isotropic response at sites from YJG5 to YJG8. F0, F1 and A1 show smaller variations, and their trend is more similar to that shown by the directional peak values of HVIP (compare Figs. 7a and 8a, b with Figs. 11a and 12a, b).

Fig. 11
figure 11

Lower resonance peak frequencies F0 (a) and relative amplitudes A0 (b) derived from directional average of HVIP values for each noise measurement station, compared with the topographic profile

Full size image
Fig. 12
figure 12

Higher resonance frequencies F1 (a) and their relative amplitudes A1 (b) derived from directional average of HVIP values for each noise measurement station, compared with the topographic profile

Full size image

To overcome the interpretative ambiguity of single-station passive seismic techniques and introduce constraints in the modelling, the ambient noise results were combined with ERT data. The electrode array length allowed covering the part of the rock avalanche deposit between the stations YJG1 and YJG6, and for this portion of the array, it was possible to use data from ERT for constraining the inversion of ellipticity curves.

The outcomes of both Wenner and dipole-dipole profiles show the deposit composed of a superficial more resistive layer and one less resistive layer below. The Wenner profile has a better vertical resolution and was used to estimate the thickness of the superficial resistive layer at each station. The thickness estimates at the two stations located at the ends of the ERT array (YJG1 and YJG6) are less reliable, because the results of inversions are affected by anomalies related to boundary effects.

The high-frequency peak marker represents the seismic response of the most superficial layer. By comparing the frequency values and the layer thickness for each station, one can note that their variations are consistent with Eq. (1), which implies a frequency decrease as the thickness increases (Table 3). The resulting estimates of Vs for the most superficial layer give indications of a relative homogeneity of stiffness at least in this part of the deposit, so that the higher resonance frequency appears mainly controlled by the deposit thickness.

Table 3 Comparison of the superficial layer thickness, based on the Wenner resistivity profile, and the peak higher frequencies. Values of mean Vs of the surface layer are estimated from Eq. (1)
Full size table

To obtain more details on the subsoil model of the rock avalanche deposit, the ellipticity curves were inverted in terms of velocity vertical profile. The inversion code requires the definition of an initial parameterization for a fixed number of layers, with the possibility of defining a more or less wide range of admissible values for layer thicknesses and velocities. After some trial tests, we adopted a 6-layer model. For the stations from YJG2 to YJG5, the thickness of the surface layer was constrained based on the ERT results, leaving broader ranges of possible values for velocity and thickness of other layers. The results obtained for these stations guided the parameterization of the subsoil velocity model for the other stations of the profile.

The output of inversions of single station is reported in the supplementary material (from Figs. S1 to S9), and a schematic cross-section based on the interpolation of vertical profiles of S-wave velocity is shown in Fig. 13. The results indicate the presence of a surficial layer with velocity of 200–400 m/s and thickness that varies between 4 and 17 m. Underneath is a layer with a velocity up to 600 m/s and lateral variations in thickness ranging from a minimum of 3 m to over 20 m. Further below, at depth varying from about 20 to 60 m, velocities larger than 800 m/s are reached. These variations could reflect irregularities of the topographic surface over which the 2008 rock avalanche deposited its material.

Fig. 13
figure 13

Characteristics of Rayleigh wave ellipticity peaks (a) at each noise measurement station of the profile traversing the landslide dam deposit (arrowheads mark orientation of directional resonance, circles denotes more isotropic site response, colours and symbols’ sizes arranged to represent frequency and amplitude of H/V peak) and the corresponding Vs cross-section (b) derived from inversion of H/V curves. The black dashed line represents the streambed profile sub-parallel to the direction of Vs cross-section

Full size image

By considering the streambed profile (Fig. 13b), it becomes evident that the layers with velocities lower than 800 m/s are those exposed along the gully walls. In situ inspections of this part of the landslide dam revealed no evidence of a superimposition of deposits from two distinct rock avalanche events. Therefore, it is likely that within the Vs range between 200 and 800 m/s, the vertical changes in velocity reflect a layering within the deposit of the last (2008) rock avalanche. The underlying layer with velocity over 800 m/s could perhaps represent the substratum of this deposit, corresponding to the material deposited by an older rock avalanche. Figure 13b indicates that the interface between the two rock avalanche deposits should crop out at the foot of the gully wall section corresponding to the location of the stations (YJG6-8). In this section, however, the base of the gully is covered by the 2008 rock avalanche material fallen from the unstable upper gully walls (cf., Figs. 2 and 3).

Discussion and conclusions

Figure 13 summarises the results of this study, in terms of site response properties (resonance frequency and directivity) and of inferred subsoil properties (S-wave velocities), derived from the analysis of ambient noise recordings carried out through the new technique of determining the Rayleigh wave ellipticity from instantaneous polarization properties (HVIPs). In comparison to the standard HVNR technique (Nakamura 1989), devised to interpret noise data on site conditions characterized by flat horizontal isotropic layering, the new technique provided more detailed and stable results, revealing the complex site resonance properties of the Yang Jia Gou rock avalanche deposit. Indeed, although the HVNR technique is able to identify the fundamental resonance frequency, which provides information on mean velocity and thickness of relatively soft deposits overlying a stiffer substratum, this technique failed in revealing higher frequency resonances related to the presence of layering in the 2008 rock avalanche deposit. Furthermore, the amplitude of H/V peaks derived from HVNR analysis, which could provide information on the velocity of substratum, appears altered by variations of the noise wavefield composition. This can include different types of waves (P, S, Love, Rayleigh) in unknown proportion, causing considerable uncertainties in inferring velocity models from the H/V curve inversion.

The HVIP technique, isolating only the contribution of Rayleigh waves to the H/V ratios, provides more stable results, which can be more reliably interpreted. Its capacity of estimating the H/V ratios on an instantaneous basis allows obtaining more details on complex resonance pattern (e.g. multiple peaks at relatively similar frequencies), which the spectral average over time windows of a few tens of seconds, adopted by the HVNR technique, could hide.

Both methods, however, suffer from interpretative ambiguities related to the fact that resonance frequencies depend on the ratio between velocities and thicknesses. Indeed, multiple velocity models can result compatible with the H/V curves by modifying values of both these parameters by the same factor. Therefore, there is the need of independent constraints for one of the two unknowns. For this purpose, we tested the combined use of electrical resistivity tomography (aimed at constraining at least the thickness of the most surficial layers) with ambient noise analysis, which has the capacity of providing information on sub-soil properties down to larger depths.

The application of this approach to the case of the Yang Jia Gou rock avalanche triggered by the 2008 Wenchuan earthquake shows that the combination of ERT and ambient noise methods can provide consistent information on geometrical and mechanical properties of the rock avalanche deposit forming the landslide dam. This information is relevant for the assessment of the debris flow hazards linked to the ongoing erosion of the landslide dam by the Yangjia gully stream. In particular, the combination of geophysical data revealed:

  1. 1.

    The presence of an overlap of layers having different proportions of solid (rock blocks) and fragmented material, as well as different degree of compaction, possibly derived from successive slope failure events or from successive deposition stages of the same rock avalanche, along with an estimate of their thickness

  2. 2.

    Lateral changes in thickness and mechanical properties of the most surficial layer characterized by the presence of more compact and resistive material, corresponding to rock blocks englobed within an incoherent matrix composed of the fragmented and sheared material; the thickness of the surficial layer decreases towards the distal portion of the landslide dam deposit; this surface layer is the main contributor of loose material to recurrent debris flows originating during the monsoon periods

  3. 3.

    Lateral variation of the overall thickness of the landslide dam deposit, possibly reflecting the substratum undulations related to the presence of the deposit of an older rock-avalanche deposit.

Although based on a limited number of measurement stations, the results indicate the applicability of ambient noise measurements acquired by three-component single-station instrumentation to infer properties of thick deposits with rough topographic surfaces. Such investigation targets pose difficulties to the use of geophysical survey techniques requiring the deployment of long linear array of sensors. The successful exploitation of ambient noise analysis methods, however, necessitates the application of advanced techniques of data processing and analysis. As a future perspective, further advancements are expected by combining different techniques of ambient noise analysis, relying on synchronized simultaneous recordings of the noise wavefield at different measurement points. This is the topic of our ongoing research.