# Physical modelling of a surface-wave survey over a laterally varying granular medium with property contrasts and velocity gradients

Paolo Bergamo, Ludovic Bodet, Laura Valentina Socco, Régis Mourgues, Vincent Tournat

► To cite this version:

Paolo Bergamo, Ludovic Bodet, Laura Valentina Socco, Régis Mourgues, Vincent Tournat. Physical modelling of a surface-wave survey over a laterally varying granular medium with property contrasts and velocity gradients. Geophysical Journal International, 2014, 197 (1), pp.233-247. 10.1093/gji/ggt521. hal-01196396

HAL Id: hal-01196396

https://hal.science/hal-01196396

Submitted on 7 Mar 2022

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Physical modelling of a surface-wave survey over a laterally varying granular medium with property contrasts and velocity gradients

Paolo Bergamo,<fnref n="1" /> Ludovic Bodet,<fnref n="2" /> Laura Valentina Socco,<fnref n="1" /> Régis Mourgues<fnref n="3" /> and Vincent Tournat<fnref n="4" />

Accepted 2013 December 23. Received 2013 December 12; in original form 2012 April 30

## SUMMARY

Laboratory experiments using laser-based ultrasonic techniques can be used to simulate seismic surveys on highly controlled small-scale physical models of the subsurface. Most of the time, such models consist in assemblies of homogeneous and consolidated materials. To enable the physical modelling of unconsolidated, heterogeneous and porous media, the use of granular materials is suggested here. We describe a simple technique to build a two-layer physical model characterized by lateral variations, strong property contrasts and velocity gradients. We use this model to address the efficiency of an innovative surface-wave processing technique developed to retrieve 2-D structures from a limited number of receivers. A step by step inversion procedure of the extracted dispersion curves yields accurate results so that the 2-D structure of the physical model is satisfactorily reconstructed. The velocity gradients within each layer are accurately retrieved as well, confirming current theoretical and experimental studies regarding guided surface acoustic modes in unconsolidated granular media.

**Key words:** Fourier analysis; Surface waves and free oscillations; Guided waves; Wave propagation; Acoustic properties.

## 1 INTRODUCTION

Non-contacting ultrasonic techniques, such as mechanical-wave generation using high frequency laser sources and particle motion measurement using laser-Doppler vibrometers, represent powerful tools to address the physics of various processes in geophysics. Their non-contacting character makes it possible to generate mechanical waves or to record particle motion at the surface of many different types of materials, without any coupling. Additionally, these measurement devices allow very fine resolutions and present high-density sampling abilities. In seismology, laboratory experiments involving lasers are chosen to study the propagation of seismic waves in random heterogeneous media, when numerical models fail to predict the actual complexity of Earth materials (Nishizawa et al. 1997; Sivaji et al. 2002; Spetzler et al. 2002; Scales & Malcolm 2003; Nishizawa & Kitagawa 2007). However such experiments must not be considered as a mere alternative to numerical modelling since both approaches often appear to be highly complementary (Scales & van Wijk 1999; van Wijk et al. 2004) and to provide insights for specific applications (van Wijk & Levshin 2004).

Physical modelling and laser experiments are similarly proposed to tackle theoretical and methodological issues related to the domain of exploration seismology, more particularly when experimental

validations of processing or inversion techniques are required. In this context, lasers are mainly used to reproduce typical field seismic acquisition set-ups at the laboratory scale (e.g. Hayashi & Nishizawa 2001; Campman et al. 2004, 2005; Bodet et al. 2005, 2009; Dewangan et al. 2006; Kaslilar 2007; Blum et al. 2011; Bretauveau et al. 2011; De Cacqueray et al. 2011). Most of the time, the controlled character of both the geometry and mechanical properties of specifically designed ‘small-scale physical models’ makes it possible to address, using real data, the efficiency, robustness or limitations of studied methods. Homogeneous and consolidated materials such as metals and thermoplastics are thus frequently used because they offer a wide range of mechanical parameters and they can be easily manufactured into different shapes and sizes. Increasing degrees of complexity can then be chosen for the models construction, in order to mimic even complex structures (Bretauveau et al. 2011). However from civil engineering to seismology, there is an obvious need to study the propagation of seismic waves in more complex and realistic media in terms of intrinsic parameters.

Granular materials [such as natural sand, granular silica or glass beads (GB)] offer an evident flexibility in terms of physical models construction and choice of parameters. Their mechanical properties are the subject of active investigations in various application fields (Valverde & Castellanos 2006; Hentschel & Page 2007) and they

<footnote n="1">¹Department of Environment, Land and Infrastructure Engineering, Politecnico di Torino, Turin, Italy. E-mail: paolo.bergamo@polito.it</footnote>

<footnote n="2">²Sorbonne Universités, UPMC Univ Paris 06, UMR 7619, METIS, F-75005, Paris, France</footnote>

<footnote n="3">³UMR CNRS 6112 LPGN, Université du Maine, Le Mans, France</footnote>

<footnote n="4">⁴UMR CNRS 6613 LAUM, Université du Maine, Le Mans, France</footnote>
are used for both physical and geological modelling purposes (e.g. Sherlock & Evans 2001; Graveleau 2008; Jacob et al. 2008; Buddensiek 2009; Bodet et al. 2010, 2011; Krawczyk et al. 2013). Bodet et al. (2010) recently addressed the ability of laser-Doppler experiments in the systematic characterization of granular materials involved in geological analogue modelling. The methodology was validated on an unconsolidated granular laboratory medium, perfectly characterized in terms of elastic parameters by Jacob et al. (2008). A mechanical source was developed and a laser-Doppler visbrometer was used to record small-scale seismic lines at the surface of the granular medium. Pressure-wave first arrival times and P-SV waves dispersion were then inverted for 1-D pressure- and shear-wave propagation velocity profiles with depth. Inferred velocity structures appeared to match previously thoroughly estimated properties of the probed medium, thus validating this approach (Bodet et al. 2010). This technique has been thoroughly addressed in terms of reproducibility of the medium preparation using micrometric GB (Bodet et al. 2011). These experiments involved single layer models, basically made up of uncompacted GB with a relatively uniform diameter size: the main aim of the present work is to test the feasibility of building granular physical models with a more complex geometry, that is, not only constituted by 2-D structures and property contrasts, but also with continuously varying elastic behaviours within layers.

In addition, we propose here to show that such granular physical model can help addressing the efficiency of surface-wave dispersion measurements in the straightforward characterization of laterally varying media. Surface-wave dispersion inversion is well known as a standard approach to infer shear-wave velocity ($V_S$) structure for global earth seismology, near-surface geophysics, geotechnical and civil engineering applications, or as a Non-Destructive Evaluation (NDE) tool (e.g. Gabriels et al. 1987; Wu & Liu 1999; Ruiz & Nagy 2004; Socco & Strobbia 2004; Socco et al. 2010). At the applied geophysics scale, the recorded wavefield is basically transformed to the frequency-wavenumber (or frequency-slowness) domain, in which dominant surface-wave (Rayleigh or Love, depending on the source and sensors) propagation modes can be relatively easily picked as dispersion curves. Dispersion curves are then inverted for a 1-D $V_S$ profile with depth. A well-known limitation of the method is its trade-off between lateral resolution and investigation depth (Gabriels et al. 1987). On one hand, the inverse problem formulation imposes the investigated medium to be assumed 1-D below the spread, which must be short enough to achieve lateral resolution. On the other hand, long spreads are required to record wavelengths long enough to increase investigation depth and to mitigate near-field effects (Bodet et al. 2005, 2009). Several methods have been proposed to overcome such limitations (Bohlen et al. 2004; Hayashi & Suzuki 2004; Nagai et al. 2005; Grandjean & Bitri 2006; Neduca 2007; Luo et al. 2009; Socco et al. 2009; Boiero & Socco 2010; Strobbia et al. 2011). Bergamo et al. (2012) recently proposed a processing technique to retrieve 2-D structures from a limited number of receivers, using a moving Gaussian window that sweeps the receiver spread. This technique was successfully tested on synthetic data sets and applied to a real case (Bergamo et al. 2012). This work tends to bridge between the tests on synthetic data and the applications to real acquisitions. We address here the efficiency and robustness of the technique on a real data set extracted from a small-scale physical model consisting of two layers separated by a dipping interface and characterized by two distinct non-linear velocity gradients with depth.

