Rayleigh’s, Stoneley’s, and Scholte’s Interface Waves in Elastic Models Using a Boundary Element Method

interface waves for three canonical models, that is, interfaces formed by vacuum-solid, solid-solid, and liquid-solid. These interfaces excited by dynamic loads cause the emergence of Rayleigh’s, Stoneley’s, and Scholte’s waves, respectively. To perform the study, the indirect boundary element method is used, which has proved to be a powerful tool for numerical modeling of problems in elastodynamics. In essence, the method expresses the diffracted wave ﬁeld of stresses, pressures, and displacements by a boundary integral, also known as single-layer representation, whose shape can be regarded as a Fredholm’s integral representation of second kind and zero order. This representation can be considered as an exempliﬁcation of Huygens’ principle, which is equivalent to Somigliana’s representation theorem. Results in frequency domain for the three types of interfaces are presented; then, using the fourier discrete transform, we derive the results in time domain, where the emergence of interface waves is highlighted.


Introduction
The study of interface waves has always attracted the interest of the scientific community because of the importance and complexity of the waves that propagate in such interfaces.

Journal of Applied Mathematics
For example, Rayleigh's waves are one of the three types of interface waves, which travel in vacuum-solid surfaces. In isotropic solids the particle motion is elliptical and retrograde, for shallow depths, with respect to the direction of propagation Rayleigh, 1 . Today, many engineering and seismology studies e.g., 2-5 are focused on understanding Rayleigh's waves. Recent research concerning Rayleigh's waves is also carried out on nondestructive testing for detecting defects.
Stoneley's waves occur at the interface between two solids 6 . The higher energy, as well as Rayleigh's waves, is present in the interface and shows an exponential decay away from the interface. Important applications around this type of interface waves can be found in 7-10 .
Scholte's waves are presented at the interface of fluid-solid media 11-13 . Similarly, most of the energy in this type of wave is presented in the interface and decays exponentially into the solid medium and fluid one. Some applications, mainly applied to seabed, can be seen in 8, 14-17 . To model realistic problems and complex geometries, numerical methods are a good option. Methods like finite element 18 , finite difference 19 , boundary element 20, 21 , spectral and pseudospectral elements 22-24 have been extensively used.
Particularly, the boundary element method BEM has been useful to deal with interface problems. For instance, based on the BEM, a coupled model was developed to investigate the dynamic interaction between an offshore pile, a porous seabed, and seawater when subjected to the pseudo-Stoneley wave along the seabed and the seawater interface 25 . They found that the maximum pore pressure of the seabed usually occurs at the region near the interfaces between the seabed and the seawater.
Numerical modeling to simulate the propagation of acoustic and elastic waves generated by a borehole source embedded in a layered medium was formulated in terms of the boundary element technique, where Green's functions were calculated by the discrete wavenumber method. Results display Stoneley's wave reflections at the bed boundaries and show the importance of the diffraction that takes place where the borehole wall intersects the layer interfaces 26 .
Characterization of surface cracks using Rayleigh's wave excitations was dealt by an indirect boundary element method. The variations of spectral ratios between the transmitted and incident waves were studied as a function of the crack depth. They were used to design an efficient procedure for the determination of crack depths 27 .
BEM formulations have been also used to study scattering of Rayleigh's wave by cavities 28 . Moreover, BEM methods were developed to study reflection and transmission of Rayleigh's surface waves by a juncture normal to free surface, between identical or different materials 29 .
Propagation of Scholte's waves in a water-filled borehole in an anisotropic solid by a time-domain boundary element method was studied in 30 , where detailed arrival identification and interpretation of acoustic and elastic waves propagating along the fluid-solid interfaces were pointed out.
In this paper, a numerical method known as the indirect boundary element method IBEM is used to study, in frequency and time domain, the behavior of three canonical models of interface giving rise to the emergence of Rayleigh's, Stoneley's, and Scholte's waves. To validate the equations used here, we included in The appendix, a comparison between the IBEM and discrete wave number DWN , for the case of a fluid-solid interface, where an initial pressure is applied in the fluid using a Hankel's function of second kind and zero order.

Indirect Boundary Element Method Formulation
For the three interface models studied in this work, the resulting state of tractions, displacements, and pressures at any point of the models can be expressed as the sum of an incident field and a diffracted one. All interface models, vacuum-solid, solid-solid, fluid-solid, are considered as the union of two half spaces or regions. The source generates, in the region where it acts, an incident field of pressures and displacements for the fluid or an incident field of tractions and displacement for solids. In the region where the source is not applied, only diffracted fields are expected.
For comparison and validation, the results obtained by this numerical technique are compared with respect to those obtained by the DWN. This comparison was performed for the case of a fluid-solid interface using a source in the fluid, which was represented by a Hankel's function of second kind and zero order; this validation is detailed in The appendix. However, for the purpose of making comparisons between the three interface models, the source is applied in the solid S 1 , as shown in Figure 1.
Therefore, for each model, the total field for the S 1 region can be expressed as t S 1 x t oS 1 x t d S 1 x for tractions, and u S 1 x u oS 1 x u d S 1 x for displacements. The regions S 2 and F will only have diffracted wave fields, because no sources are applied in such regions. The superindexes o and d stand for incident and diffracted wave field, respectively.

Integral Representations of Diffracted Fields
The diffracted wave field of displacements and tractions for the elastic solid S 1 can be expressed by means of where u i x represents the ith component of the displacement at point x, G ij x; ξ represents Green's function, which is the displacement produced in the direction i at x due to the application of a unit force in direction j at point ξ and φ j ξ is the force density in direction j at point ξ. This integral representation can be obtained from Somigliana's identity 4 . t i x is the ith component of traction, c 1 0.5, if x tends to the boundary S "from inside" the region, c 1 −0.5 if x tends to S "from outside" the región, and c 1 0 if x is not at S. T ij x; ξ is the traction Green's function, that is, the traction in the direction i at a point x, associated to the unit vector n i x , due to the application of a unit force in the direction j at ξ on S. The 2D Green functions for unbounded spaces can be obtained in 5, 31 . Diffracted wave field for the region S 2 has similar form to 3.1 .
For the region F fluid in Figure 1 c , the diffracted fields for displacements and pressures are written as Stoneley's-wave where Ψ · represents the force density for the fluid, G F · is the Green function for the fluid, and c 2 defines the region orientation and can assume a value of −0.5, 0, or 0.5 see explanation for c 1 , given above . ρ is the mass density, and ω represents the circular frequency.

Boundary Conditions
Boundary conditions for each of the models presented in Figure 1 are set as follows: Solid-solid interface

3.5
Journal of Applied Mathematics

3.6
Equation 3.4 represents the traction state on a free surface, which promotes the emergence of Rayleigh's waves. Boundary conditions 3.5 ensure continuity between two materials with different mechanical properties and exhibit the existence of Stoneley's waves. Finally, 3.6 are the appropriate boundary conditions between an acoustic medium and an elastic solid one.

Discretization Scheme
For purposes of exemplification, the discretization procedure for the equations corresponding to the interface of Figure 1 b is illustrated, such interface is related to the emergence and propagation of Stoneley's waves, which exist in the interface between two elastic solids. Then, from the equation of continuity 3.5 , it can be said that the traction and displacement states may be expressed, respectively, as Reordering these last two equations, one has

4.3
According to the integral representations 3.1 , 4.3 can be written as where t o S 1 i and u o S 1 i represent, respectively, the stress and displacement wave fields produced by the source, both applied in the region S 1 . In general, the interface between two solids may be discretized according to Figure 2. If we assume that the force densities φ i x are constant in each boundary element that forms the surfaces of the regions S 1 and S 2 and the Gaussian integration is performed or analytical integration where Green's function is singular , then 4.4 can be rewritten as where δ ij represents the Kronecker delta and ΔS n is the length of each boundary element. Equations 4.5 represent the system of Fredholm's integral equations to be solved. Once the unknowns are found, it is possible to determine the state of tractions and displacement at any point within the regions S 1 and S 2 using 3.1 , plus the incident wave field.
For the other types of interface, it is possible to follow the same discretization scheme, applying the corresponding boundary conditions. Our integral representations can be used to handle nonflat interfaces, which is the subject of our current research.

Numerical Examples
For validation purposes we refer the reader to the appendix. There, results achieved by the IBEM are compared with those obtained by the DWN method, for the case of fluid-solid interfaces. Good agreement is seen between both methods.
In this section, numerical simulations for the three canonical interface models are developed. Flat interfaces are considered for the three models. The material properties used for the analysis are shown in Table 1; these material properties were consulted from 32-34 . α and β are the compressional and shear wave velocities ms −1 , respectively, and ρ is the material density kgm −3 . Materials with α 2670 ms −1 correspond to sandstone, while those with α 4810 ms −1 correspond to limestone. Results are described in the following paragraphs. Figure 3 presents the displacement spectra, for the three interfaces models studied, for the first receiver detailed in Figure 1. The depth to which the source is applied is H 0.05 m and the horizontal distance from the source to the receiver is 20H 1.0 m. The response corresponding to Materials 1 is shown in Figures 3 a and 3 b , while that associated with the Materials 2 are graphed in Figures 3 c and 3 d . For the analysis, a frequency increment of 150 Hz was considered and a maximum frequency of 19200 Hz was reached.
Results associated with Materials 1 show clear amplifications for vacuum-solid and water-solid interfaces for both directions of displacement. At low frequencies, these two interfaces describe similar amplitudes for frequencies lower than 1000 Hz . For frequencies close to 19000 Hz, the behavior becomes asymptotic and almost negligible. For the model formed by water and sandstone, strong variations of displacement are noted for the range of 1000 to 12000 Hz, mainly for the x 3 component. As was expected, the model formed by the S 1 and S 2 shows amplitudes of displacement that are much smaller, and its behavior describes soft patterns. From the frequency of 10000 Hz, the behavior is almost negligible in both directions.
The response obtained from Materials 2 for the three models is depicted in Figures  3 c and 3 d . It is possible to appreciate that for both directions of displacement, the spectra show soft trajectories and almost negligible from the frequency of 8000 Hz. This behavior can be attributed to the great rigidity of the bottom material limestone in comparison with the top material air, wáter, or sandstone . This stiffness not only has dominion or control in the response displacements , but also, provides a certain similarity in the spectra, and therefore, on the Rayleigh's, Stoneley's and Scholte's interface waves for these materials. Figure 4 shows synthetic seismograms of displacement for the directions x 1 and x 3 left and right, resp. , measured by 25 receivers located as depicted in Figure 1. The first receiver is located at a horizontal distance of 20H 1.0 m from the source. The other receivers are located using a distance increment of 0.04 m. For each of the interface models studied is evident the emergence of their corresponding interface waves, that is, Rayleigh's, Stoneley's, and Scholte's waves, for vacuum-solid, solid-solid, and water-solid interfaces, Figures 4 a ,  4 b , and 4 c , respectively.   Table 1 are plotted a and b , while those obtained from Materials 2 are presented in c and d .
limestone, the second is close to the shear waves velocity of limestone, and the last shows a velocity of less than the shear wave velocity of sandstone and can be associated with the propagation of Stoneley's waves. In Figure 4 c two wave fronts are highlighted. The first traveling at a velocity of 2670 ms −1 and is associated with the compressional wave velocity of sandstone. The second front corresponds to the propagation of Scholte's waves whose speed is 970 ms −1 and obviously carries a significant amount of energy, mainly in its x 3 component.
In the case of interface waves for Materials 2 of Table 1, it should be emphasized that results for this material describes similar behavior to Materials 1. In general, the amplitudes of displacement, for both components x 1 and x 3 , are lower than those obtained for Materials 1. This is due to the great rigidity of limestone compared with sandstone. In the synthetic seismograms presented in Figure 5, different wave front arrivals and their corresponding speeds Scholte's waves Water-sandstone interface  Table 1 are used for the analysis.
are highlighted. It can be seen the amount of energy that leads the different waves that propagate in the interface, manifesting itself in the amplitudes of the displacement generated. In Figure 3 c , it is important to mention the influence of fluid layer, showing a wave front that propagates at a speed of 1501 ms −1 . Figures 6 a and 6 b show the time response for the entire 2D water-sandstone interface model. To this purpose, a net of 51 × 51 receivers, spaced using a distance increment of 0.04 m, is required. Column a plots the results of pressure in the fluid and displacements in the x 1 direction for the solid, while the column b plots pressures in the fluid and displacements in the x 3 direction for the solid. This numerical simulation is shown for three different times. Rayleigh's waves  Table 1 are used for the analysis.
For the time t 0.000911 s, the source has hit the solid boundary and a diffracted wave in the fluid and reflected waves in the solid can be seen, generating the emergence of P y S wave fronts. For the time t 0.001432 s, the above mentioned waves go away from the source while the presence of interface waves is clearly evident to this time. These are the Scholte's waves and are highlighted using circles in Figures 6 a and 6 b . For the time t 0.001953 s, the propagation of interface waves is very visibly and shows a delay with respect to the P and S wave fronts in the solid. Scholte's waves for this case propagate with a velocity of 970 ms −1 .

Conclusions
In this work, we expand the use of the indirect boundary element method to study the propagation of elastic waves in vacuum-solid, solid-solid, and fluid-solid interfaces. In this  numerical technique, based on Huygens' principle and the Somigliana's representation theorem, the fields of pressures, tractions, and displacements are expressed in terms of singlelayer boundary integral equations. Green's functions for tractions and displacements, for unbounded space, were used, but they are enforced to meet the proper boundary conditions that prevail for each interface model studied.
Firstly, spectra of displacements were included, and some aspects about the behavior for the two groups of materials studied were pointed out. Therefore, a fast fourier transform algorithm was applied to obtain time responses for the three interface models. In all cases, the existence and propagation of Rayleigh, Stoneley, and Scholte's waves are manifested, highlighting the important amount of energy that they transport.
The results obtained from our numerical technique were compared with the DWN; a good agreement between the different approaches was evident. Therefore, the IBEM can be considered as a good technique to model interface problems. Complex geometries solved by means of IBEM are the target of our current research.

Appendix
This section presents the validation of the indirect boundary element method for the solution of interface problems and compares results by IBEM with those obtained by the discrete wave number. For this purpose, a fluid-solid interface has been selected. In this case, a source is applied in the fluid, which is represented by a Hankel's function of second kind and zero order. The following is a brief description of the DWN method 35 .
The discrete wave number method is one of the techniques to simulate earthquake ground motions. The seismic wave radiated from a source is expressed as a wavenumber integration 29 . The main idea of the method is to represent a source as a superposition of homogenous plane waves propagating in discrete angles. As long as the medium has no inelastic damping, the denominator of the integrand becomes zero for a particular wavenumber and, consequently, the numerical integration becomes impossible. To solve this problem, a method to incorporate a complex frequency was proposed as early as the proposal of the discrete wavenumber method itself.
The incident pulse at the fluid source can be expressed as where p 0 F x is the incident pulse at the fluid, x {x 1 , x 3 }, C ω represents a scale factor for the incident pulse, H 2 0 · is Hankel's function of second kind and zero order, ω is circular frequency, c F represents the compressional wave velocity for the fluid, and r r x is the distance from the receiver to the source. k is the wavenumber, η ω 2 /c F 2 − k 2 with Im η 0. If we express k in discrete values, then we have k n nΔk and η n ω 2 /c F 2 − k 2 n with Im η n 0. If we assume that the whole pressure and displacement field in the fluid, that is, free and diffracted field, can be expressed, respectively, by A n represents the unknown coefficient to evaluate the diffracted field in the fluid and a is the distance from the source to the elastic solid.
For the solid, we assume that the potential of displacement has the form φ B n e −ik n x 3 e −iγ n x 1 −a and ψ C n e −ik n x 3 e −iυ n x 1 −a , where γ n ω 2 /α 2 − k n 2 with Im γ n 0 and υ n ω 2 /β 2 − k n 2 with Im υ n 0. α and β are the compressional and shear wave velocities, respectively.
The displacement field for the solid can be expressed as u ∂φ/∂x 1 − ∂ψ/∂x 3 and w ∂φ/∂x 3 ∂ψ/∂x 1 . The stress field can be obtained by the well-known equation: where σ ij x is stress tensor, λ and μ are Lamé's constants, ε ij is strain tensor and δ ij is Kronecker's delta. The boundary conditions to be enforced are represented by 3.6 . Once the boundary conditions have been applied, the unknown coefficients A n , B n , and C n are obtained, the whole pressure field in the fluid is finally determined by means of A.2 .
For validation purposes, water-sandstone and water-limestone interfaces are considered for material properties see Table 1 . Figure 7 shows the spectra of pressures for the models analyzed. For all cases, the initial pressure source was generated in the fluid at a distance of 0.05 m from the elastic solid boundary. The receiver is placed at a horizontal distance of 1.0 m from the source. The frequency analysis is done considering a frequency increment of