Accurate computation of leaky modes for anomalous layered models

Current research in seismic exploration or seismological engineering has been focused in shallow subsurface, which usually consists of unconsolidated sediments or weathered bedrock, or other types of anomalous medium properties characterized by high Poisson’s ratios. In such cases leaky waves may occur either as noise needed to be suppressed or as useful information to reveal the underlying seismic structure. Analysis may be undertaken in terms of leaky modes, which have been conventionally computed based on Haskell matrix method that suffers inherent shortcoming of loss of precision at high frequencies, or have been calculated in the real-number domain but with crude accuracy. An alternative method is proposed herein based on the generalized reflection/transmission coefficients, which searches the roots in the complex wavenumber domain and calculates the leaky modes naturally. This method is then tested for several typical anomalous layered models, and has been proved to be accurate and stable, and thus superior over the traditional methods. A powerful and efficient analyzing tool is thus offered and may be promisingly applied in seismic prospecting, geotechnical engineering, structural dynamics and other relevant areas.


Introduction
Earth's surface is often covered by unconsolidated sediments characterized by high Poisson's ratios [Stümpel et al. 1984, Roth et al. 1998, Gao et al. 2014], and thus the layer boundaries are commonly associated with a strong velocity contrast.Leaky waves, which can be regarded as the superposition of leaky modes [Haddon 1984], are therefore evident in the early arrivals [Roth et al. 1998, Gao et al. 2014]; in certain cases, slowly-attenuating leaky waves of phase velocities faster than S waves in the half-space may exist, which correspond to real frequency and real wavenumber pairs [García-Jerez and Sánchez-Sesma 2015].
The presence of leaky waves may either cause difficulty in correctly identifying dispersion curves of nor-mal modes so that the final inverted S-wave velocity model could deviate from the true model [Gao et al. 2014], or constructively afford useful information for characterizing P-wave velocity structure [Su and Dorman 1965] or for performing full-waveform inversions [Roth et al. 1998], and also help fill in the gaps of Rayleigh dispersion curves of normal modes for models with the half-space being a low-velocity layer [Yang and Yi 2005].
Investigations on leaky waves by computing dispersion of leaky modes has been performed by many researchers [Gilbert 1964, Su and Dorman 1965, Watson 1972, Dainty 1971, Kennett 1983], who largely adopted Haskell's method [Haskell 1953[Haskell , 1962] ] or its variants.However, this method suffers the well-known inherent numerical instability, as has been explained by Wu and Chen [2016].Failure to accurately find leaky modes cannot guarantee an effective dispersion analysis [Roth et al. 1998].On the contrary, accurate and complete computation of dispersion curves of leaky modes may bring better agreement between the forward calculation and the inverted dispersion image.In this paper, we shall show how to accurately calculate leaky modes for anomalous layered models, which are commonly encountered in exploration geophysics or geotechnical engineering.

Method
The highly efficient as well as stable and accurate method by Wu and Chen [2016], which is based on the generalized reflection/transmission [Luco andApsel 1983, Chen 1993], is slightly modified and utilized here to be able to find the complex roots of a secular function constructed according to the boundary conditions.As leaky modes are usually complex roots of a secular function, the only modification that we made here is allowing the wavenumber k to be complex, assuming the frequency ω is real.The phase velocity c of a leaky mode is thus defined as c=ω/Re(k).For the multivalued functions of vertical wavenumbers corresponding to P and S waves respectively where N is the number of the layers above the halfspace for a horizontally stratified model, we are mainly concerned with the roots defined on the (+,-) Riemann sheet [De Bremaecker 1967, Watson 1972] with Re(γ)>0 and Re(ν)<0 so that the computed phase velocities are primarily between the S-wave velocity b (N+1) and the P-wave velocity α (N+1) of the bottom layer of the model.
Therefore, we are required to find the complex roots k of the dispersion equation where are 2x2 submatrices of the normalized layer matrix for the first layer, is the downward generalized reflection coefficient associated with the lower boundary of the first layer, and containing phase factors of upgoing waves, is a submatrix of the diagonal matrix Λ 1 (z 0 ) for the first layer (see the detailed derivation and explicit expressions for these quantities in Wu et al. 2016).Note here the only one dispersion equation ( 2) is used to find the leaky modes, while the concept of the family of secular functions may have to be employed as in Wu and Chen [2016] to find all the physically existent normal modes for a model with low-velocity layers.
When a suitable estimated solution is given, the Newton-Raphson method [Press et al. 1992] is then used to simultaneously locate the real and the imaginary parts of the complex roots k to the given precision very quickly.Thus, the leaky modes at each frequency can be found efficiently and accurately.The priori estimate for the initial solution for the next frequency can be obtained from the roots found for the frequency just before.The correctness of our algorithm was verified by comparing our results with those by Watson [1972], who adopted the matrix method based on Haskell's method described by Watson [1970].

Models and Results
We then first apply our algorithm to the twolayer model with high Poisson's ratios (Table 1), which was used by Roth et al. [1998], who only showed possible regions of leaky modes on the dispersion diagram because the leaky roots were found rather crudely.Shown in Figure 1 are the normal modes (solid lines) and leaky modes (open circles) computed by our method for the model.The fundamental and many higher-order leaky modes are accurately computed with the root precision up to 1E-8 and clearly displayed on the dispersion diagram.
Next, we conduct the computation of leaky modes for the three-layer model by Gao et al. [2014] which contains a low-velocity half-space.The normal modes (represented by solid dots in Figure 2) are absent in the frequency range from ~0.6 Hz to ~22 Hz.
The leaky modes (denoted by open circles in Figure 2) act as an extension of normal modes that they naturally fill the gap of the dispersion curve composed of normal modes.Besides, many higher-order dispersion curves of leaky modes are also calculated.It should be noted that, the "gap" here lies between the dispersion curves of normal modes of the same order, and is different from that between the dispersion curves of normal modes and those of leaky modes for Love waves, a phenomenon observed by Malischewsky [1985].The former arises because the S-wave velocity of the half-space is smaller than those of the above layers, and the latter may be filled by the socalled open modes [Budden 1961, Malischewsky 1985].Lastly, we have examined the characteristic re- det  1.  3 with the velocities decreasing with depth, which was adopted by Yang and Yi [2005].The presence of a lowvelocity layer sandwiched in a model or the case that the medium under the ground becomes softer with depth is usually encountered in Rayleigh-wave exploration in geotechnical engineering.The gap of the dispersion curves of normal modes may make it difficult to appropriately determine the exploration depth of Rayleigh waves [Zhang and Lu 2003].Figure 3 shows only the first dispersion curve consisting of fundamental normal modes (solid dots) and leaky modes (open circles), with which engineers are mainly concerned, because the excitation energy of these modes would be dominant for a surface source.For such an anomalous model that the velocities decrease with depth, the computed dispersion curve behaves as a monotonically nondecreasing function with frequency.

Discussion
For the three cases under study our method produces more accurate results compared to those of the previous researchers.For models typically with high Poisson's ratios such as that shown in Table 1, leaky waves may show significant energy on the seismogram recorded at proper distances [Roth et al. 1998, Gao et al. 2014].Dispersion analysis in terms of mode theory would then afford complementary information for seismic surveying, especially for inverting the P-wave velocity structure [Oliver 1964, Su and Dorman 1965, Znak and Kashtan 2015] due to it being better sensed by leaky modes than by shear normal modes.Therefore, accurate computation of leaky modes is necessary, as done here, for it may also have great implications for high-resolution seismic surveying and full waveform inversions [Roth et al. 1998].
In exploration geophysics and geotechnical engineering, anomalous models having low-velocity layers as well as high Poisson's ratios are frequently encountered.Gao et al. [2014] used the model shown in Table 2 to illustrate that leaky waves would appear with their phase velocities being larger than the maximum S-wave velocity of the model.They searched the leaky modes for the model in the real-number domain and their results may not be expected to be accurate (cf. Figure 7e in their paper).In addition, the search for the leaky modes was not complete, and only the fundamental and some of the first higherorder modes of the dispersion curves were computed.When our algorithm has been applied to the model, multiple dispersion curves (Figure 2) are generated and the computation of these modes is also more accurate.
To fill the gap of the dispersion curves due to the presence of a low-velocity half-space in a model, Yang and Yi [2005] eschewed direct computation of the   3.
complex roots associated with leaky modes by appending beneath the half-space of the original model a layer of the same medium properties as those of the layer whose shear wave velocity is the greatest in all the layers to ensure that the (N+1)th layer is sufficiently thick.The effect is that the computation of the complex leaky-mode solutions is converted into that of the real normal-mode roots of the secular equation for the new model.Because of adding a layer to the original model and thickening the (N+1)th layer, the resultant modes are greatly increased and the time consumption accordingly larger.In our method the computation of leaky modes is directly performed in complex wavenumber domain, which agrees with the physical characteristics of the problem.The gap for the first dispersion curve of the fundamental normal modes is naturally bridged by the computed leaky modes, which then meet the normal modes at higher frequencies (see Figures 2 and 3).

Conclusions
We have extended the highly efficient algorithm based on the generalized reflection/transmission coefficients for computing normal modes to the putation of the leaky modes in complex wavenumber domain.Numerical stability and accuracy as well as computing efficiency are guaranteed.For the anomalous layered models with high Poisson's ratios or lowvelocity layers underlain by high-velocity surficial layers considered in this paper, leaky modes lying on many dispersion branches are computed with higher accuracy compared with previously published results.With this effective and efficient tool provided, various potential applications in geophysics, geotechnical engineering and structural dynamics would then be made feasible and soar in the near future.5 COMPUTATION OF LEAKY MODES Wu, B., and Chen, X. (2016)

Figure 1 .
Figure 1.Dispersion curves of normal modes (solid lines) and leaky modes (open circles) for the model shown in Table1.
soil foundation model shown in Table

Figure 3 .
Figure 3. Dispersion curves of normal modes (solid dots) and leaky modes (open circles) for the model shown in Table3.

Table 1 .
Parameters of a two-layer model with high Poisson's ratio.

Table 2 .
Parameters of a three-layer model with the half-space being a low-velocity layer.Dispersion curves of normal modes (solid dots) and leaky modes (open circles) for the model shown in Table2.

Table 3 .
Parameters of a soil foundation model with the velocities decreasing with depth.