In the following, we briefly describe the setup of the laboratory experiment, which is based on theory of wave propagation in loose granular materials and on preliminary numerical simulations. Then we present the results of acquisition and processing and we describe the inversion approach used to estimate the velocity model. We finally present and discuss the results.

## 2 PHYSICAL MODEL CONSTRUCTION AND EXPERIMENTAL SET UP

### 2.1 Relevant aspects of wave propagation in unconsolidated granular materials

Uncompacted granular packed structures at low pressures are characterized by a strong depth dependence of seismic properties distribution. As advocated by Gassmann (1951), in such materials P- and S-wave velocities show a power-law dependence on overburden pressure:

$$V_{P,S} = \gamma_{P,S}(p)^{\alpha_{P,S}}, \quad (1)$$

where p is the pressure, $\gamma_{P,S}$ is a depth-independent coefficient mainly depending on the elastic properties of the grains (e.g. Makse et al. 2004 for details about the theoretical background) and $\alpha_{P,S}$ is the power-law exponent. In a single material model, $p = gz\rho$ with g the gravity acceleration, z the depth and $\rho$ the bulk density, which can be considered depth independent. The presence of a strong velocity gradient with depth in combination with the free surface gives rise to the propagation of dispersive, guided surface acoustic modes (GSAM; Gusev et al. 2006; Aleshin et al. 2007; Jacob et al. 2008; Tournat & Gusev 2010). GSAM propagate close to the surface and according to their direction of polarization we distinguish shear horizontal (SH) modes and modes polarized in the vertical plane, named P-SV as they result from the interaction of longitudinal (P) and shear vertical (SV) waves. Bodet et al. (2010, 2011) proved the ability of surface wave method in handling P-SV modes for the characterization of near-surface continuous velocity gradients: they demonstrated the possibility to consider the recorded wavefield in the framework of elastic-wave propagation in stratified media instead of GSAM theory. Another interesting element of wave propagation in loose granular media is that low order, slower P-SV modes exhibit a power-law dispersion relation, for which the modes eigenfrequencies are proportional to $k^{1-\alpha_S}$ (with $k$ = wavenumber), and therefore phase velocities are proportional to $\lambda^{\alpha_S}$ (with $\lambda$ = wavelength; see Gusev et al. 2006; Aleshin et al. 2007; Jacob et al. 2008). If phase velocities are proportional to $\lambda^{\alpha_S}$, it is possible to estimate the ratio between the wavelengths and the depths at which the corresponding phase velocities refer to (which we call pseudo-depths): this ratio can be considered independent on the wavenumber for short-wavelength P-SV modes (as in our case: see Jacob et al. 2008) and it has a different value for every mode.

For the interpretation of our data we therefore exploited the following aspects of wave propagation in unconsolidated granular materials: (i) the presence of a definite power-law linking seismic velocities to pressure, (ii) the possibility to handle P-SV modes via the surface wave method and (iii) the proportionality, valid in low-order P-SV modes, between phase velocity and $\lambda^{\alpha_S}$.

### 2.2 Preliminary numerical simulations

Prior to the model construction, we performed a series of numerical finite element method (FEM) simulations for the calibration of the physical model: we were interested in estimating the optimal acquisition geometry, the frequency band in which we could perform
<table><caption>**Table 1.** Synthetic models. Parameters $\alpha_S$, $\gamma_S$ and $\rho$ of Gassmann relationship (eq. 1) we assumed for the materials constituting the synthetic models A and B and upper layer thickness for model B.</caption><thead><tr><th></th><th>Parameters</th><th>Model A</th><th>Model B</th></tr></thead><tbody><tr><td rowspan="4">Layer 1</td><td>$\alpha_S$</td><td>0.33</td><td>0.33</td></tr><tr><td>$\gamma_S$</td><td>5.25</td><td>5.25</td></tr><tr><td>$\rho$ (kg m<sup>-3</sup>)</td><td>1580</td><td>1580</td></tr><tr><td>Thickness (m)</td><td>—</td><td>0.05</td></tr><tr><td rowspan="3">Layer 2</td><td>$\alpha_S$</td><td>—</td><td>0.23</td></tr><tr><td>$\gamma_S$</td><td>—</td><td>19</td></tr><tr><td>$\rho$ (kg m<sup>-3</sup>)</td><td>—</td><td>1600</td></tr></tbody></table>

surface-wave analysis and the maximum investigation depth. The results we obtained were also useful for the interpretation of experimental data. We used the FEM code COMSOL Multiphysics®. In particular, we ran two numerical analyses simulating the propagation of a seismic signal (a Ricker pulse centred at 1.5 kHz, similar to the source we intended to use for laboratory experiments) into two synthetic models: the first model (in the following model A) is composed by one material only, while the second (model B) is a two-layer model with 1-D geometry. In both models we simulated the seismic properties of unconsolidated granular materials: according to the experience of Jacob *et al.* (2008) and Bodet *et al.* (2010, 2011) for wave propagation in granular materials we imposed the relationship introduced by Gassmann (1951) where P- and S-wave velocity are considered as power-law dependent on overburden pressure (eq. 1).

Table 1 reports the values of $\alpha_S$, $\gamma_S$ and bulk density ($\rho$) we assumed for the materials constituting the synthetic models A and B. Note that the properties of the upper layer of model B are the same as the material of model A and coincide with the values of $\alpha_S$, $\gamma_S$ and $\rho$ retrieved by Bodet *et al.* (2010) from experiments on a granular material similar to the ones we intended to use for our laboratory model. From the FEM simulations on model A and B we measured the vertical velocity at the top surface of the models, obtaining two synthetic seismograms, constituted by 499 traces with offset comprised between 2 and 500 mm and 1 mm spacing. We applied an f-k analysis to both seismic gathers and we obtained the two dispersion curves displayed in Figs 1(a) and (b): as we recorded the vertical component of the velocity, these curves represent the dispersion relation of P-SV modes. Both curves stretch over a wide frequency band (200–1000 Hz): the curve from the two-layer model presents greater phase velocity values with respect to the curve from the homogeneous model for frequencies <400 Hz, suggesting that the investigation depth reaches the lower layer of model B and it is therefore appropriate for the scale of our laboratory experiments. These numerical simulations also show that, in case of a single layer model, the dispersion curve fundamental mode represented in terms of pseudo-depth (which is a fraction of the wavelength) versus pseudo-$V_S$ (phase velocity times 1.1, Abbiss 1981) follows the power-law trend described by eq. (1) for $V_S$ (see Fig. 1c). This agreement is partially verified in the case of two layer model (Fig. 1d). In this case, the acoustic impedance contrast between the two layers causes a mode jump from fundamental to first higher mode at low frequencies, so that the curve

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Synthetic data. (a) Theoretical dispersion curve (grey line) for the $V_S$ profile of the synthetic model A and dispersion curve extracted from model A synthetic seismogram (black dots). (b) Theoretical dispersion curve (grey line) for the $V_S$ profile of the synthetic model B and dispersion curve extracted from model B synthetic seismogram (black dots). (c) $V_S$ profile for model A (grey line) and dispersion curve extracted from model A synthetic seismogram (black dots) represented in pseudo-depth (wavelengths/2.62) versus pseudo-$V_S$ (phase velocities times 1.1). (d) $V_S$ profile for model B (grey line) and dispersion curve extracted from model B synthetic seismogram (black dots) represented in pseudo-depth (wavelengths/2.62) versus pseudo-$V_S$ (phase velocities times 1.1).</figcaption></figure>
represented as pseudo-depth versus pseudo-Vs does not follow the trend predicted by eq. (1) for pseudo-depths greater than the interface depth (Fig. 1d). The agreement between the Vs profile and the dispersion curve fundamental mode is explained by GSAM theory: in uncompacted granular materials low order, slower P-SV modes exhibit a power-law dispersion relation, for which phase velocities are proportional to λαs (with λ = wavelength; see Section 2.1). It is therefore possible to estimate the ratio between the wavelengths and the depths at which the corresponding fundamental mode phase velocities refer to. To determine this ratio we performed a parametric analysis: the numerical dispersion curve from model A was translated into the pseudo-Vs versus pseudo-depth plane assuming 1.1 as the ratio between pseudo-Vs and the fundamental mode phase velocities (Abbiss 1981) and by successively adopting a different value for the ratio between wavelengths and pseudo-depths, chosen within the range 1.5–3. For every value of this ratio, γs and αs were estimated by linearizing eq. (1) and applying a least square approach:

$$ \begin{bmatrix} \log(\gamma_s) \\ \alpha_s \end{bmatrix} = (\mathbf{G}^T \mathbf{G})^{-1} \mathbf{G}^T \log(\mathbf{v}_s), \quad (2) $$

where vs is the vector of pseudo-Vs of the extracted dispersion curve and

$$ \mathbf{G} = \begin{bmatrix} 1 & \log(\rho g z_1) \\ \vdots & \vdots \\ 1 & \log(\rho g z_{n-1}) \\ 1 & \log(\rho g z_n) \end{bmatrix} \quad (3) $$

with g the gravity acceleration, ρ the bulk density and z₁...zₙ the pseudo-depths of the extracted dispersion curve.

To linearize eq. (1) we express the power-law in logarithmic scale where it becomes a straight line, whose slope is αs and log(γs) is the intercept. Hence, the value of the ratio used to transform wavelength into pseudo-depth does not affect the estimation of αs (Fig. 2a). The lowest error in the estimation of γs is obtained for the ratio value of 2.62 (Fig. 2b). Hereafter, we will then use this value to express the ratio between wavelengths and pseudo-depths. Once defined the wavelength/pseudo-depth ratio, the approach expressed by eq. (2) enables to estimate γs and αs from the experimental dispersion curve in a way that does not require forward modelling computation. It is worth stressing that the described analysis is valid for single material models, hence, in case of a two layer model, it

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Synthetic data, parametric analysis to determine the optimal value of the wavelength/pseudo-depth ratio. (a) Relative error in the estimation of αs as a function of the wavelength/pseudo-depth ratio. (b) Relative error in the estimation of γs as a function of the wavelength/pseudo-depth ratio.</figcaption></figure>

can be applied only to the dispersion curve data points which are relevant to the uppermost layer only and are not affected by the presence of the deeper layer.

## 2.3 Physical model construction

In this study, the medium preparation basically reproduced the technique previously proposed by Jacob et al. (2008) and Bodet et al. (2010) and originally designed to set-up an unconsolidated granular packed structure under gravity. When the experiments require great quantities of such material, it remains difficult to ensure a perfect homogeneity of the medium. This is why this technique has been recently thoroughly addressed in terms of reproducibility of the medium preparation using micrometric GB (Bodet et al. 2011). Similar GB were chosen here to build a two-layer physical model presenting both in-depth property gradients and a sloping interface: the experimental setup and the model are shown in Fig. 3 and described with detail in Fig. 4.

The GB presented two ranges of diameter: 200–300 µm (GB1) and 100–200 µm (GB2). These GB were sieved or poured into a 1000 × 800 × 300 mm wooden box, presented on Fig. 4. Compared to the experiments cited above the box dimensions were chosen large enough to limit wave-reflections from sidewalls and consequently avoid boundary effects. However, due to obvious practical limitations in terms of material quantities and model weight, the height of the box remained too small to avoid bottom reflections. The bottom of the box moreover consisted of a micrometric metallic

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Experimental set-up. In the photograph it is possible to single out the seismic source (labelled with ‘A’), the granular medium (‘B’), the wooden box containing the glass beads (‘C’) and the laser-Doppler unit (‘D’). Dimensions of the set-up are given in detail in Fig. 4.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Experimental set-up. The medium consisted of two glass beads layers, presenting two ranges of diameter: 200–300 µm (GBL1) and 100–200 µm (GBL2), prepared in a wooden box and separated by a sloping interface. The bulk porosities and densities were estimated as equal to 0.38 and 1550 kg m⁻³ for GBL1 and 0.42 and 1440 kg m⁻³ for GBL2. A laser-Doppler vibrometer was used to record seismograms of particle vertical velocity at the surface of the medium excited by a mechanical source at position (x_s, y_s). Following a step by step procedure, the laser scanned the surface and, thanks to an oscilloscope, particle normal velocities were recorded at each offset along a linear section (line), the source remaining still.</figcaption></figure>

sieve glued on a perforated plate that was originally built to inject overpressured air into the granular material and to study the influence of pore overpressure on seismic-wave propagation at very low confining pressures (Bodet *et al.* 2011). Strong bottom reflections were thus expected due to high impedance contrast between the granular material and the metallic sieve. In order to mitigate such bottom reflections and their possible influence on the data, several countermeasures were proposed by Bodet *et al.* (2011), among which the simple laying of a 4-mm-thick paper-board plate onto the sieve (depicted by the thick black line on Fig. 4).

The bottom layer of the physical model (GBL1 on Fig. 4) was built by sieving GB1 through a 400 µm sieve. This layer was compacted every centimetre by means of a flattening tool and vibrations applied to the box. Its final thickness was 205 mm. The bulk porosity ($\varphi_{GBL1}$) and density ($\rho_{GBL1}$) of the medium were estimated, by means of laboratory measurements, as equal to 0.38 and 1550 kg m⁻³, respectively. A sloping surface was then thoroughly dug by half centimetre steps between the positions x = 350 and 650 mm, as presented on Fig. 4. The slope was approximately 7.6°. The sloping interface was graded flat by stamping the surface with a flattening tool, as recommended by Buddensiek (2009). The top layer (GBL2 on Fig. 4) was finally achieved by thoroughly pouring GB2 onto GBL1. Its final thickness varied from 20 to 60 mm. We did not compact this layer in order to ensure a porosity contrast between GBL1 and GBL2. Again, the bulk porosity ($\varphi_{GBL2}$) and density ($\rho_{GBL2}$) of GBL2 were estimated, by means of laboratory measurements, as equal to 0.42 and 1440 kg m⁻³, respectively.

The physical model presents a total thickness of 225 mm and can be considered as a sequence, in the x-direction, of three different geometries: a 2-layer model with an upper layer thickness of 20 mm between the positions x = 0 and 350 mm; a 2-layer model with an

upper layer thickness varying from 20 to 60 mm between the positions x = 350 and 650 mm; a 2-layer model with an upper layer thickness of 60 mm between the positions x = 650 and 1000 mm (Fig. 4). The contrast between the two layers is ensured via the estimated porosities which compare well with previous studies (Sherlock 1999; Sherlock & Evans 2001; Graveleau 2008; Buddensiek 2009). The porosity of GB in GBL1 appears close to the random close packing limit, while it tends to the random loose packing limit in GBL2 (Valverde & Castellanos 2006). The additional complexity and interest of this physical model is that each layer presents a gravity-induced rigidity gradient that will be described in the following.

## 2.4 Data acquisition

As previously proposed for the physical modelling of surface-wave propagation at the laboratory scale (van Wijk *et al.* 2004; Campman *et al.* 2005; Bodet *et al.* 2005, 2009), laser-Doppler acquisitions were performed here to simulate a seismic survey at the surface of the physical model. The experimental set-up, adapted to granular materials by Jacob *et al.* (2008) and Bodet *et al.* (2010, 2011), basically involved a laser-Doppler vibrometer to record seismograms of particle vertical velocity at the surface of the medium which was excited by a mechanical source at position (x_s, y_s), as presented in Fig. 4. The source consisted in a 120 mm long, 3 mm diameter metallic stick attached to a low-frequency (LF) shaker and buried into the granular material with an angle of 20° from the normal to the free surface. The force source signal was sent from a waveform generator to the LF-shaker, exciting the stick. A detailed view of the source set-up and its characteristics is given in Fig. 5. The force source signal (blue lines on Fig. 5) was a Ricker pulse with
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Source set-up and characteristics. (a) The force source signal (blue line on the left inset) was sent from a waveform generator to a low-frequency (LF) shaker exciting a metal stick buried in the granular material. The laser beam was set at the zero offset position ($x_s$ on Fig. 4) to record the stick normal velocity (red line on the right inset). (b) The force source signal was a Ricker pulse with its frequency spectrum centred on 1.5 kHz (blue lines). Except for a slight shift toward low frequencies due to possible 'ringing of the stick', the signal recorded at zero offset position (red lines), differentiated to acceleration here, is in very good agreement with the original force source signal (blue lines).</figcaption></figure>

its frequency spectrum centred on 1.5 kHz (Fig. 5b). Before each
seismogram acquisition, the laser beam was set at the zero offset
position ($x_s$ on Fig. 4) to record the stick normal velocity (red lines
on Fig. 5). Despite the accuracy of the LF shaker in restoring the the-
oretical signal coming from the waveform generator (as shown on
Fig. 5b), it is experimentally difficult to avoid low frequency vibra-
tions or ‘ringing’ of the stick (due to its length) right after the Ricker
pulse. In similar experiments previously performed, several coun-
termeasures were tested, trying to mitigate these effects (Jacob *et al.*
2008; Bodet *et al.* 2010, 2011), but it appeared impossible to com-
pletely remove them. To record a seismogram, the source remained
still and the laser-Doppler vibrometer scanned by constant steps

(1, 2 or 3 mm) the surface of the granular medium. The laser spot
size on the surface was between 1 and 2 mm in diameter for each
step and was primarily determined by the quality of the received
Doppler signal (note that it represents ~10 bead diameters within
the spot diameter). Using an oscilloscope, up to 300 traces were
recorded in linear single-channel walkway mode, parallel to the
long edges of the box (see Table 2 for details about the acquisition
geometries). Each trace was stacked 50 times and the time sampling
rate was 100 kHz over 5002 samples. Each seismogram was then
300 mm in distance and 50 ms in time.

As shown on Fig. 6(a) a first seismogram of particle verti-
cal velocity was recorded at the surface of GBL1 before the
<table><caption>**Table 2.** Laboratory experiments. Acquisition parameters of collected seismograms (illustrated on Fig. 6).</caption><thead><tr><th>Line number</th><th>Source position x<sub>s</sub> (mm)</th><th>Spacing between traces (mm)</th><th>Total number of traces</th><th>Source-first receiver distance (mm)</th><th>Source-last receiver distance (mm)</th><th>Sampling rate (kHz)</th><th>Number of time samples</th></tr></thead><tbody><tr><td>1</td><td>350</td><td>3</td><td>100</td><td>3</td><td>300</td><td>100</td><td>5002</td></tr><tr><td>2</td><td>50</td><td>2</td><td>150</td><td>2</td><td>300</td><td>100</td><td>5002</td></tr><tr><td>3</td><td>370</td><td>1</td><td>300</td><td>1</td><td>300</td><td>100</td><td>5002</td></tr><tr><td>4</td><td>950</td><td>2</td><td>150</td><td>2</td><td>300</td><td>100</td><td>5002</td></tr><tr><td>5</td><td>650</td><td>1</td><td>300</td><td>1</td><td>300</td><td>100</td><td>5002</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Data acquisition. (a) A seismogram was recorded at the surface of GBL1 before the preparation of the sloping interface (Line 1). (b) Four seismograms were recorded on the two-layer physical model (Lines 2, 3, 4, 5). The lines are red when the laser scanned lines along the positive x-direction (direct shot) and green along the negative x-direction (reverse shot), the stars giving sources position.</figcaption></figure>

preparation of the sloping interface. This seismogram should allow a good characterization of the bottom layer. Once the final model was prepared, four seismograms were recorded at its surface (Fig. 6b). Only one shot was recorded above each horizontally layered (1-D) part of the model (Line 2 and Line 4 on Fig. 6b). However, the 2-D central part was probed using both direct and reverse shooting

(Line 3 and Line 5 on Fig. 6b) in order to provide an optimal data set to be processed using the technique proposed by Bergamo *et al.* (2012). The recorded lines were 300 mm long with a space sampling interval of 2 mm above the flat-interface parts of the model (Lines 2 and 4) and of 1 mm above the sloping interface (Lines 3 and 5). All recorded seismograms are shown in Fig. 7.

The recorded wavefields clearly present coherent events which we identify as P-SV dispersive modes (marked with P-SV) and fast guided modes considered as mainly longitudinal by Jacob *et al.* (2008) and therefore indicated with P (Fig. 7). Along Lines 1 and 5, it is possible to guess bottom-reflected arrivals (marked with bR in the Lines 1 and 5 insets of Fig. 7). The seismograms are corrupted by strong amplitudes typically related to the source ringing (labelled with Sr in Fig. 7). These experimental artefacts will be mitigated by deconvolution with the source signal, recorded at the zero offset position (Fig. 5b). The P-SV events identified on every seismogram can be quantitatively interpreted using typical seismic surface wave processing techniques, as proposed by Bodet *et al.* (2010). We propose to use the physical model and recorded data and to apply the technique previously developed by Bergamo *et al.* (2012) with the aim of retrieving lateral variations of $V_S$ in 2-D media.

## 3 DISPERSION CURVES EXTRACTION

Once traces were collected, a pre-processing stage was necessary before the dispersion curves extraction. First of all, a high-pass filter at 100 Hz was applied in order to suppress ambient noise that might have corrupted the acquisitions. Secondly, raw seismograms (Fig. 7) appeared affected by evident source ringing as noticed in the previous section (Fig. 5). To mitigate this effect on the seismograms, raw traces were deconvolved with the source signal (i.e. the signal recorded by setting the laser beam on the stick just before it plunges into the GB bed, see Fig. 5). The operation of deconvolution is

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Raw data. Recorded seismograms globally show that, despite strong amplitudes associated to the source ringing (Sr), the wavefields clearly present coherent events such as P- and P-SV wave trains. Along the Line 1 and the Line 5, it is possible to guess bottom-reflected arrivals (bR).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Experimental data. Seismic sections of Fig. 7 after high-pass filtering and deconvolution.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Experimental data. Dispersion curves from the five survey lines of Fig. 6. ‘P–SV’ labels P–SV modes while ‘P’ indicates a predominantly longitudinal, fast guided mode.</figcaption></figure>

carried out by convolving the raw seismograms for a filter operator which is the mathematical inverse of the source signal (Yilmaz 1988). All recorded seismograms after filtering and deconvolution are presented in Fig. 8. Once every seismogram was pre-processed, we extracted a dispersion curve from every shot by applying an *f*-k transform to the seismograms and then picking the energy maxima on the spectra: the picked dispersive events are displayed in Fig. 9. From all seismic lines three dispersive events were retrieved: the two slower events, with phase velocities below 120 m s⁻¹, are the *P*–*SV* fundamental and first higher mode (labelled with *P*–*SV* in Fig. 9), while the third event, with velocities between 100 and 180 m s⁻¹ and frequencies lying in the 700–1300 Hz range can be identified as a predominantly longitudinal, fast guided mode (indicated with *P* in Fig. 9: see Jacob *et al.* 2008). Not surprisingly, the dispersive events from Line 1 are characterized by higher phase velocity values, as the lower layer is stiffer than the upper one. As for the other four dispersion curves (Lines 2–5), *P*–*SV* fundamental mode phase velocities are quite similar from 650 Hz on: at lower frequencies the dispersion curves are affected by the lower layer and show a rapid increase of the velocity. In particular, the frequency at which this increase takes place depends on the depth of the interface between

GBL1 and GBL2 and therefore it is roughly the same (about 500 Hz) for the two dispersion curves relevant to the dipping interface part of the model (curves from Lines 3 and 5). This threshold frequency becomes smaller (around 400 Hz) for the dispersion curve from the portion of the model where the interface is deeper (dispersion curve from Line 4) and it is greater (650 Hz) for the curve referring to the part of the model where the interface is shallower (curve from Line 2). Moreover, to get a better description of the lateral variation of the physical model, we applied to the central shots (Lines 3 and 5) an algorithm developed by Bergamo *et al.* (2012) based on a spatial windowing of the traces to extract a set of local dispersion curves from a single shot. This method applies to the seismogram a spatial windowing which is based on Gaussian windows whose maxima span the survey line, thus assigning from time to time a different weight to the seismic traces. Therefore, we are able to extract a set of dispersion curves, each of them referring to a different subsurface portion. These dispersion curves are reported on Fig. 10, where they are referred to the different positions of the Gaussian window maxima along the receiver spread. Again, for every position of the Gaussian window maximum three dispersive events were identified: *P*–*SV* fundamental and first higher mode (with phase
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Experimental data: dispersion curves retrieved from seismograms of Lines 3 and 5 by applying a space-varying spatial windowing. ‘P–SV’ labels P–SV modes while ‘P’ indicates a predominantly longitudinal, fast guided mode.</figcaption></figure>

velocities lower than 120 m s⁻¹, indicated with P–SV in Fig. 10) and a predominantly longitudinal, fast guided mode (with velocities between 80 and 170 m s⁻¹ and frequencies within 700–1700 Hz, labelled with P in Fig. 10). Similarly as observed in Fig. 9 the velocity of the P–SV fundamental modes starts to differ at around 400–600 Hz, depending on the interface depth of the vertical profile at the dispersion curves reference points (i.e. the locations of the Gaussian window maxima). As shown in Figs 9 and 10 dispersion data were mainly collected in the 200–1600 Hz frequency band. The thickness of the physical models is then close to the largest wavelengths for P–SV modes (observed at frequencies close to 200 Hz). However, it is worth remarking that P–SV modes wavelengths rapidly decrease as frequency increases, down to approximately 10 mm at high frequencies (~1600 Hz). Moreover, since most of P–SV modes energy propagates between the free surface and depths close to half a wavelength (Aleshin et al. 2007), dispersion data are not affected by the base of the box, as previously ensured by Jacob et al. (2008) and Bodet et al. (2010, 2011) using similar configurations. To address the possible influence of the boundaries of the model on the acquired data, we also performed additional numerical simulations and forward modelling computations, not shown here, which confirmed the negligibility of the effects of the bounding box on processed GSAM dispersion data.

## 4 DISPERSION CURVES INVERSION

The last stage of the process involved the dispersion curves inversion. As already successfully done by Jacob *et al.* (2008) and Bodet *et al.* (2010, 2011) we adopted the relationship introduced by Gassmann (1951) where P- and S-wave velocity are considered as power-law dependent on overburden pressure (eq. 1). We therefore assumed that in the two-layer physical model the S-wave velocities follow a power-law trend which is controlled by seven parameters: the two sets of $\alpha_S$, $\gamma_S$ and $\rho$ for the upper and lower layer, and the interface depth. Within each layer, the bulk density can be assumed constant with depth even if it may slightly vary in the vicinity of the free surface, at pressures less than 75 Pa, as previously noted by Jacob *et al.* (2008). However, the contribution of such density variation to velocity change with depth is indeed negligible (Gusev

*et al.* 2006). Bulk densities for GBL1 and GBL2 are not unknown since they were estimated by means of laboratory experiments (see Section 2.3). The whole inversion process was subdivided into three successive steps:

(i) estimating parameters $\alpha_S$ and $\gamma_S$ for GBL1;

(ii) estimating parameters $\alpha_S$ and $\gamma_S$ for GBL2;

(iii) estimating the interface depth.

This step-by-step process was made possible by the assumption of the power-law dependence of $V_S$ with overburden pressure within each model layer and by the model preparation and data acquisition procedures: the process ensured an accurate estimation of the unknown parameters and it enabled to limit solution non-uniqueness.

### 4.1 Retrieving parameters $\alpha_S$ and $\gamma_S$ for the lower layer

As described in the experimental set-up and data acquisition section (Section 2.4), a seismogram was acquired on top of GBL1 before digging the slope and adding the GB for GBL2 (Line 1, Fig. 6a). A dispersion curve was extracted by picking the energy maxima on the f-k spectrum derived from the seismogram (black crosses on Fig. 9). This dispersion curve was inverted for $\alpha_S$ and $\gamma_S$ of GBL1. Being the physical model before the addition of GBL2 a single layer model, we used the dispersion curve from Line 1 in pseudo-depth versus pseudo-$V_S$ domain for the estimation of parameters $\alpha_S$ and $\gamma_S$ controlling the S-wave vertical velocity distribution in GBL1, as illustrated in the analysis of synthetic data (Section 2.2). We adopted a grid search approach (on the fundamental mode only): its results are reported in Fig. 11, where the misfit evaluated for every couple of $\alpha_S$ and $\gamma_S$ is displayed on a logarithmic scale (Fig. 11a). The couple yielding the minimum misfit value is $\alpha_S = 0.345$ and $\gamma_S = 6.5$: Fig. 11(b) reports the good agreement between the experimental points and the power-law S-wave velocity trend determined by these values of $\alpha_S$ and $\gamma_S$. On Fig. 11(a), many ($\alpha_S$, $\gamma_S$) couples exhibit low misfit values around the minimum. Actually, two different ($\alpha_S$, $\gamma_S$) couples can produce velocity profiles very close to each other (even if the corresponding power laws only theoretically cross at one depth below zero). Associated dispersion curves are consequently very close to each other too.
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Estimating $\alpha_S$ and $\gamma_S$ for GBL1. (a) Misfit surface in the $\alpha_S-\gamma_S$ domain for the inversion of the dispersion curve from Line 1: the white cross marks the position of the best fitting point; (b) comparison between the experimental dispersion curve represented in terms of pseudo-depth versus pseudo-$V_S$ and the best fitting $V_S$ profile.</figcaption></figure>

This generates equivalent misfit values according to the dispersion estimations errors, as previously noted by Bodet *et al.* (2010).

## 4.2 Retrieving parameters $\alpha_S$ and $\gamma_S$ for the upper layer

The second step of the inversion process involved the estimation of the couple of $\alpha_S$ and $\gamma_S$ controlling the power-law behaviour of the S-wave velocity in GBL2. Differently from what had been done for GBL1, values for $\alpha_S$ and $\gamma_S$ were determined for every dispersion curve extracted from the seismograms acquired on the top of the upper layer (Fig. 6b). The lower layer, as described in the experimental set-up section, was arranged into the box with special cares (i.e. GB were sieved and compacted every centimetre) and, at the end of the deposition process, a greater compaction was achieved by vibrating the wooden box. Such deposition process and compaction are quite likely to ensure material homogeneity and constant values of $\alpha_S$ and $\gamma_S$. On the contrary, the upper layer GB were simply poured, so that slight lateral heterogeneities and different S-wave velocity behaviours (and therefore $\alpha_S$ and $\gamma_S$ values) are expected over the model. Therefore we inverted for $\alpha_S$ and $\gamma_S$ of GBL2 all

the dispersion curves we derived from the seismic data acquired on top of GBL2 (Figs 9 and 10), each time retrieving a different couple of values: however, the different estimates for $\alpha_{SGBL2}$ and $\gamma_{SGBL2}$ show good consistency.

