# Combining interferometry and wave equation tomography for near surface characterization: 3-D imaging of the Harmalière alpine landslide

Giuseppe Provenzano, Grégory Bièvre, Romain Brossier, Stéphane

Garambois, Ludovic Métivier

<figure><img src="image_1.png" /></figure>

## To cite this version:

Giuseppe Provenzano, Grégory Bièvre, Romain Brossier, Stéphane Garambois, Ludovic Métivier. Combining interferometry and wave equation tomography for near surface characterization: 3-D imaging of the Harmalière alpine landslide. Geophysical Journal International, 2025, 243 (2), pp.ggaf354. 10.1093/gji/ggaf354. hal-05336544

HAL Id: hal-05336544

https://hal.science/hal-05336544v1

Submitted on 29 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Combining interferometry and wave equation tomography for near surface characterization: 3-D imaging of the Harmalière alpine landslide

Giuseppe Provenzano<sup>1</sup>, Grégory Bièvre,<sup>1</sup> Romain Brossier,<sup>1</sup> Stéphane Garambois<sup>1</sup> and
Ludovic Métivier<sup>1,2</sup>

Accepted 2025 September 3. Received 2025 September 2; in original form 2025 May 19

## SUMMARY

Harmalière is an active landslide located in a mountainous region of southern France where the presence of a thick layer of clays provides favourable conditions for the development of slowly moving landslides. However, at Harmalière, the alternation of sudden reactivation and quiet episodes suggests that specific structural and geomechanical properties control its kinematics. In order to shed light on its subsurface properties, we deployed for one month a dense network of seismic nodes within the landslide and recorded active-source and ambient noise seismic data. These data sets have been first independently processed with dedicated interferometry-based processing and inversion workflows to reconstruct *P*-wave (active sources) and *S*-wave (seismic noise) velocity models. Full wave-equation tomography is then performed to improve the reliability and resolution of the obtained elastic model by iteratively fitting the virtual gathers obtained by cross-correlation of the ambient-noise recordings. As opposed to conventional ambient-noise tomography, the approach fully accounts for topography and 3-D elastic heterogeneities. The obtained high-resolution 3-D models are then qualitatively interpreted in terms of landslide properties and geological lithologies, that can influence landslide kinematics.

**Key words:** Interferometry; Geomechanics; Seismic noise; Seismic tomography; Waveform inversion.

## 1 INTRODUCTION

Landslides pose a risk to settlements and infrastructures, causing loss of lives and significant economic damages (Mansour et al. 2011; Froude & Petley 2018). Initiation of mass-movement, its volume and dynamics depend on a combination of external forcing factors (e.g. precipitations) and inherent pre-conditioning properties of the shallow earth (10–100 m depth), namely lithology, permeability, distribution of effective stress and response to loading and shear (Jongmans et al. 2009; Renalier et al. 2012; L’Heureux et al. 2012; Cappa et al. 2014; Lacroix et al. 2020). Spatial and temporal variations of these properties are conducive to changes in geophysical observables (Whiteley et al. 2019) such as electrical resistivity, density and seismic wave velocity. In particular, seismic waves are sensitive to anomalies of shear strength and porosity that may be key to remotely identify the pre-conditioning factors and the evolution of an unstable slope (e.g. Jongmans et al. 2009; Renalier et al. 2010; Vardy et al. 2017), with a resolution in three-dimensions that depends on the acquisition design and frequency content of the energy

sources (e.g. Huang & Schuster 2014). Complementary to coring campaigns and direct measurements on sediment samples, appropriately designed geophysical imaging is a potential cost-effective non-invasive characterization tool to fill the landslide knowledge gap in three dimensions, thus providing more robust and comprehensive hazard estimates.

The Avignonet (AVI) and Harmalière (HAR) landslides, in the French Western Alps around 30 km south of Grenoble and east of the village of Sinard, are slow-moving landslides affecting areas of 2 and 3 km<fnref n="2" />, respectively (Fig. 1a). Both are characterized by a landslide-prone fine-grained Quaternary unit (the Trièves Clayes) whose thickness ranges between 0 and 250 m. This unit is made up by alternating glacio-lacustrine laminated silt and clay, and non-laminated lacustrine sediments mixed with moraines (Fig. 1c). It overlays an irregularly shaped bedrock made of Jurassic marly limestones, and it is capped by moraines. Detailed descriptions of the geological, geotechnical and geophysical settings are provided in Bièvre & Crouzet (2021), Fiolleau et al. (2021) and Jongmans et al. (2009). Despite the same geological and meteorological

<footnote n="1"><sup>1</sup>Univ. Grenoble Alpes, ISTERe, F-3800 Grenoble, France. Email: ludovic.metivier@univ-grenoble-alpes.fr</footnote>

<footnote n="2"><sup>2</sup>CNRS, LJK, F-3800 Grenoble, France</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Location of the study site and of the experiment. (a) Localization of the Harmalière (HAR) and Avignonet (AVI) landslides, and of the seismic experiments. The blue dots represent the geophone array and the co-located active sources, the red dots shows the location of the circular array used by Fiolleau *et al.* (2021); (b) zoom on the HAR landslide with mapped surface deformation and suggested transition between plastic and fluid regime; (c) Geological cross-section along the central, active part of the HAR landslide, representing the current state of knowledge of the survey area.</figcaption></figure>

setting, these two landslides exhibit different deformation mechanisms, suggesting the role of local subsurface heterogeneities. While AVI is gently sliding eastward with motions ranging between a few to some tens of cm yr⁻¹, HAR shows sudden reactivations, the most important being in March 1981 and June 2016. During this last event, Lacroix *et al.* (2018) measured motions up to around 40 m within a month. In addition, and differently from AVI, the central part of HAR varies between a plastic slide at the top and a fluid flow at the toe (Figs 1b,c), although the precise location and geometry of the transition have not been yet identified (e.g. Fiolleau *et al.* 2020).

Bièvre *et al.* (2011) assessed the pseudo-3-D lacustrine body thickness by forward-computing the horizontal-to-vertical spectral ratio (HVSR) of ambient noise recordings at sparse locations. They found that the different landsliding styles may be explained by local bedrock paleotopography, with an important bulge at the toe of AVI that could act as a buttress and prevent reactivations as observed in HAR. A similar approach, complemented by a local 1-D S-wave profile obtained by array ambient noise analysis (spatial autocorrelation, SPAC), within HAR was adopted by Fiolleau *et al.* (2021). The sparseness of the measures, along with the simplifying assumptions of the methods (vertically varying stratified medium, local 1-D parametrization), however, may be limiting factors in

capturing key subsurface complexities, such as the position and geometry of the transition between slide and flow within HAR.

Within the RESOLVE project (Gimbert *et al.* 2021), a dense 3-D array of broad-band geophone-nodes has been deployed in May-June 2021 to record continuously seismic noise for a one-month period in the Harmalière site (Voisin *et al.* 2021). In addition, active-source data have been recorded by the same array. In this study, we use both the active-source and ambient noise densely sampled vertical-component wavefields, exploiting their complementary characteristics in frequency content as well as in dominant seismic phases, to infer the 3-D *P*- and *S*-wave velocity distributions within the landslide body and below. To this purpose we follow the flowchart presented in Fig. 2.

The active data set is pre-processed by super-virtual refraction interferometry (SVRI; Bharadwaj *et al.* 2012) to improve the signal-to-noise ratio (SNR) of *P*-wave refractions and thus better constrain the 3-D *P*-wave model obtained by first-arrival traveltime tomography. A pseudo-3-D *S*-wave velocity volume is derived by inverting the fundamental mode of Rayleigh wave group velocity dispersion curves emerging from the cross-correlation of the ambient noise recordings (Shapiro & Campillo 2004). In order to overcome the assumptions of the conventional surface-wave inversion
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Summary of data processing and inversion steps for the 3-D elastic velocity model building.</figcaption></figure>

approach by fully accounting for 3-D propagation and topographic
effects, we adopt a full waveform inversion approach (Virieux et al.
2017) based on a spectral-element modelling (SEM46; Trinh et al.
2019; Cao et al. 2022) and a windowed traveltime-based objective
function (multitaper misfit function; Tape et al. 2010). The wave-
equation based tomography (WET) uses the interferometry-based
elastic macro-model as initial guess, and updates jointly P- and S-
wave velocities by iteratively fitting the cross-correlation waveform
traveltimes. This approach is comparable to the crustal (≈ 100 Km)
studies of Lu et al. (2020) and Nouibat et al. (2023), at a much
smaller (≈ 100 m) investigation scale.

