Article Highlights

  • Surface wave tomography is an effective technique to identify lateral velocity variations, also for environmental application.

  • The proposed inverse scheme for 2-D active acquisition geometries could be used in place of any multichannel technique.

  • A new procedure for the choice of the roughness damping parameter is suggested, that considers only physically consistent phase variations.

1 Introduction

Detecting shallow seismic velocity variations in space is of great interest for many fields of application, such as environmental geophysics, engineering and mining exploration. Surface wave methods allow retrieving near-surface shear wave velocity models through dispersion analysis and depth inversion. Investigation depth mainly depends on the sampled frequencies. Passive techniques based on microtremors recordings such as the beamforming method (Asten and Henstridge 1984), the spatial autocorrelation method (SPAC; Aki 1957), the extended spatial autocorrelation method (ESPAC; Ling and Okada 1993) and refraction microtremor (Re-Mi; Louie 2001; Strobbia et al. 2015) measure the very low frequency surface wave content of both natural and anthropic noise sources, and allow the investigation of depths from several hundred meters to a few kilometres (Okada and Suto 2003). On the contrary controlled-source methods generally penetrate shallower depths, from tens to hundreds meters depending on the available acquisition equipment and on the local site conditions. Sampled frequencies also have a huge impact on the achievable lateral resolution: for the same propagation length, shorter wavelengths compute more cycles compared to higher wavelengths, which results in a finer spatial sampling and a better recovery of the lateral variability, but also in a much rapid attenuation. For this reason, many authors propose to use frequency-dependent offset ranges when performing surface wave analysis with active data (Neducza 2007; Strobbia et al. 2011; Vignoli et al. 2016). Besides the frequency content, the choice of the most appropriate surface wave analysis technique could make a significant difference when trying to resolve spatial heterogeneities. This choice most often relies on existing information about the site to be investigated, either coming from geological/geotechnical data or from preliminary seismic analysis, such as the expected presence of sharp variations (e.g., faults, voids, lithological contacts, etc.). In fact, surface wave propagation in a strongly laterally heterogeneous medium leads to relevant changes in the phase behaviour. Boaga et al. (2014) show the surface-wave modes energy distribution for different subsoil conditions, including a two-layer system with an abrupt lateral discontinuity. If this geological contrast is ignored, the fundamental and higher order modes of the left and right parts of the section could be misinterpreted, bringing to a wrong subsurface reconstruction. Moreover, specific 2-D geometries such as sediment-filled valleys could bring to the generation of higher modes Rayleigh waves from the valley edges due to resonant SV-waves (Bard and Bouchon 1985). This effect is particularly relevant when performing 1-D microtremor analysis (Claprood et al. 2011). Most common surface wave analysis techniques for active seismic data fall into two different categories: (i) techniques that aim to retrieve smooth lateral variations; (ii) techniques for detecting sharp velocity variations. A brief overview of the techniques belonging to the two groups is reported below.

When lateral variations are sufficiently smooth, the path-average (PAVA) approximation (Woodhouse and Dziewonski 1984), for which a 2-D structure is approximated by its average 1-D profile, may be used locally. Many authors attempted to perform pseudo 2-D analysis by integrating the widespread multichannel analysis of surface waves (MASW; Park et al. 1999) with local 1-D approximations. The most common approach consists in performing spectral analysis on moving windows along the acquisition array, attributing the estimated dispersion curve to the centre of each window (e.g., Bohlen et al. 2004). Socco et al. (2009) and Boiero and Socco (2010) associated the same windowing procedure to a laterally constrained depth inversion, in order to force the lateral continuity of the final shear velocity section. Although windowing in the space domain is effective in recovering smooth lateral variations, it also produces a spectral degradation: when choosing the size of the sliding window, a trade-off between spectral and lateral resolution needs to be found. Luo et al. (2008) overcame this problem by using a high-resolution linear Radon transform (LRT), producing a more defined spectral image compared to standard spectral transformation such as fk or \(\omega\)p. This method is also effective in multi-modal detection and separation. Hayashi and Suzuki (2004) proposed an alternative approach to spatial windowing, where the multi-channel analysis is applied to CMP cross-correlation gathers, making the retrieved dispersion curves representative of local site conditions. Grandjean and Bitri (2006) described a similar procedure applied to multi-fold data, where data-redundancy is used to improve the S/N ratio on the computed spectra. A valid alternative to MASW is the multi-offset phase analysis (MOPA; Strobbia and Foti 2006), which performs a frequency-dependent linear regression of the phase versus offset to derive the dispersion relation. This method is robust but has the drawback of being applicable only with single-mode propagation. Recently, an extension of this method for the two-modes case has been proposed by Barone et al. (2020). The MOPA can be applied on sliding windows in order to catch smooth lateral variations; the great advantage of this technique consists in minimizing the number of selected channels without compromising the lateral resolution (Vignoli et al. 2011). However, in order to assure a robust analysis at low frequencies and maximize lateral resolution at higher frequencies, a moving window with a frequency-dependent length should be used (Vignoli et al. 2016).

In case of abrupt lateral variations such as faults, presence of buried structures, lithologic contacts etc., the 1-D PAVA approximation cannot be adopted and a different approach should be used. A common strategy requires to first locate the discontinuities, then divide the medium into sub portions of homogeneous characteristics and analyse them separately. Abrupt lateral variations cause the surface wave energy to be partly back-reflected and back-scattered at the encountered interface. For this reason, many techniques to detect sharp lateral changes are based on the analysis of seismic attenuation (Colombero et al. 2019). Nasseri-Moghaddam et al. (2005) proposed a method based on the computation of the energy-distance (ED) and logarithmic decrement (LD) parameters to detect shallow cavities. The same analysis has been improved by Bergamo and Socco (2014) and extended to multi-fold data. Other techniques are based on the separation between incident and back-scattered wavefield: Xia et al. (2007) used Rayleigh wave diffractions to detect near-surface features such as voids or faults, while Kaslilar (2007) and Liu et al. (2017) proposed to invert the scattered surface waves for imaging near-surface heterogeneities. Sudden lateral variations also produce visible changes in the surface-wave phase slope: Vignoli and Cassiani (2010) further developed the MOPA technique to search for sudden lateral velocity changes (the so-called knee points of the phase versus offset distribution); Vignoli et al. (2011) extended this procedure to multi-shot data, where the application of an fk filter and the statistical analysis of redundant data allow a much better location of the discontinuities.

The availability of a-priori information about the site of investigation, however, is not a common situation. Moreover, in the same area both smooth and sharp lateral variations could coexist, making the choice of the most suitable surface wave analysis technique challenging. In the last decade surface wave tomography (SWT), a widespread technique in global seismology, has become very attractive also for local scale applications, since it allows recovering both smooth and sharp lateral variations. Because of its simplicity, the most widely used procedure to perform SWT consists in inverting all measured travel-times between source-receiver (or receiver–receiver) couples at one specific frequency to produce a phase velocity spatial distribution. This operation is repeated for all frequencies of analysis and local dispersion curves are extracted. Then, individual dispersion curves are depth inverted and the obtained 1-D profiles are assembled to compose a 2-D or 3-D shear-wave velocity model. This is what we call the “two-step approach”. Some authors attempted to make this process less fragmented. Mohammadi et al. (2020) formulated the simultaneous tomography of all periods (STP), where traveltimes of all ray paths associated with all periods are inverted simultaneously to retrieve the phase/group velocity maps. This method takes into account the continuity of phase velocity measurements in frequency other then in space. Boschi and Ekström (2002) proposed the “one-step approach”, that consists in the direct inversion of traveltimes to retrieve a 3-D shear wave velocity distribution, without intermediate reconstruction of phase velocity maps. This method makes use of sensitivity kernels to link the Earth model parameters at depth with the phase velocity at each frequency of analysis. Other formulations of this method, based on the ray-theory (high frequency approximation), are shown in An et al. (2009), Boiero (2009) and Fang et al. (2015), while Zhou et al. (2005) use the Born (finite frequency) approximation. In the latter case, the real advantage consists in imaging lateral heterogeneities that are smaller than the characteristic surface wave wavelength.

First SWT applications used surface wave signal from known earthquakes recorded by sparse seismological stations to map global to regional phase or group velocity distributions at different periods (Trampert and Woodhouse 1996; Ekstöm et al. 1997; Ritzwoller and Levshin 1998; Boschi and Dziewonski 1999; Barmin et al. 2001). Lately, ambient noise tomography (ANT) has been developed, which uses cross-correlations of low-frequency seismic noise as input signals. While ANT has been extensively used at a regional scale (Shapiro et al. 2005; Lin et al. 2008, 2009; Kästle et al. 2018), several studies using local dense seismic networks have shown its great potential at much a smaller scale. Picozzi et al. (2009) and Pilz et al. (2012) performed a small 3-D tomography test to image a known structure at the Nauen test site in Germany, aiming to verify the suitability of this technique for engineering applications; Mordret et al. (2013) adopted ANT at the Valhall oil field cross-correlating signals recorded by the 2320 four-component sensors of the ‘Life of Field Seismic’ network, obtaining Scholte wave group velocity maps of impressive resolution below 2 Hz; Lin et al. (2013) applied eikonal tomography to more than 13.5 million cross-correlations of ambient noise recorded with a dense 3-D seismic array in Long Beach, California; Hollis et al. (2018) applied ANT to mining exploration to image old mine workings and mineralized geological bodies; Mordret et al. (2019) used four weeks long recordings from a spatially dense Nodal array to image the San Jacinto fault in southern California with unprecedented resolution. As for SWT, several attempts of using it with active seismic data instead of recorded earthquakes have been made, at very different scales. Haney and Douma (2012) applied group velocity tomography and phase velocity inversion to estimate static corrections for a dense 3-D active seismic survey at the Coronation Field (Canada). Gouédard et al. (2012) applied eikonal tomography to a high-resolution exploration dataset acquired in Oman. Socco et al. (2014) proposed the application of surface wave tomography to a 2-D active seismic dataset as an alternative to classical MASW, with the aim of imaging an Alpine fault in New Zealand. Rector et al. (2015) used surface wave tomography to image voids in an old mine working site in Nevada (USA). Duret et al. (2016) combined standard P-wave refraction and SWT to compute primary statics for a 3-D narrow-azimuth land dataset. Krohn and Routh (2017a) and Krohn and Routh (2017b) developed a tomography method that evaluates the complex slowness distribution by fitting the trace transfer function, i.e. the waveform change due to the surface waves propagation, and applied it to both 2-D and 3-D densely sampled active data, to predict and subtract the surface-wave signal from the recorded seismograms. Finally, Da Col et al. (2020) applied SWT to a mine in Eastern Finland, with the aim of imaging the rock formation carrying the mineralized ore body.

