# Sensitivity of Rayleigh wave ellipticity and implications for surface wave inversion

Sensitivity of Rayleigh wave ellipticity and implications for surface wave inversion Summary The use of Rayleigh wave ellipticity has gained increasing popularity in recent years for investigating earth structures, especially for near-surface soil characterization. In spite of its widespread application, the sensitivity of the ellipticity function to the soil structure has been rarely explored in a comprehensive and systematic manner. To this end, a new analytical method is presented for computing the sensitivity of Rayleigh wave ellipticity with respect to the structural parameters of a layered elastic half-space. This method takes advantage of the minor decomposition of the surface wave eigenproblem and is numerically stable at high frequency. This numerical procedure allowed to retrieve the sensitivity for typical near surface and crustal geological scenarios, pointing out the key parameters for ellipticity interpretation under different circumstances. On this basis, a thorough analysis is performed to assess how ellipticity data can efficiently complement surface wave dispersion information in a joint inversion algorithm. The results of synthetic and real-world examples are illustrated to analyse quantitatively the diagnostic potential of the ellipticity data with respect to the soil structure, focusing on the possible sources of misinterpretation in data inversion. Inverse theory, Computational seismology, Site effects, Surface waves and free oscillations, Wave propagation 1 INTRODUCTION In many engineering and environmental problems, the shear wave velocity of near surface materials is a key parameter because it is directly linked to the ground stiffness. More specifically, when estimating the effects of the local geology on seismic amplification (the so-called ‘site effects’), a weighted average of the shear wave velocity of the upper 30 m of soil (VS30) is used in modern seismic codes (CEN 2004) as a soil classification parameter to choose the proper design spectra. Malischewsky & Scherbaum (2004) and Tuan et al. (2011) investigated extensively the behaviour of the ellipticity function for simple two-layer models, showing the relationships between the Rayleigh fundamental mode and the S-wave resonant frequency and establishing a link between the ellipticity peak and the impedance contrast between the surface layer and the bottom half-space. The inversion of Rayleigh wave ellipticity to infer the shear wave soil profile has been used for a long time now in crustal seismology (Boore & Toksöz 1969; Sexton et al.1977) and has gained new popularity in recent years for both crustal (Yamanaka et al.1994; Tanimoto & Alvizuri 2006; Tanimoto & Rivera 2008) and shallow applications (Scherbaum et al.2003; Arai & Tokimatsu 2004; Boaga et al.2013; Hobiger et al.2013; Castellaro 2016). This renewed interest, especially for near surface characterization, is mainly due to the fact that the behaviour of the ellipticity function versus frequency can be evaluated by rapid and cost-effective single-station Horizontal versus Vertical Spectral Ratio (HVSR) measurements (Fäh et al.2001; Hobiger et al.2009). Unluckily, the ellipticity function alone is not generally sufficient to determine unambiguously the soil structure because of scaling effects and intrinsic non-uniqueness (Scherbaum et al.2003; Hobiger et al.2013). Because of these limitations, ellipticity data are generally inverted jointly along with surface wave dispersion measurements to better constrain the soil structure. Maranò et al. (2012) and Maranò et al. (2017) applied newly developed array estimation techniques to pinpoint accurately the frequency of singularities using the ellipticity angle. Although the method requires large aperture, multichannel arrays, it allows to identify Rayleigh ellipticity over a broad range of frequencies even near singularities (where the particle motion is vertically or horizontally polarized) and to retrieve information about the prograde or retrograde particle motion of Rayleigh waves. A common approach for near surface characterization is to use small-aperture arrays for dispersion measurements in conjunction with single-station HVSR measurements to estimate the ellipticity at lower frequency (Scherbaum et al.2003; Arai & Tokimatsu 2005; Hobiger et al.2013). This way, the deployment of large aperture arrays to retrieve the deep soil structure can be avoided. In near surface applications, the presence of low-velocity layers (LVLs) and high stiffness contrasts within the soil column are generally responsible for mode jumps, ‘osculation points’ and higher-mode domination over specific frequency ranges (Cercato 2009; Cercato et al.2010; Cercato 2011). Misinterpretation of these complicated features can lead to pitfalls and failures in surface wave inversion. It has been recognized that low-velocity surface layers can greatly affect the Rayleigh wave particle motion (Tsuboi & Saito 1983) and that the shape of the ellipticity function can be considered as an index of mode contamination and structural anomalies such as high stiffness contrasts (Boaga et al.2013; Castellaro 2016). Hobiger et al. (2013) made a systematic study of the estimation of shear wave velocity structures by jointly inverting Rayleigh wave ellipticity and short-array Rayleigh wave dispersion measurements. They identified two classes of data to describe the main features of the observed ellipticity: data exhibiting singularities and well defined peaks or, conversely, data showing regular and smooth variation in the ellipticity function. Their study was limited to fundamental mode data only and to normally dispersive sites (where the shear wave velocity is monotonically increasing with depth). Under these limitations, they analysed the performance of the joint inversion on several synthetic simulations and experimental sites and concluded that, for ellipticity curves with singularities, the right flank of the ellipticity peak (and not the values just around the peak) carries the most valuable information on the soil structure. On the other hand, for ellipticity curves without singularities, the whole area of the fundamental peak should be inverted if such a peak exists. From a more general point of view, the diagnostic potential of the ellipticity function to invert for the soil structure can be assessed by a quantitative sensitivity analysis, which requires the computation of the derivatives of ellipticity with respect to the layer parameters. It is generally recognized that the layer thickness and the S-wave velocity have the most significant effects on the ellipticity function (Arai & Tokimatsu 2004). In previous studies, which were mostly focused on deep crustal Earth models, the derivatives of the ellipticity function were obtained numerically (Arai & Tokimatsu 2004), or employing variational theory (Tsuboi & Saito 1983; Tanimoto & Tsuboi 2009). The first method requires extensive forward modelling of the ellipticity function and is subject to numerical errors. The second is analytical in its formulation, although it requires the integration of the eigenfunctions across each layer, which was performed numerically in the cited works. The sensitivity of Rayleigh wave ellipticity with respect to the subsoil structure has been very rarely explored for near surface characterization and shallow applications. On this basis, the main objectives of this work are: To derive a new fully analytical method for computing the sensitivity of Rayleigh wave ellipticity to structural parameters in a layered elastic half-space. To perform thorough sensitivity analyses on selected Earth models using this new formulation, to analyse quantitatively the diagnostic potential of the ellipticity data with respect to the soil structure, also focusing on the possible sources of misinterpretation in data inversion. To examine and give quantitative criteria on how ellipticity data can efficiently complement surface wave dispersion information in a joint inversion algorithm, to reduce the cost and increase the depth of investigation in common applications for seismic site characterization. This paper is divided as follows: the next section (Section 2) introduces the new analytical formulation for computing the sensitivity of the ellipticity function. Successively, a few theoretical examples of sensitivity computation will be presented in Section 3. Section 4 is devoted to the analysis of the complementarity between ellipticity and surface wave data by synthetic inversion, whereas a real world application will be eventually presented in Section 5, to point out the fruitful integration between ellipticity a and surface wave dispersion data, whereas Sections 5 and 6 are devoted to Discussion and Conclusions, respectively. 2 ANALYTICAL FORMULATION FOR COMPUTING THE SENSITIVITY OF RAYLEIGH WAVE ELLIPTICITY The earth model under consideration is a 1-D elastic layered half-space (Fig. 1). This parametric model, made up by horizontal, elastic, isotropic and homogeneous layers, is bounded above by a free surface and below by a homogeneous isotropic elastic half-space of infinite thickness. For the ith soil layer, αi is P-wave velocity, βi is S-wave velocity, ρi is (volume) density and hi is layer thickness. Figure 1. View largeDownload slide The plane-layered earth model of reference. Each layer is elastic, homogeneous and isotropic. Figure 1. View largeDownload slide The plane-layered earth model of reference. Each layer is elastic, homogeneous and isotropic. The complete set of structural parameters of the earth model in Fig. 1 is:   $$\mathbf {Q}=(\alpha _1,\ldots ,\alpha _N,\beta _1,\ldots ,\beta _N,\rho _1,\ldots ,\rho _N,h_1,\ldots ,h_{N-1}).$$ (2.1)For a given frequency ω, the modal ellipticity function ε of Rayleigh waves can be computed as the amplitude of the ratio between the horizontal (Ux) and vertical (Uz) modal eigenfunctions at the ground surface (Boore & Toksöz 1969; Takeuchi & Saito 1972):   $$\varepsilon =\left|\frac{ U_x(\omega ,z=0)}{ U_z(\omega ,z=0)}\right|.$$ (2.2)This expression holds for any mode, provided that the mode exists at the given frequency ω. In eq. (2.2), the dependence of the eigenfunctions on the layer parameters appearing in eq. (2.1) is omitted for simplicity. To solve the Rayleigh wave eigenvalue problem for such a model, matrix propagator theory is frequently used (Haskell 1964; Dunkin 1965; Wang & Herrmann 1980). A key point in this class of methods is the computation of the (4 × 4) layer stack matrix R of the layered half-space as:   $$\mathbf {R} = \mathbf {\bar{T}}\mathbf {G}_{N-1}\ldots \mathbf {G}_1$$ (2.3)where $$\mathbf {\bar{T}}$$ is the inverse of the bottom half-space matrix and Gi is the propagator matrix of layer i. When performed in a straightforward way according to eq. (2.3), the computation of matrix R is affected by numerical instability at high frequencies: this problem can be addressed by employing subdeterminants of order two of the layer matrices, which are free of growing exponential terms (Dunkin 1965; Watson 1970). In our notation, a subdeterminant (or minor determinant) of order two of matrix R is indicated as:   $$r\mid^{ij}_{kl}\, =\, r_{ik}r_{jl}-r_{il}r_{jk}.$$ (2.4)For the restatement of the Rayleigh value eigenproblem in terms of order-two subdeterminants, Dunkin (1965) introduced the following theorem: Theorem 1. Given the layer stack matrix R in eq. (2.3), a generic second order subdeterminant $$r\mid ^{ij}_{kl}$$ of R can be obtained as:   $$r\mid ^{ij}_{kl}=\bar{t}_N\mid ^{ij}_{ab}g_{N-1}\mid ^{ab}_{cd}\ldots g_{1}\mid ^{ef}_{kl}$$ (2.5)where the summed pairs of indices are to be only distinct pairs of distinct indices. The analytical expressions of the subdeterminants of $$\bar{\mathbf {T}}$$ and layer propagator matrices Gi (i = 1, …, N − 1) appearing in eq. (2.5), are reported in Appendix, where the definitions of the main parameters appearing in the matrix coefficients are also introduced. As the Dunkin's notation was slightly modified in this paper for the sake of coherence with current modern notation and previous work (Cercato 2007), a link with the parameter definition appearing in the original paper (Dunkin 1965) is explicitly given in Table A1. Both the dispersion function of Rayleigh waves FR and the Rayleigh wave ellipticity ε can be expressed in terms of subdeterminants of matrix R (Dunkin 1965; Harkrider 1970) so that:   $$F_R\left[\omega ,\mathbf {Q},c(\omega ,\mathbf {Q})\right]=r\mid ^{12}_{12}=r_{11}r_{22}-r_{12}r_{21}$$ (2.6)  $$\varepsilon \left[\omega ,\mathbf {Q},c(\omega ,\mathbf {Q})\right]=\left| -\frac{r_{11}}{r_{12}}\right|=\left|\frac{r\mid ^{12}_{14}}{r\mid ^{12}_{13}}\right| \varepsilon \left[\omega ,\mathbf {Q},c(\omega ,\mathbf {Q})\right]= -\frac{r_{11}}{r_{12}}=\frac{r\mid ^{12}_{14}}{r\mid ^{12}_{13}}$$ (2.7)where c is phase velocity, ω is (circular) frequency and Q is the set of structural parameters as defined in eq. (2.1). Applying the Dunkin's theorem to the secular equation FR in eq. (2.6), one gets the well-known dispersion equation:   $$F_R(\omega,k)\,=\,r\mid ^{12}_{12}=\bar{t}_N\mid ^{12}_{ab}g_{N-1}\mid ^{ab}_{cd}\ldots g_{1}\mid^{ef}_{12}\!.$$ (2.8)In a previous paper (Cercato 2007) it was demonstrated how the partial derivatives of modal phase velocity with respect to structural parameters can be obtained taking advantage of the Implicit function theorem (Rudin 1976) as:   $$\frac{\partial c}{\partial {q_i}}=c_{q_i}=\frac{(F_R)_{q_i}}{F_R^{\prime }}=\frac{\left(r\mid ^{12}_{12}\right)_{q_i}}{\left(r\mid ^{12}_{12}\right)^{\prime }}.$$ (2.9)Following the notation in Cercato (2007), qi is used to indicate a generic parameter of layer i, so that qi can be either the P-wave velocity αi, the S-wave velocity βi, the density ρi or the layer thickness hi. Partial differentiation of a function g with respect to qi will be indicated as $$g_{q_i}$$, whereas the prime mark (i.e. g΄) will be used to indicate partial differentiation with respect to phase velocity c. This should be kept in mind for all the expressions appearing in the following. The partial derivatives of FR are obtained using two separate recursions from the half-space to the surface. One for the structural parameter qi:   $$\frac{\partial F_R}{\partial q_i}=\left(r\mid ^{12}_{12}\right)_{q_i}=\bar{t}_N\mid ^{12}_{ab}g_{N-1}\mid ^{ab}_{cd}\ldots \left(g_i\mid ^{ef}_{gh}\right)_{q_i}\ldots g_{1}\mid ^{il}_{12}$$ (2.10)and one for phase velocity c (΄ prime mark):   \begin{eqnarray} r_{i}^{\prime }\mid ^{12}_{12}&=& r_{i+1}^{\prime }\mid ^{12}_{12}g_i\mid ^{12}_{12}+r_{i+1}\mid ^{12}_{12}g_i^{\prime }\mid ^{12}_{12} +r_{i+1}^{\prime }\mid ^{12}_{13}g_i\mid ^{13}_{12}+r_{i+1}\mid ^{12}_{13}g_i^{\prime }\mid ^{13}_{12}+ r_{i+1}^{\prime }\mid ^{12}_{14}g_i\mid ^{14}_{12} \nonumber\\ && +\, r_{i+1}\mid ^{12}_{14}g_i^{\prime }\mid ^{14}_{12} +r_{i+1}^{\prime }\mid ^{12}_{23}g_i\mid ^{23}_{12}+r_{i+1}\mid ^{12}_{23}g_i^{\prime }\mid ^{23}_{12}+ r_{i+1}^{\prime }\mid ^{12}_{24}g_i\mid ^{24}_{12}+r_{i+1}\mid ^{12}_{24}g_i^{\prime }\mid ^{24}_{12} \nonumber\\ && +\, r_{i+1}^{\prime }\mid ^{12}_{34}g_i\mid ^{34}_{12}+r_{i+1}\mid ^{12}_{34}g_i^{\prime }\mid ^{34}_{12}, \end{eqnarray} (2.11)where the number of computations can be further simplified because of certain symmetries in the layer matrix subdeterminants (Watson 1970; Cercato 2007). In fact, the implementation of the recursions in eqs (2.10) and (2.11) requires the computation of a subset of matrix R subdeterminants, namely $$r_i\mid ^{12}_{12}$$, $$r_i\mid ^{12}_{13}$$, $$r_i\mid ^{12}_{14}$$, $$r_i\mid ^{12}_{24}$$ and $$r_{i}\mid ^{12}_{34}$$, to carry on the recursion layer-by-layer, together with their partial derivatives with respect to c and qi. As will be shown in the following, once that these recursions have been performed up to the free surface, the sensitivity of the ellipticity of Rayleigh waves with respect to layer parameters can be computed with little additional effort. The only difficulty is that the ellipticity ε depends on the layer parameters not only directly, but also via the intermediate variable c, the phase velocity, that changes when a single layer parameter qi is allowed to vary. Therefore, to find the total variation of ε associated with an arbitrary perturbation of a layer parameter qi, we must add the indirect dependency on phase velocity and compute the total derivative (Rudin 1976; Riley 2006) of the ellipticity with respect to the structural parameter of interest. Using a notation which is common in multivariate calculus, the total derivative with respect to the parameter qi will be indicated as d/dqi, meaning that it is not assumed that phase velocity is held constant while qi varies, in contrast with partial differentiation ∂/∂qi. Starting from eq. (2.7) and applying the rule of derivation of quotient of functions to the total derivative, one gets:   $$\frac{{\rm d} \varepsilon }{{\rm d} q_i}=\frac{\frac{{\rm d} r\mid ^{12}_{14}}{{\rm d} q_i}r\mid ^{12}_{13} -\frac{{\rm d} r\mid ^{12}_{13}}{{\rm d} q_i}r\mid ^{12}_{14}}{(r\mid ^{12}_{13})^2}= \frac{1}{r\mid ^{12}_{13}}\left(\frac{{\rm d} r\mid ^{12}_{14}}{{\rm d} q_i}-\frac{{\rm d} r\mid ^{12}_{13}}{{\rm d} q_i}\varepsilon \right).$$ (2.12)To derive the total derivatives of the subdeterminants $$r\mid ^{12}_{13}$$ and $$r\mid ^{12}_{13}$$ appearing in eq. (2.12), one can take advantage of the rule of derivation of composite functions (chain rule) so that:   $$\frac{{\rm d} r\mid ^{12}_{pq}}{{\rm d} q_i}=\frac{\partial r\mid ^{12}_{pq}}{\partial q_i}+ \frac{\partial r\mid ^{12}_{pq}}{\partial c}\frac{\partial c}{\partial q_i}\qquad pq=13,14.$$ (2.13)This way, the total derivatives of subdeterminants $$r\mid ^{12}_{13}$$ and $$r\mid ^{12}_{13}$$ with respect to qi can be expressed in terms of partial derivatives. As described above, these partial derivatives can be computed analytically according to Cercato (2007), making use of a recursive scheme from the bottom to the top of the layer stack, employing the analytical partial derivatives of a subset of layer matrix subdeterminants of order two. The total derivative of ellipticity as stated in eq. (2.12) can be directly computed with little effort, substituting the proper values of $${\rm d} r\mid ^{12}_{pq} / {\rm d} q_i$$ (pq = 13, 14) computed from the partial derivatives according to (2.13). The proposed method is analytical and exact, in the sense that the numerical accuracy depends only on machine precision. Once that the derivatives are computed, the sensitivity of Rayleigh wave ellipticity can be obtained as the normalized absolute value of the total derivative (Arai & Tokimatsu 2004, 2005):   $$S_\varepsilon (q_i)=\left|\frac{q_i}{\varepsilon }\frac{{\rm d} \varepsilon }{{\rm d} q_i}\right|$$ (2.14)which is a pure numerical quantity. 3 THEORETICAL EXAMPLES Consider a typical crustal model like Model TaN in Table1. The Rayleigh wave ellipticity and phase velocity for the fundamental mode of Rayleigh waves are displayed versus period in Fig. 2(a). The derivatives of the ellipticity function were computed analytically using the new approach. The sensitivity kernels of the ellipticity function versus depth at 0.2 Hz (period 5 s) for the P-wave and S-wave velocities and density are shown in Fig. 2(b), to be compared to fig. 2 in Tanimoto & Tsuboi (2009). To compute the sensitivity kernels according to the newly developed approach, the model layers are divided into fictitious sublayers (500m thick in this case). The computed sensitivity values are assigned to the mean depths of the fictitious layers. Figure 2. View largeDownload slide Model TaN in Table 1 after Tanimoto & Tsuboi (2009). (a) Rayleigh wave ellipticity and phase velocity for the fundamental mode of Rayleigh waves. (b) Sensitivity kernels of the ellipticity (fundamental mode) at 0.2 Hz for shear velocity, density and compressional velocity. (c) Sensitivity of the ellipticity (fundamental mode) with respect to the layer S-wave velocity. (d) Sensitivity of the ellipticity (fundamental mode) with respect to layer thickness. Figure 2. View largeDownload slide Model TaN in Table 1 after Tanimoto & Tsuboi (2009). (a) Rayleigh wave ellipticity and phase velocity for the fundamental mode of Rayleigh waves. (b) Sensitivity kernels of the ellipticity (fundamental mode) at 0.2 Hz for shear velocity, density and compressional velocity. (c) Sensitivity of the ellipticity (fundamental mode) with respect to the layer S-wave velocity. (d) Sensitivity of the ellipticity (fundamental mode) with respect to layer thickness. Table 1. Structural parameters of the three earth models used for the theoretical examples. Layer  α (m/s)  β (m/s)  ρ (kg/m−3)  h (m)  Model TaN  Three Layered model in Tanimoto & Tsuboi (2009)  1  5500  3000  2500  15000  2  6500  3600  2850  15000  3  7800  4500  3300  ∞  Model HoA  Model A in Hobiger et al.(2013)  1  540  120  1800  5  2  900  200  1800  15  3  1440  320  1800  45  4  2810  625  1800  135  5  6250  2500  2000  ∞  Model HoB  Model B in Hobiger et al.(2013)  1  540  120  1800  5  2  900  200  1800  15  3  1440  320  1800  45  4  2430  540  1800  135  5  2520  840  2000  ∞  Layer  α (m/s)  β (m/s)  ρ (kg/m−3)  h (m)  Model TaN  Three Layered model in Tanimoto & Tsuboi (2009)  1  5500  3000  2500  15000  2  6500  3600  2850  15000  3  7800  4500  3300  ∞  Model HoA  Model A in Hobiger et al.(2013)  1  540  120  1800  5  2  900  200  1800  15  3  1440  320  1800  45  4  2810  625  1800  135  5  6250  2500  2000  ∞  Model HoB  Model B in Hobiger et al.(2013)  1  540  120  1800  5  2  900  200  1800  15  3  1440  320  1800  45  4  2430  540  1800  135  5  2520  840  2000  ∞  View Large At this frequency, the sensitivity of the ellipticity function is confined to the first layer only. If we examine Figs 2(c) and (d), reporting the trend versus period of the sensitivity of the ellipticity function with respect to layer shear velocity and layer thickness, respectively, it is self-evident that at low frequency the shear velocities are the governing parameters on the sensitivity of ε. As shown also in Tanimoto & Tsuboi (2009), for such a simple model the numerical differentiation is well in agreement with the analytical results, in the frequency range of interest for crustal seismology. In this case, the advantage of the analytical approach is mainly the speed of computation. Let us now examine the two slightly more complicated and shallower earth models, named HoA and HoB after Hobiger et al. (2013), whose structural parameters are summarized in Table1. The ellipticity and the phase velocity of the fundamental mode of Rayleigh waves are displayed in Figs 3(a) and (c) for model HOA and in Figs 3(b) and (d) for model HOB, respectively. The theoretical ellipticity curve exhibits singularities (peaks and troughs) for model HoA (Fig. 3a), whereas the theoretical ellipticity curve is not discontinuous in the mathematical sense in the case of low contrasts and smooth variations of the shear wave velocities as in model HoB (Fig. 3b). Figure 3. View largeDownload slide Model HoA and Model HoB in Table 1 after Hobiger et al. (2013). Ellipticity function and phase velocity of the fundamental mode. (a) Model HoA : ellipticity function versus frequency. (b) Model HoB: ellipticity function versus frequency. (c) Model HoA: phase velocity versus frequency. (d) Model HoB: phase velocity versus frequency. Figure 3. View largeDownload slide Model HoA and Model HoB in Table 1 after Hobiger et al. (2013). Ellipticity function and phase velocity of the fundamental mode. (a) Model HoA : ellipticity function versus frequency. (b) Model HoB: ellipticity function versus frequency. (c) Model HoA: phase velocity versus frequency. (d) Model HoB: phase velocity versus frequency. We report the results of these two models for comparison to previously published results (Hobiger et al.2013) and to illustrate and support some general findings, which have been verified by many simulations on different earth models which are not shown here for the sake of brevity. The derivatives of the ellipticity function with respect to shear wave velocity were computed analytically using the new approach and their values are displayed versus frequency in Fig. 4(a) for model HoA and in Fig. 4(b) for model HoB, respectively. The behaviour of the sensitivity of the ellipticity function with respect to the layer parameters will be commented later in this paragraph. Let us analyse first the precision of the new analytical approach when compared to numerical differentiation for these two earth models. Figs 4(c) and (d) show the absolute percentage error between the analytical values of the derivatives and their numerical approximation obtained employing a two-point central finite difference (FD) formula, with a spacing equal to 1m in computing close point solutions of different shear wave velocities. The accuracy of numerical differentiation depends on the interval step, on the precision of the phase velocity modal solutions and on the choice of the FD formula (two, three points approximations or more). The analytical results are obtained at lower computational cost since all the parameters needed for calculating the derivative—as stated in eq. (2.12)—can be computed using a single recursion through the entire layer stack. The FD formulation suffers from numerical errors at low values of the derivative and at the discontinuity points, which is a potential drawback in the Jacobian calculation as employed, for instance, in linearized inversion, especially at the high frequency range of interest for engineering and exploration seismology. For near surface characterization, the analytical approach is thus superior in both precision and computational velocity (efficiency) when compared to the standard numerical approach. Figure 4. View largeDownload slide Model HoA and Model HoB in Table 1 after Hobiger et al. (2013). (a) Model HoA: derivative of the ellipticity function with respect to layer shear wave velocity. (b) Model HoB: derivative of the ellipticity function with respect to layer shear wave velocity. (c) Model HoA: absolute percentage error between the analytical and the numerical computation of the derivatives. (d) Model HoB: absolute percentage error between the analytical and the numerical computation of the derivatives. Figure 4. View largeDownload slide Model HoA and Model HoB in Table 1 after Hobiger et al. (2013). (a) Model HoA: derivative of the ellipticity function with respect to layer shear wave velocity. (b) Model HoB: derivative of the ellipticity function with respect to layer shear wave velocity. (c) Model HoA: absolute percentage error between the analytical and the numerical computation of the derivatives. (d) Model HoB: absolute percentage error between the analytical and the numerical computation of the derivatives. The sensitivity functions with respect to all the layer parameters are shown in Fig. 5 for model HoA and in Fig. 6 for model HoB, respectively. For model HoA (Fig. 5), the ellipticity function is discontinuous and the sensitivity goes virtually to infinity at the frequency value of the peaks. Therefore, around the peaks the inverse problem is highly instable, because the sensitivity of all the model parameters increases drastically. For this reasons, the experimental data points of the ellipticity function around the main peaks should be avoided in the inversion, because little variations in one of the model parameters may be responsible of large variations in the ellipticity function, adding instability to the inversion process. Figure 5. View largeDownload slide Model HoA in Table 1 after Hobiger et al. (2013). Sensitivity analysis of the fundamental mode. (a) Derivative of the ellipticity function with respect to layer shear wave velocity. (b) Derivative of the ellipticity function with respect to layer P-wave velocity. (c) Derivative of the ellipticity function with respect to layer density. (d) Derivative of the ellipticity function with respect to layer thickness. Figure 5. View largeDownload slide Model HoA in Table 1 after Hobiger et al. (2013). Sensitivity analysis of the fundamental mode. (a) Derivative of the ellipticity function with respect to layer shear wave velocity. (b) Derivative of the ellipticity function with respect to layer P-wave velocity. (c) Derivative of the ellipticity function with respect to layer density. (d) Derivative of the ellipticity function with respect to layer thickness. Figure 6. View largeDownload slide Model HoB in Table 1 after Hobiger et al. (2013). Sensitivity analysis of the ellipticity function of the fundamental of Rayleigh waves to the layer parameters. (a) Sensitivity of the ellipticity function with respect to layer S-wave velocity. (b) Sensitivity of the ellipticity function with respect to layer P-wave velocity. (c) Sensitivity of the ellipticity function with respect to layer density. (d) Sensitivity of the ellipticity function with respect to layer thickness. Figure 6. View largeDownload slide Model HoB in Table 1 after Hobiger et al. (2013). Sensitivity analysis of the ellipticity function of the fundamental of Rayleigh waves to the layer parameters. (a) Sensitivity of the ellipticity function with respect to layer S-wave velocity. (b) Sensitivity of the ellipticity function with respect to layer P-wave velocity. (c) Sensitivity of the ellipticity function with respect to layer density. (d) Sensitivity of the ellipticity function with respect to layer thickness. On the other hand, the sensitivity functions shown for model HoB (Fig. 6) depict a very different situation, with finite values of the sensitivities in the entire frequency range and extremely low values for the layer P-wave velocity α and density ρ. The sensitivity values are not evenly distributed along frequency, and only for a few cases a well-defined frequency range, where the sensitivity of a certain layer parameter is dominating the other values, is clearly notable, such as, for instance, in the 08–1.5 Hz frequency band for the shear wave velocity sensitivity of layer 4 in Fig. 6a. For both models, the sensitivities indicate that shear wave velocity and layer thickness play a major role as far as ellipticity inversion is concerned, although in the case of discontinuous ellipticity functions it is necessary to discard the ellipticity values around the peaks to maintain a well-posed inverse problem. As stated in the Introduction, among the objectives of this work, a relevant topic is to evaluate how ellipticity data can efficiently complement surface wave dispersion information. To this end, the comparison of the sensitivities of ε and c with respect to the layer parameters (especially the more relevant to the inversion process such as shear wave velocity β) may indicate quantitatively the potential of a coupled inversion when compared to the uncoupled approach. In Fig. 7 we report the sensitivities of phase velocity with respect to shear wave velocity and layer thickness. We already showed that the sensitivity of ε with respect to β, displayed in Figs 5(a) and 6(a) for model HoA and HoB, respectively, is concentrated for most layer parameters in the same frequency range. This is also true, as noted above, for model HoB. On the other hand, for both earth models, the phase velocity sensitivity (Figs 7 a and b) is evenly distributed on the layer shear velocities as far as frequency is concerned. This means, in principle, that the dispersion data provide enough sensitivity to resolve the S-wave velocity of each layer within specific frequency ranges, whereas the ellipticity function—being a fractional composite quantity according to the definition in eq. (2.2)—has a more complex sensitivity distribution along frequency. The frequency distribution of the sensitivities is evidently complicated by the presence of high stiffness contrasts as in model HoA, for both c and ε. If we compare Figs 7(a) and (b), it is self-evident how the high stiffness contrast in model HoA caused by the increased values of β for layers 4 and 5 (when compared to model HoB as in Table1) is responsible for heavy changes in the frequency distribution of the sensitivities in the 0.5–3 Hz frequency band. The same considerations are still valid for the ellipticity function (Fig. 5 a and Fig. 6a, for instance), although it should be noted that the ellipticity function is less influenced by the value of the shear wave velocity of the bottom half-space when there is a remarkable increase between the velocities of the finite-thickness layers as between layers 3 and 4 in model HoA. An analogous behaviour is shown for the sensitivities with respect to layer thickness (Figs 5d and 7c for model HoA and Figs 6d and 7d for model HoB). Figure 7. View largeDownload slide Model HoA and Model HoB in Table 1 after Hobiger et al. (2013). Sensitivity of the phase velocity for the fundamental mode. (a) Model HoA: sensitivity of the phase velocity with respect to layer S-wave velocity. (b) Model HoB: sensitivity of the phase velocity with respect to layer S-wave velocity. (c) Model HoA: sensitivity of the phase velocity with respect to layer thickness. (d) Model HoB: sensitivity of the phase velocity with respect to layer thickness. Figure 7. View largeDownload slide Model HoA and Model HoB in Table 1 after Hobiger et al. (2013). Sensitivity of the phase velocity for the fundamental mode. (a) Model HoA: sensitivity of the phase velocity with respect to layer S-wave velocity. (b) Model HoB: sensitivity of the phase velocity with respect to layer S-wave velocity. (c) Model HoA: sensitivity of the phase velocity with respect to layer thickness. (d) Model HoB: sensitivity of the phase velocity with respect to layer thickness. 4 SURFACE WAVE INVERSION AND ELLIPTICITY: SYNTHETIC INVERSION EXAMPLES On the basis of the considerations above, we illustrate in this paragraph a few synthetic inversions based on the theoretical data computed for the models in Table1. The idea is to show the impact of the introduction of ellipticity data into surface wave modal inversion of phase velocity. The Very Fast Simulated Annealing (VFSA) code (Ingber 1989; Sen & Stoffa 2013) is used to invert the data. The implementation described in Cercato (2011) has been adapted to perform the joint inversion of phase velocity and ellipticity in the VFSA inversion scheme. In the VFSA algorithm, the random models are drawn from a Cauchy like distribution which is, in turn, function of temperature (Ingber 1989). The acceptance criterion of a random model is the key of the method: random models that lower the misfit are always accepted, whereas if a random model increases the error function, it is accepted according to a probability criterion which is also temperature dependent. This is the key feature of the SA algorithm to assure that the solution is not trapped into local minima of the error function (Sen & Stoffa 2013). The energy function is defined as the relative RMS (Root mean square) proportional error:   $$E_m=\sqrt{\frac{1}{N}\sum _{j=1}^N\left|\frac{d^{{\rm obs}}_j-d^{{\rm pre}}_j}{d^{{\rm pre}}_j}\right|^2}.$$ (4.1) In eq. (4.1), for each frequency fj of the N data points, the observed data points $$d^{{\rm obs}}_j$$ and the predicted data points $$d^{{\rm pre}}_j$$ can be either phase velocity or ellipticity values. The uncertainties on the inverted solutions are obtained by importance sampling (Sen & Stoffa 1996; Bhattacharya et al.2003; Sen & Stoffa 2013) where the VFSA algorithm is used as a maximum a posteriori (MAP ) estimator of the posterior probability distribution (PPD). Sen & Stoffa (1996) showed that the models sampled by multiple VFSA runs provide a good estimate, although underrated, of uncertainty which, theoretically, should be obtained using a more numerically intensive Gibbs’ sampling procedure at the constant temperature of T = 1. We employed 20 different VFSA runs (with different starting models) for each inversion in order to guarantee the convergence of the to guarantee a faster convergence than the PPD. The number of iteration for each independent run is set to 10 000 (at least), so that, overall a complete data inversion consists of at least 200 000 sampled models. No misfit threshold is set to reject models, so that occasionally random models with large errors can be included in the accepted population depending on the size of the model parameter's space. The posterior covariance analysis on the accepted models allows the interpreter to appraise the uncertainty estimation in terms of standard deviations on the layer-parameter unknowns and their mutual correlation and interdependence via the posterior correlation matrix, as illustrated in Cercato (2011). A comparison of different methods for estimating the PPD via importance sampling is reported in Sen & Stoffa (2013), chap.8, and is well beyond the scope of this paper. In the case of synthetic inversion, the use of posterior uncertainty assessment applied to synthetic data inversion is intended to quantify the impact of ellipticity information on the solution appraisal and therefore on the overall inversion uncertainty. Among the numerous synthetic inversions on different earth models, only a couple of examples are reported, in order to illustrate the most relevant findings related to the sensitivity analysis of ellipticity, while preserving the readability of this work. In the synthetic inversions, the model space is such that each model parameter (namely shear wave velocities and thicknesses) is allowed to vary between one half and two times of the real value, mimic a quite wider choice compared with usual assumptions (Hobiger et al.2013). I chose not to add noise to the data to be consistent with previous work (Hobiger et al.2013) and to avoid the mix of sensitivity and uncertainty issues in our comparative analysis. In a first synthetic inversion considering model HoA, the fundamental mode phase velocity data is inverted in the frequency band 2–30 Hz, which is the same adopted by Hobiger et al. (2013) in their synthetic inversions. The results are shown in Fig. 8. According to Fig. 5(a), at 2Hz the sensitivity of ε with respect to shear wave velocity is extremely low, so that the velocity of the bottom half-space cannot be accurately retrieved from dispersion data only. This leads to high uncertainty in the inverted bottom-half-space velocity as in Fig. 8(d), as well as in the appraisal of the thickness of the fourth layer (Fig. 8e). Figure 8. View largeDownload slide Modal VFSA synthetic inversion for model HoA in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The solid black line is the true model. (b) Observed (true) and predicted dispersion curves. (c) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (d) Uncertainty estimation for layer shear velocities. (e) Uncertainty estimation for layer thicknesses. Figure 8. View largeDownload slide Modal VFSA synthetic inversion for model HoA in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The solid black line is the true model. (b) Observed (true) and predicted dispersion curves. (c) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (d) Uncertainty estimation for layer shear velocities. (e) Uncertainty estimation for layer thicknesses. On a successive run, I inverted the same phase velocity data along with the ellipticity curve in the range 0.9–1.7 Hz, as reported in one of the synthetic tests in Hobiger et al. (2013, fig. 3 and table 2). The ellipticity data are located between the two discontinuities of the ellipticity function, on the ‘right flank’ of the fundamental-mode curve with respect to the main (low-frequency) peak. Looking again at the sensitivity curves of the ellipticity function in Figs 5(a) and (d), it can be noted that in the frequency range picked for the inversion, only the fourth layer velocity and the third and the fourth layer thicknesses exhibit considerable sensitivities. The results of this joint inversion are shown in Fig. 9. As expected from the preliminary examination of the sensitivities functions, the comparison between the standard deviations of shear velocity and thickness (Figs 8d and e) obtained from the phase-velocity inversion with the ones (Figs 9e and f) of the joint inversion, confirms that the uncertainty on the shear velocity of the bottom half-space is virtually not influenced by the introduction of the ellipticity data. It is also worth noting that both the uncertainties on the retrieved velocity and thickness of the fourth layer are remarkably reduced by the introduction of ellipticity data, confirming how ellipticity data can fruitfully complement phase velocity information when the sensitivity of the ellipticity function is finite and not in the proximity of a discontinuity in the ellipticity function. In addition, the actual behaviour of the ε near the discontinuity points (peaks of the ellipticity functions) is badly retrieved in practice from experimental data, with high uncertainty and possible misidentification of the real H/V spectral ratio values, so that the insertion of highly uncertain and highly sensitive data into the inversion process may be the source of relevant inversion pitfalls in real data interpretation. Figure 9. View largeDownload slide Joint inversion of Rayleigh wave ellipticity and phase velocity of the fundamental mode for model HoA in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The black line is the true model. (b) Observed (true) and predicted ellipticity functions. (c) Observed (true) and predicted dispersion curves. (d) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (e) Uncertainty estimation for layer shear velocities. (f) Uncertainty estimation for layer thicknesses. Figure 9. View largeDownload slide Joint inversion of Rayleigh wave ellipticity and phase velocity of the fundamental mode for model HoA in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The black line is the true model. (b) Observed (true) and predicted ellipticity functions. (c) Observed (true) and predicted dispersion curves. (d) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (e) Uncertainty estimation for layer shear velocities. (f) Uncertainty estimation for layer thicknesses. Fig. 10 shows the results of phase-velocity modal inversion of the fundamental mode in the frequency range 2–30 Hz for model HoB, aimed again at reproducing some of the results in Hobiger et al. (2013). In this case, from Fig. 6(d) it can be noted how the sensitivity of data with respect to the shear wave velocity of the bottom half-space is, although not particularly high, fairly more consistent at the lower values of the frequency range (2–3 Hz) when compared to the sensitivity of the same parameter for model HoA (Fig. 5a). Consequently, the inverted model of the modal (phase-velocity) inversion is fairly better resolved if compared to the HoA inversion in Fig. 8, particularly as far as the deeper part of the profile is concerned. Figure 10. View largeDownload slide Modal VFSA synthetic inversion for model HoB in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The solid black line is the true model. (b) Observed (true) and predicted dispersion curves. (c) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (d) Uncertainty estimation for layer shear velocities. (e) Uncertainty estimation for layer thicknesses. Figure 10. View largeDownload slide Modal VFSA synthetic inversion for model HoB in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The solid black line is the true model. (b) Observed (true) and predicted dispersion curves. (c) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (d) Uncertainty estimation for layer shear velocities. (e) Uncertainty estimation for layer thicknesses. A joint inversion for model HoB was then performed by merging the fundamental mode dispersion curve of the previous inversion with the true ellipticity data in the range 0.75–5.2 Hz, again after Hobiger et al. (2013). The results of this inversion are shown in Fig. 11. The improvement of the retrieved profiles is self-evident in both the similarity to the true structure of the best solution (minimum misfit profile in Fig. 11a) and in the solution appraisal (Figs 11e and f), although it can be noted that the introduction of the ellipticity ε tends to increase the posterior correlation among the layer parameters (Fig. 11e), suggesting possible scaling effects if the ellipticity is used alone to invert for the earth structure. In this particular case, the lower frequency-range of the ellipticity data are, according to Figs 6(a) and (d), quite sensitive to the model parameters of the deeper part of the model profile. Therefore, their introduction into the inversion process contributes to improve the quality of the reconstruction of the inverted model and its reliability, according to the posterior (although applied to theoretical data) analysis. Figure 11. View largeDownload slide Joint inversion of Rayleigh wave ellipticity and phase velocity of the fundamental mode for model HoA in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The black line is the true model. (b) Observed (true) and predicted ellipticity functions. (c) Observed (true) and predicted dispersion curves. (d) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (e) Uncertainty estimation for layer shear velocities. (f) Uncertainty estimation for layer thicknesses. Figure 11. View largeDownload slide Joint inversion of Rayleigh wave ellipticity and phase velocity of the fundamental mode for model HoA in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The black line is the true model. (b) Observed (true) and predicted ellipticity functions. (c) Observed (true) and predicted dispersion curves. (d) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (e) Uncertainty estimation for layer shear velocities. (f) Uncertainty estimation for layer thicknesses. These results offer compelling evidence that the key to evaluate the efficiency of ellipticity measurements to complement surface wave dispersion data lies on the resolving power of the ellipticity data with respect to the layers parameters, or, in other words, on how data are sensitive to the inversion unknowns. The absence of singularities in the ellipticity function should not be considered per se a warranty of efficiency of a joint inversion procedure. Nevertheless, the presence of discontinuity in the ellipticity function may be the index of structural complexity and high stiffness contrasts, which may be difficult to be resolved by inversion. Regarding the ‘gap’ that should be allowed between high-frequency surface wave dispersion data and lower frequency ellipticity data, Hobiger et al. (2013) stated that there is no simple rule for indicating the lowest acceptable value of this quantity, especially for complex ground structures. This is consistent with the fact that the sensitivity pattern strongly depends on the stiffness contrasts within the earth model and that the lack of sensitivity with regard to certain model parameters is not predictable with a simple rule, particularly in the case of singularities in the ellipticity function. The ellipticity data on the ‘left flank’ of the ellipticity curve are generally used to constrain the peak frequency, according to Hobiger et al. (2013), where it is suggested to neglect ellipticity values larger than 3.5. Data below the main peak may carry sensitive information which can be conveniently included into the inversion process, provided one avoids data too close to the frequency of the main peak. Nevertheless, the residual sensitivity at the left flank depends strongly on the impedance contrast within the soil column and generally carries little information about the shear velocity of the bottom half-space, despite of what is commonly thought. 5 REAL WORLD EXAMPLE We present in this section a real case study near the city of Gubbio (PG, Central Italy). The sedimentary basin underlying the Gubbio city and the surrounding area is characterized by significant site effects and has been the subject of extended seismological studies (Bindi et al.2009). The geologic bedrock in the centre of the Gubbio basin has been estimated to be deeper than 400 m. The study site is located at the Western edge of the basin, where the geologic bedrock is estimated to be shallower. This very site was investigated for near surface soil characterization by using active and passive seismic investigations in Cercato et al. (2010). In that study, the inverted earth models given by Rayleigh wave dispersion data alone did not allow the authors to match the fundamental resonant frequency determined by HVSR single station measurements. In this paper, we try to invert jointly the surface seismic surveys located in Fig. 12, consisting of a linear active array, a 2-D L-shaped passive array and HVSR measurements. The location of a 40m deep downhole, whose results are shown in Cercato et al. (2010), is also reported in the same figure. The characteristics of the seismic investigations are described in Table 2. Figure 12. View largeDownload slide Gubbio (PG): location of the seismic surveys. The dashed line indicates the 108 channel array for active surface wave prospecting, whereas the black solid line with black circles represents the 2-D L-shaped passive array with 24 low-frequency receivers. The grey circle is the 3C single station for microtremor measurements (HVSR). Figure 12. View largeDownload slide Gubbio (PG): location of the seismic surveys. The dashed line indicates the 108 channel array for active surface wave prospecting, whereas the black solid line with black circles represents the 2-D L-shaped passive array with 24 low-frequency receivers. The grey circle is the 3C single station for microtremor measurements (HVSR). Table 2. Acquisition parameters of the seismic investigations in Gubbio (PG), Central Italy. Survey  Receiver type  Geometry  Source type  Sampling  SW active  Geospace 4.5 Hz Vertical  1-D - linear  Shotgun  0.250 ms  SW passive  Lennartz 0.2 Hz Vertical  2-D L-shaped  –  8 ms  HVSR  Lennartz 0.2 Hz 3C  single station  –  8 ms  Survey  Receiver type  Geometry  Source type  Sampling  SW active  Geospace 4.5 Hz Vertical  1-D - linear  Shotgun  0.250 ms  SW passive  Lennartz 0.2 Hz Vertical  2-D L-shaped  –  8 ms  HVSR  Lennartz 0.2 Hz 3C  single station  –  8 ms  View Large The idea is to show the improvement in the subsoil reconstruction by the joint inversion of Rayleigh wave dispersion and ellipticity data, when compared to the results of conventional surface wave modal inversion (Cercato et al.2010). The HVSR measurements results are shown in Fig. 13(a). In standard HVSR processing, an anti-trigger algorithm is applied to the data to remove strong transients: the selected time windows are tapered before computing the Fourier amplitude spectrum for each selected window. These spectra are smoothed and then, for each window, the vertical component is divided by the quadratic mean of the two horizontal components to provide the HVSR of the window. The final HVSRs are obtained by averaging the results computed for each window. The result of such a standard processing is shown in Fig. 13(a) (grey solid line). In principle, the HVSR match the ellipticity of Rayleigh waves only if the noise wavefield is composed entirely by Rayleigh waves. In reality, the seismic noise wavefield contains body, Love and Rayleigh waves so that the HVSR can be overestimated in terms of ellipticity, particularly if the effect of Love waves on the horizontal component is not taken into account. Unfortunately, the proportion between Rayleigh and Love waves in the seismic noise wavefield is site and source dependent, so that their mutual proportion cannot be assumed a priori. To overcome these limitations, the RayDec (Rayleigh wave ellipticity by using the random decrement technique) software (Hobiger et al.2009) was used to estimate more accurately the experimental ellipticity. This method is capable of suppressing Love and body waves efficiently from the spectral ratios, providing an ellipticity estimation which is closer to theory than conventional HVSR results. The ellipticity curve obtained by RayDec processing is shown in Fig. 13(a) (black solid line). As expected, the estimated ellipticity curve exhibits lower amplitudes when compared to the HVSR curve, mainly because of the presence of Love waves and body waves on the horizontal components of seismic noise. Figure 13. View largeDownload slide Gubbio (PG): experimental data. (a) Microtremor results. The grey solid line is the result of the standard HVSR interpretation. The black solid line is the experimental ellipticity curve obtained by the RayDec method (Hobiger et al.2009), with the dashed thin lines indicating the standard errors. The white open circles are the data points used in the ellipticity inversion. (b) Phase velocity dispersion curves of the active (grey filled circles) and passive (black filled circles) seismic array. Figure 13. View largeDownload slide Gubbio (PG): experimental data. (a) Microtremor results. The grey solid line is the result of the standard HVSR interpretation. The black solid line is the experimental ellipticity curve obtained by the RayDec method (Hobiger et al.2009), with the dashed thin lines indicating the standard errors. The white open circles are the data points used in the ellipticity inversion. (b) Phase velocity dispersion curves of the active (grey filled circles) and passive (black filled circles) seismic array. The surface wave dispersion data obtained from the linear active array and the 2-D L-shaped array (Fig. 12) are merged into the fundamental-mode dispersion curve displayed in Fig. 13(b). In a first attempt, the active dispersion data (grey solid points in Fig. 13b) and the ellipticity data (open circles in Fig. 13a) are jointly inverted using the VFSA algorithm. In a second inversion run, we also added to the inversion the passive dispersion data (black solid points in Fig. 13b) obtained by the 2-D L-shaped array. The ellipticity data are selected excluding the large amplitude values around the main peak, following the indications emerged in the sensitivity evaluations on synthetic models. A few ellipticity points on the ‘left flank’ of the ellipticity curve are included to constrain the peak frequency (Hobiger et al.2013). Data below the main peak may carry sensitive information which can be conveniently included into the inversion process, provided one avoids data too close to the frequency of the main peak. Nevertheless, the residual sensitivity at the left flank depends strongly on the impedance contrast within the soil column and generally carries little information about the shear velocity of the bottom half-space, despite of what is commonly thought. Other data choices are clearly possible, including for instance the high-frequency part of the ellipticity function in the inversion problem. We do not report the results with alternative data selection for brevity and because they do not add much to our considerations, being the main focus of the proposed joint inversion the possibility to resolve the deep part of the soil velocity profile adding ellipticity data to the surface wave dispersion measurements. In fact, the reported inversions are meant to illustrate the capability of ellipticity data to resolve deep structures and to estimate the effects of the frequency gap between ellipticity and dispersion data on the inversion results. The same model space (which is shown in Table 3) is selected for both inversions. Table 3. Model parameter space for the joint inversion at the Gubbio (PG) test site. Layer  α (m s−1)  β (m s−1)  ρ (kg m−3)  h (m)  1  1200  100–200  2000  5–15  2  1630  150–400  2000  5–15  3  1650  250–500  2000  10–30  4  2000  250–1000  2000  15–45  5  2500  400–1000  2100  20–60  6  2600  600–1600  2200  30–90  7  2700  800–2000  2300  ∞  Layer  α (m s−1)  β (m s−1)  ρ (kg m−3)  h (m)  1  1200  100–200  2000  5–15  2  1630  150–400  2000  5–15  3  1650  250–500  2000  10–30  4  2000  250–1000  2000  15–45  5  2500  400–1000  2100  20–60  6  2600  600–1600  2200  30–90  7  2700  800–2000  2300  ∞  View Large The results of the first inversion run are shown in Fig. 14. The layer structure is well resolved down to about 150 m of depth, although the shear wave velocities of the bottom three layers exhibits quite large posterior uncertainty variations (Fig. 14e) when compared to the near surface layers. The results of the second inversion run, encompassing the 2-D-array passive dispersion measurements, are shown in Fig. 15. There is satisfactory agreement between the results of the two inversions (Fig. 14a and Fig. 15a). The minimum misfit models are consistent in the layer structure as well as in the velocity values. If we look at the solution appraisal (Figs 14d–f and 15d–f), it is quite evident that the passive dispersion data from the passive array contribute to a better resolution of the lower part of the profile, especially for the intermediate layers (from layer 4 to layer 6). This concurs well with the sensitivity analysis shown in Fig. 16, where we compare the sensitivities of both ellipticity and phase velocity computed for the two minimum misfit models of the above inversions (namely the red solid line in Figs 14a and 15a, respectively). In the frequency range 1.5–4 Hz, the ‘gap’ between the ellipticity and the active dispersion data in Fig. 14, the shear velocities of the intermediate layers in the minimum misfit model exhibit large sensitivities (Figs 16c and e) and the same considerations fits well to the sensitivities of the minimum misfit model of the ‘complete data’ inversion (Figs 16d and f). The frequency gap clearly affects the resolving power of the data with respect to both the thicknesses and the shear wave velocities of the intermediate layers, so that this lack of sensitivity of the data in the first inversion is responsible of an increased posterior uncertainty of the layers above the half-space. Nevertheless, the inverted soil structure is well consistent between the two inversions, indicating the effectiveness of incorporating ellipticity data along with phase-velocity dispersion values in a joint inversion algorithm. Figure 14. View largeDownload slide Gubbio (PG), inversion #1. Joint inversion of active (linear array) surface wave dispersion data and the ellipticity data shown in Fig. 13. (a) Mean model and minimum misfit model obtained by very fast simulated annealing inversion. (b) Observed and predicted dispersion curves. (c) Standard correlation matrix after very fast simulated annealing inversion. The first seven layer parameters correspond to layer shear velocities while the last six (8–13) correspond to layer thicknesses. (d) Uncertainty estimation for layer shear velocities after very fast simulated annealing inversion. (e) Uncertainty estimation for layer thicknesses after very fast simulated annealing inversion. Figure 14. View largeDownload slide Gubbio (PG), inversion #1. Joint inversion of active (linear array) surface wave dispersion data and the ellipticity data shown in Fig. 13. (a) Mean model and minimum misfit model obtained by very fast simulated annealing inversion. (b) Observed and predicted dispersion curves. (c) Standard correlation matrix after very fast simulated annealing inversion. The first seven layer parameters correspond to layer shear velocities while the last six (8–13) correspond to layer thicknesses. (d) Uncertainty estimation for layer shear velocities after very fast simulated annealing inversion. (e) Uncertainty estimation for layer thicknesses after very fast simulated annealing inversion. Figure 15. View largeDownload slide Gubbio (PG), inversion #2. Joint inversion of active (linear array) and passive (2-D L-shaped array) surface wave dispersion data and the ellipticity data shown in Fig. 13. (a) Mean model and minimum misfit model obtained by very fast simulated annealing inversion. (b) Observed and predicted dispersion curves. (c) Standard correlation matrix after very fast simulated annealing inversion. The first seven layer parameters correspond to layer shear velocities while the last six (8–13) correspond to layer thicknesses. d) Uncertainty estimation for layer shear velocities after very fast simulated annealing inversion. (e) Uncertainty estimation for layer thicknesses after very fast simulated annealing inversion. Figure 15. View largeDownload slide Gubbio (PG), inversion #2. Joint inversion of active (linear array) and passive (2-D L-shaped array) surface wave dispersion data and the ellipticity data shown in Fig. 13. (a) Mean model and minimum misfit model obtained by very fast simulated annealing inversion. (b) Observed and predicted dispersion curves. (c) Standard correlation matrix after very fast simulated annealing inversion. The first seven layer parameters correspond to layer shear velocities while the last six (8–13) correspond to layer thicknesses. d) Uncertainty estimation for layer shear velocities after very fast simulated annealing inversion. (e) Uncertainty estimation for layer thicknesses after very fast simulated annealing inversion. Figure 16. View largeDownload slide Gubbio (PG). Sensitivity analysis for the minimum misfit model (MM) of the inversions. MMA is the minimum misfit model in Fig. 14(a) and MMB is the minimum misfit model in Fig. 15(a). (a) MMA : phase velocity and ellipticity. The experimental data are the open white circles. (b) Model MMB: phase velocity and ellipticity. The experimental data are the open white circles. Results appearing in panel (a) for model MMA are redrawn here for comparison. (c) Model MMA: sensitivity of the phase velocity with respect to layer S-wave velocity. (d) Model MMB: sensitivity of the phase velocity with respect to layer S-wave velocity. (e) Model MMA: sensitivity of the ellipticity function with respect to layer S-wave velocity. (f) Model MMB: sensitivity of the ellipticity function with respect to layer S-wave velocity. Figure 16. View largeDownload slide Gubbio (PG). Sensitivity analysis for the minimum misfit model (MM) of the inversions. MMA is the minimum misfit model in Fig. 14(a) and MMB is the minimum misfit model in Fig. 15(a). (a) MMA : phase velocity and ellipticity. The experimental data are the open white circles. (b) Model MMB: phase velocity and ellipticity. The experimental data are the open white circles. Results appearing in panel (a) for model MMA are redrawn here for comparison. (c) Model MMA: sensitivity of the phase velocity with respect to layer S-wave velocity. (d) Model MMB: sensitivity of the phase velocity with respect to layer S-wave velocity. (e) Model MMA: sensitivity of the ellipticity function with respect to layer S-wave velocity. (f) Model MMB: sensitivity of the ellipticity function with respect to layer S-wave velocity. It is worth noting that our active array is quite large, so that the frequency gap between active and passive dispersion data is confined between 1.5 Hz and −4 Hz which is a gap that does not prevent to resolve the entire profile down to the half-space. The possibility to retrieve moderate-to-low frequency data employing active arrays depends on the array geometry, on the seismic source and on the soil structure (and therefore it is definitely site-dependent), so the frequency range of the dispersion data may be hardly estimated a priori. Inclusion of larger high-frequency parts of the experimental ellipticity may compensate for such a lack of frequency on the active data, although the frequency content of the seismic noise does not generally allow the interpreter to constrain the upper portion of the soil profile and the sensitivity of high-frequency ellipticity-data to the surface layers is generally low when compared to the sensitivity of phase velocity, as shown in the synthetic examples in Section 3. 6 DISCUSSION In some recent papers (Boaga et al.2013; Castellaro 2016) the possibility of using the HVSRs and the ellipticity function as an index of mode superposition and structural complexity (the presence of LVL, high stiffness contrasts, etc.) was analysed suggesting an active intervention of the interpreter (mainly derived from visual inspection and experience) to constrain the inversion process. The sensitivity analysis introduced in this paper constitute a general approach for understanding the relation and the diagnostic potential of the coupled inversion of Rayleigh wave ellipticity and surface wave dispersion, so that the ellipticity information can be efficiently introduced into the inversion process, limiting the interpreter's bias and the application of subjective choices to the inversion approach, stressing generality and reproducibility of the inversion workflow. In the synthetic inversion presented for models HoA and HoB after Hobiger et al. (2013) as well as in the real case study of Gubbio (PG), we adopted the RMS relative error as the Energy function for the VFSA algorithm, so that the inversion algorithm is not modified between synthetic and real applications. In a previous paper Cercato (2011), I analysed the possibility to employ compound energy functions with different misfit terms and constrains for surface wave inversion. In particular, smoothing terms and minimum norm weights were introduced to drive the inversion process to geologically meaningful (although biased) solutions. In principle, this can be also effective for the joint inversion of Rayleigh wave dispersion and ellipticity. In addition to that, a common approach is to use the experimental uncertainty to weight the data in the misfit (energy) function of the inversion process, so that the contribution of those data affected by large uncertainty is reduced. Being this paper focused on the sensitivity of the ellipticity function with respect to the layer parameters, I chose not to introduce the effect of weights and uncertainty in our inversion applications. From a general point of view, dispersion and ellipticity data may exhibit very different experimental uncertainty. Likewise, the uncertainty affecting active or passive surface wave dispersion data may be very different, so that the optimal choice of the energy function, being data-driven, may not be straightforward and of general validity. The propagation of data uncertainty into the inversion process may be the subject of further research and could not be thoroughly analysed in this work, in order to maintain its focus on the sensitivity analysis of the ellipticity function more than its experimental uncertainty, although these aspect are certainly related. Another issue which affects the inversion process is the number of data points to be inverted and the frequency sampling employed. We maintained the same logarithmic sampling for the ellipticity and the dispersion curves in each case study. In the author's opinion, this is a reasonable and ‘fair’ choice to avoid the misestimation of one kind of data on another, or to oversample certain frequency ranges. Due to the redundant type of acquisition measurements, the surface layers in the real case study of Gubbio (PG) might be considered as oversampled, but this has a very limited impact in practice on the seismic response at the site. We also noted that the introduction of ellipticity data has the effect of increasing both the posterior correlation between layers as well as the overall misfit of the inverted models. This may be due to the fact that ellipticity is inherently responsible of scaling effects due to its definition as a displacement ratio, whose high sensitivity with respect to layer parameters may be responsible of more difficulties when fitting the experimental data. In the author's opinion, at this stage of the research these aspects cannot be generalized directly, being clearly data-dependent, so that the picture is still incomplete regarding this aspects, which may be the subject of future work. In near surface seismic characterization, the main focus is on seismic velocities, especially shear wave velocity, so that our derivation of the Rayleigh wave eigenproblem employs the seismic velocities as independent variables defining the layer parameters (eq. 2.1) in the computation of the secular function (eq. 2.6) and of the ellipticity function (eq. 2.7). Alternative choices can be made, since the stress-strain relation expressed by the generalized Hooke's law for an isotropic elastic material can be written using different couples of elastic constants, such as, for example, Young's modulus and Poisson's ratio. In particular, Poisson's ratio has been described as an important parameter to define the ellipticity shape for simple models (Malischewsky & Scherbaum 2004; Tuan et al.2011), and an analysis of the sensitivity of the ellipticity function with respect to the Poisson's ratio may be useful in many applications and may be the subject of future studies. As described in Section 1, an alternative array-based approach based on the ellipticity angle (Maranò et al.2012, 2017) has been recently developed to invert for the subsoil structure. The representation through the ellipticity angle is capable of handling the singularities near the peaks and toughs of the ellipticity function, which is a distinctive advantage over the conventional HVSR representation. This method is based on multicomponent array estimation and it is not directly applicable to single-station measurements. As far as ellipticity inversion is concerned, in this paper the main focus is to investigate, through a sensitivity analysis, the effectiveness of the coupled inversion between single-station HVSR (ellipticity) data and small-aperture vertical-receiver dispersion data. For this reason, the ellipticity angle has little application in our case, although the method is very promising to improve the current state-of-the-art of ellipticity inversion. 7 CONCLUSIONS In this paper, a new method was presented for computing the sensitivity of the Rayleigh wave ellipticity with respect to the structural parameters of a layered elastic half-space. This method is analytical, it is based on the implicit function theorem and on the partition of the Rayleigh wave eigenproblem into subdeterminants of order two. This method is fast and accurate when compared to numerical differentiation, especially at high frequencies and it is fully analytical in its formulation. The new approach was tested to highlight the link between the sensitivity of the ellipticity function and its shape, as well as the complementarity of Rayleigh wave dispersion and ellipticity data when inverted to infer the subsoil structure. We found indeed that being the ellipticity a ratio between two displacements, its sensitivity to certain layer parameters may be excessively high within certain frequency ranges, particularly in the presence of singularities in the ellipticity function. This localized high sensitivity is not of any practical use in the inversion, because, as for model HoA, may be the source of potential pitfalls in the inversion process. This substantiates previous findings in the literature (Hobiger et al.2013), suggesting to avoid the inclusion of the ellipticity values around the main peak into the inversion process. When compared to the phase-velocity of surface waves, the sensitivity of the ellipticity function exhibits a more complex behaviour and it is not well distributed along frequencies. Its amplitude is generally higher when compared to phase velocity and may be excessively sensitive to perturbations of certain layer parameters within restricted frequency ranges. This behaviour strongly substantiates previous findings (Scherbaum et al.2003) that ellipticity alone is affected by inherent non-uniqueness due to its definition as a ratio which implies inevitably scaling effects. The complementarity between surface wave dispersion and ellipticity data was explored by synthetic examples and a real case-study: when used in conjunction with surface wave data, the ellipticity of Rayleigh waves may add valuable information to the inversion process provided, that the interpreter avoids including ellipticity data around the main peaks of the ellipticity function, in order to maintain the inversion process stable. Understanding the sensitivity of ellipticity data to the soil structure is a key point to integrate meaningfully surface wave dispersion data with ellipticity. To avoid the deployment of large aperture arrays, the use of active surface wave data and dispersion data from conventional linear arrays of relatively small aperture may be fruitful, provided that the gap between the low-frequency ellipticity data and the high-frequency surface wave data does not cause a lack of sensitivity for large portion of the soil structure, causing the impossibility to resolve certain layer parameters by the inversion process. ACKNOWLEDGEMENTS Manuel Hobiger of ETH Zurich is thanked for providing the original RayDec code to extract the ellipticity curve from microtremor data and for accurately reviewing this manuscript. The Gubbio data were acquired with colleagues Giuliano Milana (INGV, Rome), Giuseppe Di Giulio (INGV, Rome), Fabrizio Cara (INGV, Rome) and Salomon Hailemikael (ENEA, Rome) within the research projects DPC-S4 supported by the Italian Civil Protection Agency (Dipartimento di Protezione Civile, DPC) and the Italian Institute of Geophysics (Istituto Nazionale di Geofisica e Vulcanologia INGV). REFERENCES Arai H., Tokimatsu K., 2004. S-wave velocity profiling by inversion of microtremor H/V spectrum, Bull. seism. Soc. Am.  94 53– 63. Google Scholar CrossRef Search ADS   Arai H., Tokimatsu K., 2005. S-wave velocity profiling by joint inversion of microtremor dispersion curve and horizontal-to-vertical (H/V) spectrum, Bull. seism. Soc. Am.  95 1766– 1778. Google Scholar CrossRef Search ADS   Bhattacharya B.B., Shalivahan Sen M.K., 2003. Use of VFSA for resolution, sensitivity and uncertainty analysis in 1D DC resistivity and IP inversion, Geophys. Prospect.  51 393– 408. Google Scholar CrossRef Search ADS   Bindi D. et al.  , 2007. Site amplifications observed in the Gubbio Basin, Central Italy: hints for lateral propagation effects, Bull. seism. Soc. Am.  99 741– 760. Google Scholar CrossRef Search ADS   Boaga J., Cassiani G., Strobbia C.L., Vignoli G., 2013. Mode misidentification in Rayleigh waves: ellipticity as a cause and a cure, Geophysics  78 EN17– EN28. Google Scholar CrossRef Search ADS   Boore D.M., Toksöz N., 1969. Rayleigh wave particle motion and crustal structure, Bull. seism. Soc. Am.  59 331– 346. Castellaro S., 2016. The complementarity of H/V and dispersion curves, Geophysics  81 T323– T338. Google Scholar CrossRef Search ADS   CEN (European Committee for Standardization), 2004. EN 1998-1:2004. Eurocode 8: Design of Structures for Earthquake Resistance—Part 1: General Rules, Seismic Actions and Rules for Buildings . CEN, Brussels. Cercato M., 2007. Computation of partial derivative of Rayleigh-wave phase velocity using second order subdeterminants, Geophys. J. Int.  170 217– 238. Google Scholar CrossRef Search ADS   Cercato M., 2009. Addressing non-uniqueness in linearized multichannel surface wave inversion, Geophys. Prospect.  57 27– 47. Google Scholar CrossRef Search ADS   Cercato M., 2011. Global surface wave inversion with model constraints, Geophys. Prospect.  59 210– 226. Google Scholar CrossRef Search ADS   Cercato M., Cara F., Cardarelli E., di Filippo G., Di Giulio G., Milana G., 2010. Shear-wave velocity profiling at sites with high stiffness contrasts: a comparison between invasive and non-invasive methods, Near Surf. Geophys.  8 75– 94. Dunkin J.W., 1965. Computation of modal solutions in layered, elastic media at high frequencies, Bull. seism. Soc. Am.  55 335– 358. Fäh D., Kind F., Giardini D., 2001. A theoretical investigation of average H/V ratios, Geophys. J. Int.  145 535– 549. Google Scholar CrossRef Search ADS   Harkrider D.G., 1970. Surface waves in multilayered elastic media. Part II. Higher mode spectra and spectral ratio from point sources in plane layered earth models, Bull. seism. Soc. Am.  60 1937– 1987. Haskell N.A., 1964. Radiation pattern of surface waves from point sources in multi-layered medium, Bull. seism. Soc. Am.  54 377– 393. Hobiger M., Bard P.-Y., Cornou C., Le Bihan N., 2009. Single station determination of Rayleigh wave ellipticity by using the random decrement technique (RayDec), Geophys. Res. Lett.  36 L14303, doi:10.1029/2009GL038863. Google Scholar CrossRef Search ADS   Hobiger M. et al.  , 2013. Ground structure imaging by inversions of Rayleigh wave ellipticity: sensitivity analysis and application to European strong-motion sites, Geophys. J. Int.  192 207– 229. Google Scholar CrossRef Search ADS   Ingber L., 1989. Very fast simulated re-annealing, Math. Comput. Modelling.  12 967– 973. Google Scholar CrossRef Search ADS   Malischewsky P.G., Scherbaum F., 2004. Love's formula and H/V-ratio (ellipticity) of Rayleigh waves, Wave Motion  40 57– 67. Google Scholar CrossRef Search ADS   Maranò S., Reller C., Loeliger H.A., Fäh D., 2012. Seismic waves estimation and wavefield decomposition: application to ambient vibrations, Geophys. J. Int.  191 175– 188. Google Scholar CrossRef Search ADS   Maranò S., Hobiger M., Fäh D., 2017. Retrieval of Rayleigh wave ellipticity from ambient vibration recordings, Geophys. J. Int.  209 334– 352. Riley K.F., Hobson M.P., Bence S.J., 2006. Mathematical Methods for Physics and Engineering  3rd edn, Cambridge Univ. Press. Google Scholar CrossRef Search ADS   Rudin W., 1976. Principle of Mathematical Analysis  3rd edn, McGraw-Hill. Scherbaum F., Hinzen K.G., Ohrnberger M., 2003. Determination of shallow shear wave velocity profiles in the Cologne, Germany area using ambient vibrations, Geophys. J. Int.  152 597– 612. Google Scholar CrossRef Search ADS   Sen M.K., Stoffa P.L., 1996. Bayesian inference, Gibbs’ sampler and uncertainty estimation in geophysical inversion, Geophys. Prospect.  44 313– 350. Google Scholar CrossRef Search ADS   Sen M.K., Stoffa P.L., 2013. Global Optimization Methods in Geophysical Inversion  2nd edn, Cambridge Univ. Press. Google Scholar CrossRef Search ADS   Sexton J.L., Rudman A.J., Mead J., 1977. Ellipticity of Rayleigh waves recorded in the Midwest, Bull. seism. Soc. Am.  67 369– 382. Takeuchi H., Saito M., 1972. Surface waves and Earth oscillations, in Methods in Computational Physics, Vol. 11 - Surface Waves and Earth Oscillations  pp. 217– 295, ed. Bolt B.A., Academic Press. Tanimoto T., Alvizuri C., 2006. Inversion of the HZ ratio of microseism for S-wave velocity in the Crust, Geophys. J. Int.  165 323– 335. Google Scholar CrossRef Search ADS   Tanimoto T., Rivera L., 2008. The ZH ratio method for long-period seismic data: sensitivity kernels and observational techniques, Geophys. J. Int.  172 187– 198. Google Scholar CrossRef Search ADS   Tanimoto T., Tsuboi S., 2009. Variational principle for Rayleigh wave ellipticity, Geophys. J. Int.  179 1658– 1668. Google Scholar CrossRef Search ADS   Tuan T.T., Scherbaum F., Malischewsky P.G., 2011. On the relationship of peaks and troughs of the ellipticity (H/V) of Rayleigh waves and the transmission response of single layer over half-space models, Geophys. J. Int.  184 793– 800. Google Scholar CrossRef Search ADS   Tsuboi S., Saito M., 1983. Partial derivatives of Rayleigh wave particle motion, J. Phys. Earth  31 103– 116. Google Scholar CrossRef Search ADS   Wang C.Y., Herrmann R.B.H., 1980. A numerical study of P-, SV-, and SH-wave generation in a plane layered medium, Bull. seism. Soc. Am. , 70 1015– 1036. Watson T.H., 1970. A note on fast computation of Rayleigh wave dispersion in the multilayered elastic half-space, Bull. seism. Soc. Am.  60 161– 166. Yamanaka H., Takemura M., Ishida H., Niwa M., 1994. Characteristics of long-period microtremors and their applicability in exploration of deep sedimentary layers, Bull. seism. Soc. Am.  84 1831– 1841. APPENDIX: A REVIEW OF DUNKIN'S THEORY For a finite layer in the layered half-space, dropping the layer index, we define the following quantities: If we recall that, in each homogeneous elastic layer, the vertical wavenumbers are defined as:   \begin{eqnarray} \nu _\alpha &=& {\left(k^2-\frac{\omega ^2}{\alpha ^2}\right)^\frac{1}{2}} \end{eqnarray} (A.1)  \begin{eqnarray} \nu _\beta &=& {\left(k^2-\frac{\omega ^2}{\beta ^2}\right)^\frac{1}{2}} \end{eqnarray} (A.2)with α and β being respectively the P-wave and the S-wave layer velocities, introducing:   \begin{eqnarray} m_\alpha =\sqrt{\left|k^2-\frac{\omega ^2}{\alpha ^2}\right|}=\left\lbrace \begin{array}{ll}\sqrt{\frac{\omega ^2}{c^2}-\frac{\omega ^2}{\alpha ^2}}\qquad & \qquad \mathrm{if} \qquad c<\alpha \\ \sqrt{\frac{\omega ^2}{\alpha ^2}-\frac{\omega ^2}{c^2}}\qquad & \qquad \mathrm{if} \qquad c>\alpha \\ \end{array}\right. \end{eqnarray} (A.3)  \begin{eqnarray} m_\beta =\sqrt{\left|k^2-\frac{\omega ^2}{\beta ^2}\right|}=\left\lbrace \begin{array}{ll}\sqrt{\frac{\omega ^2}{c^2}-\frac{\omega ^2}{\beta ^2}}\qquad & \qquad \mathrm{if} \qquad c<\beta \\ \sqrt{\frac{\omega ^2}{\beta ^2}-\frac{\omega ^2}{c^2}}\qquad & \qquad \mathrm{if} \qquad c>\beta \\ \end{array}\right. \end{eqnarray} (A.4)we can derive expression of the vertical wavenumbers which are real valued or pure imaginary. In the bottom half-space, due to the radiation condition:   \begin{eqnarray} \nu _\alpha = m_\alpha \end{eqnarray} (A.5)  \begin{eqnarray} \nu _\beta = m_\beta \end{eqnarray} (A.6)In the other layers:   \begin{eqnarray} \nu _\alpha =\left\lbrace \begin{array}{ll}m_\alpha \qquad & \qquad \mathrm{if} \qquad c<\alpha \\ i{\,\,}m_\alpha \qquad & \qquad \mathrm{if} \qquad c>\alpha \\ \end{array}\right. \end{eqnarray} (A.7)  \begin{eqnarray} \nu _\beta =\left\lbrace \begin{array}{ll}m_\beta \qquad & \qquad \mathrm{if} \qquad c<\beta \\ i{\,\,}m_\beta \qquad & \qquad \mathrm{if} \qquad c>\beta \\ \end{array}\right. \end{eqnarray} (A.8)Following Dunkin (1965), we define:   \begin{eqnarray} S_\alpha = \frac{k}{\nu _\alpha }\sinh (\nu _\alpha h)=\left\lbrace \begin{array}{ll}\frac{k}{m_\alpha } \sinh (m_\alpha h)\qquad & \qquad \mathrm{if} \qquad c<\alpha \\ \frac{k}{m_\alpha } \sin (m_\alpha h)& \qquad \mathrm{if} \qquad c>\alpha \\ \end{array}\right. \end{eqnarray} (A.9)  \begin{eqnarray} S_\beta = \frac{k}{\nu _\beta }\sinh (\nu _\beta h)=\left\lbrace \begin{array}{ll}\frac{k}{m_\beta } \sinh (m_\beta h)\qquad & \qquad \mathrm{if} \qquad c<\beta \\ \frac{k}{m_\beta } \sin (m_\beta h)& \qquad \mathrm{if} \qquad c>\beta \\ \end{array}\right. \end{eqnarray} (A.10)  \begin{eqnarray} C_\alpha = \cosh (\nu _\alpha h)=\left\lbrace \begin{array}{ll}\cosh (m_\alpha h)\qquad & \qquad \mathrm{if} \qquad c<\alpha \\ \cos (m_\alpha h)\qquad & \qquad \mathrm{if} \qquad c>\alpha \\ \end{array}\right. \end{eqnarray} (A.11)  \begin{eqnarray} C_\beta = \cosh (\nu _\beta h)=\left\lbrace \begin{array}{ll}\cosh (m_\beta h)\qquad & \qquad \mathrm{if} \qquad c<\beta \\ \cos (m_\beta h)\qquad & \qquad \mathrm{if} \qquad c>\beta \\ \end{array}\right. \end{eqnarray} (A.12)  \begin{eqnarray} \gamma =\frac{2\beta ^2}{c^2} \end{eqnarray} (A.13)  \begin{eqnarray} \bar{\nu }_\alpha =\frac{\nu _\alpha }{k} \end{eqnarray} (A.14)  \begin{eqnarray} \bar{\nu }_\beta =\frac{\nu _\beta }{k} \end{eqnarray} (A.15)  \begin{eqnarray} l = 2k^2 - \frac{\omega ^2}{\beta ^2}. \end{eqnarray} (A.16)Relationships between the parameters appearing in our notation and the ones in the original paper by Dunkin are summarized in Table A1. Table A1. Table linking our notation with Dunkin's. The D subscript is intended to indicate the parameter as defined in the original paper. ξD = k  sD = iω  λD = iωk−1  γD = −γ  hD = να  kD = νβ  $$\bar{h}_D=\bar{\nu }_\alpha$$  $$\bar{k}_D=\bar{\nu }_\beta$$  CHD = Cα  CKD = Cβ  SHD = Sα  SKD = Sβ  $$l_D=2k^2-\frac{\omega ^2}{\beta ^2}=k^2+\nu _\beta ^2=l$$  ξD = k  sD = iω  λD = iωk−1  γD = −γ  hD = να  kD = νβ  $$\bar{h}_D=\bar{\nu }_\alpha$$  $$\bar{k}_D=\bar{\nu }_\beta$$  CHD = Cα  CKD = Cβ  SHD = Sα  SKD = Sβ  $$l_D=2k^2-\frac{\omega ^2}{\beta ^2}=k^2+\nu _\beta ^2=l$$  View Large We hereby report the second-order subdeterminants of matrix $$\mathbf {\bar{T}}$$ (inverse of the bottom layer matrix) and of layer matrices G which are of use in forward calculations to correct for a couple of misprints in the original version appearing in Cercato (2007):   $$\bar{t}\mid ^{12}_{12}\,= \,\frac{\beta ^4}{4{\,\,}\omega ^4}\left(\frac{l^2}{\nu _\alpha \nu _\beta }-\frac{4\omega ^2}{c^2}\right)$$ (A.17)  $$\bar{t}\mid ^{12}_{13}\,=\, -\frac{1}{4\rho\, \omega ^2\nu _\beta }$$ (A.18)  $$\bar{t}\mid ^{12}_{14}\,=\, \bar{t}\mid ^{12}_{23} =i\frac{\beta ^2}{4{\,\,}\rho {\,\,}\omega ^3c}\left(\frac{l}{\nu _\alpha \nu _\beta }-2\right)$$ (A.19)  $$\bar{t}\mid ^{12}_{24}\,=\,\frac{1}{4\rho {\,\,}\omega ^2{\,\,}\nu _\alpha }$$ (A.20)  $$\bar{t}\mid ^{12}_{34}\,=\,\frac{1}{4\rho ^2{\,\,}\omega ^4} \left(\frac{\omega ^2}{c^2{\,\,}\nu _\alpha \nu _\beta }-1\right)\!.$$ (A.21)The second-order layer-matrix subdeterminants are eqs (A29)–(A33) to eq. (A48) in Cercato (2007):   $$g\mid ^{12}_{12}\,=\,g\mid ^{34}_{34}\,=\,2\gamma (1-\gamma )+(2\gamma ^2-2\gamma +1)C_\alpha C_\beta -\left[(1-\gamma )^2+\gamma ^2\bar{\nu }_\alpha ^2\bar{\nu }_\beta ^2\right]S_\alpha S_\beta$$ (A.22)  $$g\mid ^{12}_{13}\,=\,g\mid ^{24}_{34}\,=\, \frac{1}{\rho \omega c}\left(C_\alpha S_\beta -\bar{\nu }_\alpha ^2S_\alpha C_\beta \right)$$ (A.23)  $$g\mid ^{12}_{14}\,=\,g\mid ^{12}_{23}\,=\,g\mid ^{14}_{34}\,=\,g\mid ^{23}_{34}\,=\, i{\,\,}\frac{1}{\rho \omega c}\left[(1-2\gamma )(1-C_\alpha C_\beta )+ \left(1-\gamma -\gamma \bar{\nu }_\alpha ^2\bar{\nu }_\beta ^2\right) S_\alpha S_\beta \right]$$ (A.24)  $$g\mid ^{12}_{24}\,=\,g\mid ^{13}_{34}\,=\, \frac{1}{\rho \omega c}\left(\bar{\nu }_\beta ^2 C_\alpha S_\beta -S_\alpha C_\beta \right)$$ (A.25)  $$g\mid ^{12}_{34}\,=\,-\frac{1}{\rho ^2\omega ^2c^2}\left[2(1-C_\alpha C_\beta )+\left(1+\bar{\nu }_\alpha ^2 \bar{\nu }_\beta ^2\right)S_\alpha S_\beta \right]$$ (A.26)  $$g\mid ^{13}_{12}\,=\, g\mid ^{34}_{24}\,=\,\rho \omega c\left[ \gamma ^2 \bar{\nu }_\beta ^2 {\,\,}C_\alpha {\,\,}S_\beta -(1-\gamma )^2S_\alpha C_\beta \right]$$ (A.27)  $$g\mid ^{13}_{13}\,=\, g\mid ^{24}_{24}\,=\,C_\alpha {\,\,}C_\beta$$ (A.28)  $$g\mid ^{13}_{14}\,=\, g\mid ^{13}_{23}\,=\, g\mid ^{14}_{24}\,=\, g\mid ^{23}_{24}\,=\,i{\,\,}\left[(1-\gamma ) S_\alpha C_\beta + \gamma {\,\,}\bar{\nu }_\beta ^2{\,\,}C_\alpha S_\beta \right]$$ (A.29)  $$g\mid ^{13}_{24}\,=\, -\bar{\nu }_\beta ^2{\,\,}S_\alpha S_\beta$$ (A.30)  $$\begin{array}{l}g\mid ^{14}_{12}\,=\, g\mid ^{23}_{12}\,=\,g\mid ^{34}_{14}\,=\,g\mid ^{34}_{23}= \\ \qquad \qquad = i{\,\,}\rho \omega c\left\lbrace (3\gamma ^2-2\gamma ^3-\gamma )(1-C_\alpha C_\beta )+\left[ (1-\gamma )^3 - \gamma ^3\bar{\nu }_\alpha ^2\bar{\nu }_\beta ^2\right]S_\alpha S_\beta \right\rbrace \end{array}$$ (A.31)  $$g\mid ^{14}_{13}\,=\,g\mid ^{23}_{13}\,=\,g\mid ^{24}_{14}\,=\,g\mid ^{24}_{23}= -i{\,\,}\left[(1-\gamma ){\,\,}C_\alpha {\,\,} S_\beta +\gamma \bar{\nu }_\alpha ^2{\,\,} S_\alpha C_\beta \right]$$ (A.32)  $$g\mid ^{14}_{14}=g\mid ^{23}_{23}\,=\, 1-2{\,\,}\gamma {\,\,}(1-\gamma ){\,\,}(1-C_\alpha C_\beta )+ \left[(1-\gamma )^2 + \gamma ^2{\,\,}\bar{\nu }_\alpha ^2 {\,\,}\bar{\nu }_\beta ^2\right]{\,\,}S_\alpha S_\beta$$ (A.33)  $$g\mid ^{14}_{23}\,=\, g\mid ^{23}_{14}\,=\,g\mid ^{14}_{14}-1$$ (A.34)  $$g\mid ^{24}_{12}=g\mid ^{34}_{13}= \rho \omega c{\,\,}\left[ (1-\gamma )^2{\,\,}C_\alpha S_\beta - \gamma ^2\bar{\nu }_\alpha ^2{\,\,}S_\alpha C_\beta \right]$$ (A.35)  $$g\mid ^{24}_{13}\,=\, -\bar{\nu }_\alpha ^2{\,\,}S_\alpha S_\beta$$ (A.36)  $$g\mid ^{34}_{12}\,=\, -\rho ^2\omega ^2 c^2\left\lbrace 2\gamma ^2(1-\gamma )^2(1-C_\alpha C_\beta )+[(1-\gamma )^4+\gamma ^4\bar{\nu }_\alpha ^2\bar{\nu }_\beta ^2]S_\alpha S_\beta \right\rbrace\!.$$ (A.37) © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Sensitivity of Rayleigh wave ellipticity and implications for surface wave inversion