The resulting 3-D *P*-wave velocity and *S*-wave velocity models
serve as a basis for a preliminary geomechanical interpretation of
the Harmalière landslide. A bulge-like structure is highlighted in
the central part of the landslide, which could explain the specific
reactivation pattern of the Harmalière landslide compared with the
neighbouring Avignonet landslide. After a reactivation at the head-
scarp, the eroded material accumulates uphill behind this bulge-like
structure, before a specific stress threshold is reached, and the land-
slide transitions from an elasto-plastic slide to a visco-plastic flow
downhill.

The paper is organized as follows. In Section 2, we outline the data
characteristics and pre-processing. In Section 3, the interferometry-
based methodologies used to obtain the elastic ($V_p - V_s$) macro-
model are presented. In Section 4, the wave equation tomogra-
phy method is introduced, and its application to the ambient noise
recordings starting from the interferometry-based model is pre-
sented. In Section 5, a preliminary geomechanical interpretation of
the landslide is introduced, based on the reconstructed 3-D $V_p$ and
$V_s$ models and the $V_p/V_s$ ratio. A discussion is then proposed in Sec-
tion 6, before conclusions and perspectives are given in Section 7.
A data availability statement is added in the next section.

## 2 DATA

An array of 100 three-component 5-Hz geophones (Magseis Fair-field ZLand 3C) was deployed in the target area according to a grid with a nominal 70 m spacing (Fig. 1a) to record continuously natural seismicity for a 33-d period (Voisin et al. 2021). This was complemented by the acquisition of an active data set using 100 vertical hammer-strike sources co-located with the receivers, with
nominal offsets ranging between 70 and 1000 m. In this section, the data characteristics and the dedicated processing steps are outlined. The passive data is publicly available (see data availability statement).

## 2.1 Active source processing by super-virtual interferometry

In the active-source data set, a total of 94 usable sources and 99 active receivers have been selected based on their noise energy. The three hammer-strikes performed per location (Fig. 3a) are processed before stacking to maximize the improvement in signal-to-noise (SNR) ratio: (1) a data-driven identification of the trigger times compensates for the lack of source-receiver synchronization, by means of a constrained auto-picker; (2) a multichannel shaping deconvolution makes up for the non-repeatable source signature of the hammer-strikes, by matching the waveform shape of each channel to a reference source wavelet; (3) common-source gathers are stacked by scaling each trace by the inverse of its average power, in order to attenuate the impact of amplitude variations.

After stacking, a frequency decomposition of an example common-source gather is instructive of the spectral diversity of the data set (Figs 3b,c,d). Most of the recorded energy sits at frequencies between 20 and 40 Hz (Fig. 3c), and is carried by dispersive slow arrivals. A large proportion of the source energy in land seismic surveys is indeed expected to be converted to multi-scattered guided waves and surface waves in the very shallow, slow and heterogeneous subsurface (Pugin & Yilmaz 2019; Stork 2020; Yilmaz et al. 2022). Little or no coherent active-source energy is instead observable below 10 Hz (Fig. 3b), due to the characteristics of the vertical hammer-strike source. Identification of body waves refracted arrivals is possible between 40 and 90 Hz (Fig. 3d), but only up to ≈400 m offset, beyond which the SNR deteriorates below an acceptable level for first-arrival traveltime picking.

In order to extend the usable aperture of the P-wave refraction data set, super-virtual refraction interferometry is applied (SVRI; Bharadwaj et al. 2012; Hanafy & Al-Hagan 2012; Place & Malehmir 2016). SVRI exploits the redundant sampling of refracting paths in multichannel surveys to enhance the SNR of critical refraction. This is performed in two steps:

(i) Datuming: each pair of common-receiver trace gather ($R(x_a), R(x_b)$) is cross-correlated, and the contributions from all the stationary sources ($s$) are stacked to obtain a virtual refraction transfer function (Green's function) between $x_a$ and $x_b$. In this work, we use a normalized cross-correlation function (cross-coherence) in the frequency domain

$$G(x_b, x_a, \omega) = \sum_s \frac{R^*(s, x_a, \omega) R(s, x_b, \omega)}{\|R(s, x_a, \omega)\| \|R(s, x_b, \omega)\|}, \quad (1)$$

where $\omega$ is the angular frequency and the operator * indicates complex conjugation. The obtained Green's function $G(x_b, x_a)$, in the time domain, has its maximum at the relative refraction delay time between the two stations.

(ii) De-datuming (or modelling): for each common-source trace gather, the super-virtual trace at position $x_b$ is obtained by stacking the convolutions between each trace recorded at position $x$ and the virtual transfer function $G$ between $x$ and $x_b$ obtained in (eq. 1).

$$R_{SVRI}(s, x_b, \omega) = \sum_x R(s, x, \omega) G(x_b, x, \omega), \quad (2)$$

After inverse Fourier transformation, the super-virtual trace $R_{SVRI}(s, x_b, t)$ contains the enhanced refraction at the absolute source-to-receiver refraction traveltime.

Prior to SVRI, the 3-D data set is windowed in a time-offset region where refractions can be plausibly recorded, based on a prior knowledge of the expected P-wave model. In order to ensure SVRI is applied to traces for which refractions share portions of the same stationary phase path (Wapenaar et al. 2010a), eqs (1) and (2) are applied iteratively to narrow corridors of 'in-line' sources and receiver gathers, until the full azimuth range is covered. This yields a redundant coverage whereby the same station pairs may fall into more than one corridor, and thus undergoes cascading cross-correlations and convolutions. Cross-coherence is used instead of cross-correlation in order to reduce the pulse broadening resulting from the iterative cross-correlations in the datuming step (Lu & Chavez-Perez 2020), while phase-weighted stacking is used in both steps to enhance noise attenuation (Schimmel et al. 2011).

Traveltime picks are automatically determined by recursive STA/LTA on both the SVRI and the raw data set. Non-reciprocal times (i.e. with discrepancies larger than 9 per cent) and trigger amplitudes lower than a threshold are then discarded. Based on these criteria, the number of valid picks is 50 per cent larger in the SVRI data set than in the raw one, being in particular richer in offsets between 400 and 700 m. In Fig. 4 we present the picks performed on one active shot-gather, before and after the SVRI processing. These picks appear as blue dots, the rejected ones circled in red, the accepted ones circled in green. The picks are superimposed to the window envelope of each trace of the gathers, again, before and after SVRI. The enhancement of the pick reliability thanks to SVRI is striking.

## 2.2 Ambient seismic noise

The lack of coherent energy below 10 Hz limits the penetration of the active-source surface waves to the top 5–10 m for realistic shallow subsurface velocities (e.g. Zhou et al. 2011). The ambient noise recordings, on the other hand, have a broad frequency band limited only by the instrument transfer function, whose corner frequency is 5 Hz with an acceptable response down to 1 Hz. These data can therefore be exploited to fill the low-frequency active-source spectral gap by extracting via interferometry the coherent surface waves contained in the noise field.

The ambient noise recordings have thus been used to extract Empirical Green's functions (EGF) by cross-correlating independent pair of stations, under the assumption of a stationary and uniform ambient noise field (Shapiro & Campillo 2004; Campillo 2006).

Data pre-processing follows a rather conventional workflow (Bensen et al. 2007). The noise traces are decimated to a sampling rate of 15 Hz, in order to compress its size while respecting Nyquist-Shannon criterium for a maximum frequency of 7.5 Hz. Within this band, a comb filter is applied. In three frequency sub-bands, namely 1–3, 3–5 and 5–7 Hz, the time-series are normalized by the value of the envelope (module of the Hilbert transform), thus equalizing the contributions of weak and strong signal and attenuating the impact of short high-amplitude transients (e.g. micro-earthquakes, slide-quakes).

Each 24-hr long segment is divided into 15 min-long windows, on which cross-coherence is calculated in a range of time-lags between (–7,7) s. Cross-coherence is performed in the frequency domain,
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) Raw active-source example gather. Note the three strikes at quasi-regular intervals between 20 and 50s listening time, before stacking. Spectral analysis after stacking (b–d): gathers below 10 Hz (b), between 10 and 40 Hz (c) and above 40 Hz (d).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Example processed and windowed envelope of an active-source recording before (a) and after (b) SVRI and first-arrival picking for *V*<sub>p</sub> traveltime tomography. Data is shown in greyscale, white corresponding to high amplitudes. Acceptable picks in green, rejected ones in red. Time-axis static-shifted by 0.2 s to help auto-picking.</figcaption></figure>