In this study we want to apply a linear inversion scheme derived from SWT to a 2-D small scale controlled source dataset acquired for environmental purposes. The seismic line crosses longitudinally an old landfill located 15 km far from Vienna, Austria, a site now completely integrated in an urban centre. The purpose of the study is to characterize the shape of the landfill body, with a specific focus on the urban waste deposit, which causes leachate and gas generation problems. To this aim, we set up a simplified inverse scheme adopting the assumption of straight rays between source and receiver, i.e. ignoring the second (not sampled) spatial dimension. Phase differences between adjacent receivers are used as observations for the least-squares minimization problem, which is solved for each individual frequency. Special attention has been put on the choice of the regularization constraints. Finally, local dispersion curves are depth inverted and the obtained profiles assembled to compose a quasi-2-D shear wave velocity model. The obtained results are in agreement with different electrical and electromagnetic geophysical measurements performed along the same transect.

2 Surface Wave Tomography

The tomographic problem for one single frequency can be expressed as a minimization problem in a least square sense. Observations (travel times between two points) are linked to the model parameters (phase velocities or slownesses) through the forward computation, which can be expressed as follows:

$$\begin{aligned} t=\int _{{\rm PATH}}{\frac{1}{v(\mathbf{x })}{\rm d}l}=\int _{{\rm PATH}}{p(\mathbf{x }){\rm d}l} \end{aligned}$$
(1)

where \(v(\mathbf{x })\) and \(p(\mathbf{x })\) are, respectively, the phase velocity and the phase slowness at the position \(\mathbf{x } = (x_1, x_2, x_3)\), and dl is the infinitesimal segment of the travel path over which the integration is performed. In a 3-dimensional space the problem defined in Eq. (1) is highly nonlinear, because the wave path depends on the phase velocity distribution in space, which is by definition unknown. For this reason, the forward problem is commonly reformulated in terms of traveltime perturbations:

$$\begin{aligned} \delta t=t-t_0 =\int _{{\rm PATH}}{p(\mathbf{x }){\rm d}l}-\int _{{\rm PATH}}{p_0(\mathbf{x }){\rm d}l} \end{aligned}$$
(2)

where \(p_0\) an \(t_0\) refer to the phase slownesses and the travel times computed using a reference model, which we assume being not too far from the real one. Under this condition, the integration path can be assumed unchanged for the reference and the real model, therefore:

$$\begin{aligned} \delta t = \int _{{\rm PATH}}{\delta p(\mathbf{x }){\rm d}l} \end{aligned}$$
(3)

We can express Eq. (3) in a matrix form:

$$\begin{aligned} \mathbf{A } \mathbf{m } = \mathbf{d } \end{aligned}$$
(4)

Where m is the model vector, containing the n unknown model parameters, d is the data vector, containing the m observations and A is the \(m \times n\) matrix whose elements, Aij, contain the length of the ray path in each jth model element for the ith observation.

The least-squares solution to Eq. (4) is:

$$\begin{aligned} \mathbf{m } = (\mathbf{A }^{\rm T} \mathbf{A })^{-1} \mathbf{A }^{\rm T} \mathbf{d } \end{aligned}$$
(5)

where \((\mathbf{A }^{\rm T} \mathbf{A })^{-1}\) is the pseudo-inverse of A. However, in most tomographic problems A is likely to be sparse: as a consequence, \(\mathbf{A }^{\rm T} \mathbf{A }\) is close to be singular (i.e. not invertible). The tomographic problem, as is, is ill-posed and a regularization is needed to stabilize the solution.

Many sorts of regularizations may be imposed, such as the “norm damping” and the “roughness damping” (also called “smoothing”). The first consists in minimizing the norm of the solution (vector m), the latter in minimizing the norm of its gradient. Boschi and Dziewonski (1999) demonstrated how norm damping could produce artefacts in less sampled areas, where the solution is forced to resemble the reference model. Roughness damping, whose only effect is to smooth the solution, is therefore preferable.

The regularized tomographic problem can be written as:

$$\begin{aligned} \left[ \begin{array}{c} \mathbf{A } \\ \alpha \mathbf{G } \end{array} \right] \mathbf{m } = \left[ \begin{array}{c} \mathbf{d } \\ 0 \end{array} \right] \end{aligned}$$
(6)

