In the present communication, scattering of elastic waves in fluid-layered solid interfaces is studied. The indirect boundary element method is used to deal with this wave propagation phenomenon in 2D fluid-layered solid models. The source is represented by Hankel’s function of second kind and this is always applied in the fluid. Our method is an approximate boundary integral technique which is based upon an integral representation for scattered elastic waves using single-layer boundary sources. This approach is typically called indirect because the sources’ strengths are calculated as an intermediate step. In addition, this formulation is regarded as a realization of Huygens’ principle. The results are presented in frequency and time domains. Various aspects related to the different wave types that emerge from this kind of problems are emphasized. A near interface pulse generates changes in the pressure field and can be registered by receivers located in the fluid. In order to show the accuracy of our method, we validated the results with those obtained by the discrete wave number applied to a fluid-solid interface joining two half-spaces, one fluid and the other an elastic solid.

1. 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 [812]. 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 [1418]. Attenuation and dispersion in interface waves due to the presence of fractures were studied in [11, 1922]. 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 fluid-solid 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.

2. 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 whereis incident pulse at the fluid,,is scale coefficient for the incident pulse,is Hankel’s function of second kind and zero order,is circular frequency,is compressional wave velocity in the fluid,is the distance from the receiver to the source (incident pulse),is the wavenumber, andwith Im. If we expressin discrete values, then we haveandwith,wavenumber increment, andis 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 whereis density of the fluid,is the distance between the source and the interface, andrepresents the diffraction coefficient for the fluid to be found.

For the solid, we assume that potentials of displacement have the formand, wherewithandwith.andare the compressional and shear wave velocities, respectively. Similarly,andrepresent the diffraction coefficients, for the solid, to be found.

The displacement field for the solid is expressed asand. The stress field can be obtained from the well-known equation: whereis stress tensor,andare Lamé’s constants,is strain tensor, andis Kronecker’s delta.

The boundary conditions to be applied to the interface are for displacements, and for tractions.

Once the boundary conditions were applied, the unknown coefficients,, andare 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

3. 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 well-known 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 whereis force density for the fluid,is Green’s function for the fluid, andand defines the region orientation (see explanation for, 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 domainwith 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: whereisth component of the displacement at point,is Green’s function, which represent the displacement produced in the directionatdue to the application of a unitary force in directionat point, andis the force density in the directionat point. The productis the force distribution at the surface(the subscriptsare 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 ifis 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, whenis equal toon the surface. From a limiting process established on equilibrium considerations, around an internal vicinity of the boundary, one can write, foron, whereos theth component of tractions,iftends to the boundary“from inside” the region,iftends to“from outside” the region, orifis not at.is the traction Green’s function, that is, the traction in the directionat a point, linked to the unit vector, due to the application of a unitary force in the directionaton. The 2D Green’s functions for infinite spaces can be seen in [32, 33].

3.1. 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 fluid-solid 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 regionsand(boundariesand). These are defined as follows:

3.2. Discretization

Here, the discretization of (22) to (26) is presented. Considering that force densitiesandare 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 (and) 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 [3538]. 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 fluid-layered solids and provides the following aspects.(a)It contains the complete formulation applied to fluid-layered 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.(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.

4. 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 casesandshown in Table 1.

As shown in Figure 1, one receiver located at the distances m and m was considered. In Figure 3, the spectrum of pressure for such receiver is shown for casesand. 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 and asterisks for case. 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 (to) solved by IBEM are displayed. For the casesand, the thickness of the layer is m and m, for both cases. For these layered models, h refers to the layer thickness (e.g., sandstone is the layer for case, while granite is the layer for case). Casesandcorrespond 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) and circles (case), both shown in Figures 4(a) and 4(b). In Figure 4(a), the response for caseis also included (for m and 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. In other words, the material that constitutes the half-space (see regionin Figure 2(b)) controls this phenomenon (i.e., granite for caseand sandstone for case). At high frequencies the behavior is controlled by the shallow material (see regionin 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 caseis delimited by the simple models (casesand 4). At low frequencies, its response is very close to that obtained by case. In contrast, at high frequencies, such response is near to case(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,, and(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 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 m from the source, and the rest of the receivers were placed using a distance increment of 0.05 m. The distance m (see Figure 2) for all cases. For case(fluid-sandstone-granite interface), two layer thicknesses were considered, which are m (Figure 5(c)) and m (Figure 5(d)).

In Figure 5(a), the wave propagation phenomenon for case(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 withand the Scholte’s interface wave is illustrated using. The super indexrepresents “Sandstone,” whileis for “Fluid.” Borejko [34] also found this kind of waves by means of theoretical and experimental studies. For this case, the velocities measured werems−1,ms−1, and 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 ms−1, ms−1, and ms−1, respectively. The super indexrefers 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), two layer thicknesses were studied, as mentioned above. In this case, the influence of granite is evident. In Figure 5(c) (for m), the direct wave represented byand the four fronts associated with Scholte’s wave velocity () are clearly observed. Moreover, two wave fronts that travel at a velocity of ms−1 are identified. This velocity coincides with(shear wave velocity for sandstone, also known aswave velocity). The subindexstands forwave velocity. In Figure 5(d), these same wave fronts appear, but less interactions are evidently appreciated due to the considered layer thickness ( m).

5. 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.


Thanks are given to G. Sánchez, E. Plata, E. Chávez, and L. A. Cruz of Unidad de Servicios de Información (USI) of Instituto de Ingeniería, UNAM; their help was crucial to locate useful references. Partial supports from DGAPA-UNAM, Project IN121709, SENER-CONACYT Project 128376, and from the Instituto Mexicano del Petróleo are greatly appreciated. This work has been partial supported by the National Polytechnic Institute of México under project SIP-20130622.