where a power-spectrum normalization is added with respect to con-
ventional cross-correlation, thus effectively implementing for each
data chunk a spectral whitening. Finally, the 14-s long (–7,7) cross-
coherence functions are phase-weighted stacked in order to enhance
the contribution of coherent waveforms (Schimmel *et al.* 2011), and

a virtual-source gather is obtained for each station (Fig. 5). For sim-
plicity, hereafter these are referred to as cross-correlation traces.

Note that the ambient noise EGFs derivation by cross-correlation
is equivalent to the first step of SVRI (eq. 1) where, instead of a
known in-line source position, if the noise is uniformly distributed,
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Example cross-coherence noise function for negative and positive time-lags, sorted in increasing offset with respect to the reference station. Dispersive surface waves emerge from the ambient noise interferometry with higher wavespeeds at lower frequencies. 1.5–2.5 Hz frequency band (top), 3–5 Hz frequency band (bottom).</figcaption></figure>

only the stationary phases contribute constructively to reconstructing the Green's function between two stations (Wapenaar et al. 2010a).

A frequency-time analysis (FTAN, e.g. Bensen et al. 2007) is applied to the vertical-to-vertical (ZZ) cross-coherence functions to obtain the time-delay of the surface waves envelope between each station pair. If the cross-correlation functions are reliable estimates of the EGF, those can be interpreted as group slownesses once divided by the interstation distances. Since no significant-energy overtones have been observed, they are attributed to the Rayleigh wave fundamental mode.

However, the anisotropy of the ambient noise field may bias the group velocity estimates with spurious apparent velocities corresponding to dominant arrivals outside the stationary phase zone

(Wapenaar et al. 2010a). Although dispersion can be qualitatively observed (Fig. 5), the example EGF shows a strong asymmetry between positive and negative time-lags in particular at higher frequencies, which is indicative of an anisotropic noise sources distribution. This is probably due to the impact of anisotropic noise sources in the nearby settlements and transit of motor vehicles.

A weak symmetry criterion, accounting only for the similarity between positive and negative time-lags is here used as a measure of the reliability of the kinematic estimates. For each frequency, the positive and negative time-lags are picked independently via a constrained automatic procedure. If the discrepancy between the positive and negative apparent group velocities is larger than 9 per cent of their average, the value is discarded; otherwise, the average
is kept as a reliable measure for that frequency and distance. A similar approach is used in Nouibat et al. (2022).

A final statistical analysis is performed to discard velocity measurements that appear as outliers of an estimated Gaussian distribution, that is, whose discrepancy with the mean is larger than two-times the standard deviation. The selected velocity data points as well as the corresponding full waveform cross-correlation traces, are going to be used as data sets for the following inversion steps.

# 3 BUILDING THE ELASTIC MACRO-MODEL

## 3.1 P-wave traveltime tomography

3-D P-wave first arrival traveltime tomography is performed using a 3-D ray-based linearized inversion approach (TOMO-TV; Latorre et al. 2004). The starting model is determined by grid search over an ensemble of 1-D profiles within a range of sensible velocities for the expected sediment types (Fig. 1). A relatively strong regularization parameter ($\alpha = 4$) for the LSQR inversion is obtained by L-curve analysis. Anisotropic Laplacian smoothing is used, with correlation lengths ten times larger in the x and y horizontal directions than in the vertical direction z, to reflect an expected stronger vertical variability. After 50 iterations, the inversion converges attaining an overall final misfit decrease of 60 per cent.

## 3.2 S-wave ambient noise tomography

Unlike 3-D body-wave tomography, the conventional surface-wave inversion workflow produces pseudo-3-D results in a two-step workflow (e.g. Shapiro et al. 2005a; Nouibat et al. 2022), referred to as ambient noise tomography (ANT): (1) 2-D frequency-dependent velocity maps inversion outputs local velocity-frequency relationships for each grid-point (dispersion curves); (2) 1-D depth inversion of the dispersion curves obtained in step 1 yields a wavelength-to-depth mapping for each grid point.

The local dispersion curves obtained for each independent station pair are regionalized by performing a 2-D traveltime tomography on the time-delay data set (e.g. Shapiro & Campillo 2004; Shapiro et al. 2005b; Nouibat et al. 2022) in a homogeneous starting model (300 m s⁻¹). A nonlinear tomographic approach (TOMO-TV, e.g. in Latorre et al. 2004) with ray tracing at each iteration is used as opposed to a more conventional straight ray approximation (e.g. Karimpour et al. 2022). Thirty iterations of LSQR (Paige & Saunders 1982) inversion are performed for each frequency band, resulting in a densely sampled dispersion data set in a 10 m grid volume.

The uneven ray coverage resulting from the cross-correlation data selection is apparent in particular in the southernmost part of the target and at higher frequencies (Fig. 6), where noisier data are observed, probably due to poorer station coupling. A derivative-weighted-sum approach is used to attribute a standard deviation value to each velocity-frequency pair, inversely proportional to the ray density. A poor ray coverage corresponds to a large uncertainty ($\sigma = \pm 50$ per cent of the estimated slowness), whereas the maximum ray coverage corresponds to a high-confidence slowness value.

For each grid point of the frequency-dependent tomographic maps (Fig. 6), the group velocity dispersion curve, along with the associated standard-deviation, is interpreted as that of the fundamental mode of Rayleigh waves and inverted for a local 4-layers S-wave depth model using a Neighbourhood algorithm (NA; Cambridge

1999) in the parallel version implemented in the Geopsy-Dinver software package (Rickwood & Cambridge 2006; Wathelet 2008).

Prior to inversion, each dispersion curve is cleaned from outliers by a 3-point median filter and resampled with 50 values between 1 and 6 Hz. For each layer, the unknowns are the S-wave velocity and the thickness, thus resulting in a highly nonlinear inverse problem with 8 degrees of freedom. The 4-layers parametrization has been chosen heuristically by assessing the inversion quality on the average dispersion curve. As shown by Sajeva et al. (2017), NA is a well suited approach to locate the global minimum of non-convex objective functions with this dimensionality.

NA first generates 120 independent models within a uniform a-priori distribution; then, a total of 6000 models are generated by iteratively resampling the parameter space in the Voronoi cells of the best 80 models, thus progressively focusing the exploration of the objective function around the optimal solution. Finally the best model is drawn from the distribution for each location (Fig. 7), shifted to the local elevation and used to build the best-fitting pseudo-3-D S-wave model (Fig. 8a).

## 3.3 Resulting P-wave and S-wave velocity models

In Fig. 8, the ANT $V_s$ (a) and tomographic $V_p$ (b) volumes are presented. It is worth pointing out the consistency between the shallow low velocity structures in the southern part of the NS slice, even if the two models have been obtained independently; on the other hand, dissimilar trends may be the impact of 3-D changes in $V_p/V_s$ ratios in the shallow subsurface.

# 4 WAVE-EQUATION TOMOGRAPHY

## 4.1 Waveform modelling and inversion

The obtained pseudo-3D ANT $V_s$ model suffers from two main limitations: (1) it does not account for the impact of topography on surface wave propagation (Snieder 1986; Li et al. 2019), and (2) each inverted depth profile relies on a local 1-D approximation, thus neglecting the impact of the important lateral heterogeneities in the shallow subsurface (e.g. Schäfer et al. 2013).