where G is the gradient operator and \(\alpha\) is an arbitrary constant determining the weight of the roughness damping condition over data fitting. Details about how to construct the G matrix in the simplified 2-D case will be given in the following section. Solution to Eq. (6) is given by:

$$\begin{aligned} \mathbf{m } = ( \mathbf{A }^{\rm T} \mathbf{A } + \alpha ^2 \mathbf{G }^{\rm T} \mathbf{G } )^{-1} \mathbf{A }^{\rm T} \mathbf{d } \end{aligned}$$
(7)

Equation (7) corresponds with the solution of an inverse problem with a first-order Tikhonov regularization (Tikhonov and Arsenin 1977), expressed as:

$$\begin{aligned} \left\| \mathbf{A } \mathbf{m } - \mathbf{d } \right\| + \lambda \left\| \mathbf{G } \mathbf{m } \right\| = min \end{aligned}$$
(8)

For this to be true, the Tikhonov parameter \(\lambda\) must equal \(\alpha ^2\).

In case data uncertainty is known, Eq. (7) can be replaced by its weighted formulation:

$$\begin{aligned} \mathbf{m } = ( \mathbf{A }^{\rm T} \mathbf{W } \mathbf{A } + \alpha ^2 \mathbf{G }^{\rm T} \mathbf{G } )^{-1} \mathbf{A }^{\rm T} \mathbf{W } \mathbf{d } \end{aligned}$$
(9)

where W is a \(m \times m\) diagonal matrix containing the inverse of data variances.

Equation (7) or Eq. (9) may be solved with a direct computation in case of a small scale problem. However, as the size of A matrix grows, inverting the term within parenthesis requires important computing power. Many alternative techniques, most of them making use of A matrix sparseness, are used to solved this kind of large scale minimization problems, such as Conjugate Gradient (Hestenes and Stiefel 1952; Inoue et al. 1990), LSQR (Paige and Saunders 1982) and Subspace Inversion Scheme (Kennett et al. 1988). Fully nonlinear methods such as Monte Carlo schemes (Bodin and Sambridge 2009) are also used, allowing a wider model space sampling, but they are normally highly computationally expensive.

Concerning the choice of the \(\alpha\) constant: this is generally done either by selecting the smoothest model that fits data within its noise range (“Chi-Square” criteria) or by finding the optimum compromise between data fitting and model smoothing (“L-shaped curve” criteria). While the first approach is rigorous, but requires a reliable estimate of data uncertainties, the second is qualitative. In the sequel we will propose a method that, in the specific case of an active 2-D acquisition, operates a frequency-dependent selection of the best regularization based on the physical properties of surface waves: for each frequency, the optimum \(\alpha\) constant is chosen such as only lateral variations that can be physically resolved are included in the analysis. Data uncertainties will be used to weight the solution.

3 The 2-D Active Case

All concepts described in the previous section are general and applicable to any kind of surface wave tomographic problem. Here we focus on the specific case of a 2-D active seismic acquisition, where only one spatial dimension is sampled (sources and receivers are aligned). 2-D active acquisitions usually present a very regular geometry, with fix inter-receiver and inter-source spacing (see Fig. 1a). As a consequence, data uniformly cover each portion of the line.

First of all, a basic assumption has to be made: we assume straight rays between each source-receiver couple. For this reason, the forward problem can be expressed in terms of absolute traveltimes (Eq. 1), and no reference model is needed. In particular, by multiplying both terms of Eq. (1) by \(2 \pi f\), with f being frequency, we obtain:

$$\begin{aligned} \varDelta \phi =\int _{{\rm PATH}}{-k(\mathbf{x }){\rm d}l} \end{aligned}$$
(10)

where \(\phi\) is the phase and k is the wavenumber. Equation (10) is a linear expression, valid for one single frequency and one mode of propagation. Such a definition of the forward problem can be translated into matrix notation (Eq. 4), with vector d containing the m observed phase differences between couples of receivers and vector m the n unknown wavenumber values along the acquisition profile. A convenient choice is to measure phase differences between adjacent receivers: by doing so, the length of all raypaths is constant and equal to the inter-receiver distance. As a consequence, matrix A has only one non-zero element per raw, which makes \(\mathbf{A }^{\rm T} \mathbf{A }\) diagonal. Moreover, due to the uniform data coverage over the investigated line, diagonal elements of \(\mathbf{A }^{\rm T} \mathbf{A }\) have similar values, therefore \((\mathbf{A }^{\rm T} \mathbf{A })^{-1}\) gives a numerically stable solution. For all reasons above, regularizing the problem is not necessary.

However, field data are always affected by a certain amount of noise, which propagates through the inverse problem to its solution: when multiple repetitions of the same shot are available, the variance of the observed phase differences can be computed and used to weight the least-squares solution (Eq. 9). This can significantly improve the final model reliability, but it does not guarantee its smoothness and for this reason regularizing the problem is often useful. Equation (6) includes the roughness damping regularization constrain within the forward calculation. For the 2-D case matrix G represents the first-order differential operator that is applied to the model vector m. G is a \(n-1 \times n\) matrix, with the ith and ith + 1 column elements of each ith rows being − 1 and 1. As for the choice of the \(\alpha\) parameter, this is a critical issue: with \(\alpha\) too low, the final solution is dominated by the influence of data noise, while with \(\alpha\) too high, small scale lateral variations, which are our main target, are averaged out together. A real balance between these two extreme situations must be found.

Due to the reduced size of the 2-D problem, an exact solution to Eq. (7) or Eq. (9) can be easily provided. Once vector m is found, phase velocities are derived from wavenumbers as:

$$\begin{aligned} c =\frac{2 \pi f}{k} \end{aligned}$$
(11)

It is easy to notice how the inverse problem described here detaches from the original meaning of tomography, both for the unique direction sampled and for the local character of the phase measures. For this reason, in the sequel we will refer to it as “tomography-derived” linear inversion of surface wave data.

As an example, a simple 2-D acquisition geometry has been simulated (Fig. 1a), with 14 source positions and 48 channels. Source spacing and receiver spacing are, respectively, 8 m and 2 m. Then, a horizontal phase velocity distribution for one frequency and one mode has been generated, with 47 constant velocity cells (each cell stays between two adjacent receivers). Synthetic phase differences for all adjacent receivers and all shots were computed. The inverse problem was solved for different levels of Gaussian noise and different values of \(\alpha\), both using the weighted (Eq. 9) and non-weighted (Eq. 7) formulations. Weights were computed as the inverse of the squared data errors. Results are resumed in Fig. 1b. The reader can observe how, provided a reliable estimate of errors, the weighted formulation gives a more stable result in presence of noise. The reader will also notice the smoothing effect given by the imposed regularization: while it assures a better lateral continuity it can overdamp small scale lateral variations. For this reason, the \(\alpha\) parameter must be carefully chosen.

Fig. 1
figure 1

a 2-D Acquisition geometry of the synthetic test. b Results of the tomography-derived linear inversion with different percentages of Gaussian noise in the data and different values of \(\alpha\): the black line represents the true phase velocity distribution, the red and blue lines show the result of the weighted and non-weighted inversions, respectively

4 Real Data Example

In this section we present a field data example where the importance of detecting lateral velocity variations is central. The “Heferlbach” landfill (Fig. 2) is located in Mannsworth, 15 km south of Vienna (Austria). Between 1965 and 1974 this site has been filled with both urban and construction waste, for a total estimated volume of 240.000 m\(^3\) (top soil included) and an average height of 3.7 m (Valtl 2005). Several samples collected in 2001 (Valtl 2005) and 2012 (Brandstätter et al. 2013) permitted to punctually reconstruct the waste distribution in depth. However, due to the scarce documentation about the filling period, the precise geometry of the landfill is not known. With the time the landfill underwent both subsidence and leachate/biogas generation problems, reason why it is still studied and monitored. Nowadays, the landfill area hosts green spaces, including a children playground and several cultivated fields, and is adjacent to many residential buildings and a sport field.

In 2019 we carried out a geophysical campaign in order to characterize the South-Est ending of the landfill: special attention was put on the identification of the lateral limit of the urban waste deposit. To this end, the same line was acquired both with electromagnetic induction (EMI), electrical resistivity tomography (ERT) and seismic methods: while the first two are commonly used in landfills, being sensitive to the presence of pollutants, the third method is rarely used in this kind of environment. However, a significant change in the mechanical response of the medium is expected between original sediments, urban waste and construction waste. Moreover, the concurrent use of different geophysical techniques permits to retrieve a more complete picture of the investigated medium, allowing an easier data interpretation. In this study we will mainly focus on the surface wave analysis of seismic data, and only briefly describe the results from EMI and ERT techniques for comparison and interpretation purposes.

Fig. 2
figure 2

On the left: Presumed approximate perimeter and area of the “Heferlbach” landfill, Mannsworth. The 2019 geophysical campaign interested the Southern-Eastern part of the landfill (dotted-line rectangle). On the right: Satellite image of the prospection area, with the 272 m long seismic acquisition line in yellow and the extension of the three acquisition profiles in cyan, green and blue. The approximate landfill margins are displayed in red

4.1 Seismic Measurements

Seismic data were acquired along a 272 m long line longitudinally crossing the South-East part of the landfill. Due to logistic limitations, the line was split into three different profiles: while the first and second profiles are next to each other, the second and third profiles partially overlap for a length of 12 m (Fig. 2). Seismic acquisition was performed independently for each profile, using the same acquisition geometry shown in the previous chapter (Fig. 1a). Records were acquired with 48 vertical geophones with a central frequency of 4.5 Hz. Each shot position has been energized four to five times with a 8 kg sledgehammer, in order to increase the S/N ratio.