, Volume 213 (1) – Apr 1, 2018
22 pages

/lp/ou_press/sensitivity-of-rayleigh-wave-ellipticity-and-implications-for-surface-1qq2saIhCo
Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx558
Publisher site
See Article on Publisher Site

### Abstract

Summary The use of Rayleigh wave ellipticity has gained increasing popularity in recent years for investigating earth structures, especially for near-surface soil characterization. In spite of its widespread application, the sensitivity of the ellipticity function to the soil structure has been rarely explored in a comprehensive and systematic manner. To this end, a new analytical method is presented for computing the sensitivity of Rayleigh wave ellipticity with respect to the structural parameters of a layered elastic half-space. This method takes advantage of the minor decomposition of the surface wave eigenproblem and is numerically stable at high frequency. This numerical procedure allowed to retrieve the sensitivity for typical near surface and crustal geological scenarios, pointing out the key parameters for ellipticity interpretation under different circumstances. On this basis, a thorough analysis is performed to assess how ellipticity data can efficiently complement surface wave dispersion information in a joint inversion algorithm. The results of synthetic and real-world examples are illustrated to analyse quantitatively the diagnostic potential of the ellipticity data with respect to the soil structure, focusing on the possible sources of misinterpretation in data inversion. Inverse theory, Computational seismology, Site effects, Surface waves and free oscillations, Wave propagation 1 INTRODUCTION In many engineering and environmental problems, the shear wave velocity of near surface materials is a key parameter because it is directly linked to the ground stiffness. More specifically, when estimating the effects of the local geology on seismic amplification (the so-called ‘site effects’), a weighted average of the shear wave velocity of the upper 30 m of soil (VS30) is used in modern seismic codes (CEN 2004) as a soil classification parameter to choose the proper design spectra. Malischewsky & Scherbaum (2004) and Tuan et al. (2011) investigated extensively the behaviour of the ellipticity function for simple two-layer models, showing the relationships between the Rayleigh fundamental mode and the S-wave resonant frequency and establishing a link between the ellipticity peak and the impedance contrast between the surface layer and the bottom half-space. The inversion of Rayleigh wave ellipticity to infer the shear wave soil profile has been used for a long time now in crustal seismology (Boore & Toksöz 1969; Sexton et al.1977) and has gained new popularity in recent years for both crustal (Yamanaka et al.1994; Tanimoto & Alvizuri 2006; Tanimoto & Rivera 2008) and shallow applications (Scherbaum et al.2003; Arai & Tokimatsu 2004; Boaga et al.2013; Hobiger et al.2013; Castellaro 2016). This renewed interest, especially for near surface characterization, is mainly due to the fact that the behaviour of the ellipticity function versus frequency can be evaluated by rapid and cost-effective single-station Horizontal versus Vertical Spectral Ratio (HVSR) measurements (Fäh et al.2001; Hobiger et al.2009). Unluckily, the ellipticity function alone is not generally sufficient to determine unambiguously the soil structure because of scaling effects and intrinsic non-uniqueness (Scherbaum et al.2003; Hobiger et al.2013). Because of these limitations, ellipticity data are generally inverted jointly along with surface wave dispersion measurements to better constrain the soil structure. Maranò et al. (2012) and Maranò et al. (2017) applied newly developed array estimation techniques to pinpoint accurately the frequency of singularities using the ellipticity angle. Although the method requires large aperture, multichannel arrays, it allows to identify Rayleigh ellipticity over a broad range of frequencies even near singularities (where the particle motion is vertically or horizontally polarized) and to retrieve information about the prograde or retrograde particle motion of Rayleigh waves. A common approach for near surface characterization is to use small-aperture arrays for dispersion measurements in conjunction with single-station HVSR measurements to estimate the ellipticity at lower frequency (Scherbaum et al.2003; Arai & Tokimatsu 2005; Hobiger et al.2013). This way, the deployment of large aperture arrays to retrieve the deep soil structure can be avoided. In near surface applications, the presence of low-velocity layers (LVLs) and high stiffness contrasts within the soil column are generally responsible for mode jumps, ‘osculation points’ and higher-mode domination over specific frequency ranges (Cercato 2009; Cercato et al.2010; Cercato 2011). Misinterpretation of these complicated features can lead to pitfalls and failures in surface wave inversion. It has been recognized that low-velocity surface layers can greatly affect the Rayleigh wave particle motion (Tsuboi & Saito 1983) and that the shape of the ellipticity function can be considered as an index of mode contamination and structural anomalies such as high stiffness contrasts (Boaga et al.2013; Castellaro 2016). Hobiger et al. (2013) made a systematic study of the estimation of shear wave velocity structures by jointly inverting Rayleigh wave ellipticity and short-array Rayleigh wave dispersion measurements. They identified two classes of data to describe the main features of the observed ellipticity: data exhibiting singularities and well defined peaks or, conversely, data showing regular and smooth variation in the ellipticity function. Their study was limited to fundamental mode data only and to normally dispersive sites (where the shear wave velocity is monotonically increasing with depth). Under these limitations, they analysed the performance of the joint inversion on several synthetic simulations and experimental sites and concluded that, for ellipticity curves with singularities, the right flank of the ellipticity peak (and not the values just around the peak) carries the most valuable information on the soil structure. On the other hand, for ellipticity curves without singularities, the whole area of the fundamental peak should be inverted if such a peak exists. From a more general point of view, the diagnostic potential of the ellipticity function to invert for the soil structure can be assessed by a quantitative sensitivity analysis, which requires the computation of the derivatives of ellipticity with respect to the layer parameters. It is generally recognized that the layer thickness and the S-wave velocity have the most significant effects on the ellipticity function (Arai & Tokimatsu 2004). In previous studies, which were mostly focused on deep crustal Earth models, the derivatives of the ellipticity function were obtained numerically (Arai & Tokimatsu 2004), or employing variational theory (Tsuboi & Saito 1983; Tanimoto & Tsuboi 2009). The first method requires extensive forward modelling of the ellipticity function and is subject to numerical errors. The second is analytical in its formulation, although it requires the integration of the eigenfunctions across each layer, which was performed numerically in the cited works. The sensitivity of Rayleigh wave ellipticity with respect to the subsoil structure has been very rarely explored for near surface characterization and shallow applications. On this basis, the main objectives of this work are: To derive a new fully analytical method for computing the sensitivity of Rayleigh wave ellipticity to structural parameters in a layered elastic half-space. To perform thorough sensitivity analyses on selected Earth models using this new formulation, to analyse quantitatively the diagnostic potential of the ellipticity data with respect to the soil structure, also focusing on the possible sources of misinterpretation in data inversion. To examine and give quantitative criteria on how ellipticity data can efficiently complement surface wave dispersion information in a joint inversion algorithm, to reduce the cost and increase the depth of investigation in common applications for seismic site characterization. This paper is divided as follows: the next section (Section 2) introduces the new analytical formulation for computing the sensitivity of the ellipticity function. Successively, a few theoretical examples of sensitivity computation will be presented in Section 3. Section 4 is devoted to the analysis of the complementarity between ellipticity and surface wave data by synthetic inversion, whereas a real world application will be eventually presented in Section 5, to point out the fruitful integration between ellipticity a and surface wave dispersion data, whereas Sections 5 and 6 are devoted to Discussion and Conclusions, respectively. 2 ANALYTICAL FORMULATION FOR COMPUTING THE SENSITIVITY OF RAYLEIGH WAVE ELLIPTICITY The earth model under consideration is a 1-D elastic layered half-space (Fig. 1). This parametric model, made up by horizontal, elastic, isotropic and homogeneous layers, is bounded above by a free surface and below by a homogeneous isotropic elastic half-space of infinite thickness. For the ith soil layer, αi is P-wave velocity, βi is S-wave velocity, ρi is (volume) density and hi is layer thickness. Figure 1. View largeDownload slide The plane-layered earth model of reference. Each layer is elastic, homogeneous and isotropic. Figure 1. View largeDownload slide The plane-layered earth model of reference. Each layer is elastic, homogeneous and isotropic. The complete set of structural parameters of the earth model in Fig. 1 is:   $$\mathbf {Q}=(\alpha _1,\ldots ,\alpha _N,\beta _1,\ldots ,\beta _N,\rho _1,\ldots ,\rho _N,h_1,\ldots ,h_{N-1}).$$ (2.1)For a given frequency ω, the modal ellipticity function ε of Rayleigh waves can be computed as the amplitude of the ratio between the horizontal (Ux) and vertical (Uz) modal eigenfunctions at the ground surface (Boore & Toksöz 1969; Takeuchi & Saito 1972):   $$\varepsilon =\left|\frac{ U_x(\omega ,z=0)}{ U_z(\omega ,z=0)}\right|.$$ (2.2)This expression holds for any mode, provided that the mode exists at the given frequency ω. In eq. (2.2), the dependence of the eigenfunctions on the layer parameters appearing in eq. (2.1) is omitted for simplicity. To solve the Rayleigh wave eigenvalue problem for such a model, matrix propagator theory is frequently used (Haskell 1964; Dunkin 1965; Wang & Herrmann 1980). A key point in this class of methods is the computation of the (4 × 4) layer stack matrix R of the layered half-space as:   $$\mathbf {R} = \mathbf {\bar{T}}\mathbf {G}_{N-1}\ldots \mathbf {G}_1$$ (2.3)where $$\mathbf {\bar{T}}$$ is the inverse of the bottom half-space matrix and Gi is the propagator matrix of layer i. When performed in a straightforward way according to eq. (2.3), the computation of matrix R is affected by numerical instability at high frequencies: this problem can be addressed by employing subdeterminants of order two of the layer matrices, which are free of growing exponential terms (Dunkin 1965; Watson 1970). In our notation, a subdeterminant (or minor determinant) of order two of matrix R is indicated as:   $$r\mid^{ij}_{kl}\, =\, r_{ik}r_{jl}-r_{il}r_{jk}.$$ (2.4)For the restatement of the Rayleigh value eigenproblem in terms of order-two subdeterminants, Dunkin (1965) introduced the following theorem: Theorem 1. Given the layer stack matrix R in eq. (2.3), a generic second order subdeterminant $$r\mid ^{ij}_{kl}$$ of R can be obtained as:   $$r\mid ^{ij}_{kl}=\bar{t}_N\mid ^{ij}_{ab}g_{N-1}\mid ^{ab}_{cd}\ldots g_{1}\mid ^{ef}_{kl}$$ (2.5)where the summed pairs of indices are to be only distinct pairs of distinct indices. The analytical expressions of the subdeterminants of $$\bar{\mathbf {T}}$$ and layer propagator matrices Gi (i = 1, …, N − 1) appearing in eq. (2.5), are reported in Appendix, where the definitions of the main parameters appearing in the matrix coefficients are also introduced. As the Dunkin's notation was slightly modified in this paper for the sake of coherence with current modern notation and previous work (Cercato 2007), a link with the parameter definition appearing in the original paper (Dunkin 1965) is explicitly given in Table A1. Both the dispersion function of Rayleigh waves FR and the Rayleigh wave ellipticity ε can be expressed in terms of subdeterminants of matrix R (Dunkin 1965; Harkrider 1970) so that:   $$F_R\left[\omega ,\mathbf {Q},c(\omega ,\mathbf {Q})\right]=r\mid ^{12}_{12}=r_{11}r_{22}-r_{12}r_{21}$$ (2.6)  $$\varepsilon \left[\omega ,\mathbf {Q},c(\omega ,\mathbf {Q})\right]=\left| -\frac{r_{11}}{r_{12}}\right|=\left|\frac{r\mid ^{12}_{14}}{r\mid ^{12}_{13}}\right| \varepsilon \left[\omega ,\mathbf {Q},c(\omega ,\mathbf {Q})\right]= -\frac{r_{11}}{r_{12}}=\frac{r\mid ^{12}_{14}}{r\mid ^{12}_{13}}$$ (2.7)where c is phase velocity, ω is (circular) frequency and Q is the set of structural parameters as defined in eq. (2.1). Applying the Dunkin's theorem to the secular equation FR in eq. (2.6), one gets the well-known dispersion equation:   $$F_R(\omega,k)\,=\,r\mid ^{12}_{12}=\bar{t}_N\mid ^{12}_{ab}g_{N-1}\mid ^{ab}_{cd}\ldots g_{1}\mid^{ef}_{12}\!.$$ (2.8)In a previous paper (Cercato 2007) it was demonstrated how the partial derivatives of modal phase velocity with respect to structural parameters can be obtained taking advantage of the Implicit function theorem (Rudin 1976) as:   $$\frac{\partial c}{\partial {q_i}}=c_{q_i}=\frac{(F_R)_{q_i}}{F_R^{\prime }}=\frac{\left(r\mid ^{12}_{12}\right)_{q_i}}{\left(r\mid ^{12}_{12}\right)^{\prime }}.$$ (2.9)Following the notation in Cercato (2007), qi is used to indicate a generic parameter of layer i, so that qi can be either the P-wave velocity αi, the S-wave velocity βi, the density ρi or the layer thickness hi. Partial differentiation of a function g with respect to qi will be indicated as $$g_{q_i}$$, whereas the prime mark (i.e. g΄) will be used to indicate partial differentiation with respect to phase velocity c. This should be kept in mind for all the expressions appearing in the following. The partial derivatives of FR are obtained using two separate recursions from the half-space to the surface. One for the structural parameter qi:   $$\frac{\partial F_R}{\partial q_i}=\left(r\mid ^{12}_{12}\right)_{q_i}=\bar{t}_N\mid ^{12}_{ab}g_{N-1}\mid ^{ab}_{cd}\ldots \left(g_i\mid ^{ef}_{gh}\right)_{q_i}\ldots g_{1}\mid ^{il}_{12}$$ (2.10)and one for phase velocity c (΄ prime mark):   \begin{eqnarray} r_{i}^{\prime }\mid ^{12}_{12}&=& r_{i+1}^{\prime }\mid ^{12}_{12}g_i\mid ^{12}_{12}+r_{i+1}\mid ^{12}_{12}g_i^{\prime }\mid ^{12}_{12} +r_{i+1}^{\prime }\mid ^{12}_{13}g_i\mid ^{13}_{12}+r_{i+1}\mid ^{12}_{13}g_i^{\prime }\mid ^{13}_{12}+ r_{i+1}^{\prime }\mid ^{12}_{14}g_i\mid ^{14}_{12} \nonumber\\ && +\, r_{i+1}\mid ^{12}_{14}g_i^{\prime }\mid ^{14}_{12} +r_{i+1}^{\prime }\mid ^{12}_{23}g_i\mid ^{23}_{12}+r_{i+1}\mid ^{12}_{23}g_i^{\prime }\mid ^{23}_{12}+ r_{i+1}^{\prime }\mid ^{12}_{24}g_i\mid ^{24}_{12}+r_{i+1}\mid ^{12}_{24}g_i^{\prime }\mid ^{24}_{12} \nonumber\\ && +\, r_{i+1}^{\prime }\mid ^{12}_{34}g_i\mid ^{34}_{12}+r_{i+1}\mid ^{12}_{34}g_i^{\prime }\mid ^{34}_{12}, \end{eqnarray} (2.11)where the number of computations can be further simplified because of certain symmetries in the layer matrix subdeterminants (Watson 1970; Cercato 2007). In fact, the implementation of the recursions in eqs (2.10) and (2.11) requires the computation of a subset of matrix R subdeterminants, namely $$r_i\mid ^{12}_{12}$$, $$r_i\mid ^{12}_{13}$$, $$r_i\mid ^{12}_{14}$$, $$r_i\mid ^{12}_{24}$$ and $$r_{i}\mid ^{12}_{34}$$, to carry on the recursion layer-by-layer, together with their partial derivatives with respect to c and qi. As will be shown in the following, once that these recursions have been performed up to the free surface, the sensitivity of the ellipticity of Rayleigh waves with respect to layer parameters can be computed with little additional effort. The only difficulty is that the ellipticity ε depends on the layer parameters not only directly, but also via the intermediate variable c, the phase velocity, that changes when a single layer parameter qi is allowed to vary. Therefore, to find the total variation of ε associated with an arbitrary perturbation of a layer parameter qi, we must add the indirect dependency on phase velocity and compute the total derivative (Rudin 1976; Riley 2006) of the ellipticity with respect to the structural parameter of interest. Using a notation which is common in multivariate calculus, the total derivative with respect to the parameter qi will be indicated as d/dqi, meaning that it is not assumed that phase velocity is held constant while qi varies, in contrast with partial differentiation ∂/∂qi. Starting from eq. (2.7) and applying the rule of derivation of quotient of functions to the total derivative, one gets:   $$\frac{{\rm d} \varepsilon }{{\rm d} q_i}=\frac{\frac{{\rm d} r\mid ^{12}_{14}}{{\rm d} q_i}r\mid ^{12}_{13} -\frac{{\rm d} r\mid ^{12}_{13}}{{\rm d} q_i}r\mid ^{12}_{14}}{(r\mid ^{12}_{13})^2}= \frac{1}{r\mid ^{12}_{13}}\left(\frac{{\rm d} r\mid ^{12}_{14}}{{\rm d} q_i}-\frac{{\rm d} r\mid ^{12}_{13}}{{\rm d} q_i}\varepsilon \right).$$ (2.12)To derive the total derivatives of the subdeterminants $$r\mid ^{12}_{13}$$ and $$r\mid ^{12}_{13}$$ appearing in eq. (2.12), one can take advantage of the rule of derivation of composite functions (chain rule) so that:   $$\frac{{\rm d} r\mid ^{12}_{pq}}{{\rm d} q_i}=\frac{\partial r\mid ^{12}_{pq}}{\partial q_i}+ \frac{\partial r\mid ^{12}_{pq}}{\partial c}\frac{\partial c}{\partial q_i}\qquad pq=13,14.$$ (2.13)This way, the total derivatives of subdeterminants $$r\mid ^{12}_{13}$$ and $$r\mid ^{12}_{13}$$ with respect to qi can be expressed in terms of partial derivatives. As described above, these partial derivatives can be computed analytically according to Cercato (2007), making use of a recursive scheme from the bottom to the top of the layer stack, employing the analytical partial derivatives of a subset of layer matrix subdeterminants of order two. The total derivative of ellipticity as stated in eq. (2.12) can be directly computed with little effort, substituting the proper values of $${\rm d} r\mid ^{12}_{pq} / {\rm d} q_i$$ (pq = 13, 14) computed from the partial derivatives according to (2.13). The proposed method is analytical and exact, in the sense that the numerical accuracy depends only on machine precision. Once that the derivatives are computed, the sensitivity of Rayleigh wave ellipticity can be obtained as the normalized absolute value of the total derivative (Arai & Tokimatsu 2004, 2005):   $$S_\varepsilon (q_i)=\left|\frac{q_i}{\varepsilon }\frac{{\rm d} \varepsilon }{{\rm d} q_i}\right|$$ (2.14)which is a pure numerical quantity. 3 THEORETICAL EXAMPLES Consider a typical crustal model like Model TaN in Table1. The Rayleigh wave ellipticity and phase velocity for the fundamental mode of Rayleigh waves are displayed versus period in Fig. 2(a). The derivatives of the ellipticity function were computed analytically using the new approach. The sensitivity kernels of the ellipticity function versus depth at 0.2 Hz (period 5 s) for the P-wave and S-wave velocities and density are shown in Fig. 2(b), to be compared to fig. 2 in Tanimoto & Tsuboi (2009). To compute the sensitivity kernels according to the newly developed approach, the model layers are divided into fictitious sublayers (500m thick in this case). The computed sensitivity values are assigned to the mean depths of the fictitious layers. Figure 2. View largeDownload slide Model TaN in Table 1 after Tanimoto & Tsuboi (2009). (a) Rayleigh wave ellipticity and phase velocity for the fundamental mode of Rayleigh waves. (b) Sensitivity kernels of the ellipticity (fundamental mode) at 0.2 Hz for shear velocity, density and compressional velocity. (c) Sensitivity of the ellipticity (fundamental mode) with respect to the layer S-wave velocity. (d) Sensitivity of the ellipticity (fundamental mode) with respect to layer thickness. Figure 2. View largeDownload slide Model TaN in Table 1 after Tanimoto & Tsuboi (2009). (a) Rayleigh wave ellipticity and phase velocity for the fundamental mode of Rayleigh waves. (b) Sensitivity kernels of the ellipticity (fundamental mode) at 0.2 Hz for shear velocity, density and compressional velocity. (c) Sensitivity of the ellipticity (fundamental mode) with respect to the layer S-wave velocity. (d) Sensitivity of the ellipticity (fundamental mode) with respect to layer thickness. Table 1. Structural parameters of the three earth models used for the theoretical examples. Layer  α (m/s)  β (m/s)  ρ (kg/m−3)  h (m)  Model TaN  Three Layered model in Tanimoto & Tsuboi (2009)  1  5500  3000  2500  15000  2  6500  3600  2850  15000  3  7800  4500  3300  ∞  Model HoA  Model A in Hobiger et al.(2013)  1  540  120  1800  5  2  900  200  1800  15  3  1440  320  1800  45  4  2810  625  1800  135  5  6250  2500  2000  ∞  Model HoB  Model B in Hobiger et al.(2013)  1  540  120  1800  5  2  900  200  1800  15  3  1440  320  1800  45  4  2430  540  1800  135  5  2520  840  2000  ∞  Layer  α (m/s)  β (m/s)  ρ (kg/m−3)  h (m)  Model TaN  Three Layered model in Tanimoto & Tsuboi (2009)  1  5500  3000  2500  15000  2  6500  3600  2850  15000  3  7800  4500  3300  ∞  Model HoA  Model A in Hobiger et al.(2013)  1  540  120  1800  5  2  900  200  1800  15  3  1440  320  1800  45  4  2810  625  1800  135  5  6250  2500  2000  ∞  Model HoB  Model B in Hobiger et al.(2013)  1  540  120  1800  5  2  900  200  1800  15  3  1440  320  1800  45  4  2430  540  1800  135  5  2520  840  2000  ∞  View Large At this frequency, the sensitivity of the ellipticity function is confined to the first layer only. If we examine Figs 2(c) and (d), reporting the trend versus period of the sensitivity of the ellipticity function with respect to layer shear velocity and layer thickness, respectively, it is self-evident that at low frequency the shear velocities are the governing parameters on the sensitivity of ε. As shown also in Tanimoto & Tsuboi (2009), for such a simple model the numerical differentiation is well in agreement with the analytical results, in the frequency range of interest for crustal seismology. In this case, the advantage of the analytical approach is mainly the speed of computation. Let us now examine the two slightly more complicated and shallower earth models, named HoA and HoB after Hobiger et al. (2013), whose structural parameters are summarized in Table1. The ellipticity and the phase velocity of the fundamental mode of Rayleigh waves are displayed in Figs 3(a) and (c) for model HOA and in Figs 3(b) and (d) for model HOB, respectively. The theoretical ellipticity curve exhibits singularities (peaks and troughs) for model HoA (Fig. 3a), whereas the theoretical ellipticity curve is not discontinuous in the mathematical sense in the case of low contrasts and smooth variations of the shear wave velocities as in model HoB (Fig. 3b). Figure 3. View largeDownload slide Model HoA and Model HoB in Table 1 after Hobiger et al. (2013). Ellipticity function and phase velocity of the fundamental mode. (a) Model HoA : ellipticity function versus frequency. (b) Model HoB: ellipticity function versus frequency. (c) Model HoA: phase velocity versus frequency. (d) Model HoB: phase velocity versus frequency. Figure 3. View largeDownload slide Model HoA and Model HoB in Table 1 after Hobiger et al. (2013). Ellipticity function and phase velocity of the fundamental mode. (a) Model HoA : ellipticity function versus frequency. (b) Model HoB: ellipticity function versus frequency. (c) Model HoA: phase velocity versus frequency. (d) Model HoB: phase velocity versus frequency. We report the results of these two models for comparison to previously published results (Hobiger et al.2013) and to illustrate and support some general findings, which have been verified by many simulations on different earth models which are not shown here for the sake of brevity. The derivatives of the ellipticity function with respect to shear wave velocity were computed analytically using the new approach and their values are displayed versus frequency in Fig. 4(a) for model HoA and in Fig. 4(b) for model HoB, respectively. The behaviour of the sensitivity of the ellipticity function with respect to the layer parameters will be commented later in this paragraph. Let us analyse first the precision of the new analytical approach when compared to numerical differentiation for these two earth models. Figs 4(c) and (d) show the absolute percentage error between the analytical values of the derivatives and their numerical approximation obtained employing a two-point central finite difference (FD) formula, with a spacing equal to 1m in computing close point solutions of different shear wave velocities. The accuracy of numerical differentiation depends on the interval step, on the precision of the phase velocity modal solutions and on the choice of the FD formula (two, three points approximations or more). The analytical results are obtained at lower computational cost since all the parameters needed for calculating the derivative—as stated in eq. (2.12)—can be computed using a single recursion through the entire layer stack. The FD formulation suffers from numerical errors at low values of the derivative and at the discontinuity points, which is a potential drawback in the Jacobian calculation as employed, for instance, in linearized inversion, especially at the high frequency range of interest for engineering and exploration seismology. For near surface characterization, the analytical approach is thus superior in both precision and computational velocity (efficiency) when compared to the standard numerical approach. Figure 4. View largeDownload slide Model HoA and Model HoB in Table 1 after Hobiger et al. (2013). (a) Model HoA: derivative of the ellipticity function with respect to layer shear wave velocity. (b) Model HoB: derivative of the ellipticity function with respect to layer shear wave velocity. (c) Model HoA: absolute percentage error between the analytical and the numerical computation of the derivatives. (d) Model HoB: absolute percentage error between the analytical and the numerical computation of the derivatives. Figure 4. View largeDownload slide Model HoA and Model HoB in Table 1 after Hobiger et al. (2013). (a) Model HoA: derivative of the ellipticity function with respect to layer shear wave velocity. (b) Model HoB: derivative of the ellipticity function with respect to layer shear wave velocity. (c) Model HoA: absolute percentage error between the analytical and the numerical computation of the derivatives. (d) Model HoB: absolute percentage error between the analytical and the numerical computation of the derivatives. The sensitivity functions with respect to all the layer parameters are shown in Fig. 5 for model HoA and in Fig. 6 for model HoB, respectively. For model HoA (Fig. 5), the ellipticity function is discontinuous and the sensitivity goes virtually to infinity at the frequency value of the peaks. Therefore, around the peaks the inverse problem is highly instable, because the sensitivity of all the model parameters increases drastically. For this reasons, the experimental data points of the ellipticity function around the main peaks should be avoided in the inversion, because little variations in one of the model parameters may be responsible of large variations in the ellipticity function, adding instability to the inversion process. Figure 5. View largeDownload slide Model HoA in Table 1 after Hobiger et al. (2013). Sensitivity analysis of the fundamental mode. (a) Derivative of the ellipticity function with respect to layer shear wave velocity. (b) Derivative of the ellipticity function with respect to layer P-wave velocity. (c) Derivative of the ellipticity function with respect to layer density. (d) Derivative of the ellipticity function with respect to layer thickness. Figure 5. View largeDownload slide Model HoA in Table 1 after Hobiger et al. (2013). Sensitivity analysis of the fundamental mode. (a) Derivative of the ellipticity function with respect to layer shear wave velocity. (b) Derivative of the ellipticity function with respect to layer P-wave velocity. (c) Derivative of the ellipticity function with respect to layer density. (d) Derivative of the ellipticity function with respect to layer thickness. Figure 6. View largeDownload slide Model HoB in Table 1 after Hobiger et al. (2013). Sensitivity analysis of the ellipticity function of the fundamental of Rayleigh waves to the layer parameters. (a) Sensitivity of the ellipticity function with respect to layer S-wave velocity. (b) Sensitivity of the ellipticity function with respect to layer P-wave velocity. (c) Sensitivity of the ellipticity function with respect to layer density. (d) Sensitivity of the ellipticity function with respect to layer thickness. Figure 6. View largeDownload slide Model HoB in Table 1 after Hobiger et al. (2013). Sensitivity analysis of the ellipticity function of the fundamental of Rayleigh waves to the layer parameters. (a) Sensitivity of the ellipticity function with respect to layer S-wave velocity. (b) Sensitivity of the ellipticity function with respect to layer P-wave velocity. (c) Sensitivity of the ellipticity function with respect to layer density. (d) Sensitivity of the ellipticity function with respect to layer thickness. On the other hand, the sensitivity functions shown for model HoB (Fig. 6) depict a very different situation, with finite values of the sensitivities in the entire frequency range and extremely low values for the layer P-wave velocity α and density ρ. The sensitivity values are not evenly distributed along frequency, and only for a few cases a well-defined frequency range, where the sensitivity of a certain layer parameter is dominating the other values, is clearly notable, such as, for instance, in the 08–1.5 Hz frequency band for the shear wave velocity sensitivity of layer 4 in Fig. 6a. For both models, the sensitivities indicate that shear wave velocity and layer thickness play a major role as far as ellipticity inversion is concerned, although in the case of discontinuous ellipticity functions it is necessary to discard the ellipticity values around the peaks to maintain a well-posed inverse problem. As stated in the Introduction, among the objectives of this work, a relevant topic is to evaluate how ellipticity data can efficiently complement surface wave dispersion information. To this end, the comparison of the sensitivities of ε and c with respect to the layer parameters (especially the more relevant to the inversion process such as shear wave velocity β) may indicate quantitatively the potential of a coupled inversion when compared to the uncoupled approach. In Fig. 7 we report the sensitivities of phase velocity with respect to shear wave velocity and layer thickness. We already showed that the sensitivity of ε with respect to β, displayed in Figs 5(a) and 6(a) for model HoA and HoB, respectively, is concentrated for most layer parameters in the same frequency range. This is also true, as noted above, for model HoB. On the other hand, for both earth models, the phase velocity sensitivity (Figs 7 a and b) is evenly distributed on the layer shear velocities as far as frequency is concerned. This means, in principle, that the dispersion data provide enough sensitivity to resolve the S-wave velocity of each layer within specific frequency ranges, whereas the ellipticity function—being a fractional composite quantity according to the definition in eq. (2.2)—has a more complex sensitivity distribution along frequency. The frequency distribution of the sensitivities is evidently complicated by the presence of high stiffness contrasts as in model HoA, for both c and ε. If we compare Figs 7(a) and (b), it is self-evident how the high stiffness contrast in model HoA caused by the increased values of β for layers 4 and 5 (when compared to model HoB as in Table1) is responsible for heavy changes in the frequency distribution of the sensitivities in the 0.5–3 Hz frequency band. The same considerations are still valid for the ellipticity function (Fig. 5 a and Fig. 6a, for instance), although it should be noted that the ellipticity function is less influenced by the value of the shear wave velocity of the bottom half-space when there is a remarkable increase between the velocities of the finite-thickness layers as between layers 3 and 4 in model HoA. An analogous behaviour is shown for the sensitivities with respect to layer thickness (Figs 5d and 7c for model HoA and Figs 6d and 7d for model HoB). Figure 7. View largeDownload slide Model HoA and Model HoB in Table 1 after Hobiger et al. (2013). Sensitivity of the phase velocity for the fundamental mode. (a) Model HoA: sensitivity of the phase velocity with respect to layer S-wave velocity. (b) Model HoB: sensitivity of the phase velocity with respect to layer S-wave velocity. (c) Model HoA: sensitivity of the phase velocity with respect to layer thickness. (d) Model HoB: sensitivity of the phase velocity with respect to layer thickness. Figure 7. View largeDownload slide Model HoA and Model HoB in Table 1 after Hobiger et al. (2013). Sensitivity of the phase velocity for the fundamental mode. (a) Model HoA: sensitivity of the phase velocity with respect to layer S-wave velocity. (b) Model HoB: sensitivity of the phase velocity with respect to layer S-wave velocity. (c) Model HoA: sensitivity of the phase velocity with respect to layer thickness. (d) Model HoB: sensitivity of the phase velocity with respect to layer thickness. 4 SURFACE WAVE INVERSION AND ELLIPTICITY: SYNTHETIC INVERSION EXAMPLES On the basis of the considerations above, we illustrate in this paragraph a few synthetic inversions based on the theoretical data computed for the models in Table1. The idea is to show the impact of the introduction of ellipticity data into surface wave modal inversion of phase velocity. The Very Fast Simulated Annealing (VFSA) code (Ingber 1989; Sen & Stoffa 2013) is used to invert the data. The implementation described in Cercato (2011) has been adapted to perform the joint inversion of phase velocity and ellipticity in the VFSA inversion scheme. In the VFSA algorithm, the random models are drawn from a Cauchy like distribution which is, in turn, function of temperature (Ingber 1989). The acceptance criterion of a random model is the key of the method: random models that lower the misfit are always accepted, whereas if a random model increases the error function, it is accepted according to a probability criterion which is also temperature dependent. This is the key feature of the SA algorithm to assure that the solution is not trapped into local minima of the error function (Sen & Stoffa 2013). The energy function is defined as the relative RMS (Root mean square) proportional error:   $$E_m=\sqrt{\frac{1}{N}\sum _{j=1}^N\left|\frac{d^{{\rm obs}}_j-d^{{\rm pre}}_j}{d^{{\rm pre}}_j}\right|^2}.$$ (4.1) In eq. (4.1), for each frequency fj of the N data points, the observed data points $$d^{{\rm obs}}_j$$ and the predicted data points $$d^{{\rm pre}}_j$$ can be either phase velocity or ellipticity values. The uncertainties on the inverted solutions are obtained by importance sampling (Sen & Stoffa 1996; Bhattacharya et al.2003; Sen & Stoffa 2013) where the VFSA algorithm is used as a maximum a posteriori (MAP ) estimator of the posterior probability distribution (PPD). Sen & Stoffa (1996) showed that the models sampled by multiple VFSA runs provide a good estimate, although underrated, of uncertainty which, theoretically, should be obtained using a more numerically intensive Gibbs’ sampling procedure at the constant temperature of T = 1. We employed 20 different VFSA runs (with different starting models) for each inversion in order to guarantee the convergence of the to guarantee a faster convergence than the PPD. The number of iteration for each independent run is set to 10 000 (at least), so that, overall a complete data inversion consists of at least 200 000 sampled models. No misfit threshold is set to reject models, so that occasionally random models with large errors can be included in the accepted population depending on the size of the model parameter's space. The posterior covariance analysis on the accepted models allows the interpreter to appraise the uncertainty estimation in terms of standard deviations on the layer-parameter unknowns and their mutual correlation and interdependence via the posterior correlation matrix, as illustrated in Cercato (2011). A comparison of different methods for estimating the PPD via importance sampling is reported in Sen & Stoffa (2013), chap.8, and is well beyond the scope of this paper. In the case of synthetic inversion, the use of posterior uncertainty assessment applied to synthetic data inversion is intended to quantify the impact of ellipticity information on the solution appraisal and therefore on the overall inversion uncertainty. Among the numerous synthetic inversions on different earth models, only a couple of examples are reported, in order to illustrate the most relevant findings related to the sensitivity analysis of ellipticity, while preserving the readability of this work. In the synthetic inversions, the model space is such that each model parameter (namely shear wave velocities and thicknesses) is allowed to vary between one half and two times of the real value, mimic a quite wider choice compared with usual assumptions (Hobiger et al.2013). I chose not to add noise to the data to be consistent with previous work (Hobiger et al.2013) and to avoid the mix of sensitivity and uncertainty issues in our comparative analysis. In a first synthetic inversion considering model HoA, the fundamental mode phase velocity data is inverted in the frequency band 2–30 Hz, which is the same adopted by Hobiger et al. (2013) in their synthetic inversions. The results are shown in Fig. 8. According to Fig. 5(a), at 2Hz the sensitivity of ε with respect to shear wave velocity is extremely low, so that the velocity of the bottom half-space cannot be accurately retrieved from dispersion data only. This leads to high uncertainty in the inverted bottom-half-space velocity as in Fig. 8(d), as well as in the appraisal of the thickness of the fourth layer (Fig. 8e). Figure 8. View largeDownload slide Modal VFSA synthetic inversion for model HoA in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The solid black line is the true model. (b) Observed (true) and predicted dispersion curves. (c) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (d) Uncertainty estimation for layer shear velocities. (e) Uncertainty estimation for layer thicknesses. Figure 8. View largeDownload slide Modal VFSA synthetic inversion for model HoA in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The solid black line is the true model. (b) Observed (true) and predicted dispersion curves. (c) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (d) Uncertainty estimation for layer shear velocities. (e) Uncertainty estimation for layer thicknesses. On a successive run, I inverted the same phase velocity data along with the ellipticity curve in the range 0.9–1.7 Hz, as reported in one of the synthetic tests in Hobiger et al. (2013, fig. 3 and table 2). The ellipticity data are located between the two discontinuities of the ellipticity function, on the ‘right flank’ of the fundamental-mode curve with respect to the main (low-frequency) peak. Looking again at the sensitivity curves of the ellipticity function in Figs 5(a) and (d), it can be noted that in the frequency range picked for the inversion, only the fourth layer velocity and the third and the fourth layer thicknesses exhibit considerable sensitivities. The results of this joint inversion are shown in Fig. 9. As expected from the preliminary examination of the sensitivities functions, the comparison between the standard deviations of shear velocity and thickness (Figs 8d and e) obtained from the phase-velocity inversion with the ones (Figs 9e and f) of the joint inversion, confirms that the uncertainty on the shear velocity of the bottom half-space is virtually not influenced by the introduction of the ellipticity data. It is also worth noting that both the uncertainties on the retrieved velocity and thickness of the fourth layer are remarkably reduced by the introduction of ellipticity data, confirming how ellipticity data can fruitfully complement phase velocity information when the sensitivity of the ellipticity function is finite and not in the proximity of a discontinuity in the ellipticity function. In addition, the actual behaviour of the ε near the discontinuity points (peaks of the ellipticity functions) is badly retrieved in practice from experimental data, with high uncertainty and possible misidentification of the real H/V spectral ratio values, so that the insertion of highly uncertain and highly sensitive data into the inversion process may be the source of relevant inversion pitfalls in real data interpretation. Figure 9. View largeDownload slide Joint inversion of Rayleigh wave ellipticity and phase velocity of the fundamental mode for model HoA in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The black line is the true model. (b) Observed (true) and predicted ellipticity functions. (c) Observed (true) and predicted dispersion curves. (d) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (e) Uncertainty estimation for layer shear velocities. (f) Uncertainty estimation for layer thicknesses. Figure 9. View largeDownload slide Joint inversion of Rayleigh wave ellipticity and phase velocity of the fundamental mode for model HoA in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The black line is the true model. (b) Observed (true) and predicted ellipticity functions. (c) Observed (true) and predicted dispersion curves. (d) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (e) Uncertainty estimation for layer shear velocities. (f) Uncertainty estimation for layer thicknesses. Fig. 10 shows the results of phase-velocity modal inversion of the fundamental mode in the frequency range 2–30 Hz for model HoB, aimed again at reproducing some of the results in Hobiger et al. (2013). In this case, from Fig. 6(d) it can be noted how the sensitivity of data with respect to the shear wave velocity of the bottom half-space is, although not particularly high, fairly more consistent at the lower values of the frequency range (2–3 Hz) when compared to the sensitivity of the same parameter for model HoA (Fig. 5a). Consequently, the inverted model of the modal (phase-velocity) inversion is fairly better resolved if compared to the HoA inversion in Fig. 8, particularly as far as the deeper part of the profile is concerned. Figure 10. View largeDownload slide Modal VFSA synthetic inversion for model HoB in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The solid black line is the true model. (b) Observed (true) and predicted dispersion curves. (c) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (d) Uncertainty estimation for layer shear velocities. (e) Uncertainty estimation for layer thicknesses. Figure 10. View largeDownload slide Modal VFSA synthetic inversion for model HoB in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The solid black line is the true model. (b) Observed (true) and predicted dispersion curves. (c) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (d) Uncertainty estimation for layer shear velocities. (e) Uncertainty estimation for layer thicknesses. A joint inversion for model HoB was then performed by merging the fundamental mode dispersion curve of the previous inversion with the true ellipticity data in the range 0.75–5.2 Hz, again after Hobiger et al. (2013). The results of this inversion are shown in Fig. 11. The improvement of the retrieved profiles is self-evident in both the similarity to the true structure of the best solution (minimum misfit profile in Fig. 11a) and in the solution appraisal (Figs 11e and f), although it can be noted that the introduction of the ellipticity ε tends to increase the posterior correlation among the layer parameters (Fig. 11e), suggesting possible scaling effects if the ellipticity is used alone to invert for the earth structure. In this particular case, the lower frequency-range of the ellipticity data are, according to Figs 6(a) and (d), quite sensitive to the model parameters of the deeper part of the model profile. Therefore, their introduction into the inversion process contributes to improve the quality of the reconstruction of the inverted model and its reliability, according to the posterior (although applied to theoretical data) analysis. Figure 11. View largeDownload slide Joint inversion of Rayleigh wave ellipticity and phase velocity of the fundamental mode for model HoA in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The black line is the true model. (b) Observed (true) and predicted ellipticity functions. (c) Observed (true) and predicted dispersion curves. (d) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (e) Uncertainty estimation for layer shear velocities. (f) Uncertainty estimation for layer thicknesses. Figure 11. View largeDownload slide Joint inversion of Rayleigh wave ellipticity and phase velocity of the fundamental mode for model HoA in Table 1 after Hobiger et al. (2013). (a) Mean model and minimum misfit model. The black line is the true model. (b) Observed (true) and predicted ellipticity functions. (c) Observed (true) and predicted dispersion curves. (d) Posterior correlation matrix. The first five layer parameters correspond to layer shear velocities whereas the last four correspond to layer thicknesses. (e) Uncertainty estimation for layer shear velocities. (f) Uncertainty estimation for layer thicknesses. These results offer compelling evidence that the key to evaluate the efficiency of ellipticity measurements to complement surface wave dispersion data lies on the resolving power of the ellipticity data with respect to the layers parameters, or, in other words, on how data are sensitive to the inversion unknowns. The absence of singularities in the ellipticity function should not be considered per se a warranty of efficiency of a joint inversion procedure. Nevertheless, the presence of discontinuity in the ellipticity function may be the index of structural complexity and high stiffness contrasts, which may be difficult to be resolved by inversion. Regarding the ‘gap’ that should be allowed between high-frequency surface wave dispersion data and lower frequency ellipticity data, Hobiger et al. (2013) stated that there is no simple rule for indicating the lowest acceptable value of this quantity, especially for complex ground structures. This is consistent with the fact that the sensitivity pattern strongly depends on the stiffness contrasts within the earth model and that the lack of sensitivity with regard to certain model parameters is not predictable with a simple rule, particularly in the case of singularities in the ellipticity function. The ellipticity data on the ‘left flank’ of the ellipticity curve are generally used to constrain the peak frequency, according to Hobiger et al. (2013), where it is suggested to neglect ellipticity values larger than 3.5. Data below the main peak may carry sensitive information which can be conveniently included into the inversion process, provided one avoids data too close to the frequency of the main peak. Nevertheless, the residual sensitivity at the left flank depends strongly on the impedance contrast within the soil column and generally carries little information about the shear velocity of the bottom half-space, despite of what is commonly thought. 5 REAL WORLD EXAMPLE We present in this section a real case study near the city of Gubbio (PG, Central Italy). The sedimentary basin underlying the Gubbio city and the surrounding area is characterized by significant site effects and has been the subject of extended seismological studies (Bindi et al.2009). The geologic bedrock in the centre of the Gubbio basin has been estimated to be deeper than 400 m. The study site is located at the Western edge of the basin, where the geologic bedrock is estimated to be shallower. This very site was investigated for near surface soil characterization by using active and passive seismic investigations in Cercato et al. (2010). In that study, the inverted earth models given by Rayleigh wave dispersion data alone did not allow the authors to match the fundamental resonant frequency determined by HVSR single station measurements. In this paper, we try to invert jointly the surface seismic surveys located in Fig. 12, consisting of a linear active array, a 2-D L-shaped passive array and HVSR measurements. The location of a 40m deep downhole, whose results are shown in Cercato et al. (2010), is also reported in the same figure. The characteristics of the seismic investigations are described in Table 2. Figure 12. View largeDownload slide Gubbio (PG): location of the seismic surveys. The dashed line indicates the 108 channel array for active surface wave prospecting, whereas the black solid line with black circles represents the 2-D L-shaped passive array with 24 low-frequency receivers. The grey circle is the 3C single station for microtremor measurements (HVSR). Figure 12. View largeDownload slide Gubbio (PG): location of the seismic surveys. The dashed line indicates the 108 channel array for active surface wave prospecting, whereas the black solid line with black circles represents the 2-D L-shaped passive array with 24 low-frequency receivers. The grey circle is the 3C single station for microtremor measurements (HVSR). Table 2. Acquisition parameters of the seismic investigations in Gubbio (PG), Central Italy. Survey  Receiver type  Geometry  Source type  Sampling  SW active  Geospace 4.5 Hz Vertical  1-D - linear  Shotgun  0.250 ms  SW passive  Lennartz 0.2 Hz Vertical  2-D L-shaped  –  8 ms  HVSR  Lennartz 0.2 Hz 3C  single station  –  8 ms  Survey  Receiver type  Geometry  Source type  Sampling  SW active  Geospace 4.5 Hz Vertical  1-D - linear  Shotgun  0.250 ms  SW passive  Lennartz 0.2 Hz Vertical  2-D L-shaped  –  8 ms  HVSR  Lennartz 0.2 Hz 3C  single station  –  8 ms  View Large The idea is to show the improvement in the subsoil reconstruction by the joint inversion of Rayleigh wave dispersion and ellipticity data, when compared to the results of conventional surface wave modal inversion (Cercato et al.2010). The HVSR measurements results are shown in Fig. 13(a). In standard HVSR processing, an anti-trigger algorithm is applied to the data to remove strong transients: the selected time windows are tapered before computing the Fourier amplitude spectrum for each selected window. These spectra are smoothed and then, for each window, the vertical component is divided by the quadratic mean of the two horizontal components to provide the HVSR of the window. The final HVSRs are obtained by averaging the results computed for each window. The result of such a standard processing is shown in Fig. 13(a) (grey solid line). In principle, the HVSR match the ellipticity of Rayleigh waves only if the noise wavefield is composed entirely by Rayleigh waves. In reality, the seismic noise wavefield contains body, Love and Rayleigh waves so that the HVSR can be overestimated in terms of ellipticity, particularly if the effect of Love waves on the horizontal component is not taken into account. Unfortunately, the proportion between Rayleigh and Love waves in the seismic noise wavefield is site and source dependent, so that their mutual proportion cannot be assumed a priori. To overcome these limitations, the RayDec (Rayleigh wave ellipticity by using the random decrement technique) software (Hobiger et al.2009) was used to estimate more accurately the experimental ellipticity. This method is capable of suppressing Love and body waves efficiently from the spectral ratios, providing an ellipticity estimation which is closer to theory than conventional HVSR results. The ellipticity curve obtained by RayDec processing is shown in Fig. 13(a) (black solid line). As expected, the estimated ellipticity curve exhibits lower amplitudes when compared to the HVSR curve, mainly because of the presence of Love waves and body waves on the horizontal components of seismic noise. Figure 13. View largeDownload slide Gubbio (PG): experimental data. (a) Microtremor results. The grey solid line is the result of the standard HVSR interpretation. The black solid line is the experimental ellipticity curve obtained by the RayDec method (Hobiger et al.2009), with the dashed thin lines indicating the standard errors. The white open circles are the data points used in the ellipticity inversion. (b) Phase velocity dispersion curves of the active (grey filled circles) and passive (black filled circles) seismic array. Figure 13. View largeDownload slide Gubbio (PG): experimental data. (a) Microtremor results. The grey solid line is the result of the standard HVSR interpretation. The black solid line is the experimental ellipticity curve obtained by the RayDec method (Hobiger et al.2009), with the dashed thin lines indicating the standard errors. The white open circles are the data points used in the ellipticity inversion. (b) Phase velocity dispersion curves of the active (grey filled circles) and passive (black filled circles) seismic array. The surface wave dispersion data obtained from the linear active array and the 2-D L-shaped array (Fig. 12) are merged into the fundamental-mode dispersion curve displayed in Fig. 13(b). In a first attempt, the active dispersion data (grey solid points in Fig. 13b) and the ellipticity data (open circles in Fig. 13a) are jointly inverted using the VFSA algorithm. In a second inversion run, we also added to the inversion the passive dispersion data (black solid points in Fig. 13b) obtained by the 2-D L-shaped array. The ellipticity data are selected excluding the large amplitude values around the main peak, following the indications emerged in the sensitivity evaluations on synthetic models. A few ellipticity points on the ‘left flank’ of the ellipticity curve are included to constrain the peak frequency (Hobiger et al.2013). Data below the main peak may carry sensitive information which can be conveniently included into the inversion process, provided one avoids data too close to the frequency of the main peak. Nevertheless, the residual sensitivity at the left flank depends strongly on the impedance contrast within the soil column and generally carries little information about the shear velocity of the bottom half-space, despite of what is commonly thought. Other data choices are clearly possible, including for instance the high-frequency part of the ellipticity function in the inversion problem. We do not report the results with alternative data selection for brevity and because they do not add much to our considerations, being the main focus of the proposed joint inversion the possibility to resolve the deep part of the soil velocity profile adding ellipticity data to the surface wave dispersion measurements. In fact, the reported inversions are meant to illustrate the capability of ellipticity data to resolve deep structures and to estimate the effects of the frequency gap between ellipticity and dispersion data on the inversion results. The same model space (which is shown in Table 3) is selected for both inversions. Table 3. Model parameter space for the joint inversion at the Gubbio (PG) test site. Layer  α (m s−1)  β (m s−1)  ρ (kg m−3)  h (m)  1  1200  100–200  2000  5–15  2  1630  150–400  2000  5–15  3  1650  250–500  2000  10–30  4  2000  250–1000  2000  15–45  5  2500  400–1000  2100  20–60  6  2600  600–1600  2200  30–90  7  2700  800–2000  2300  ∞  Layer  α (m s−1)  β (m s−1)  ρ (kg m−3)  h (m)  1  1200  100–200  2000  5–15  2  1630  150–400  2000  5–15  3  1650  250–500  2000  10–30  4  2000  250–1000  2000  15–45  5  2500  400–1000  2100  20–60  6  2600  600–1600  2200  30–90  7  2700  800–2000  2300  ∞  View Large The results of the first inversion run are shown in Fig. 14. The layer structure is well resolved down to about 150 m of depth, although the shear wave velocities of the bottom three layers exhibits quite large posterior uncertainty variations (Fig. 14e) when compared to the near surface layers. The results of the second inversion run, encompassing the 2-D-array passive dispersion measurements, are shown in Fig. 15. There is satisfactory agreement between the results of the two inversions (Fig. 14a and Fig. 15a). The minimum misfit models are consistent in the layer structure as well as in the velocity values. If we look at the solution appraisal (Figs 14d–f and 15d–f), it is quite evident that the passive dispersion data from the passive array contribute to a better resolution of the lower part of the profile, especially for the intermediate layers (from layer 4 to layer 6). This concurs well with the sensitivity analysis shown in Fig. 16, where we compare the sensitivities of both ellipticity and phase velocity computed for the two minimum misfit models of the above inversions (namely the red solid line in Figs 14a and 15a, respectively). In the frequency range 1.5–4 Hz, the ‘gap’ between the ellipticity and the active dispersion data in Fig. 14, the shear velocities of the intermediate layers in the minimum misfit model exhibit large sensitivities (Figs 16c and e) and the same considerations fits well to the sensitivities of the minimum misfit model of the ‘complete data’ inversion (Figs 16d and f). The frequency gap clearly affects the resolving power of the data with respect to both the thicknesses and the shear wave velocities of the intermediate layers, so that this lack of sensitivity of the data in the first inversion is responsible of an increased posterior uncertainty of the layers above the half-space. Nevertheless, the inverted soil structure is well consistent between the two inversions, indicating the effectiveness of incorporating ellipticity data along with phase-velocity dispersion values in a joint inversion algorithm. Figure 14. View largeDownload slide Gubbio (PG), inversion #1. Joint inversion of active (linear array) surface wave dispersion data and the ellipticity data shown in Fig. 13. (a) Mean model and minimum misfit model obtained by very fast simulated annealing inversion. (b) Observed and predicted dispersion curves. (c) Standard correlation matrix after very fast simulated annealing inversion. The first seven layer parameters correspond to layer shear velocities while the last six (8–13) correspond to layer thicknesses. (d) Uncertainty estimation for layer shear velocities after very fast simulated annealing inversion. (e) Uncertainty estimation for layer thicknesses after very fast simulated annealing inversion. Figure 14. View largeDownload slide Gubbio (PG), inversion #1. Joint inversion of active (linear array) surface wave dispersion data and the ellipticity data shown in Fig. 13. (a) Mean model and minimum misfit model obtained by very fast simulated annealing inversion. (b) Observed and predicted dispersion curves. (c) Standard correlation matrix after very fast simulated annealing inversion. The first seven layer parameters correspond to layer shear velocities while the last six (8–13) correspond to layer thicknesses. (d) Uncertainty estimation for layer shear velocities after very fast simulated annealing inversion. (e) Uncertainty estimation for layer thicknesses after very fast simulated annealing inversion. Figure 15. View largeDownload slide Gubbio (PG), inversion #2. Joint inversion of active (linear array) and passive (2-D L-shaped array) surface wave dispersion data and the ellipticity data shown in Fig. 13. (a) Mean model and minimum misfit model obtained by very fast simulated annealing inversion. (b) Observed and predicted dispersion curves. (c) Standard correlation matrix after very fast simulated annealing inversion. The first seven layer parameters correspond to layer shear velocities while the last six (8–13) correspond to layer thicknesses. d) Uncertainty estimation for layer shear velocities after very fast simulated annealing inversion. (e) Uncertainty estimation for layer thicknesses after very fast simulated annealing inversion. Figure 15. View largeDownload slide Gubbio (PG), inversion #2. Joint inversion of active (linear array) and passive (2-D L-shaped array) surface wave dispersion data and the ellipticity data shown in Fig. 13. (a) Mean model and minimum misfit model obtained by very fast simulated annealing inversion. (b) Observed and predicted dispersion curves. (c) Standard correlation matrix after very fast simulated annealing inversion. The first seven layer parameters correspond to layer shear velocities while the last six (8–13) correspond to layer thicknesses. d) Uncertainty estimation for layer shear velocities after very fast simulated annealing inversion. (e) Uncertainty estimation for layer thicknesses after very fast simulated annealing inversion. Figure 16. View largeDownload slide Gubbio (PG). Sensitivity analysis for the minimum misfit model (MM) of the inversions. MMA is the minimum misfit model in Fig. 14(a) and MMB is the minimum misfit model in Fig. 15(a). (a) MMA : phase velocity and ellipticity. The experimental data are the open white circles. (b) Model MMB: phase velocity and ellipticity. The experimental data are the open white circles. Results appearing in panel (a) for model MMA are redrawn here for comparison. (c) Model MMA: sensitivity of the phase velocity with respect to layer S-wave velocity. (d) Model MMB: sensitivity of the phase velocity with respect to layer S-wave velocity. (e) Model MMA: sensitivity of the ellipticity function with respect to layer S-wave velocity. (f) Model MMB: sensitivity of the ellipticity function with respect to layer S-wave velocity. Figure 16. View largeDownload slide Gubbio (PG). Sensitivity analysis for the minimum misfit model (MM) of the inversions. MMA is the minimum misfit model in Fig. 14(a) and MMB is the minimum misfit model in Fig. 15(a). (a) MMA : phase velocity and ellipticity. The experimental data are the open white circles. (b) Model MMB: phase velocity and ellipticity. The experimental data are the open white circles. Results appearing in panel (a) for model MMA are redrawn here for comparison. (c) Model MMA: sensitivity of the phase velocity with respect to layer S-wave velocity. (d) Model MMB: sensitivity of the phase velocity with respect to layer S-wave velocity. (e) Model MMA: sensitivity of the ellipticity function with respect to layer S-wave velocity. (f) Model MMB: sensitivity of the ellipticity function with respect to layer S-wave velocity. It is worth noting that our active array is quite large, so that the frequency gap between active and passive dispersion data is confined between 1.5 Hz and −4 Hz which is a gap that does not prevent to resolve the entire profile down to the half-space. The possibility to retrieve moderate-to-low frequency data employing active arrays depends on the array geometry, on the seismic source and on the soil structure (and therefore it is definitely site-dependent), so the frequency range of the dispersion data may be hardly estimated a priori. Inclusion of larger high-frequency parts of the experimental ellipticity may compensate for such a lack of frequency on the active data, although the frequency content of the seismic noise does not generally allow the interpreter to constrain the upper portion of the soil profile and the sensitivity of high-frequency ellipticity-data to the surface layers is generally low when compared to the sensitivity of phase velocity, as shown in the synthetic examples in Section 3. 6 DISCUSSION In some recent papers (Boaga et al.2013; Castellaro 2016) the possibility of using the HVSRs and the ellipticity function as an index of mode superposition and structural complexity (the presence of LVL, high stiffness contrasts, etc.) was analysed suggesting an active intervention of the interpreter (mainly derived from visual inspection and experience) to constrain the inversion process. The sensitivity analysis introduced in this paper constitute a general approach for understanding the relation and the diagnostic potential of the coupled inversion of Rayleigh wave ellipticity and surface wave dispersion, so that the ellipticity information can be efficiently introduced into the inversion process, limiting the interpreter's bias and the application of subjective choices to the inversion approach, stressing generality and reproducibility of the inversion workflow. In the synthetic inversion presented for models HoA and HoB after Hobiger et al. (2013) as well as in the real case study of Gubbio (PG), we adopted the RMS relative error as the Energy function for the VFSA algorithm, so that the inversion algorithm is not modified between synthetic and real applications. In a previous paper Cercato (2011), I analysed the possibility to employ compound energy functions with different misfit terms and constrains for surface wave inversion. In particular, smoothing terms and minimum norm weights were introduced to drive the inversion process to geologically meaningful (although biased) solutions. In principle, this can be also effective for the joint inversion of Rayleigh wave dispersion and ellipticity. In addition to that, a common approach is to use the experimental uncertainty to weight the data in the misfit (energy) function of the inversion process, so that the contribution of those data affected by large uncertainty is reduced. Being this paper focused on the sensitivity of the ellipticity function with respect to the layer parameters, I chose not to introduce the effect of weights and uncertainty in our inversion applications. From a general point of view, dispersion and ellipticity data may exhibit very different experimental uncertainty. Likewise, the uncertainty affecting active or passive surface wave dispersion data may be very different, so that the optimal choice of the energy function, being data-driven, may not be straightforward and of general validity. The propagation of data uncertainty into the inversion process may be the subject of further research and could not be thoroughly analysed in this work, in order to maintain its focus on the sensitivity analysis of the ellipticity function more than its experimental uncertainty, although these aspect are certainly related. Another issue which affects the inversion process is the number of data points to be inverted and the frequency sampling employed. We maintained the same logarithmic sampling for the ellipticity and the dispersion curves in each case study. In the author's opinion, this is a reasonable and ‘fair’ choice to avoid the misestimation of one kind of data on another, or to oversample certain frequency ranges. Due to the redundant type of acquisition measurements, the surface layers in the real case study of Gubbio (PG) might be considered as oversampled, but this has a very limited impact in practice on the seismic response at the site. We also noted that the introduction of ellipticity data has the effect of increasing both the posterior correlation between layers as well as the overall misfit of the inverted models. This may be due to the fact that ellipticity is inherently responsible of scaling effects due to its definition as a displacement ratio, whose high sensitivity with respect to layer parameters may be responsible of more difficulties when fitting the experimental data. In the author's opinion, at this stage of the research these aspects cannot be generalized directly, being clearly data-dependent, so that the picture is still incomplete regarding this aspects, which may be the subject of future work. In near surface seismic characterization, the main focus is on seismic velocities, especially shear wave velocity, so that our derivation of the Rayleigh wave eigenproblem employs the seismic velocities as independent variables defining the layer parameters (eq. 2.1) in the computation of the secular function (eq. 2.6) and of the ellipticity function (eq. 2.7). Alternative choices can be made, since the stress-strain relation expressed by the generalized Hooke's law for an isotropic elastic material can be written using different couples of elastic constants, such as, for example, Young's modulus and Poisson's ratio. In particular, Poisson's ratio has been described as an important parameter to define the ellipticity shape for simple models (Malischewsky & Scherbaum 2004; Tuan et al.2011), and an analysis of the sensitivity of the ellipticity function with respect to the Poisson's ratio may be useful in many applications and may be the subject of future studies. As described in Section 1, an alternative array-based approach based on the ellipticity angle (Maranò et al.2012, 2017) has been recently developed to invert for the subsoil structure. The representation through the ellipticity angle is capable of handling the singularities near the peaks and toughs of the ellipticity function, which is a distinctive advantage over the conventional HVSR representation. This method is based on multicomponent array estimation and it is not directly applicable to single-station measurements. As far as ellipticity inversion is concerned, in this paper the main focus is to investigate, through a sensitivity analysis, the effectiveness of the coupled inversion between single-station HVSR (ellipticity) data and small-aperture vertical-receiver dispersion data. For this reason, the ellipticity angle has little application in our case, although the method is very promising to improve the current state-of-the-art of ellipticity inversion. 7 CONCLUSIONS In this paper, a new method was presented for computing the sensitivity of the Rayleigh wave ellipticity with respect to the structural parameters of a layered elastic half-space. This method is analytical, it is based on the implicit function theorem and on the partition of the Rayleigh wave eigenproblem into subdeterminants of order two. This method is fast and accurate when compared to numerical differentiation, especially at high frequencies and it is fully analytical in its formulation. The new approach was tested to highlight the link between the sensitivity of the ellipticity function and its shape, as well as the complementarity of Rayleigh wave dispersion and ellipticity data when inverted to infer the subsoil structure. We found indeed that being the ellipticity a ratio between two displacements, its sensitivity to certain layer parameters may be excessively high within certain frequency ranges, particularly in the presence of singularities in the ellipticity function. This localized high sensitivity is not of any practical use in the inversion, because, as for model HoA, may be the source of potential pitfalls in the inversion process. This substantiates previous findings in the literature (Hobiger et al.2013), suggesting to avoid the inclusion of the ellipticity values around the main peak into the inversion process. When compared to the phase-velocity of surface waves, the sensitivity of the ellipticity function exhibits a more complex behaviour and it is not well distributed along frequencies. Its amplitude is generally higher when compared to phase velocity and may be excessively sensitive to perturbations of certain layer parameters within restricted frequency ranges. This behaviour strongly substantiates previous findings (Scherbaum et al.2003) that ellipticity alone is affected by inherent non-uniqueness due to its definition as a ratio which implies inevitably scaling effects. The complementarity between surface wave dispersion and ellipticity data was explored by synthetic examples and a real case-study: when used in conjunction with surface wave data, the ellipticity of Rayleigh waves may add valuable information to the inversion process provided, that the interpreter avoids including ellipticity data around the main peaks of the ellipticity function, in order to maintain the inversion process stable. Understanding the sensitivity of ellipticity data to the soil structure is a key point to integrate meaningfully surface wave dispersion data with ellipticity. To avoid the deployment of large aperture arrays, the use of active surface wave data and dispersion data from conventional linear arrays of relatively small aperture may be fruitful, provided that the gap between the low-frequency ellipticity data and the high-frequency surface wave data does not cause a lack of sensitivity for large portion of the soil structure, causing the impossibility to resolve certain layer parameters by the inversion process. ACKNOWLEDGEMENTS Manuel Hobiger of ETH Zurich is thanked for providing the original RayDec code to extract the ellipticity curve from microtremor data and for accurately reviewing this manuscript. The Gubbio data were acquired with colleagues Giuliano Milana (INGV, Rome), Giuseppe Di Giulio (INGV, Rome), Fabrizio Cara (INGV, Rome) and Salomon Hailemikael (ENEA, Rome) within the research projects DPC-S4 supported by the Italian Civil Protection Agency (Dipartimento di Protezione Civile, DPC) and the Italian Institute of Geophysics (Istituto Nazionale di Geofisica e Vulcanologia INGV). REFERENCES Arai H., Tokimatsu K., 2004. S-wave velocity profiling by inversion of microtremor H/V spectrum, Bull. seism. Soc. Am.  94 53– 63. Google Scholar CrossRef Search ADS   Arai H., Tokimatsu K., 2005. S-wave velocity profiling by joint inversion of microtremor dispersion curve and horizontal-to-vertical (H/V) spectrum, Bull. seism. Soc. Am.  95 1766– 1778. Google Scholar CrossRef Search ADS   Bhattacharya B.B., Shalivahan Sen M.K., 2003. Use of VFSA for resolution, sensitivity and uncertainty analysis in 1D DC resistivity and IP inversion, Geophys. Prospect.  51 393– 408. Google Scholar CrossRef Search ADS   Bindi D. et al.  , 2007. Site amplifications observed in the Gubbio Basin, Central Italy: hints for lateral propagation effects, Bull. seism. Soc. Am.  99 741– 760. Google Scholar CrossRef Search ADS   Boaga J., Cassiani G., Strobbia C.L., Vignoli G., 2013. Mode misidentification in Rayleigh waves: ellipticity as a cause and a cure, Geophysics  78 EN17– EN28. Google Scholar CrossRef Search ADS   Boore D.M., Toksöz N., 1969. Rayleigh wave particle motion and crustal structure, Bull. seism. Soc. Am.  59 331– 346. Castellaro S., 2016. The complementarity of H/V and dispersion curves, Geophysics  81 T323– T338. Google Scholar CrossRef Search ADS   CEN (European Committee for Standardization), 2004. EN 1998-1:2004. Eurocode 8: Design of Structures for Earthquake Resistance—Part 1: General Rules, Seismic Actions and Rules for Buildings . CEN, Brussels. Cercato M., 2007. Computation of partial derivative of Rayleigh-wave phase velocity using second order subdeterminants, Geophys. J. Int.  170 217– 238. Google Scholar CrossRef Search ADS   Cercato M., 2009. Addressing non-uniqueness in linearized multichannel surface wave inversion, Geophys. Prospect.  57 27– 47. Google Scholar CrossRef Search ADS   Cercato M., 2011. Global surface wave inversion with model constraints, Geophys. Prospect.  59 210– 226. Google Scholar CrossRef Search ADS   Cercato M., Cara F., Cardarelli E., di Filippo G., Di Giulio G., Milana G., 2010. Shear-wave velocity profiling at sites with high stiffness contrasts: a comparison between invasive and non-invasive methods, Near Surf. Geophys.  8 75– 94. Dunkin J.W., 1965. Computation of modal solutions in layered, elastic media at high frequencies, Bull. seism. Soc. Am.  55 335– 358. Fäh D., Kind F., Giardini D., 2001. A theoretical investigation of average H/V ratios, Geophys. J. Int.  145 535– 549. Google Scholar CrossRef Search ADS   Harkrider D.G., 1970. Surface waves in multilayered elastic media. Part II. Higher mode spectra and spectral ratio from point sources in plane layered earth models, Bull. seism. Soc. Am.  60 1937– 1987. Haskell N.A., 1964. Radiation pattern of surface waves from point sources in multi-layered medium, Bull. seism. Soc. Am.  54 377– 393. Hobiger M., Bard P.-Y., Cornou C., Le Bihan N., 2009. Single station determination of Rayleigh wave ellipticity by using the random decrement technique (RayDec), Geophys. Res. Lett.  36 L14303, doi:10.1029/2009GL038863. Google Scholar CrossRef Search ADS   Hobiger M. et al.  , 2013. Ground structure imaging by inversions of Rayleigh wave ellipticity: sensitivity analysis and application to European strong-motion sites, Geophys. J. Int.  192 207– 229. Google Scholar CrossRef Search ADS   Ingber L., 1989. Very fast simulated re-annealing, Math. Comput. Modelling.  12 967– 973. Google Scholar CrossRef Search ADS   Malischewsky P.G., Scherbaum F., 2004. Love's formula and H/V-ratio (ellipticity) of Rayleigh waves, Wave Motion  40 57– 67. Google Scholar CrossRef Search ADS   Maranò S., Reller C., Loeliger H.A., Fäh D., 2012. Seismic waves estimation and wavefield decomposition: application to ambient vibrations, Geophys. J. Int.  191 175– 188. Google Scholar CrossRef Search ADS   Maranò S., Hobiger M., Fäh D., 2017. Retrieval of Rayleigh wave ellipticity from ambient vibration recordings, Geophys. J. Int.  209 334– 352. Riley K.F., Hobson M.P., Bence S.J., 2006. Mathematical Methods for Physics and Engineering  3rd edn, Cambridge Univ. Press. Google Scholar CrossRef Search ADS   Rudin W., 1976. Principle of Mathematical Analysis  3rd edn, McGraw-Hill. Scherbaum F., Hinzen K.G., Ohrnberger M., 2003. Determination of shallow shear wave velocity profiles in the Cologne, Germany area using ambient vibrations, Geophys. J. Int.  152 597– 612. Google Scholar CrossRef Search ADS   Sen M.K., Stoffa P.L., 1996. Bayesian inference, Gibbs’ sampler and uncertainty estimation in geophysical inversion, Geophys. Prospect.  44 313– 350. Google Scholar CrossRef Search ADS   Sen M.K., Stoffa P.L., 2013. Global Optimization Methods in Geophysical Inversion  2nd edn, Cambridge Univ. Press. Google Scholar CrossRef Search ADS   Sexton J.L., Rudman A.J., Mead J., 1977. Ellipticity of Rayleigh waves recorded in the Midwest, Bull. seism. Soc. Am.  67 369– 382. Takeuchi H., Saito M., 1972. Surface waves and Earth oscillations, in Methods in Computational Physics, Vol. 11 - Surface Waves and Earth Oscillations  pp. 217– 295, ed. Bolt B.A., Academic Press. Tanimoto T., Alvizuri C., 2006. Inversion of the HZ ratio of microseism for S-wave velocity in the Crust, Geophys. J. Int.  165 323– 335. Google Scholar CrossRef Search ADS   Tanimoto T., Rivera L., 2008. The ZH ratio method for long-period seismic data: sensitivity kernels and observational techniques, Geophys. J. Int.  172 187– 198. Google Scholar CrossRef Search ADS   Tanimoto T., Tsuboi S., 2009. Variational principle for Rayleigh wave ellipticity, Geophys. J. Int.  179 1658– 1668. Google Scholar CrossRef Search ADS   Tuan T.T., Scherbaum F., Malischewsky P.G., 2011. On the relationship of peaks and troughs of the ellipticity (H/V) of Rayleigh waves and the transmission response of single layer over half-space models, Geophys. J. Int.  184 793– 800. Google Scholar CrossRef Search ADS   Tsuboi S., Saito M., 1983. Partial derivatives of Rayleigh wave particle motion, J. Phys. Earth  31 103– 116. Google Scholar CrossRef Search ADS   Wang C.Y., Herrmann R.B.H., 1980. A numerical study of P-, SV-, and SH-wave generation in a plane layered medium, Bull. seism. Soc. Am. , 70 1015– 1036. Watson T.H., 1970. A note on fast computation of Rayleigh wave dispersion in the multilayered elastic half-space, Bull. seism. Soc. Am.  60 161– 166. Yamanaka H., Takemura M., Ishida H., Niwa M., 1994. Characteristics of long-period microtremors and their applicability in exploration of deep sedimentary layers, Bull. seism. Soc. Am.  84 1831– 1841. APPENDIX: A REVIEW OF DUNKIN'S THEORY For a finite layer in the layered half-space, dropping the layer index, we define the following quantities: If we recall that, in each homogeneous elastic layer, the vertical wavenumbers are defined as:   \begin{eqnarray} \nu _\alpha &=& {\left(k^2-\frac{\omega ^2}{\alpha ^2}\right)^\frac{1}{2}} \end{eqnarray} (A.1)  \begin{eqnarray} \nu _\beta &=& {\left(k^2-\frac{\omega ^2}{\beta ^2}\right)^\frac{1}{2}} \end{eqnarray} (A.2)with α and β being respectively the P-wave and the S-wave layer velocities, introducing:   \begin{eqnarray} m_\alpha =\sqrt{\left|k^2-\frac{\omega ^2}{\alpha ^2}\right|}=\left\lbrace \begin{array}{ll}\sqrt{\frac{\omega ^2}{c^2}-\frac{\omega ^2}{\alpha ^2}}\qquad & \qquad \mathrm{if} \qquad c<\alpha \\ \sqrt{\frac{\omega ^2}{\alpha ^2}-\frac{\omega ^2}{c^2}}\qquad & \qquad \mathrm{if} \qquad c>\alpha \\ \end{array}\right. \end{eqnarray} (A.3)  \begin{eqnarray} m_\beta =\sqrt{\left|k^2-\frac{\omega ^2}{\beta ^2}\right|}=\left\lbrace \begin{array}{ll}\sqrt{\frac{\omega ^2}{c^2}-\frac{\omega ^2}{\beta ^2}}\qquad & \qquad \mathrm{if} \qquad c<\beta \\ \sqrt{\frac{\omega ^2}{\beta ^2}-\frac{\omega ^2}{c^2}}\qquad & \qquad \mathrm{if} \qquad c>\beta \\ \end{array}\right. \end{eqnarray} (A.4)we can derive expression of the vertical wavenumbers which are real valued or pure imaginary. In the bottom half-space, due to the radiation condition:   \begin{eqnarray} \nu _\alpha = m_\alpha \end{eqnarray} (A.5)  \begin{eqnarray} \nu _\beta = m_\beta \end{eqnarray} (A.6)In the other layers:   \begin{eqnarray} \nu _\alpha =\left\lbrace \begin{array}{ll}m_\alpha \qquad & \qquad \mathrm{if} \qquad c<\alpha \\ i{\,\,}m_\alpha \qquad & \qquad \mathrm{if} \qquad c>\alpha \\ \end{array}\right. \end{eqnarray} (A.7)  \begin{eqnarray} \nu _\beta =\left\lbrace \begin{array}{ll}m_\beta \qquad & \qquad \mathrm{if} \qquad c<\beta \\ i{\,\,}m_\beta \qquad & \qquad \mathrm{if} \qquad c>\beta \\ \end{array}\right. \end{eqnarray} (A.8)Following Dunkin (1965), we define:   \begin{eqnarray} S_\alpha = \frac{k}{\nu _\alpha }\sinh (\nu _\alpha h)=\left\lbrace \begin{array}{ll}\frac{k}{m_\alpha } \sinh (m_\alpha h)\qquad & \qquad \mathrm{if} \qquad c<\alpha \\ \frac{k}{m_\alpha } \sin (m_\alpha h)& \qquad \mathrm{if} \qquad c>\alpha \\ \end{array}\right. \end{eqnarray} (A.9)  \begin{eqnarray} S_\beta = \frac{k}{\nu _\beta }\sinh (\nu _\beta h)=\left\lbrace \begin{array}{ll}\frac{k}{m_\beta } \sinh (m_\beta h)\qquad & \qquad \mathrm{if} \qquad c<\beta \\ \frac{k}{m_\beta } \sin (m_\beta h)& \qquad \mathrm{if} \qquad c>\beta \\ \end{array}\right. \end{eqnarray} (A.10)  \begin{eqnarray} C_\alpha = \cosh (\nu _\alpha h)=\left\lbrace \begin{array}{ll}\cosh (m_\alpha h)\qquad & \qquad \mathrm{if} \qquad c<\alpha \\ \cos (m_\alpha h)\qquad & \qquad \mathrm{if} \qquad c>\alpha \\ \end{array}\right. \end{eqnarray} (A.11)  \begin{eqnarray} C_\beta = \cosh (\nu _\beta h)=\left\lbrace \begin{array}{ll}\cosh (m_\beta h)\qquad & \qquad \mathrm{if} \qquad c<\beta \\ \cos (m_\beta h)\qquad & \qquad \mathrm{if} \qquad c>\beta \\ \end{array}\right. \end{eqnarray} (A.12)  \begin{eqnarray} \gamma =\frac{2\beta ^2}{c^2} \end{eqnarray} (A.13)  \begin{eqnarray} \bar{\nu }_\alpha =\frac{\nu _\alpha }{k} \end{eqnarray} (A.14)  \begin{eqnarray} \bar{\nu }_\beta =\frac{\nu _\beta }{k} \end{eqnarray} (A.15)  \begin{eqnarray} l = 2k^2 - \frac{\omega ^2}{\beta ^2}. \end{eqnarray} (A.16)Relationships between the parameters appearing in our notation and the ones in the original paper by Dunkin are summarized in Table A1. Table A1. Table linking our notation with Dunkin's. The D subscript is intended to indicate the parameter as defined in the original paper. ξD = k  sD = iω  λD = iωk−1  γD = −γ  hD = να  kD = νβ  $$\bar{h}_D=\bar{\nu }_\alpha$$  $$\bar{k}_D=\bar{\nu }_\beta$$  CHD = Cα  CKD = Cβ  SHD = Sα  SKD = Sβ  $$l_D=2k^2-\frac{\omega ^2}{\beta ^2}=k^2+\nu _\beta ^2=l$$  ξD = k  sD = iω  λD = iωk−1  γD = −γ  hD = να  kD = νβ  $$\bar{h}_D=\bar{\nu }_\alpha$$  $$\bar{k}_D=\bar{\nu }_\beta$$  CHD = Cα  CKD = Cβ  SHD = Sα  SKD = Sβ  $$l_D=2k^2-\frac{\omega ^2}{\beta ^2}=k^2+\nu _\beta ^2=l$$  View Large We hereby report the second-order subdeterminants of matrix $$\mathbf {\bar{T}}$$ (inverse of the bottom layer matrix) and of layer matrices G which are of use in forward calculations to correct for a couple of misprints in the original version appearing in Cercato (2007):   $$\bar{t}\mid ^{12}_{12}\,= \,\frac{\beta ^4}{4{\,\,}\omega ^4}\left(\frac{l^2}{\nu _\alpha \nu _\beta }-\frac{4\omega ^2}{c^2}\right)$$ (A.17)  $$\bar{t}\mid ^{12}_{13}\,=\, -\frac{1}{4\rho\, \omega ^2\nu _\beta }$$ (A.18)  $$\bar{t}\mid ^{12}_{14}\,=\, \bar{t}\mid ^{12}_{23} =i\frac{\beta ^2}{4{\,\,}\rho {\,\,}\omega ^3c}\left(\frac{l}{\nu _\alpha \nu _\beta }-2\right)$$ (A.19)  $$\bar{t}\mid ^{12}_{24}\,=\,\frac{1}{4\rho {\,\,}\omega ^2{\,\,}\nu _\alpha }$$ (A.20)  $$\bar{t}\mid ^{12}_{34}\,=\,\frac{1}{4\rho ^2{\,\,}\omega ^4} \left(\frac{\omega ^2}{c^2{\,\,}\nu _\alpha \nu _\beta }-1\right)\!.$$ (A.21)The second-order layer-matrix subdeterminants are eqs (A29)–(A33) to eq. (A48) in Cercato (2007):   $$g\mid ^{12}_{12}\,=\,g\mid ^{34}_{34}\,=\,2\gamma (1-\gamma )+(2\gamma ^2-2\gamma +1)C_\alpha C_\beta -\left[(1-\gamma )^2+\gamma ^2\bar{\nu }_\alpha ^2\bar{\nu }_\beta ^2\right]S_\alpha S_\beta$$ (A.22)  $$g\mid ^{12}_{13}\,=\,g\mid ^{24}_{34}\,=\, \frac{1}{\rho \omega c}\left(C_\alpha S_\beta -\bar{\nu }_\alpha ^2S_\alpha C_\beta \right)$$ (A.23)  $$g\mid ^{12}_{14}\,=\,g\mid ^{12}_{23}\,=\,g\mid ^{14}_{34}\,=\,g\mid ^{23}_{34}\,=\, i{\,\,}\frac{1}{\rho \omega c}\left[(1-2\gamma )(1-C_\alpha C_\beta )+ \left(1-\gamma -\gamma \bar{\nu }_\alpha ^2\bar{\nu }_\beta ^2\right) S_\alpha S_\beta \right]$$ (A.24)  $$g\mid ^{12}_{24}\,=\,g\mid ^{13}_{34}\,=\, \frac{1}{\rho \omega c}\left(\bar{\nu }_\beta ^2 C_\alpha S_\beta -S_\alpha C_\beta \right)$$ (A.25)  $$g\mid ^{12}_{34}\,=\,-\frac{1}{\rho ^2\omega ^2c^2}\left[2(1-C_\alpha C_\beta )+\left(1+\bar{\nu }_\alpha ^2 \bar{\nu }_\beta ^2\right)S_\alpha S_\beta \right]$$ (A.26)  $$g\mid ^{13}_{12}\,=\, g\mid ^{34}_{24}\,=\,\rho \omega c\left[ \gamma ^2 \bar{\nu }_\beta ^2 {\,\,}C_\alpha {\,\,}S_\beta -(1-\gamma )^2S_\alpha C_\beta \right]$$ (A.27)  $$g\mid ^{13}_{13}\,=\, g\mid ^{24}_{24}\,=\,C_\alpha {\,\,}C_\beta$$ (A.28)  $$g\mid ^{13}_{14}\,=\, g\mid ^{13}_{23}\,=\, g\mid ^{14}_{24}\,=\, g\mid ^{23}_{24}\,=\,i{\,\,}\left[(1-\gamma ) S_\alpha C_\beta + \gamma {\,\,}\bar{\nu }_\beta ^2{\,\,}C_\alpha S_\beta \right]$$ (A.29)  $$g\mid ^{13}_{24}\,=\, -\bar{\nu }_\beta ^2{\,\,}S_\alpha S_\beta$$ (A.30)  $$\begin{array}{l}g\mid ^{14}_{12}\,=\, g\mid ^{23}_{12}\,=\,g\mid ^{34}_{14}\,=\,g\mid ^{34}_{23}= \\ \qquad \qquad = i{\,\,}\rho \omega c\left\lbrace (3\gamma ^2-2\gamma ^3-\gamma )(1-C_\alpha C_\beta )+\left[ (1-\gamma )^3 - \gamma ^3\bar{\nu }_\alpha ^2\bar{\nu }_\beta ^2\right]S_\alpha S_\beta \right\rbrace \end{array}$$ (A.31)  $$g\mid ^{14}_{13}\,=\,g\mid ^{23}_{13}\,=\,g\mid ^{24}_{14}\,=\,g\mid ^{24}_{23}= -i{\,\,}\left[(1-\gamma ){\,\,}C_\alpha {\,\,} S_\beta +\gamma \bar{\nu }_\alpha ^2{\,\,} S_\alpha C_\beta \right]$$ (A.32)  $$g\mid ^{14}_{14}=g\mid ^{23}_{23}\,=\, 1-2{\,\,}\gamma {\,\,}(1-\gamma ){\,\,}(1-C_\alpha C_\beta )+ \left[(1-\gamma )^2 + \gamma ^2{\,\,}\bar{\nu }_\alpha ^2 {\,\,}\bar{\nu }_\beta ^2\right]{\,\,}S_\alpha S_\beta$$ (A.33)  $$g\mid ^{14}_{23}\,=\, g\mid ^{23}_{14}\,=\,g\mid ^{14}_{14}-1$$ (A.34)  $$g\mid ^{24}_{12}=g\mid ^{34}_{13}= \rho \omega c{\,\,}\left[ (1-\gamma )^2{\,\,}C_\alpha S_\beta - \gamma ^2\bar{\nu }_\alpha ^2{\,\,}S_\alpha C_\beta \right]$$ (A.35)  $$g\mid ^{24}_{13}\,=\, -\bar{\nu }_\alpha ^2{\,\,}S_\alpha S_\beta$$ (A.36)  $$g\mid ^{34}_{12}\,=\, -\rho ^2\omega ^2 c^2\left\lbrace 2\gamma ^2(1-\gamma )^2(1-C_\alpha C_\beta )+[(1-\gamma )^4+\gamma ^4\bar{\nu }_\alpha ^2\bar{\nu }_\beta ^2]S_\alpha S_\beta \right\rbrace\!.$$ (A.37) © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations