Diffraction of Elastic Waves in Fluid-Layered Solid Interfaces by an Integral Formulation

1 Escuela Superior de Ingenieŕıa Mecánica y Eléctrica, Instituto Politécnico Nacional, Unidad Profesional Adolfo López Mateos s/n, 07738 México, DF, Mexico 2 Instituto Mexicano del Petróleo, Eje Central Lázaro Cárdenas 152, Gustavo A Madero, 07730 México, DF, Mexico 3 Escuela Superior de Ingenieŕıa Mecánica y Eléctrica, Instituto Politécnico Nacional, Unidad Profesional Azcapotzalco, 02250 México, DF, Mexico


Introduction
In many areas of physics, the study of fluid-solid interfaces has always attracted interest.For example, important developments to study the dynamic behavior of an oceanic layer over an elastic solid by means of analytical solutions can be seen in the original works of Biot [1] and Ewing et al. [2], where the attention to Stoneley and Rayleigh waves was paid.Other applications have been aimed to understand the behavior of interface waves being focused to ocean bottom [3,4].Carcione and Helle [5] studied the physics related to these interfaces in a variety of seabed mechanical properties, from soft sediments to crustal rocks.Analytical results to show the appearance of Rayleigh waves in oceanic ambient excited by deep earthquakes were presented, for instance, in [6,7].
Attenuation and dispersion of interface waves were investigated in multilayered cases [8][9][10][11][12].The inverse problem of determining the mechanical parameters of layered media in contact with fluids by measuring the variation of the pressure fields in the fluid was published by Zein et al. [13].Studies applied to porous media have evidenced the huge influence of porosity in wave propagation, particularly, when the medium is partially saturated [14][15][16][17][18]. Attenuation and dispersion in interface waves due to the presence of fractures were studied in [11,[19][20][21][22].The use of Green's functions for layered acoustic and elastic formations in 3D was applied in Tadeu and António [23,24], which can be used in numerical modeling.
In the field of numerical methods such as finite element and finite difference methods, the solution of fluid-solid interfaces was carried out, for instance, in Zienkiewicz and Bettess [25] and Thomas et al. [26], respectively.In addition, spectral element and pseudospectral formulations have revealed to be accurate methods for modeling realistic geometries (see, e.g., Carcione et al. [27]).On the other hand, Komatitsch et al. [28] developed the spectral element method to deal with more complex problems and numerous advantages over classical approaches were remarked.
In this communication, the indirect boundary element method (IBEM) to study the interactions between acoustic and elastic waves, near a fluid-and elastic-layered solid interface, is applied.Monopole point source, characterized by Hankel's function of second kind, is employed to produce an initial pressure wave in the fluid.This formulation could be considered as a numerical realization of Huygens' principle, in which the diffracted waves are built at the boundary from which they are radiated.Mathematically speaking, this is completely equivalent to the well-known Somigliana representation theorem.The accuracy of our results for a fluidsolid interface is verified with respect to those obtained by means of the discrete wave number (DWN).Observations previously described by other authors are highlighted.In the following sections, the formulations of both DWN and IBEM applied to fluid-solid interfaces are detailed.

Formulation of the Problem by Means of the Discrete Wave Number
The DWN method is one of the techniques to simulate seismic motion.The seismic waves radiated from a source could be expressed as an integral in the wavenumber domain.Moreover, the source is represented as a superposition of homogenous plane waves propagating in discrete angles (see, e.g., Bouchon and Aki [29]).In this method, when denominator of the integrand becomes zero for a particular wavenumber, then the numerical integration could be a difficult task.To overcome this problem, an approach to include complex frequency was suggested as early as the proposal of the DWN method itself.In what follows, a brief description of this method applied to fluid-solid interfaces is shown.
The incident pulse at the fluid, as shown in Figure 1, can be given as where  0  (x) is incident pulse at the fluid, x = { 1 ,  3 }, () is scale coefficient for the incident pulse,  (2)  0 (⋅) is Hankel's function of second kind and zero order,  is circular frequency,   is compressional wave velocity in the fluid,  = (x) is the distance from the receiver to the source (incident pulse),  is the wavenumber, and  = √ ( 2 /  2 ) −  2 with Im  < 0. If we express  in discrete values, then we have   = Δ and   = √ ( 2 /  2 ) −  2   with Im   < 0, Δ = wavenumber increment, and  is the imaginary unit.
If we consider that the entire pressure and displacement fields in the fluid are expressed as the sum of the free and diffracted fields, then one has where   is density of the fluid,  is the distance between the source and the interface, and   represents the diffraction coefficient for the fluid to be found.
The boundary conditions to be applied to the interface are for displacements, and for tractions.
Once the boundary conditions were applied, the unknown coefficients   ,   , and   are obtained and then the whole pressure field in the fluid is finally calculated by means of (2).The system of equation to be solved is given by

Formulation of the Problem by Means of the Indirect Boundary Element Method
For the IBEM, the incident pulse (source) in the fluid is also represented by ( 1) and applied as shown in Figure 2(a).Assume that the equation that governs the wave motion in the fluid is given by If we consider that stresses in the fluid can be linked to the pressure generated by the incident pulse, then one can express this as Then, the displacement field can be represented by the wellknown form: To express the diffracted wave field (for pressure and displacement, resp.) in the fluid due to the source impacting the elastic medium, we propose the use of the following integral representation: where where Ψ(⋅) is force density for the fluid,   (⋅) is Green's function for the fluid, and  1 = −0.5 and defines the region orientation (see explanation for  2 , given below).
The complete pressure and displacement field in the fluid, besides free and diffracted ones, can be written, respectively, as Since the source is only applied in the fluid, only diffracted waves appear in the solid and they can be established as follows.
Consider a domain  with a boundary .If the domain is occupied by an elastic solid, the displacement field under harmonic excitation can be expressed, neglecting body forces, by means of the single-layer boundary integral equation: where   (x) is th component of the displacement at point x,   (x; ) is Green's function, which represent the displacement produced in the direction  at x due to the application of a unitary force in direction  at point , and   () is the force density in the direction  at point .
The product   ()  is the force distribution at the surface  (the subscripts ,  are limited to be 1 or 3).The subscript in the differential shows the variable over which the integration is done.This integral equation can be obtained from Somigliana's representation [30].Furthermore, it was demonstrated that if   () is continuous along , then in that case, the displacement field is continuous across  [31].
The integral representation ( 16) permits the calculation of stresses and tractions by using the direct application of Hooke's law and Cauchy's equation, respectively, except for singularities on the boundary, that is, when x is equal to  on the surface .From a limiting process established on equilibrium considerations, around an internal vicinity of the boundary, one can write, for x on , where   (x) os the th component of tractions,  2 = 0.5 if x tends to the boundary  "from inside" the region,  2 = −0.5ifx tends to  "from outside" the region, or  2 = 0 if x is not at .   (x; ) is the traction Green's function, that is, the traction in the direction  at a point x, linked to the unit vector   (x), due to the application of a unitary force in the direction  at  on .The 2D Green's functions for infinite spaces can be seen in [32,33].

Boundary Conditions.
From the configuration depicted in Figure 2(b), it is convenient to partition the domain in three regions (, , and ), for which proper boundary conditions should be established.These conditions for fluidsolid interfaces can be written as follows.
For the fluid-solid interface, For the continuous solid interface, Expressing (18) as a function of the diffracted field ( 16) for the solid, and incident and diffracted fields (( 10) and ( 12), resp.) for the fluid, one obtains: The traction free condition ( 19) is expressed from its integral form (17), resulting in The Equation ( 20) can be expressed by means of ( 17), (1), and (11), and then one has Equations ( 21) express the continuity condition that must exist between the interface of regions  and  (boundaries  2  and ).These are defined as follows:

Discretization.
Here, the discretization of ( 22) to ( 26) is presented.Considering that force densities (x) and Ψ(x) are constant on each boundary element that forms the surfaces of regions , , and  (see Figure 2(b)) and Gaussian integration (or analytical integration, where Green's function is singular) is performed, ( 22) is rewritten as, Equation ( 23) leads to Equation ( 24) can be written as Equations ( 25) and ( 26) lead to where Equations ( 27), ( 30), (32), and (34) form a system of integral equations to be solved.Once the force densities ((x) and Ψ(x)) were obtained, the whole displacement and pressure fields in the fluid are found by means of ( 14) and ( 15), respectively.For the solid, the complete displacement and traction fields are calculated applying ( 16) and ( 17), respectively.In the following section, we verify the accuracy of the IBEM and DWN for several cases applied to an interface that joins two half-spaces, one fluid and the other solid.Moreover, we present results for layered models using IBEM.Additional details on the discretization scheme can be consulted in [35][36][37][38].Work related to fluid-solid interfaces using IBEM was also presented in [39], in which general aspects of wave propagation in fluid-solid media were pointed out.The novelty of the present paper is related to the application of the IBEM to model wave propagation in layered solid-fluid systems.In this sense, this paper deals specifically with fluidlayered solids and provides the following aspects.(e) Finally, the discussions focused on the diffracted waves that appear due to the presence of the layer.
These results are shown in the following sections.

Validation and Application of the IBEM and DWN Methods
To verify the accuracy of both formulations (IBEM and DWN), we considered various models with fluid-solid interfaces.Borejko [34] developed theoretical and experimental studies in order to show the emergence of interface waves.The material properties for his models are described in Table 1.
For comparison purposes, we chose an interface model joining two half-spaces, one fluid and the other an elastic solid, for cases 1 and 2 shown in Table 1.
As shown in Figure 1, one receiver located at the distances  = 0.05 m and  = 1.00 m was considered.In Figure 3, the spectrum of pressure for such receiver is shown for cases 1 and 2. Results from IBEM are plotted with solid and dashed lines to represent Limestone and Pitch, respectively.Calculations from DWN are depicted with circles for case 1 and asterisks for case 2. Good agreement can be appreciated between IBEM and DWN results.It is clear from this figure that the responses for both materials vary significantly and can be associated with the relative value of the shear wave velocity of the solid in comparison with the fluid wave velocity [34].
For shear wave velocities higher than the compressional wave velocity for the fluid, the spectrum of pressures shows more simple patterns in comparison with the opposite case, where some resonant peaks are observed.
In Figure 4, spectra of pressures for four cases (3 to 6) solved by IBEM are displayed.For the cases 5 and 6, the thickness of the layer is ℎ = 0.05 m and ℎ = 0.10 m, for both cases.For these layered models, h refers to the layer thickness (e.g., sandstone is the layer for case 5, while granite is the layer for case 6).Cases 3 and 4 correspond to simple models (let us refer to them as simple models), in other words, an interface that joins two half-spaces, one fluid and the other solid.These cases are plotted using dotted lines (case 3) and circles (case 4), both shown in Figures 4(a) and 4(b).In Figure 4(a), the response for case 5 is also included (for ℎ = 0.05 m and ℎ = 0.10 m).
An interesting fact that emerges from the results observed in Figure 4(a) is that the behavior shown by the layered model (sandstone-granite) is delimited by the simple models (granite and sandstone).Besides, at low frequencies the layered model behavior shows a clear tendency to case 3.In other words, the material that constitutes the half-space (see region  in Figure 2(b)) controls this phenomenon (i.e., granite for case 5 and sandstone for case 6).At high frequencies the behavior is controlled by the shallow material (see region  in Figure 2(b)) due to shorter wavelengths, and then the material that constitutes the half-space has no influence in the response.
Analogously, the behavior shown by case 6 is delimited by the simple models (cases 3 and 4).At low frequencies, its response is very close to that obtained by case 4. In contrast, at high frequencies, such response is near to case 3 (granite).For both figures (Figures 4(a) and 4(b)), there is a transition zone, which depends on the layer thickness (ℎ).
In Figure 5, pressure fields in time domain are shown for cases 3, 4, and 5 (see Table 1).To this end, a fast Fourier transform (FFT) algorithm was applied using a Ricker wavelet as source time function (i.e., the second derivative of a Gaussian) with a characteristic period of 2.604 × 10 −5 sec; this source is used in all analyses presented here.All models were analyzed for frequency increments of 150 Hz, reaching a final frequency of 19200 Hz.By means of the Fourier transformation, it is possible to observe different kinds of waves that emerge.
For all cases, 25 receivers were located in the fluid.The first one was placed in a distance of  = 1.0 m from the  source, and the rest of the receivers were placed using a distance increment of 0.05 m.The distance  = 0.05 m (see Figure 2) for all cases.For case 5 (fluid-sandstone-granite interface), two layer thicknesses were considered, which are ℎ = 0.20 m (Figure 5(c)) and ℎ = 0.40 m (Figure 5(d)).
In Figure 5(a), the wave propagation phenomenon for case 4 (fluid-sandstone interface) is shown.Here, it is possible to observe the influence of the  (compressional wave velocity of sandstone, also known as -wave velocity) represented as    ; the direct wave that travels in the fluid and that is perceived by the receivers is shown with    and the Scholte's interface wave is illustrated using    .The super index  represents "Sandstone, " while  is for "Fluid." Borejko [34] also found this kind of waves by means of theoretical and experimental studies.For this case, the velocities measured were    ≈ 2600 ms −1 ,    = 1500 ms −1 , and    = 937.5 ms −1 .In Figure 5(b), the wave fronts that emerge from the fluid-granite interface interactions are depicted.Here, it is possible to identify wave velocities associated with pseudo Rayleigh, direct-and Scholte's waves, which propagate at    = 3076.9ms −1 ,    = 1500 ms −1 , and    = 1500 ms −1 , respectively.The super index  refers to "Granite." Our results for these last two cases agree with those obtained by Borejko [34].It is important to mention that Scholte's wave travels at a velocity close to the direct wave in the fluid and, then, only one wave front is seen.This was also reported by Borejko [34].
For the interface model (case 5), two layer thicknesses were studied, as mentioned above.In this case, the influence of granite is evident.In Figure 5(c) (for ℎ = 0.20 m), the direct wave represented by    and the four fronts associated with Scholte's wave velocity (   ) are clearly observed.Moreover, two wave fronts that travel at a velocity of    = 1090 ms −1 are identified.This velocity coincides with  (shear wave velocity for sandstone, also known as  wave velocity).The subindex  stands for  wave velocity.In Figure 5(d), these same wave fronts appear, but less interactions are evidently appreciated due to the considered layer thickness (ℎ = 0.40 m).

Conclusions
The indirect boundary element method to study the wave propagation phenomenon in 2D fluid-layered solid interfaces was used.This indirect formulation can give to the analyst a deep physical insight on the generated diffracted waves because it is closer to the physical reality and can be regarded as a realization of Huygens' principle.In any event, it is fully equivalent to the classical Somigliana representation theorem.In order to verify the accuracy, we tested our method by comparing results with the analytical solution known as discrete wave number.A near interface pulse generates scattered waves that can be registered by receivers located in the fluid.Results were presented in frequency and time domains, where some aspects related to different wave types that emerge from this kind of problems were pointed out.The results between IBEM and DWN describe the same physics and are, for engineering purposes, quite approximated.

Figure 1 :
Figure 1: Fluid-solid interface solved by means of the DWN technique.

Figure 2 :
Figure 2: (a) Fluid-layered solid interface excited by a source in the fluid; (b) mesh used to solve the problem with boundary element.
(a) It contains the complete formulation applied to fluidlayered solid interfaces.Moreover, the paper describes the precise boundary conditions and discretization scheme for fluid-layered solid media.(b) It contains pressure spectra that clearly show the influence of a fluid-layered medium, according to its layer thickness.(c) It presents synthetic seismograms of pressures for several interfaces like fluid-sandstone-granite interfaces and others.(d) It pointed out the kind of interface waves that emerge in each case.

Figure 3 :
Figure3: Results from IBEM and DWN are shown.Spectrum of pressures by IBEM, for cases 1 and 2, is plotted with solid (Limestone) and dashed (pitch) lines, while those obtained by DWN are displayed using asterisk (limestone) and circle (pitch) lines.Here, good agreement between both methods is also observed.

Figure 4 :
Figure 4: Spectrum of pressures by IBEM for cases 3 to 6 is shown.Cases 3 and 4 are plotted with dotted (granite) and circle (sandstone) lines, respectively.Moreover, in (a) case 5, sandstone-granite layered model for ℎ = 0.05 m (solid line) and for ℎ = 0.10 m (dashed line) is plotted, and (b) case 6, granite-sandstone layered model for ℎ = 0.05 m (solid line) and for ℎ = 0.10 m (dashed line) is displayed.