Tomography-derived linear inversion and MOPA were applied to the 2-D seismic data collected from the three profiles, and their results were compared. Both techniques require single-mode propagation, therefore we first performed a spectral analysis to identify the Rayleigh wave fundamental mode and remove possible higher-order modes contamination. To this aim, all seismograms recorded at the same shot position were stacked together. Each stacked shot gather has been split into forward (positive offsets) and reverse (negative offsets) traces and transformed into the fk domain, where a filter was manually designed. The filter was applied to all individual (pre-stack) gathers. Filtered fk spectra were then back-transformed into the xt domain (Fig. 3). By performing this operation, we could appreciate how the fundamental mode always carries most of the energy, and that its velocities smoothly change along the line. We could also identify a suitable frequency interval for the surface wave analysis: although for some shots fundamental mode energy spans from 10 to 50 Hz (Fig. 3), for others higher frequencies are strongly attenuated, so that our analysis has been limited to the 10–25 Hz range. Finally, a simple fft was applied to each filtered shot gather to extract the phase of the signal, which was subsequently unwrapped in offset. As for tomography-derived analysis: the linear inverse problem for the 2-D active case, described in the previous section, was solved for each frequency within the spectral range of interest. Since multiple repetitions of the same shots were available, statistical data errors could be evaluated. In particular, for each shot position average phase differences between couples of receivers and their standard deviations were extracted: the first were used to fill the data vector d, the latter to compose the weight matrix W. Two inversion steps were performed: in the first step, the forward problem did not involve any regularization, allowing obtaining a first rough 2-D phase velocity model. This preliminary model was averaged in space to extract one single (average) dispersion curve, from which the frequency-wavelength relation was inferred. In the second step the inverse problem, including the roughness damping regularization constrain, was solved. The near-offset region was excluded, corresponding to offsets lower than half wavelength. The optimum \(\alpha\) parameter was chosen through a minimization process. Different \(\alpha\) values were tested and, for each of them, the modelled phase differences were computed as:

$$\begin{aligned} \mathbf{d }_{{\rm model}} = \mathbf{A } \mathbf{m } \end{aligned}$$
(12)

Then, observed phase differences were spatially smoothed: a moving average with a window length equal to the wavenumber for that frequency was used. This operation is needed to remove unrealistic lateral variations caused by noise: in fact, lateral resolution of surface waves is strictly related to the wavenumber, so that smaller measured variations do not reflect true changes in the medium properties. Finally, we defined a Cost function that corresponds to the L1 distance between the vector of the modelled phase differences, \(\mathbf{d }_{{\rm model}}\), and the vector of the smoothed phase differences, \(\mathbf{d }_{{\rm smooth}}\). In formula:

$$\begin{aligned} {\rm Cost} = \left| \mathbf{d }_{{\rm model}} - \mathbf{d }_{{\rm smooth}} \right| = \sum _{i=1}^m \left| d \phi _{i {\rm model}} - d \phi _{i {\rm smooth}} \right| \end{aligned}$$
(13)

By minimizing this Cost function, we find the optimum value of \(\alpha\), let us call it \(\alpha _{opt}\), which is therefore frequency-dependent. Figure 4 shows the relation described by Eq. (13) for one frequency (f = 14.65 Hz) and, for three points of this curve, the observed, smoothed and modelled phase differences for one shot and all active receivers: it can be noticed how the \(\alpha\) value corresponding to the minimum represents the best compromise between data (including high frequency noise) overfitting and model oversmoothing. The \(\alpha _\mathrm{opt}\) values for all frequencies are shown in Fig. 5: a general trend of the curve is visible, which implies higher smoothing for lower frequencies. The result of the tomography-derived linear inversion for all frequencies between 10 and 25 Hz is shown in Fig. 6a. The pseudo-section, which collects all dispersion curves belonging to different spatial positions of the line, highlights a clear phase velocity change starting around \(X = 180\) m, which indicates a smooth transition from loose/soft materials to more compact/resistant ones.

Fig. 3
figure 3

a Raw shot gather, corresponding to the source position \(X = 278\) m. b Fk spectrum of the (flipped) raw gather in (a). c Filtered fk spectrum in (b). d Shot gather after reverse transformation of the fk spectrum in (c)

Fig. 4
figure 4

On the left: Cost function to be minimized to find the optimum \(\alpha\) value. On the right: Raw, smoothed and modelled phase differences for one shot and all active receivers. Case B represents the trade-off between noise fitting (CASE A) and oversmoothing (CASE C)

Fig. 5
figure 5

Optimum \(\alpha\) values found through the proposed minimization procedure

Fig. 6
figure 6

Rayleigh-wave phase velocity pseudo-sections from a tomography-derived linear inversion and b MOPA

As for the MOPA analysis: the method performs a weighted linear regression of phase versus offset in each analysis window, and extracts the average slope (wavenumber). In formula:

$$\begin{aligned} \left\| \boldsymbol{ \varPhi} - \mathbf{G } \mathbf{M } \right\| = {\rm min} \end{aligned}$$
(14)

where \(\varPhi\) is the phase vector, G is the data kernel matrix, containing the offset information, and M is the vector of the unknown polynomial coefficients, including the wavenumber. The weighted least-squares solution to Eq. (14) is:

$$\begin{aligned} \mathbf{M } = ( \mathbf{G }^{\rm T} \mathbf{W } \mathbf{G } )^{-1} \mathbf{G }^{\rm T} \mathbf{W } \mathbf \varPhi \end{aligned}$$
(15)