As already shown in Section 2.2, in case of a two-layer model, the experimental dispersion curve fundamental mode represented in terms of pseudo-depth versus pseudo-$V_S$ follows the exponential trend described by eq. (1) down to the interface depth (Fig. 1d). The basic idea to estimate $\alpha_{SGBL2}$ and $\gamma_{SGBL2}$ was then to use the experimental curve fundamental mode points whose pseudo-depths is smaller than the interface depth between GBL1 and GBL2. These experimental curve points in fact correspond to surface wave frequency components which travel within subsurface portions which are not thick enough to meet the influence of GBL1: therefore these points would have the same phase velocities of points with the same frequency from a dispersion curve extracted from data acquired over a homogeneous model constituted by GBL2 only. However, the interface depth between GBL1 and GBL2 is not precisely known, as GBL2 smaller GB might have leaked down among larger GBL1 GB (thus turning the interface into a transition zone): moreover,
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Retrieving parameters $\alpha_S$ and $\gamma_S$ for GBL2 for the dispersion curve at $x = 0.462$ m. (a) Experimental dispersion curve represented in pseudo-$V_S$ versus pseudo-depth: the horizontal dashed line marks the border between the curve points relevant to GBL2 only (black dots) and the points belonging to the portion of the curve influenced by the presence of GBL1 (grey dots). The grey line indicates the theoretical $V_S$ profile obtained from the chosen value of $\alpha_S$ and $\gamma_S$. Note that to estimate $\alpha_{SGBL2}$ and $\gamma_{SGBL2}$ only black points were used. (b) Misfit as a function of pseudo-depth: every dot indicates the misfit between the experimental curve points whose pseudo-depths are equal or less than the dot pseudo-depth and the theoretical $V_S$ profile computed by inverting the considered experimental points for $\alpha_S$ and $\gamma_S$. (c) Estimated values of $\alpha_S$ versus pseudo-depth. Every dot represents the value of $\alpha_S$ obtained by inverting the experimental curve points with pseudo-depth equal or less than the pseudo-depth of the dot. The black circle marks the final result for $\alpha_S$ for the considered dispersion curve. (d) Same as (c) but for $\gamma_S$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Estimated values of $\alpha_S$ and $\gamma_S$ for GBL2. Values are displayed at the $x$ coordinate of the reference point of the relevant dispersion curve.</figcaption></figure>

interface depth is also one of the parameters to be retrieved through inversion. As it is not possible to identify a priori the dispersion curve points relevant to GBL2 only, to estimate $\alpha_{SGBL2}$ and $\gamma_{SGBL2}$ we adopted the following approach for each dispersion curve:

(i) Fundamental mode dispersion curve is represented in terms of pseudo-depth versus pseudo-$V_S$: the experimental points are sorted with increasing values of pseudo-depth ($z$) (Fig. 12a).

(ii) A ‘for’ loop over the $n$ points of the experimental curve is implemented. For every $i$th iteration (with $i$ ranging from 2 to $n$) only the first $i$ points are considered (i.e. the points up to the $i$th $z$-value): parameters $\alpha_{Si}$ and $\gamma_{Si}$ are computed with a least-squares approach (see eqs 2 and 3), minimizing the misfit $m_i$ between

$v_{S,theor} = \gamma_{Si}[\rho g z(1 : i)]^{\alpha_{Si}}$ and the experimental pseudo-$V_S$ values $v_{S,exp}(1 : i)$. At the end of every iteration, $\alpha_{Si}$, $\gamma_{Si}$ and the relevant misfit $m_i$ are represented at depth $z_i$ (see Figs 12b–d).

(iii) Eventually, the couple of $\alpha_S$ and $\gamma_S$ describing the $V_S$ behaviour of GBL2 is chosen. Up to a certain depth, in fact, the misfit is low and almost constant: beyond a certain depth threshold, it suddenly increases, because eq. (1) is no more able to model the behaviour of S-wave velocity with depth because of the influence of the deeper layer (in Fig. 12b this threshold is met at a depth of around 35 mm). Therefore, we chose the $\alpha_S$ and $\gamma_S$ values relevant to the greatest depth that still provides low misfit. Moreover, it has been noticed that, not surprisingly, this depth threshold roughly coincides with the expected interface depth between GBL2 and GBL1.
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Estimation of the interface depth for the experimental dispersion curve at $x = 0.462$ m. (a) Comparison between the experimental curve and the theoretical curves corresponding to the 10 best-fitting $V_S$ profiles (colour scale ranges from yellow = best profile to blue = 10th best profile). (b) Best fitting $V_S$ profiles.</figcaption></figure>

Fig. 13 represents the couples of $\alpha_S$ and $\gamma_S$ for GBL2 that were retrieved from all dispersion curves. $\alpha_S$ and $\gamma_S$ values are located at the corresponding dispersion curve position. The estimated $\alpha_S$ and $\gamma_S$ values (lying within 0.31–0.35 and 4.3–6, respectively) do not have a large variability and they are in agreement with the results previously achieved by Jacob *et al.* (2008) and Bodet *et al.* (2010, 2011) on similar physical models. Their slight variability is also satisfying since the different ($\alpha_S, \gamma_S$) couples lead to close $V_S$ values. Significantly, the power-law parameters retrieved for the upper layer are generally lower than the ones estimated for the lower layer ($\alpha_{SGBL1} = 0.345$ and $\gamma_{SGBL1} = 6.5$) and, in any case, they yield slower S-wave velocity profiles. In particular, $\alpha_{SGBL1}$ roughly coincides with the upper limit of the values for $\alpha_{SGBL2}$ and $\gamma_{SGBL1}$ is greater than any value retrieved for $\gamma_{SGBL2}$: moreover, relatively high values of $\alpha_{SGBL2}$ are associated to relatively low values for $\gamma_{SGBL2}$ and vice versa, so that all retrieved couples of $\alpha_{SGBL2}$ and $\gamma_{SGBL2}$ ensure lower $V_S$ values with respect to the ones determined by $\alpha_{SGBL1}$ and $\gamma_{SGBL1}$. The GBL1 layer, in fact, was vibrated in order to get a good compaction degree, and greater compaction implies a stiffer medium and therefore higher S-wave velocities.

### 4.3 Interface depth estimation

In the previous two steps we have characterized GBL1 and GBL2 in terms of S-wave velocity propagation by estimating, respectively $\alpha_{SGBL1}$, $\gamma_{SGBL1}$ and $\alpha_{SGBL2}$, $\gamma_{SGBL2}$ which express the dependence of $V_S$ on overburden pressure (see eq. 1). $\alpha_{SGBL1}$, $\gamma_{SGBL1}$ and $\alpha_{SGBL2}$, $\gamma_{SGBL2}$ were retrieved separately but can actually be combined to define the $V_S$ profile in our two-layer model: once reconstructed the overburden pressure $p$ as a function of depth, $\alpha_{SGBL1}$ and $\gamma_{SGBL1}$ define the trend of $V_S$ in GBL1 while $\alpha_{SGBL2}$ and $\gamma_{SGBL2}$ define the behaviour of $V_S$ within GBL2. However, to reconstruct the overburden pressure, it is necessary to estimate the depth of the interface between GBL1 and GBL2, being $\rho_{GBL1}$ and $\rho_{GBL2}$ known (see Section 2.3). To invert the experimental dispersion curves for the interface depth a grid search approach was chosen, made up of the following steps and individually applied to each dispersion curve:

(i) The interface depth is moved by 1 mm steps within the range 3–90 mm.

(ii) The S-wave velocity profile determined by parameters $\alpha_S$, $\gamma_S$ and $\rho$ for the upper and lower layer (considering the total load related
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Inversion results. 1-D $V_s$ profiles obtained by combining the estimated values of $\alpha_s$ and $\gamma_s$ for the two layers and the interface depth. $V_s$ profiles are located at the maxima of the Gaussian windows used to extract the corresponding dispersion curves.</figcaption></figure>

to both upper and lower layer portion) and by the current interface depth is discretized into 1-mm-thick layers, so that a theoretical dispersion curve can be computed through the Haskell and Thomson approach (Thomson 1950; Haskell 1953 and 1964).

(iii) The misfit is evaluated by considering both the P-SV fundamental and first higher mode of the experimental dispersion curve, to fully exploit the information contained in the experimental data.

(iv) The results are plotted according to the misfit and the interface depth yielding the minimum misfit is chosen.

Fig. 14(a) represents the experimental dispersion curve located at $x = 0.462$ m (Fig. 10) compared with the theoretical dispersion curves corresponding to the 10 best fitting interface depths while the corresponding $V_s$ profiles are depicted in Fig. 14(b): note that the estimated interface depths are quite close to 35 mm, which is the expected value for the interface depth at that location. Fig. 15 portrays the 1-D $V_s$ profiles that were obtained by combining the values of $\alpha_s$ and $\gamma_s$ for the two layers and of the interface depths that were retrieved for every inverted dispersion curves: as far as the interface depth estimation is concerned, the agreement with the expected values is good, as the maximum discrepancy between estimated and expected interface does not exceed 1 cm. Indeed, the global result is satisfactory, as the lateral variation is correctly retrieved and the S-wave velocity behaviour is coherent among all profiles.

## 5 CONCLUSIONS

This work shows the feasibility of building sophisticated physical models using micrometric GB and adopting different granulometries. We were able to create a stack of two GB's layers, presenting non-linear velocity gradients and different stiffness. Thanks to a simple but thorough step by step deposition procedure, these two layers were separated by a dipping interface. This provided a physical model, perfectly controlled in terms of geometry and very

well constrained in terms of elastic parameters thanks to previous theoretical and experimental results (Tournat & Gusev 2010). Based on an experimental set-up previously developed by Jacob *et al.* (2008) and improved by Bodet *et al.* (2010, 2011), we performed a small-scale seismic survey at the surface of the physical models using a laser-Doppler vibrometer and a mechanical source. We applied the algorithm by Bergamo *et al.* (2012) to extract several local surface wave dispersion curves from the acquired seismograms. We retrieved the S-wave velocity profiles along the recording line by means of a step by step inversion procedure. The 2-D structure of the analogue model was satisfactorily reconstructed, and the estimated elastic properties of the probed medium are in agreement with the results of previous theoretical and experimental studies (Aleshin *et al.* 2007; Gusev & Tournat 2008; Jacob *et al.* 2008).

Our results confirmed the ability of the algorithm by Bergamo *et al.* (2012) in extracting a dense set of reliable local dispersion curves whose inversion allow the model lateral variations to be reliably reconstructed. The *a priori* knowledge of a power-law variation of seismic properties with depth and the existence of an interface have contributed to obtain accurate results. A further step of the research work would be to attempt an inversion process inverting for all parameters at the same time, in order to evaluate the effect of equivalence phenomena and to test the effectiveness and reliability of a more ‘blind’ inversion procedure. Joint analysis of P-wave first arrival times and surface wave dispersion data, or multimodal dispersion inversion for P- and S-wave velocities could be addressed as well, as recently suggested by Boiero *et al.* (2009, 2011, 2013) and Maraschini *et al.* (2010).

The present study additionally confirmed the great potential of granular materials in terms of physical models construction and choice of parameters for the physical modelling of mechanical-wave propagation in various application fields (e.g. acoustics, exploration seismic), to address a wide range of theoretical and practical issues (e.g. heterogeneous media, complex structures, pore fluids, etc.). Krawczyk *et al.* (2013), for instance, recently pointed out that
Bodet *et al.* (2010) were only able to provide a local estimation
analogue models elastic properties. This study shows that thanks to
advanced processing approaches (Bergamo *et al.* 2012), the analy-
sis of guided and surface waves should make it possible to provide
structural images of the models. The experimental set-up presented
here should eventually provide an interesting tool for the imaging
of dry analogue models in geology, as suggested by Bodet *et al.*
(2010).

ACKNOWLEDGEMENTS

The authors acknowledge financial support from the Italian Ministry of Education, University and Research for funding the Ph.D. grant of Paolo Bergamo, and from Région Pays de la Loire, France.

REFERENCES

Abbiss, C.P., 1981. Shear wave measurements of the elasticity of the ground, *Geotechnique*, **31**, 91–104.

Aleshin, V., Gusev, V. & Tournat, V., 2007. Acoustic modes propagating along the free surface of granular media, *J. acoust. Soc. Am.*, **121**, 2600–2611.

Bergamo, P., Boiero, D. & Socco, L.V., 2012. Retrieving 2D structures from surface wave data by means of a space-varying spatial windowing, *Geophysics*, **77**, EN39–EN51.

Blum, T.E., van Wijk, K., Snieder, R. & Willis, M.E., 2011. Theory and laboratory experiments of elastic wave scattering by dry planar fractures, *J. geophys. Res.* **116**, B08218, doi:10.1029/2011JB008295.

Bodet, L., Abraham, O. & Clorennec, D., 2009. Near-offsets effects on Rayleigh-wave dispersion measurements: physical modelling, *J. appl. Geophys.*, **68**, 5–103.

Bodet, L., Dhemaied, A., Mourgues, R., Tournat, V. & Rejiba, F., 2011. Laser-Doppler acoustic probing of granular media with in-depth property gradient and varying pore pressures, in *Proceedings of the International Congress on Ultrasonics (ICU)*, University of Gdansk, 5–8 September, Poland.

Bodet, L., Jacob, X., Tournat, V., Mourgues, R. & Gusev, V., 2010. Elasticity profile of an unconsolidated granular medium inferred from guided waves: toward acoustic monitoring of analogue models, *Tectonophysics*, **476**, 99–104.

Bodet, L., van Wijk, K., Bitri, A., Abraham, O., Côte, Ph., Grandjean, G. & Leparoux, D., 2005. Surface-wave inversion limitations from laser-Doppler physical modeling, *J. Environ. Eng. Geophys.*, **10**, 151–162.

Bohlen, T., Kugler, S., Klein, G. & Theilen, F., 2004. 1.5-D inversion of lateral variation of Scholte wave dispersion, *Geophysics* **69**, 330–344.

Boiero, D. & Socco, L.V., 2010. Retrieving lateral variations from surface wave dispersion curves analysis, *Geophys. Prospect.*, **58**, 977–976.

Boiero, D., Calzoni, C. & Socco, L.V., 2011. Joint inversion of surface-wave dispersion and P-wave refraction data for laterally varying layered models, in *Proceedings of the 73rd EAGE Conference*, 23–26 May, Vienna, Austria.

Boiero, D., Maraschini, M. & Socco, L.V., 2009. P and S wave velocity model retrieved by multimodal surface-wave analysis, in *Proceedings of the 71st EAGE Conference*, 8–11 June, Amsterdam, The Netherlands.

Boiero, D., Wiarda, E. & Vermeer, P., 2013. Surface- and guided-wave inversion for near-surface modeling in land and shallow marine seismic data, *Leading Edge*, **32**(6) 638–646

Bretaudeau, F., Leparoux, D., Durand, O. & Abraham, O., 2011. Small-scale modeling of onshore seismic experiment: a tool to validate numerical modeling and seismic imaging methods, *Geophysics*, **76**, T101–T112.

Buddensiek, M.L., 2009. Seismic imaging of sandbox models, *PhD thesis*, Freien Universität Berlin, Berlin, Germany.

Campman, X.H., van Wijk, K., Scales, J.A. & Herman, G.C., 2005. Imaging and suppressing near-receiver scattered surface waves, *Geophysics*, **70**, V21–V29.

Campman, X.H., van Wijk, K., Riyanti, C.D., Scales, J.A. & Herman, G.C., 2004. Imaging scattered seismic surface waves, *Near Surf. Geophys.*, **2** 223–230.

De Cacqueray, B., Roux, Ph., Campillo, M., Catheline, S. & Boue, P., 2011. Elastic-wave identification and extraction through array processing: an experimental investigation at the laboratory scale, *J. appl. Geophys.*, **74**, 81–88.