In order to overcome these shortcomings, 3-D $V_p - V_s$ wave-equation tomography (WET) is applied on the ambient noise data, using the ANT $V_s$ model and the SVRI tomographic $V_p$ model as initial guess. WET iteratively updates the elastic subsurface model (m) by minimizing the misfit (F) between the calculated (or synthetic, $d_{cal}$) and observed ($d_{obs}$) waveforms. The latter are obtained as the convolution of the ZZ particle-velocity cross-correlation functions, considered to be an approximation of the Green's functions between vertical-force sources and vertical component receivers (Campillo & Paul 2003; Wapenaar et al. 2010b), with a source time function (source wavelet, w) taken as a band-pass filtered version of a Dirac impulse function in the frequency-band of the cross-correlation data. For a pair of stations (a, b), we have

$$ \begin{cases} d_{\text{obs}}(x_a, x_b, t) = G_{\text{obs}}(x_a, x_b, t) \star w(t), \\ d_{\text{cal}}(x_a, x_b, t)[m] = G_{\text{cal}}(x_a, x_b, t)[m] \star w(t), \end{cases} \quad (3) $$

where $G_{\text{obs}}$ are the estimated Green's functions from ambient noise interferometry, x indicates the location, t the time dependence and ★ the convolution operator, while $G_{\text{cal}}$ are numerical Green's function computed using the spectral element full waveform modelling and inversion code SEM46 (Trinh et al. 2019; Cao et al. 2022).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Regionalization of Δt(f) by traveltime tomography gives 2-D Group velocity (GV) maps [$V_g(f, x, y)$]. Examples at 1.8 and 3.2 Hz of group velocity values (left hand side) and corresponding ray coverage (right hand side). In the ray maps, blue triangles and red dots represent respectively sources and receivers positions.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Example of Rayleigh group velocity dispersion curve inversion. Profile extracted from the central part of the survey area. Observed dispersion data in black dots. Misfit values are normalized to standard deviation of the slowness data. (Fig. 6).</figcaption></figure>

WET solves the following minimization problem on the misfit function F for a chosen parametrization of the subsurface model m

$$ \min_{m} \sum_{x_a} \sum_{x_b} F(d_{\text{cal}}(x_a, x_b)[m], d_{\text{obs}}(x_a, x_b)). \quad (4) $$

The objective function F is the multitaper traveltime (MTT) distance introduced by Tape et al. (2010) as the sum of the squares of the frequency-dependent phase-traveltime residuals. As opposed to a time-domain cross-correlation objective function (Luo & Schuster 1991), MTT accounts for frequency-dependent traveltime residuals, extracted from robust multitaper spectral estimates of observed and calculated data.

To compute the MTT objective function, following Tape et al. (2010), we first calculate the multitapered least-squares transfer function T between dcal and obs using N tapers

$$ T(\omega) = \frac{\sum_{k=1}^{N} \hat{d}_{\text{obs}}^{k}(\omega) \hat{d}_{\text{cal}}^{k}(\omega)^{*}}{\sum_{k=1}^{N} \hat{d}_{\text{cal}}^{k}(\omega) \hat{d}_{\text{cal}}^{k}(\omega)^{*}}, \quad (5) $$

where d̂_obs and d̂_cal are respectively the frequency-domain observed and calculated data multiplied by the kfh taper ξ

$$ \left\{ \begin{array}{l} \hat{d}_{\text{obs}}^{k}(\omega) = \mathcal{F}(\xi_k(t)d_{\text{obs}}(t)), \\ \hat{d}_{\text{cal}}^{k}(\omega) = \mathcal{F}(\xi_k(t)d_{\text{cal}}(t)), \end{array} \right. \quad (6) $$

with F the Fourier transform operator.

In this work, N = 3 orthogonal tapers (ξ) are used. The summation over orthogonal tapers aims at obtaining a robust spectral

estimate of a signal, as the mean of (simulated) independent realizations, namely one for each taper. The phase differences are then computed as the inverse tangent of the ratio between the imaginary and real part of the multitapered transfer function

$$ \Delta\phi(\omega) = \tan^{-1}\left(\frac{\Re(T(\omega))}{\Im(T(\omega))}\right). \quad (7) $$

From Δφ, the frequency-dependent traveltime Δt(ω) is defined as

$$ \Delta t(\omega) = -\frac{\Delta\phi(\omega)}{\omega}. \quad (8) $$

The MTT objective function comparing observed dobs and calculated data dcal is finally defined as

$$ F(d_{\text{cal}}[m], d_{\text{obs}}) = \sum_{\omega} \Delta t^2(\omega). \quad (9) $$

The spectral element modelling code we use to compute the Green's function (Gcal) uses a structured mesh with vertically deformed hexahedral elements in order to account for topography. The latter is extracted and resampled from a high resolution (1 m) LiDaR (Light Detection and Ranging) airborne survey. An adaptive mesh is built with element size varying as a function of the local wavelength (λ) in the starting model. Since wavelengths in the order of few hundreds metres are propagated at the lower end of the spectrum (≈ 1.5 Hz), sponge layers with a width comparable to the model size are introduced to reduce the impact of off-boundaries reflections everywhere except at the free-surface.

At each iteration, the gradient of the objective function is computed via the adjoint state method (Plessix 2006) as the zero-lag cross-convolution of the incident and adjoint fields. Similarly to the cross-correlation time-lag objective function (Luo & Schuster
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Pseudo-3-D S-wave ANT volume (a) and 3-D $V_p$ SVRI tomography model (b).</figcaption></figure>

1991), but in a frequency-dependent approach, the adjoint source of MTT depends only on the first-order derivatives of the calculated (multitapered) data, weighted by the phase-traveltime residuals.

We update $V_s$ and $V_p$ simultaneously, to account for the non-zero impact of $V_p$ on Rayleigh waves in this relatively high frequency band. The multiparameter inversion uses a quasi-Newton (l-BFGS) optimization (Nocedal & Wright 2006). Combined with a source-wavefield illumination compensation pre-conditioner, l-BFGS aims at approximating the effect of the inverse Hessian on the model update direction (Métivier & Brossier 2016). Regularization is implemented by means of a double Bessel filter applied to the gradient to approximate a Laplace smoothing (Trinh et al. 2017), with coherent-lengths equal to 0.1λ.

## 4.2 Data pre-conditioning and preliminary analysis

As in the regional seismology application of Nouibat et al. (2023), the objective function is computed in a windowed portion of the traces centred around the peak of the synthetic data amplitude, with length equal to the wavelet duration ± the maximum expected time-lag. In this work, windows with a poor waveform similarity, namely with cross-correlation lags and value respectively above and below a threshold or large amplitude discrepancies are rejected (as in, e.g. *FlexWin*, Maggi et al. 2009).

An offset-dependent weight-function is designed to penalize the contribution of data outliers, for example, waveforms outside a realistic slowness range in the frequency band of interest, based on the characteristic of the average phase dispersion curve.
In defining the maximum usable frequency for WET, the spatial sampling of the wavefield ought to be considered. With a nominal node spacing of 70 m, wavelengths shorter than 140 m would be spatially aliased according to the Nyquist-Shannon criterium. Although Brenders & Pratt (2007) show that full waveform tomography may tolerate cross-line sampling intervals coarser than Nyquist's, Nuber et al. (2017) points out that sampling should not be relaxed beyond one minimum shear wavelength. Considering a minimum $V_s$ in the order of 200 m s⁻¹, we low-pass filter the cross-correlation waveforms below 3.5 Hz, thus limiting the minimum wavelength in the order of the node spacing and avoiding major aliasing artefacts.

Pre-processing consists of

(i) selection of virtual gathers with at least 10 valid traces (as per the symmetry criteria used in ANT within the band of interest);

(ii) median filter to remove zero-slowness phases;

(iii) convolution with a band-limited source wavelet obtained by minimum-phase filtering a Kronecker delta with corner frequencies (−3 dB) equal to 1.3 and 3.0 Hz.

In Fig. 9, the $V_s$ depth sensitivities are shown for the fundamental mode Rayleigh wave phase velocity at the two corner frequencies of the inverted band and its peak frequency. Despite being computed on an average 1-D starting $V_p - V_s$ model, these functions contain useful insights about the expected investigation depths. As well as constraining the top 30 m of the subsurface, model updates down to 100 m depths should be expected in this band. The $V_p$ sensitivity kernel (not shown here) has similar shape, but is two orders of magnitude smaller.

## 4.3 Results

WET converges after nine iterations to the $V_s$ model shown in Fig. 10, reducing the MT-traveltime misfit to around 45 per cent of the initial one. This corresponds to an increased waveform alignment for most traces at the measurement windows, although the relatively high misfit plateau reflects the non-reproducibility of the full cross-correlation waveforms at some locations (Fig. 11).

The final $V_s$ model shows significant broad-band updates with respect to the starting ANT model both in the shallow and deep subsurface, with higher values at depths between 40 and 70 m, as well as a sharper low-velocity anomaly in the southernmost part of the model (Figs 12 and 13). On the other hand, the update in $V_p$ (not shown here) is minor, due to the small sensitivity of Rayleigh wave kinematics to P-wave velocity.

The depth slices shown in Figs 13(e)–(h) show that the WET $V_s$ model contains both long-wavelength and high-resolution updates with respect to the initial ANT model (Figs 13a–d). In addition, the near-vertical artefacts of the pseudo-3-D starting model are attenuated (Figs 12b,d versus Figs 12a,c), thanks to the true-3-D modelling and the appropriate representation of the slope topography.

# 5 GEOMECHANICAL INTERPRETATION

The seismic imaging methods used in this study yield an estimate of the 3-D distributions of two seismic parameters, $V_p$ and $V_s$, from which $V_p/V_s$ ratio is obtained a-posteriori. In this section, we jointly interpret them to provide the preliminary geomechanical interpretation displayed in Fig. 14.

The estimated 3-D WET $V_p$ model shows values between 800 and 1300 m s⁻¹ in the top 10 m depth below the ground surface up-hill (between 600 and 800 m along cross-section, Fig. 14a), which

suggests the presence of unsaturated clays. The relatively coarse sensors spacing prevented the picking of direct body-wave travel-times, hence limiting the resolution of the shallowmost $V_p$ model (Figs 14a, b). In the central part of the model, between 250–300 and 400–600 m, the estimated $V_p$ remains around 1300 m s⁻¹ down to 50 m depth, Fig. 14a), which might be indicative of an extension of the unsaturated zone to larger depths. This is deeper than what is observed at the rear of the Harmalière headscarp and in the Avignonet landslide, where the water table is generally found a few metres below the surface (Jongmans et al. 2009; Bièvre et al. 2016). The observation of this apparent deep water table remains to be confirmed with a higher resolution model. Below, P-wave velocity values increase quickly to 1500–1800 m s⁻¹, possibly marking the transition to a saturated clayey horizon. P-wave velocity values reaching 2000 m s⁻¹ are observed underneath this saturated horizon, in the lower half part of the model. The upper bound of 2000 m s⁻¹ for P-wave velocity imposed in the model as a bound constraint in the inversion process (to limit the computing cost of the spectral element modelling) prevents from clearly interpreting these layers. They can be associated with either firm, saturated clays, fluvial terraces made of locally cemented pebbles to boulders or the carbonate substratum, the latter two being both visible at outcrops and reported in the literature (Monjuvent 1973; Bièvre et al. 2011).