where W is the weight matrix.

In order to get the phase vector and the weight matrix, the mean unwrapped phases for repeated shots and their standard deviations were computed. This differs from the previously illustrated tomography-derived approach, where the mean and standard deviations were computed on phase differences between receivers.

Following the procedure in Vignoli et al. (2016), we estimated phase velocities on sliding windows of frequency-dependent width. In order to find the optimum compromise between lateral resolution and robustness, the width of the window for one frequency was set as close as possible to its wavelength. Similarly, the excluded near-offset region was set to half wavelength from the source position. In order to be consistent with the previous analysis, the same frequency-wavelength relation derived from the tomography-derived inversion was used. We attributed the result of each linear regression (a wavenumber value, then converted to phase velocity) to the centre of each window. For this reason, this method does not sample regions close to the profile ends, particularly for lower frequencies. MOPA has been run for each shot position independently: as a result, as many phase velocity pseudo-sections as shot positions were produced, although all of them present missing data due to the near-offset exclusion. A final complete pseudo-section has been obtained by averaging all phase velocity measurements corresponding to the same frequency and position (Fig. 6b). There is a good correspondence between the two pseudo-sections of Fig. 6. The most glaring difference is that the tomography-derived result interests the whole acquisition line, while the result from MOPA appears discontinuous both in space and frequency. While this difference does not significantly affect a direct interpretation of the result in the frequency-offset domain, it could have a huge impact on the final depth inverted section, due to the lack of low frequencies in some positions. Other minor differences are visible: the MOPA section appears to be slightly less reliable due to the presence of sudden unrealistic lateral velocity changes.

The pseudo-section output by the tomography-derived analysis has been depth-inverted with the Dinver package from the Geopsy team (Wathelet 2008), that uses the Neighbourhood algorithm (Sambridge 1999) to solve the inverse problem. Dinver requires a set of initial parameters to define the solution space. We defined a five-layer model, with increasing layer thickness with depth, plus the half-space. In each layer, S-wave velocity, P-wave velocity and density are constant, and their value moves within pre-defined ranges (Table 1). We inverted each dispersion curve independently with no lateral constrain. Figure 7 shows the results of the inversion of two dispersion curves at different positions: on the right, the 1000 minimum misfit models generated by the algorithm are plotted together with the final smooth model, which represents the average of the best 100 models; on the left, the observed dispersion curve is shown, together with the 1000 minimum misfit modelled curves. Although the two starting curves present different characteristics, a good convergence is obtained in both cases, especially for the first 8 meters depth.

The complete inverted section is shown in Fig. 8a: we can observe a shallow low velocity body (Vs \(\sim\) 100 m/s) presenting a very regular shape, with a lower limit at 4 m depth and a smooth but well defined lateral change around \(X = 180\) m. These results are confirmed by a preliminary MASW analysis on the same data (Carrera 2019b). Note that the tomography-derived inversion allows a much finer lateral resolution than standard MASW, so that sharp lateral variations can be effectively identified and correctly located in space.

Table 1 Depth inversion parameters defining the solution space
Fig. 7
figure 7

Result of the inversion of two dispersion curves at different positions. On the left: observed dispersion curve (light grey) and 1000 best modelled dispersion curves. On the right: Final Vs model (light grey) and 1000 best generated models. The colour of both modelled dispersion curves and generated Vs model indicates the relative phase velocity misfit, averaged over all frequencies

Fig. 8
figure 8

a Quasi-2-D shear wave velocity section obtained through independent inversions of local dispersion curves. b Conductivity section from 2-D inversion of ERT data

4.2 ERT and EM Measurements

Figure 8b shows the section obtained by the electrical resistivity tomography (ERT) acquisition, to be compared with seismic results. ERT data were acquired along the same line of seismic using a roll-along scheme of four profiles with 72 electrodes 1 m-spaced, with 12–13 electrodes overlap between the lines. This allowed covering 250 m of ERT over the 272 m of the seismic line. We used the multiple gradient (MG) protocol to raise the S/N ratio on a Syscal Pro resistivimeter. The obtained pseudo-section has been inverted with CRTomo, a smoothness-constraint algorithm by Kemna (2000). The resistivity section shows a clear transition from low to high resistivity values around \(X = 180\) m, in good agreement with the previously detected seismic discontinuity. Moreover, in the left part of the section a boundary between more resistive materials on the top and a more conductive medium on the bottom is visible between 2 and 3 m depth. By crossing information from seismic and electrical measures we can infer the presence of a 4 m-thick (including 1 m top soil according to Brandstätter et al. 2013) urban waste deposit in the northern-western area, extending southeastward until \(X = 180\) m, position at which a smooth change in the material properties is verified. This interpretation is consistent with a recent study on the same area using Induced Polarization (IP) data (Flores-Orozco et al. 2020). In the bottom part of the urban waste deposit and underneath it, the more conductive zone could indicate the former presence of organic leachate, which tends to flow downwards, and/or a higher level of groundwater accumulation due to the weak compaction of the material. As for the area between \(X = 180\) m and \(X = 272\) m, the analysis of samples collected in 2001 suggests the presence of mainly construction waste material (Valtl 2005), which is compatible with the encountered shear wave velocities above 200 m/s.

Finally, Fig. 9 shows the results of the electromagnetic induction survey (EMI) collected along the same line. EMI data were acquired using a CMD-Explorer (www.gfinstruments.cz), a multi-depth probe which works with a single frequency of 10 KHz and three different receiver coil spacing (respectively, 1.48 m, 2.82 m and 4.49 m, see Table 2 for details). The probe operates at carrying height with transmitting coils either parallel or perpendicular to the ground, called, respectively, horizontal mode (HMD or High mode) and vertical mode (VMD or Low mode). High and Low modes have different penetration depths, corresponding to the nominal depths of maximum sensitivity of the instrument, and both were tested in the study site. The maximum nominal penetration depth for CMD-Explorer is of 6.7 m, obtained with the maximum coil separation and High mode (see Table 2). The EM equipment retrieves in real time the soil apparent conductivity (\(\sigma _a\)) by the quadrature values ratio between the primary generated magnetic field (Hp) and the secondary recorded one (Hs):

$$\begin{aligned} \sigma _a=\frac{4}{\omega \mu _0 s^2}\frac{H_{\rm s}}{H_{\rm p}} \end{aligned}$$
(16)

being \(\omega\) is the angular frequency, s the inter-coil spacing and \(\mu _0\) the magnetic permeability of the vacuum (4\(\pi\) 10−7 NA−2). EMI data in the studied site were collected in continuous with a sampling period of 0.5 s. Results show again a gradual decrease of the apparent conductivity southeastward starting from \(X = 180\) m, where values below 30 mS/m are observed. These values are in agreement with ERT data (Fig. 8b) and further validate the interpretation of both seismic and ERT results.

Table 2 Technical specifications of the multi-coil CMD Explorer FDEM (www.gfinstruments.cz)
Fig. 9
figure 9

Apparent conductivity values from the EMI acquisition. EMI data were collected with a GF-Explorer probe using both possible coil orientations (Low and High modes, see Table 2). Results from the probe coil 3 (coil separation of 4.49 m and transmitter frequency of 10 kHz) collected in the a low-VMD mode (3.3 m nominal depth) and in the b high-HMD mode (6.7 m nominal depth). The white dotted line indicates the beginning of the lateral variation identified through the surface wave analysis. The red line defines the presumed perimeter of the landfill

5 Discussion

This work demonstrates how a frequency-dependent tomography-derived linear inversion of 2-D active seismic data can be effective in detecting lateral velocity variations, even when the acquisition line is segmented with scarce (or no) overlap between the lines. The results obtained are consistent with the ones given by the MOPA technique, which is by far considered one the most effective techniques to detect lateral velocity variations from surface wave active data. The great advantage of the tomography-derived linear inversion consists in obtaining a velocity section covering the entire acquisition transect, since observations are phase differences between couples of adjacent receivers. For this reason, no overlap between consecutive profiles is required when acquiring long lines. A possible common disadvantage of the proposed approach and of the MOPA method is the need of isolating one single mode of propagation (usually, the fundamental mode), which implies spectra computation and filtering. This operation could be challenging in the case of abrupt lateral velocity variations, where a rapid change of the energy distribution among the different modes and of their velocities could bring to mode misinterpretation.

The choice of the \(\alpha\) parameter according to the minimization procedure proposed here overcomes the problems related to a good and reliable noise estimation for the Chi-Square test and to the subjectivity of the L-curve criterion. The optimum \(\alpha\) value gives the best match between modelled and smoothed phase differences, where the smoothing is operated in a physically consistent fashion.

Results are in very good agreement with previously collected electrical and electromagnetic data, highlighting the presence of a lateral discontinuity. The application of surface wave analysis to environmental shallow problems such landfill characterization is scarcely adopted. However, the proposed approach shows a great potential in differentiating different types of waste by their mechanical properties.

6 Conclusions

We presented a case study where tomography-derived linear inversion of surface wave data is applied to a small-scale 2-D seismic survey for the characterization of a landfill. The novelty of our work consists in (i) the use of surface wave methods to image a waste body and its geometry, (ii) the formulation of a simplified linear problem for the 2-D case and (iii) a new procedure to choose the optimum regularization based on physical properties of wave propagation.

Results obtained with this procedure were first compared with those obtained using a MOPA technique: the two pseudo-sections are very similar but the tomography-derived method produces a more continuous information and less artefacts, proving to be a very robust technique.

The quasi-2-D shear wave velocity section, which is a combination of all 1-D inverted profiles, clearly highlights the presence of a lower velocity body in the left/upper part of the section extending until \(X = 180\) m and a depth of 4 m, with velocities compatible with compacted urban wastes (between 90 and 120 m/s). A smooth transition to higher velocity materials is visible from \(X = 180\) m, where urban waste has probably been replaced with construction waste (shear wave velocities above 200 m/s).

Seismic results are consistent with both ERT and EMI data acquired along the same line, and their joined observation validates the interpretation given with the proposed surface wave analysis.