Dewangan, P., Tsvankin, I., Batzle, M., van Wijk, K. & Haney, M., 2006. PS-wave moveout inversion for tilted TI media: a physical-modeling study, *Geophysics*, **71**, D135–D143.

Gabriels, P., Snieder, R. & Nolet, G., 1987. In situ measurements of shear-wave velocity in sediments with higher mode Rayleigh waves, *Geophys. Prospect.*, **35**, 187–196.

Gassmann, F., 1951. Elastic waves through a packing of spheres, *Geophysics*, **16**, 673–685.

Grandjean, G. & Bitri, A., 2006. 2M-SASW: multifold multichannel seismic inversion of local dispersion of Rayleigh waves in laterally heterogeneous subsurfaces: application to the Super-Sauze earthflow, France, Near Surf. Geophys., **4**, 367–375.

Graveleau, F., 2008. Interactions tectonique, érosion, sédimentation dans les avant-pays de chaînes : modélisation analogique et étude des piémonts de l'Est du Tian Shan (Asie centrale). *PhD thesis*, Université Montpellier II, Montpellier, France.

Gusev, V. & Tournat, V., 2008. How acoustic waves are guided in buried subsurface channels in unconsolidated granular media, *Phys. Rev., E* **78**, 036602.

Gusev, V., Aleshin, V. & Tournat, V., 2006. Acoustic waves in an elastic channel near the free surface of granular media, *Phys. Rev. Lett.*, **96**, 214301.

Haskell, N., 1953. The dispersion of surface waves on multilayered media, *Bull. seism. Soc. Am.*, **43**, 17–34.

Haskell, N., 1964. Radiation pattern of surface waves from point sources in a multilayered medium, *Bull. seism. Soc. Am.*, **54**, 377–393.

Hayashi, K. & Nishizawa, O., 2001. Laboratory studies of surface waves using a laser-Doppler vibrometer, in *Proceedings of the 5th Int. Symp. on Recent Advances in Exploration Geophysics*, Kyoto, pp. 1–8.

Hayashi, K. & Suzuki, H., 2004. CMP cross-correlation analysis of multi-channel surface-wave data, *Explor. Geophys.*, **35**, 7–13.

Hentschel, M.L. & Page, N.W., 2007. Elastic properties of powders during compaction. Part 1: pseudo-isotropic moduli, *J. Mater. Sci.*, **42**, 1261–1268.

Jacob, X., Aleshin, V., Tournat, V., Leclaire, P., Lauriks, W. & Gusev, V.E., 2008. Acoustic probing of the jamming transition in an unconsolidated granular medium, *Phys. Rev. Lett.*, **100**, 158003.

Kaslilar, A., 2007. Inverse scattering of surface waves: imaging of near-surface heterogeneity, *Geophys. J. Int.*, **171**, 352–367.

Krawczyk, C.M., Buddensiek, M.-L., Oncken, O. & Kukowski, N., 2013. Seismic imaging of sandbox experiments-laboratory hardware setup and first reflection seismic sections, *Solid Earth*, **4**, 93–104.

Luo, Y., Xia, J., Miller, R.D., Xu, Y., Liu, J. & Liu, Q., 2009. Rayleigh-wave mode separation by high-resolution linear Radon transform, *Geophys. J. Int.*, **179**, 254–264.

Makse, H.A., Gland, N., Johnson, D.L. & Schwartz, L.M., 2004. Granular packings—nonlinear elasticity, sound propagation, and collective relaxation dynamics, *Phys. Rev. E*, **70**, 061302.

Maraschini, M., Ernst, F., Foti, S. & Socco, L.V., 2010. A new misfit function for multimodal inversion of surface waves, *Geophysics*, **75**(4), G31–G43.

Nagai, K., O'Neill, A., Sanada, Y. & Ashida, Y., 2005. Genetic algorithm inversion of Rayleigh wave dispersion from CMPCC gathers over a shallow fault model, *J. Environ. Eng. Geophys.*, **10**, 275–286.

Neduca, B., 2007. Stacking of surface waves, *Geophysics*, **72**, V51–V58.

Nishizawa, O. & Kitagawa, G., 2007. An experimental study of phase angle fluctuation in seismic waves in random heterogeneous media: time-series analysis based on multivariate AR model, *Geophys. J. Int.*, **169**, 149–160.

Nishizawa, O., Satoh, T., Lei, X. & Kuwahara, Y., 1997. Laboratory studies of seismic wave propagation in inhomogeneous media using a laser Doppler vibrometer, *Bull. seism. Soc. Am.*, **87**, 809–823.
Ruiz, A. & Nagy, P.B., 2004. Laser-ultrasonic surface wave dispersion measurements on surface-treated metals, *Ultrasonics*, **42**, 665-669.

Scales, J.A. & Malcolm, A.E., 2003. Laser characterization of ultrasonic wave propagation in random media, *Phys. Rev. E* **67**, 046618.

Scales, J.A. & van Wijk, K., 1999. Multiple scattering attenuation and anisotropy of ultrasonic surface waves, *Appl. Phys. Lett.*, **74**, 3899-3901.

Sherlock, D.H., 1999. Seismic imaging of sandbox models, *PhD thesis*, Curtin University of Technology, Perth, Australia.

Sherlock, D.H. & Evans, B.J., 2001. The development of seismic reflection sandbox modelling, *AAPG Bull.*, **85**, 1645-1659.

Sivaji, C., Nishizawa, O., Kitagawa, G. & Fukushima, Y., 2002. A physical-model study of the statistics of seismic waveform fluctuations in random heterogeneous media, *Geophys. J. Int.*, **148**, 575-595.

Socco, L.V. & Strobbia, C., 2004. Surface-wave method for near surface characterization: a tutorial, *Near Surf. Geophys.*, **2**, 165-185.

Socco, L.V., Foti, S. & Boiero, D., 2010. Surface-wave analysis for building near-surface velocity models—established approaches and new perspectives, *Geophysics*, **75**, A83–A102.

Socco, L.V., Boiero, D., Foti, S. & Wisén, R., 2009. Laterally constrained inversion of ground roll from seismic reflection records, *Geophysics*, **74**, G35–G45.

Spetzler, J., Sivaji, C., Nishizawa, O. & Fukushima, Y., 2002. A test of ray theory and scattering theory based on a laboratory experiment using

ultrasonic waves and numerical simulation by finite-difference method, *Geophys. J. Int.*, **148**, 165-178.

Strobbia, C., Laake, A., Vermeer, P. & Glushchenko, A., 2011. Surface waves: use them then lose them. Surface-wave analysis, inversion and attenuation in land reflection seismic surveying, *Near Surf. Geophys.*, **9**, 503–514.

Thomson, W.T., 1950. Transmission of elastic waves through a stratified solid medium, *J. appl. Phys.*, **21**, 89-93.

Tournat, V. & Gusev, V., 2010. Acoustics of uncosolidated “model” granular media: an overview of recent results and several open problems, *Acta Acust. Acust.*, **96**, 208–224.

Valverde, J.M. & Castellanos, A., 2006. Random loose packing of cohesive granular materials, *Europhys. Lett.*, **75**, 985-991.

van Wijk, K. & Levshin, A.L., 2004. Surface wave dispersion from small vertical scatterers, *Geophys. Res. Lett.*, **31**, doi:10.1029/2004GL021007.

van Wijk, K., Komatitsch, D., Scales, J.A. & Tromp, J., 2004. Analysis of strong scattering at the micro-scale, *J. acoust. Soc. Am.*, **115**, 1006–1011.

Wu, T.-T. & Liu, Y.H., 1999. Inverse determinations of thickness and elastic properties of a bonding layer using laser-generated surface waves, *Ultrasonics*, **37**, 23-30.

Yilmaz, O., 1988. *Seismic Data Processing*, 2nd edn, Vol. 2, pp. 82–153, Society of Exploration Geophysicists.