Cross-sections of the $V_s$ WET model (Figs 14c, d) display a low velocity layer in the near-surface (lower than 200 m s⁻¹), whose thickness decreases from around 30 m at the top to 10 m at the toe of the slope (Fig. 14c). At the top, the low-to-high velocity transition depth is in good agreement with a proximal 1-D $V_s$ profile obtained 75 m up-hill of the section (Fiolleau et al. 2021) (Fig. 14c), although the $V_s$ values in the WET 3-D model are lower.

The corresponding interface can be interpreted as the base of the landslide, in analogy with the nearby Avignonet landslide, where similar 1-D and pseudo 3-D ambient noise tomography $V_s$ models (Jongmans et al. 2009; Renalier et al. 2010) are in good agreement with results from borehole (Bièvre & Crouzet 2021) and inclinometer data (Jongmans et al. 2009).

Due to the upper frequency limit of 6 Hz used in this work, the $V_s$ variability is not well constrained in the landslide body itself. Thus, it is not possible to detect an expected $V_s$ decrease between the slide (uphill) and the flow (downhill) with the necessary resolution. However, below this horizon, interpretable heterogeneities can be observed and interpreted. Uphill, from 500 m to the end of the longitudinal cross-section, $V_s$ values ranging between 450 and 500 m s⁻¹ may correspond to non-disrupted clays. The velocity of this layer is slightly lower than the 580 m s⁻¹ value reported by Fiolleau et al. (2021) (Fig. 14c). This could originate from the different setups between the experiments (1-D versus 3-D, single-mode versus multimodal inversion) and/or non-negligible lateral variations between the two investigation sites. Between 300 and 500 m, a ‘bump-like’ feature is characterized by higher $V_s$ values (larger than 450 m s⁻¹ and up to 620 m s⁻¹). Downhill, $V_s$ values decrease to values typical of remoulded clays, in the range of 300–350 m s⁻¹. The cross-section in Fig. 14(d) highlights a lateral transition from those relatively low values to the highest estimated $V_s$ values in the model, that can be attributed to fluvial terrace sediments.

$V_p/V_s$ ratio cross-sections are presented in Figs 14(e) and (f). A high variability is observed within the landslide body, with ratios ranging between 5 and 10. This may reflect the high uncertainty of $V_p$ and $V_s$ estimates within this layer. Below the landslide body, the $V_p/V_s$ ratio decreases with depth, from around 5.5–6 to 3.5. It is notable that, in the central part of the longitudinal cross-section (Fig. 14e), no $V_p/V_s$ ratio values greater than 5 are observed,
<figure><img src="image_11.png" /><figcaption>**Figure 9.** $V_s$ depth sensitivity kernel of Rayleigh wave phase velocities at the corner frequencies of the WET band (1.3,3.0 Hz) and its central frequency (2.1 Hz) [Computed using *disba* by Luu (2021)].</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** 3-D Spectral element mesh and vertical slices of the final WET $V_s$ volume.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Example gathers trace-by-trace waveform data fit (red = final, black = observed, blue = starting) at the measurement window of the MTT objective function. Traces are sorted in ascending absolute offset, but offset interval is irregular reflecting the 3-D acquisition.</figcaption></figure>

highlighting a sharp transition from 5 to 10 in the slide to less than 5 at greater depth. On the other hand, the ‘bump-like’ feature observed in the reconstructed $V_s$ model, and interpreted as a fluvial terrace, does not exhibit any $V_p/V_s$ difference with respect to the firm clays at the same elevation uphill. Finally, at depth, in the central part of the transverse cross-section (distance between 200 and 300 m in Fig. 14f), a slight increase in the $V_p/V_s$ ratio, could correspond to the carbonate bedrock.

The presence of the rigid bump between 300 and 500 m along the longitudinal cross-section coincides with the location of the slide-to-flow transition (Figs 14a and b). This structure looks similar to bedrock bumps affecting glacial flows (section 7.2 in Cuffey & Paterson 2010), where lower displacement (resp. higher) is observed uphill (resp. downhill) this morphological feature. The small bump observed here could act as a mechanical buttress controlling the deformation regime of the landslide body. After a reactivation at the headscarp, the eroded material is transported downslope and accumulates immediately uphill the buttress. A possible scenario is that once the mass of accumulated material reaches a specific threshold, depending on its water content, the shear stress exceeds a critical threshold that favour the transition between an elasto-plastic (slide) and a visco-plastic (flow) behaviour (Carrière et al. 2018; Fiolleau et al. 2021).

## 6 DISCUSSION

The starting model has been built relying, both for active-source ($V_p$) and ambient-noise data ($V_s$), on interferometric techniques to extract coherent information from apparently non-interpretable data sets. SVRI on active-source data refractions is key in extending the usable offset range closer to the nominal one (950 m) offset. In practice, this was true up to a distance of 750–800 m, beyond which the SNR of the input data is too low for SVRI to be effective. In 3-D, this should in the future account for frequency-dependent end-fire lobes (Wapenaar et al. 2010a) to select only stationary refraction paths from the full survey. The approach adopted in this paper is more pragmatically the sequential selection of inline portions, that is, corresponding to very narrow lobes, which is reasonable for the high frequencies used here (larger than 40 Hz).

The frequency of the ambient-noise empirical GF reconstruction is limited to 7 Hz, and confidently usable below 6 Hz. Even at lower frequencies, the number of reliable traces, for example, respecting a weak (traveltime) symmetry requirement, is in the order of 5 per cent. This reflects the erratic and anisotropic nature of the noise field at those (relatively) high frequencies, which could be compensated for with longer stacking periods. The lower frequency limit is instead determined by the instrument response and, more importantly, by the maximum available offsets (Bensen et al. 2007), that is, 1000 m in this work, which implies high uncertainties for wavelengths larger than 500 m. Wider apertures may grant deeper investigation depths, thus further constraining the bedrock shear properties, which may be important, for example, for seismic site effects estimation.

Despite the theoretical similarity between wave-equation tomography and full waveform inversion (e.g. FWI; Virieux et al. 2017), the resolution of the former is inherently lower than conventional FWI, primarily as a consequence of the traveltime-based objective function (Luo & Schuster 1991). Furthermore, the theoretical $\lambda/2$ depth resolution of FWI is attained, according to the diffraction tomography theory (Devaney 1984), thanks to vertical incidence phases (Huand & Schuster 2014), as opposed to the horizontally
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Vertical slices of the S-wave velocity models cutting through the middle of the volume (position shown in Fig. 10). (a,c) ANT (starting) model; (b,d) WET (final) model.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** *Vs* depth slices. Depth, below the free-surface, increases from left to right: 14 m (a,e), 24 m (b,f), 32 m (c,g) and 64 m (d,h). Panels (a,b,c,d): ANT pseudo-3-D volume. Panels (e,f,g,h): WET 3-D volume. The white lines correspond to the cross-sections presented in the interpretation Fig. 14. The position of these cross-sections is also shown in the overview Fig. 1(b).</figcaption></figure>

propagating surface waves. This effectively limits the expected resolution of surface wave WET to the size of the first Fresnel zone between stations pairs. However, as demonstrated by Fu *et al.* (2017), sub-wavelength resolution imaging of shallow anomalies can be

obtained by inversion of back-scattered surface waves if the target anomaly is located in the near-field, that is, at a depth smaller than $\lambda/3$. We have not investigated whether or not this is the case in our ambient noise waveform inversion. Higher resolution misfit
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Inline and cross-line views of the WET models for *V*<sub>p</sub> (a,b), *V*<sub>s</sub> (c,d) and for the *V*<sub>p</sub>/*V*<sub>s</sub> ratio (e,f), overlaid with a structural interpretation. Panels (b), (d), (f) are cross-sections transverse to the inline cross-sections (a), (c), (e), respectively. The relative position of these cross-line views one with respect to the other is indicated by panels (a–f) on the sections. Note that the position of these inline and cross-line sections is indicated in white line in the depth slices Fig. 13 and in blue lines in the overview Fig. 1(b).</figcaption></figure>

function accounting for amplitude discrepancies (GSOT; e.g. Métivier *et al.* 2019) and (AWI; Warner & Guasch 2016) and more importantly least-squares FWI, may instead suffer from the inaccurate frequency-dependent amplitudes of the empirical Green's function derived from ambient noise cross-correlations, as well as from the difficulty to model the highly scattering nature of the shallow subsurface (e.g. Yilmaz *et al.* 2022).

It is also worth pointing out that the ambient noise waveform inversion, as implemented in, for example, Nouibat *et al.* (2023) and in this work, does not require an estimation of the source time-space parameters, since the source time function is assumed to be known (eq. 3).

Further investigations on this data set should attempt to account for attenuation, since the impact of visco-elasticity and multiscattering is expected to be non-negligible in the imaging of shallow subsurface targets (Groos *et al.* 2014), in particular when attempting to fit higher frequencies data.

Even if only $V_s$ is significantly updated during the multiparameter WET, the accuracy of the starting $V_p/V_s$ ratio values has an important impact in determining the depth sensitivity of Rayleigh waves phase velocities. Therefore, an accurate $V_p$, such as the one obtained by P-wave traveltime tomography, is key to obtain a reliable $V_s$ from Rayleigh waves inversion. When P-wave traveltime data is not available, Poisson's ratio should be inferred from surface wave data (e.g. Socco & Comina 2017). Although they may have posed additional constraints on the P-wave model building, the body-wave active-source waveforms have not been included in the WET data set because of their high frequency content (> 40 Hz), which makes their modelling and inversion computationally not attractive.

Useful confidence intervals may be obtained for the ANT+SVRI (starting) model within an affordable computing effort, for example, by bootstrapping in traveltime tomography (e.g. Vanorio *et al.* 2005), or by inferring them from the posterior distribution of models obtained by Monte-Carlo based inversion of dispersion curves (e.g. Nouibat *et al.* 2022). Propagating the starting model confidence interval to WET should account for the nonlinearity of the waveform-based inversion, whereby changes in the starting model may lead to cycle-skipping (Virieux *et al.* 2017). Along with a significant computing cost and the high-dimensionality of the problem, this makes uncertainty estimation an open research topic for full-waveform methods (Gebraad *et al.* 2020; Zhang & Curtis 2020; Hoffmann *et al.* 2024). In this study, we may preliminarily argue, based on our analysis on frequency-dependent group velocity tomography, that uncertainties are going to be larger in the southern part of the model, where illumination is weaker due to the higher trace rejection rate.

Regarding the geomechanical interpretation of the obtained results, we want to stress that both reconstructed P-wave and S-wave velocity models have low resolution in the first 10 m below the surface, which might explain the important (and unexpected) variability of the $V_p/V_s$ ratio in the landslide body zone. This is inherent to the methods and the data used to reconstruct these models. Also, we note that the detection of remoulded clays downslope was not performed before. This might question if this horizon is actually made of clays, which calls for complementary investigations using other geophysical imaging methods (for instance electrical resistivity or impedance). The distinction between clays, terraces and bedrock downhill could also be achieved using electrical
resistivity which has already illustrate its ability to distinguish these units (Jongmans *et al.* 2009; Bièvre *et al.* 2016).

## 7 CONCLUSIONS

Taking advantage of a dense seismic node acquisition, a high resolution 3-D elastic model of the Harmalière landslide has been obtained combining ambient noise and low-power, low-impact active sources. In order to make the most of the 3-D node data sets, we have applied interferometric techniques for the reconstruction of reasonable *P*-wave and *S*-wave velocity macromodels, which are then updated using a multiparameter full-waveform traveltime tomography.

The wave-equation based elastic model represents a significant improvement in accuracy and resolution with respect to the current state of knowledge of the study area, as well as with respect to conventional ambient-noise tomography, thanks to the true-3-D approach accounting for the complex topography of the slope. The obtained results provide new insights on the shallow geology of the site, the physical properties distribution within the landslide body and the structural controlling factors of the landslide dynamics, namely the transition between plastic and fluid deformation regimes and the sudden re-activations.

Future work should explore the attenuation properties within the slope and investigate the added value of the 3-component data set for the imaging of the very shallow subsurface, for example, by inverting for Love waves. These obtained seismic tomography models could be jointly interpreted in the future with complementary geological/geophysical data sets, both legacy and novel (e.g. electrical resistivity and geodesy), in order to provide a fully non-invasive 3-D characterization of the landslide body and better constrain the moisture content and stress state within the slope, ultimately providing a robust baseline for monitoring,

## ACKNOWLEDGMENTS

This study was partially funded by the SEISCPEO consortium (http://seiscipe2.osug.fr), sponsored by Aker BP, DUG, EXXONMOBIL, GEOLINKS, JGI, PETROBRAS, SHEARWATER, SHELL, SINOPEC, TOTALENERGIES and VIRIDIEN. This study was granted access to the HPC resources provided by the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities, the HPC resources of Cray Marketing Partner Network (https://partners.cray.com), and those of CINES/IDRIS/TGCC under the allocation 046091 made by GENCI. We used Seismic Unix, Python, Paraview and Matlab for data processing and figure production. The authors would like to thank all past and present SEISCPEO members, Jean Virieux for providing TOMO-TV, Pierre Boué and the Waves team at ISTerre, Grenoble for the useful discussion about ambient noise interferometry.

## DATA AVAILABILITY

The node data recorded on the Harmalière landslide used in this study is publicly available at the following link https://www.fdsn.org/networks/detail/ZN.2021/. Intermediate results (pre-processed data, first-arrival pickings, tomography models) are available upon reasonable request to the contact author. The full waveform inversion methodology used in this study is implemented in the SEM46 FORTRAN90 code developed within the framework of the SEISCPEO

project, funded by industrial partners. Under the contracts gathering the partners of the SEISCPEO project, it has been agreed that this code is not open-source. However the SEM46 code can be (and is currently) shared with academic researchers. The interested readers are encouraged to contact LM (ludovic.metivier@univ-grenoble-alpes.fr) and RB (romain.brossier@univ-grenoble-alpes.fr).

## REFERENCES

Bensen, G.D., Ritzwoller, M.H., Barmin, M.P., Levshin, A.L., Lin, F., Moschetti, M.P., Shapiro, N.M. & Yang, Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, *J. geophys. Int.*, **169**(3), 1239–1260.

Bharadwaj, P., Schuster, G., Mallinson, I. & Dai, W., 2012. Theory of supervirtual refraction interferometry, *J. geophys. Int.*, **188**, 263–273.

Bièvre, G., Kniess, U., Jongmans, D., Pathier, E., Schwartz, S., van Westen, C.J., Villemin, T. & Zumbo, V., 2011. Paleotopographic control of landslides in lacustrine deposits (Trièves Plateau, French Western Alps), *Geomorphology*, **125**(1), 214–224.

Bièvre, G., Jongmans, D., Goutaland, D., Pathier, E. & Zumbo, V., 2016. Geophysical characterization of the lithological control on the kinematic pattern in a large clayey landslide (Avignonet, French Alps), *Landslides*, **13**(3), 423–436.

Bièvre, G. & Crouzet, C., 2021. Multi-proxy analysis of boreholes in re-molded quaternary paraglacial deposits (avignonet landslide, french western alps), *Engineering Geology*, **286**, 106073.

Brenders, A.J. & Pratt, R.G., 2007. Efficient waveform tomography for lithospheric imaging: implications for realistic 2D acquisition geometries and low frequency data, *J. geophys. Int.*, **168**, 152–170.

Campillo, M., 2006. Phase and correlation in 'random' seismic fields and the reconstruction of the green function, *Pure appl. Geophys.*, **163**, 475–502.

Campillo, M. & Paul, A., 2003. Long Range Correlations in the Diffuse Seismic Coda, *Science*, **299**, 547.

Cao, J., Brossier, R., Górszczyk, A., Métivier, L. & Virieux, J., 2022. 3D multi-parameter full-waveform inversion for ocean-bottom seismic data using an efficient fluid-solid coupled spectral-element solver, *J. geophys. Int.*, **229**(1), 671–703.

Cappa, F., Guglielmi, Y., Viseur, S. & Garambois, S., 2014. Deep fluids can facilitate rupture of slow-moving giant landslides as a result of stress transfer and frictional weakening, *Geophys. Res. Lett.*, **41**(1), 61–66.

Carrière, S. *et al.*, 2018. Rheological properties of clayey soils originating from flow-like landslides, *Landslides*, **15**(8), 1615–1630.

Cuffey, K.M. & Paterson, W. S.B., 2010. *The Physics of Glaciers*, 4th edn, Butterworth-Heinemann/Elsevier.

Devaney, A., 1984. Geophysical diffraction tomography, *IEEE Trans. Geosci. Remote Sens.*, **GE-22**(1), 3–13.

Fiolleau, S., Jongmans, D., Bièvre, G., Chambon, G., Baillet, L. & Vial, B., 2020. Seismic characterization of a clay-block rupture in Harmalière landslide, French Western Alps, *J. geophys. Int.*, **221**, 1777–1788.

Fiolleau, S., Jongmans, D., Bièvre, G., Chambon, G., Lacroix, P., Helmstetter, A., Wathelet, M. & Demierre, M., 2021. Multi-method investigation of mass transfer mechanisms in a retrogressive clayey landslide (Harmalière, French Alps), *Landslides*, **8**, 1981–2000.

Rroude, M.J. & Petley, D.N., 2018. Global fatal landslide occurrence from 2004 to 2016, *Nat. Hazards Earth Syst. Sci.*, **18**(8), 2161–2181.

Fu, L., Liu, Z. & Schuster, G., 2017. Superresolution near-field imaging with surface waves, *J. geophys. Int.*, **212**(2), 1111–1122.

Gebraad, L., Boehm, C. & Fichtner, A., 2020. Bayesian elastic full-waveform inversion using hamiltonian monte carlo, *J. geophys. Res.: Solid Earth*, **125**(3), e2019JB018428.

Gimbert, F. *et al.*, 2021. A multi-physics experiment with a temporary dense seismic array on the Argentière glacier, French Alps: The resolve project, *Seismol. Res. Lett.*, **92**(2A), 1185–1201.

Groos, L., Schäfer, M., Forbriger, T. & Bohlen, T., 2014. The role of attenuation in 2D full-waveform inversion of shallow-seismic body and Rayleigh waves, *Geophysics*, **79**(6), R247–R261.
Hanafy, S. & Al-Hagan, O., 2012. Super-virtual refraction interferometry: an engineering field data example, *Near Surf. Geophys.*, **10**, 1-7.

Hoffmann, A., Brossier, R., Métivier, L. & Tarayoun, A., 2024. Local uncertainty quantification for 3-D time-domain full-waveform inversion with ensemble Kalman filters: application to a North Sea OBC data set, *J. geophys. Int.*, **237**(3), 1353-1383.

Huang, Y. & Schuster, G.T., 2014. Resolution limits for wave equation imaging, *J. Appl. Geophys.*, **107**, 137-148.

Huang, Y. & Schuster, G., 2014. Resolution limits for wave equation imaging, *J. Appl. Geophys.*, **107**, 137-148.

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.

Karimpour, M., Slob, E. & Socco, L.V., 2022. Comparison of straight-ray and curved-ray surface wave tomography approaches in near-surface studies, *Solid Earth*, **13**(10), 1569-1583.

Lacroix, P., Bièvre, G., Pathier, E., Kniess, U. & Jongmans, D., 2018. Use of sentinel-2 images for the detection of precursory motions before landslide failures, *Remote Sens. Environ.*, **215**, 507-516.

Lacroix, P., Handwerger, A.L. & Bièvre, G., 2020. Life and death of slow-moving landslides, *Nat. Rev. Earth Environ.*, **1**, 404-419.

Latorre, D., Virieux, J., Monfret, T., Monteiller, V., Vanorio, T., Got, J.L. & Lyon-Caen, H., 2004. A new seismic tomography of Aigion area (Gulf of Corinth, Greece) from the 1991 data set, *J. geophys. Int.*, **159**, 1013-1031.

Li, J., Lin, F.-C., Allam, A., Ben-Zion, Y., Liu, Z. & Schuster, G., 2019. Wave equation dispersion inversion of surface waves recorded on irregular topography, *J. geophys. Int.* **217**(1), 346-360.

Lu, K. & Chavez-Perez, S., 2020. 3d super-virtual refraction interferometry, *J. geophys. Int.*, **85**, Q1-Q10.

Lu, Y., Stehly, L., Brossier, R., Paul, A. & Group, A.W., 2020. Imaging Alpine crust using ambient noise wave-equation tomography, *J. geophys. Int.*, **222**(1), 69-85.

Luo, Y. & Schuster, G.T., 1991. Wave-equation traveltime inversion, *Geophysics*, **56**(5), 645-653.

Luu, K., 2021. disba: Numba-accelerated computation of surface wave dispersion. https://zenodo.org/records/14534395.

L'Heureux, J.-S. et al., 2012. Identification of Weak Layers and Their Role for the Stability of Slopes at Finneidfjord, Northern Norway, pp. 321-330, Springer Netherlands.

Maggi, A., Tape, C., Chen, M., Chao, D. & Tromp, J., 2009. An automated time-window selection algorithm for seismic tomography, *J. geophys. Int.*, **178**, 257-281.

Mansour, M.F., Morgenstern, N.R. & Martin, C.D., 2011. Expected damage from displacement of slow-moving slides, *Landslides*, **8**(1), 117-131.

Métivier, L. & Brossier, R., 2016. The seiscipe optimization toolbox: A large-scale nonlinear optimization library based on reverse communication, *Geophysics*, **81**(2), F11-F25.

Métivier, L., Brossier, R., Mérigot, Q. & Oudet, E., 2019. A graph space optimal transport distance as a generalization of LP distances: application to a seismic imaging inverse problem, *Inverse Probl.*, **35**(8), 085001.

Monjuvent, G., 1973. La transfluence Durance-Isère. Essai de synthèse du Quaternaire du bassin du Drac (Alpes françaises), *Géol. Alpine*, **49**, 57-118.

Nocedal, J. & Wright, S.J., 2006. *Numerical Optimization*, 2nd edn, Springer.

Nouibat, A. et al., 2022. Lithospheric transdimensional ambient-noise tomography of w-Europe: implications for crustal-scale geometry of the w-Alps, *J. geophys. Int.*, **229**(2), 862-879.

Nouibat, A., Brossier, R., Stehly, L., Cao, J., Paul, A. & Cifalps Team and AlpArray Working Group, 2023. Ambient-noise wave-equation tomography of the Alps and Ligurian-Provence Basin, *J. geophys. Res.: Solid Earth*, **128**(10), e2023JB026776.

Nuber, A., Manukyan, E. & Maurer, H., 2017. Optimizing measurement geometry for seismic near-surface full waveform inversion, *J. geophys. Int.*, **210**(3), 1909-1921.

Paige, C.C. & Saunders, M.A., 1982. ALGORITHM 583 LSQR : Sparse linear equations and least squares problems, *ACM Trans. Math. Softw.*, **8**(2), 195-209.

Place, J. & Malehmir, A., 2016. Using supervirtual first arrivals in controlled-source hardrock seismic imaging—well worth the effort, *J. geophys. Int.*, **206**, 716-730.

Plessix, R.E., 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications, *J. geophys. Int.*, **167**(2), 495-503.

Pugin, A. & Yilmaz, O., 2019. Optimum source-receiver orientations to capture pp, ps, sp, and ss reflected wave modes, *Leading Edge*, **38**(1), 45-52.

Renalier, F., Jongmans, D., Campillo, M. & Bard, P., 2010. Shear wave velocity imaging of the avignonet landslide (france) using ambient noise cross correlation, *J. geophys. Res.: Earth Surface*, **115**(F3), 2009JF001538.

Renalier, F., Bièvre, G., Jongmans, D., Campillo, M. & Bard, P.Y., 2012. 24. Clayey Landslide Investigations Using Active and Passive Vs measurements in *Advances in Near-surface Seismology and Ground-penetrating Radar*, pp. 397-413, Society of Exploration Geophysicists.

Rickwood, P. & Sambridge, M., 2006. Efficient parallel inversion using the neighbourhood algorithm, *Geochem. Geophys. Geosyst.*, **7**(11).

Sajeva, A., Aleardi, M., Galuzzi, B., Stucchi, E., Spadavecchia, E. & Mazzotti, A., 2017. Comparing the performances of four stochastic optimisation methods using analytic objective functions, 1d elastic full-waveform inversion, and residual static computation, *Geophys. Prospect.*, **65**, 322-346.

Sambridge, M.S., 1999. Geophysical inversion with a neighbourhood algorithm - I. Searching a parameter space, *J. geophys. Int.*, **138**, 479-494.

Schäfer, M., Groos, L., Forbriger, T. & Bohlen, T., 2013. 2D full waveform inversion of recorded shallow seismic Rayleigh waves on a significantly 2D structure, in *Proc. 19th European Meeting of Environmental and Engineering Geophysics, Expanded Abstracts, Bochum, Germany*, EAGE, B Houten, the Netherlands.

Schimmel, M., Stutzmann, E. & Gallart, J., 2011. Using instantaneous phase coherence for signal extraction from ambient noise data at a local to a global scale, *J. geophys. Int.*, **184**(1), 494-506.

Shapiro, N. & Campillo, M., 2004. Emergence of broadband Rayleigh Waves from correlations of the ambient seismic noise, *Geophys. Res. Lett.*, **31**.

Shapiro, N.M., Campillo, M., Stehly, L. & Ritzwoller, M.H., 2005a. High-resolution surface-wave tomography from ambient seismic noise, *Science*, **307**(5715), 1615-1618.

Shapiro, N.M., Campillo, M., Stehly, L. & Ritzwoller, M.H., 2005b. High-resolution surface-wave tomography from ambient seismic noise, *Science*, **307**(5715), 1615-1618.

Snieder, R., 1986. The influence of topography on the propagation and scattering of surface waves, *Phys. Earth planet. Inter.*, **44**(3), 226-241.

Socco, L.V. & Comina, C., 2017. Time-average velocity estimation through surface-wave analysis: Part 2—P-wave velocity, *Geophysics*, **82**(3), U61-U73.

Stork, C., 2020. How does the thin near surface of the earth produce 10-100 times more noise on land seismic data than on marine data?, *First Break*, **38**, 67-75.

Tape, C., Liu, Q., Maggi, A. & Tromp, J., 2010. Seismic tomography of the southern California crust based on spectral-element and adjoint methods, *J. geophys. Int.*, **180**, 433-462.

Trinh, P.T., Brossier, R., Métivier, L., Virieux, J. & Wellington, P., 2017. Bessel smoothing filter for spectral element mesh, *J. geophys. Int.*, **209**(3), 1489-1512.

Trinh, P.T., Brossier, R., Métivier, L., Tavard, L. & Virieux, J., 2019. Efficient 3D time-domain elastic and viscoelastic Full Waveform Inversion using a spectral-element method on flexible Cartesian-based mesh, *Geophysics*, **84**(1), R75-R97.

Vanorio, T., Virieux, J., Capuano, P. & Russo, G., 2005. Three-dimensional seismic tomography from *P* wave and *S* wave microearthquake travel times and rock physics characterization of the Campi Flegrei Caldera, *J. geophys. Res.*, **110**, B03201. doi: 10.1029/2004JB003102.

Vardy, M.E., Vanneste, M., Henstock, T.J., Clare, M.A., Forsberg, C.F. & Provenzano, G., 2017. State-of-the-art remote characterization of shallow marine sediments: the road to a fully integrated solution, *Near Surf. Geophys.*, **15**(4), 387-402.
Virieux, J., Asnaashari, A., Brossier, R., Métivier, L., Ribodetti, A. & Zhou, W., 2017. An introduction to Full Waveform Inversion, in *Encyclopedia of Exploration Geophysics*, pp. R1–1-R1-40, edsGrechka & Wapenaar K., V., Society of Exploration Geophysics.

Voisin, C., Garambois, S., Roux, P., Bièvre, G. & RESIF, 2021. *Resolve-harmaliere, temporary experiment for generate high-resolution seismic imaging and monitoring of the harmaliere landslide, france (resif-sismob).*

Wapenaar, K., Draganov, D., Snieder, R., Campman, X. & Verdel, A., 2010a. Tutorial on seismic interferometry: Part 1—basic principles and applications, *Geophysics*, **75**(5), A195–A209.

Wapenaar, K., Slob, E., Snieder, R. & Curtis, A., 2010b. Tutorial on seismic interferometry: Part 2—underlying theory and new advances, *Geophysics*, **75**(5), A211–A227.

Warner, M. & Guasch, L., 2016. Adaptive waveform inversion: theory, *Geophysics*, **81**(6), R429–R445.

Wathelet, M., 2008. An improved neighborhood algorithm: parameter conditions and dynamic scaling, *Geophys. Res. Lett.*, **35**, L09301.

Whiteley, J.S., Chambers, J.E., Uhlemann, S., Wilkinson, P.B. & Kendall, J.M., 2019. Geophysical monitoring of moisture-induced landslides: A review, *Rev. Geophys.*, **57**, 106–145.

Yilmaz, O., Gao, K., Delic, M., Xia, J., Huang, L., Jodeiri, H. & Pugin, A., 2022. A reality check on full-wave inversion applied to land seismic data for near-surface modeling, *Leading Edge*, **41**(1), 40–46.

Zhang, X. & Curtis, A., 2020. Variational full-waveform inversion, *J. geophys. Int.*, **222**(1), 406–411.

Zhou, Y., Liu, Q. & Tromp, J., 2011. Surface wave sensitivity: mode summation versus adjoint SEM, *J. geophys. Int.*, **187**, 1560–1576.