On wavemodes at the interface of a fluid and a fluid-saturated poroelastic solid

On wavemodes at the interface of a fluid and a fluid-saturated

poroelastic solid

K. N. van Dalen,a兲 G. G. Drijkoningen, and D. M. J. Smeulders

Department of Geotechnology, Faculty of Civil Engineering and Geosciences, Delft University of Technology, Stevinweg 1, 2628 CN Delft, The Netherlands

共Received 21 April 2009; revised 12 January 2010; accepted 12 January 2010兲

Pseudo interface waves can exist at the interface of a fluid and a fluid-saturated poroelastic solid. These waves are typically related to the pseudo-Rayleigh pole and the pseudo-Stoneley pole in the complex slowness plane. It is found that each of these two poles can contribute共as a residue兲 to a full transient wave motion when the corresponding Fourier integral is computed on the principal Riemann sheet. This contradicts the generally accepted explanation that a pseudo interface wave originates from a pole on a nonprincipal Riemann sheet. It is also shown that part of the physical properties of a pseudo interface wave can be captured by loop integrals along the branch cuts in the complex slowness plane. Moreover, it is observed that the pseudo-Stoneley pole is not always present on the principal Riemann sheet depending also on frequency rather than on the contrast in material parameters only. Finally, it is shown that two additional zeroes of the poroelastic Stoneley dispersion equation, which are comparable with the P¯ -poles known in nonporous elastic solids, do have physical significance due to their residue contributions to a full point-force response. © 2010 Acoustical Society of America. 关DOI: 10.1121/1.3308473兴

PACS number共s兲: 43.35.Pt, 43.20.Jr, 43.20.Gp, 43.20.Bi 关JM兴 Pages: 2240–2251

I. INTRODUCTION

Interface waves such as Rayleigh and Stoneley waves are often used to investigate materials. One can think of ap-plications in ultrasonic testing of structures, borehole logging in geotechnical and reservoir engineering, and surface seis-mics in geophysics, see, e.g., Refs.1–4. In the case of porous materials, interface waves carry information on elastic prop-erties but also on propprop-erties such as porosity, permeability, and fluid mobility.1Rosenbaum5found that, compared to all other surface and body wavemodes, the Stoneley-type wave that travels along the open-pore interface of a fluid and a porous medium, carries the best measure of permeability.

Several theoretical studies were performed on interface waves that propagate along the boundary of a porous me-dium. These studies were carried out in the context of Biot’s theory for wave propagation in fluid-saturated poroelastic solids. Deresiewicz6showed the existence of a Rayleigh-type wave that propagates along the free surface of a poroelastic half-space and analyzed the frequency-dependent phase ve-locity and attenuation.

For a fluid/poroelastic-medium configuration, Rosen-baum5 predicted the existence of the pseudo-Rayleigh 共pR兲 and the pseudo-Stoneley共pSt兲 wave. The latter was explic-itly named as such by Feng and Johnson,7,8 since a pseudo interface wave has part of its energy leaking into slower bulk modes as it propagates along the interface. Feng and Johnson7 also showed the existence of another interface wavemode, the nonleaky true interface wave. It was found that the existence of this wave depends on whether or not the pores are open for pore fluid to flow across the interface.

Feng and Johnson8 derived Green’s functions 共impulse re-sponses兲 for high-frequency Biot theory that confirmed the existence of the three different waves.

Experimental evidence was found for all three types of interface wavemodes, see, e.g., Refs.9–11.

Feng and Johnson7 argued that other zeroes of the po-roelastic Stoneley dispersion equation have no physical sig-nificance as pseudo interface modes. The corresponding propagation velocities would be larger than that of shear waves, which is not realistic.

In order to obtain the characteristics of the interface wavemodes, Feng and Johnson7 used the zeroes of the non-viscid poroelastic Stoneley dispersion equation in the com-plex plane to obtain the propagation velocities and attenua-tions. Gubaidullin et al.12 went a step further and analyzed the frequency dependence of the interface wavemodes by incorporating the viscous loss mechanism of Johnson et al.13 They also used the zeroes of the dispersion equation to de-rive the characteristics of the interface waves. The same ap-proach was adopted by Edelman and Wilmanski,14Albers,15 and Markov.16In most of the papers, specific restrictions for the involved square roots 共i.e., their Riemann sheets兲 are given.

The generally accepted explanation for a pseudo inter-face wave is that it originates from a zero that forms a pole singularity on another Riemann sheet than the so-called “principal” sheet. It affects the behavior of the integrand on the principal Riemann sheet by causing a local maximum in the integrand.17In case the pole lies close to the real axis, it might have a contribution to the Green’s function.

In a series of publications, Allard et al.11,18,19studied the propagation of interface waves along the boundaries of po-roelastic and nonporous elastic media. In the case of an air/ air-saturated poroelastic-solid configuration, they found that

a兲_{Author to whom correspondence should be addressed. Electronic mail:}

taking the residue of the pseudo-Rayleigh pole is sufficient to describe the entire pseudo-Rayleigh waveform.19 For the water/water-saturated poroelastic-solid configuration, they found that the pseudo-Stoneley pole residue describes the entire waveform of the pseudo-Stoneley wave.11 However, for the water/elastic-solid configuration, they found that the pseudo-Rayleigh waveform is strongly affected by the loop integrals along the branch cuts.

In summary, taking just the location of the zeroes of the dispersion equation rather than computing the full transient response is a very fast way to predict the kinematic proper-ties of pseudo interface waves, but the question arises if these predictions are always complete.

Therefore, in this paper we analyze the three-dimensional transient wave propagation due to a point force applied at the interface of a fluid and a fluid-saturated po-roelastic solid. The aims are as follows.

共1兲 To investigate if a zero of the poroelastic Stoneley dis-persion equation indeed yields the pertinent physical properties of the corresponding pseudo interface wave-mode. This is done by quantitative comparison between the residues of specific poles and the full transient re-sponse.

共2兲 To verify if a pseudo interface wave indeed necessarily originates from a pole on a nonprincipal Riemann sheet. 共3兲 To verify the physical significance of additional zeroes of the poroelastic Stoneley dispersion equation that are not related to pseudo interface waves.7

The paper is organized as follows. In Sec. II, we present the model to analyze the fluid/poroelastic-medium configu-ration. Subsequently, in Sec. III, the derivation of Green’s function is summarized. The implementation of the numeri-cal integration is discussed in Sec. IV. We discuss the results in Sec. V. The conclusions are given in Sec. VI.

II. MODEL

To study the transient wave propagation in a fluid/ poroelastic-medium configuration, we consider a configura-tion that consists of a fluid half-space on top of a fluid-saturated poroelastic half-space. A vertical point force F共t兲 is applied at the interface关see Fig.1共a兲; Fig.1共b兲is referred to later兴. Both half-spaces are considered to be homogeneous and isotropic. The configuration is similar to the one applied by Gubaidullin et al.12but extended to three dimensions.

The behavior in the lower half-space 共x3⬎0兲 is

gov-erned by the well-known Biot equations of motion for a fluid-saturated poroelastic solid that were extensively dis-cussed in this journal, see, e.g., Refs. 20and21. Following Biot’s theory, we assume that for long wavelength distur-bances with respect to the characteristic pore scale, average local displacements can be defined for the solid 共frame兲 u共x,t兲=共u₁, u₂, u₃兲T _{and the fluid U共x,t兲=共U}

1, U2, U3兲T.

Considering a cube of unit size of bulk material共porosity␾兲, the forces per unit bulk area applied to that part of the cube faces occupied by the solid are denoted by ␶ij. They are constituted by both fluid pressure pf and intergranular stresses␴ijaccording to

␶ij= −␴ij−共1 −␾兲pf␦ij, 共1兲

where ␦ij is the Kronecker delta. The total normal tension force per unit bulk area applied to the fluid faces of the unit cube, denoted by ␶, is constituted by pf only,

␶= −␾pf. 共2兲

Here, ␴ij and pf are defined positive in compression and, consequently,␶ijand␶are positive in tension, see also Ref. 12. In the case of isotropic materials, the stress-strain rela-tions for the solid and the fluid can be written as

␶ij= G共⳵iuj+⳵jui兲 + A⳵kuk␦ij+ Q⳵kUk␦ij, 共3兲

␶= Q⳵kuk+ R⳵kUk, 共4兲

where Einstein’s summation convention for repeated indices is applied, and⳵j=⳵/⳵xj. A, Q, and R are generalized elastic constants that can be related via Gedanken experiments to porosity, grain bulk modulus Ks, fluid bulk modulus Kf, bulk modulus of porous drained frame Kb, and shear modulus G of both drained frame and total composite.22,23The physical background of Eqs.共3兲and共4兲is discussed in more detail in Ref. 20.

The equations of motion are found from combination of momentum conservation and the stress-strain relations, Eqs. 共3兲 and共4兲, and can be written as20,21

␳11⳵t 2 u +␳12⳵t 2 U + bⴱ⳵t共u − U兲 = P ⵜ ⵜ · u − G ⵜ ⫻ ⵜ ⫻ u + Qⵜ ⵜ · U, 共5兲 ␳12⳵t2u +␳22⳵t2U − bⴱ⳵t共u − U兲 = Q ⵜ ⵜ · u + R ⵜ ⵜ · U, 共6兲 where the asterisk denotes convolution, P = A + 2G, and the effective densities are defined as

␳11=共1 −␾兲␳s−␳12,

␳22=␾␳f−␳12,

␳12= −共␣⬁− 1兲␾␳f, 共7兲

where the tortuosity ␣_⬁ⱖ1, and hence ␳₁₂ⱕ0. Solid and fluid densities are denoted as ␳s and ␳f, respectively. The

(b) Fluid-saturated poroelastic solid Fluid 3 x 1 x ( ) F t Fluid-saturated poroelastic solid Fluid 3 x 1 x ( ) F t (a) x3 1 x pSt P2pSt F P2 P1 S pR FP1 FpR SP1 3 x 1 x pSt P2pSt F P2 P1 S pR FP1 FpR SP1

FIG. 1. 共a兲 Point force F共t兲 applied at the interface of a fluid-saturated poroelastic half-space and a fluid half-space. Both half-spaces are homoge-neous and isotropic.共b兲 Schematic snapshot of the full response with sepa-rate arrivals: fast compressional共P1兲 wave, slow compressional 共P2兲 wave,

shear共S兲 wave, fluid 共F兲 wave, Rayleigh 共pR兲 wave, and pseudo-Stoneley共pSt兲 wave. The double-mode symbols 共e.g., SP1兲 indicate lateral

waves共¯兲. The first symbol denotes the wavemode of the specific arrival; the second denotes the one from which it is radiated. Here, the F-wave velocity is assumed higher than the P2-wave velocity. For clarity, we

linear time-convolution operator b was formulated in the fre-quency domain as the viscous correction factor by Johnson et

al.,13 according to

bˆ共␻兲 = b₀

共

1 +1₂iM␻/␻c

兲

1/2_{, 共8兲}

where the viscous damping factor b0=␾2␩/k0. Here, the

dy-namic fluid viscosity is denoted by␩ and k0 represents the

zero-frequency Darcy permeability. The shape factor M is usually taken equal to 1.24 The rollover frequency, which represents the transition from low-frequency viscosity-dominated to high-frequency intertia-viscosity-dominated behavior, is defined as␻c=␩␾/共␣⬁␳fk0兲.

The behavior of the upper 共fluid兲 half-space 共x3⬍0兲 is

governed by the acoustic wave equation

␳F⳵t2pF= KFⵜ2pF, 共9兲

where KFand␳Fdenote the bulk modulus and density of the fluid, respectively, and pF denotes the fluid pressure.

We assume that the behavior at the interface is governed by conventional open-pore conditions, i.e., by continuity of volume flux and fluid pressure, and vanishing intergranular vertical and shear stresses. The force is applied to the solid. The open-pore boundary is a realistic choice to model the fluid/poroelastic-medium interface,1 and a limiting case of the situation where a finite surface flow impedance is con-sidered, see e.g., Refs.12,25, and26. It implies that the true interface wave is absent in the response.7,8 Hence, in the limit of x3→0, the following conditions should be satisfied

共1 −␾兲u3+␾U3− UF,3= 0, 共10兲

pf− pF= 0, 共11兲

␴13= 0, 共12兲

␴23= 0, 共13兲

␴33= F共t兲␦共x1兲␦共x2兲, 共14兲

where␦共¯兲 denotes the Dirac delta function, and UF,3 de-notes the vertical particle displacement in the upper half-space. The fact that the intergranular stress ␴33 is zero does

not imply that the total solid stress␶33vanishes, see Eq.共1兲.

The medium is considered to be at rest at tⱕ0. At infi-nite distance from the source, the motions are bounded. III. GREEN’S FUNCTIONS

In this section, we summarize the derivation of the Green’s functions 共impulse responses兲 as described by the solution to the set of governing equations, Eqs.共5兲,共6兲, and 共9兲–共14兲. The main part of the derivation is given in Appen-dices A and B and we refer to them where necessary.

In order to analyze the response in the plane-wave do-main, the Fourier transform is applied over time according to

uˆ共x,␻兲 =

冕

−⬁ ⬁

u共x,t兲exp共− i␻t兲dt, 共15兲

where␻denotes radial frequency. It is assumed that u共x,t兲 is real valued and hence, it is sufficient to consider␻ⱖ0.

Fol-lowing Aki and Richards,4 the Fourier transform over hori-zontal spatial coordinates can be defined as

u ˜共p,x₃,␻兲 =

冕

−⬁ ⬁

冕

−⬁ ⬁ uˆ共x,␻兲exp共i␻p · r兲dx₁dx₂, 共16兲 where p =共p1, p2兲T is the horizontal slowness vector and r

=共x₁, x₂兲T _{is the horizontal space vector. The transforms are} applied similarly to the other field quantities. The hat refers to the共x,␻兲-domain and the tilde to the 共p,x₃,␻兲-domain.

The response in the共p,x₃,␻兲-domain is described by the physical quantities u˜i and p˜f in the lower half-space, col-lected in the vector w˜ =共u˜₁, u˜₂, u˜₃, −␾˜pf兲T, and by p˜F in the upper half-space, see Eqs.共A1兲,共A4兲, and共A6兲. The expres-sions for the response are obtained using Helmholtz decom-position of the equations of motion and substitution of the general solutions into the boundary conditions. This gives a set of equations that is solved analytically共see Appendix A兲. The response can be written in terms of Green’s func-tions according to w˜ = g˜+_{Fˆ =} n˜ + ⌬St Fˆ , 共17兲 p ˜F= g˜−Fˆ = n ˜− ⌬St Fˆ , 共18兲

where g˜+_{and g˜}−_{are the Green’s functions in the lower and}

upper media, respectively, n˜+ _{and n˜}− _{are the corresponding}

numerators, and Fˆ is the Fourier transform of the force sig-nature. From Eqs. 共A1兲,共A4兲, and 共A6兲, it follows that g˜+

consists of a superposition of all possible body modes: the fast共P1兲 and slow 共P2兲 compressional waves, and the

verti-cally polarized shear 共SV兲 wave. The horizontally polarized shear 共SH兲 mode is not excited by the vertical force. The Green’s function g˜− _{only contains the fluid 共F兲}

compres-sional mode. Both Green’s functions have the “poroelastic Stoneley-wave denominator”⌬St=⌬St共p,␻兲 that is associated with interface waves along the fluid/poroelastic-medium in-terface, which is very similar to the “Scholte-wave denomi-nator” for a fluid/elastic-solid interface.27 Here, p =共p₁2 + p₂2兲1/2denotes the magnitude of the horizontal slowness.

The body-wave slownesses sj, j =兵P1, P2, F , S其, are

de-fined in Appendix A 共TableIII兲. The corresponding vertical slownesses are defined as qj=共sj

2

− p2兲1/2, where Im共qj兲ⱕ0 due to Sommerfeld’s radiation condition.

To find the Green’s functions in the 共x,␻兲-domain, the inverse Fourier transform is applied according to

gˆ+₌ ␻ 2 共2␲兲2

冕

−⬁ ⬁

冕

−⬁ ⬁ _n˜+_共p,x 3,␻兲 ⌬St exp共− i␻p · r兲dp₁dp₂, 共19兲 where ␻ⱖ0. We only show the derivation of gˆ+_{, but the}

expressions for gˆ−_{are obtained by simply replacing n˜}+_{by n˜}−_.

When cylindrical coordinates are introduced, Eq.共19兲can be written as

gˆ+₌ ␻ 2 4␲

冕

_−⬁ ⬁ _˜n+_共p,⳵ ␣,x3,␻兲 ⌬St H₀共2兲共␻pr兲pdp, 共20兲

where r =共x₁2+ x₂2兲1/2 and in which the horizontal derivatives ⳵␣, ␣=兵1,2其, are applied to the Hankel function H0

共2兲_共¯兲,

see Eqs.共B3兲 and共B4兲 共Appendix B兲.

Now we change the real-axis integral into a contour in-tegral in the complex p-plane. The idea is that by integration in the complex plane, contributions from loop integrals and from pole residues can be distinguished. We choose branch cuts along the hyperbolic lines3 Im共qj兲=0. In this way Im共qj兲ⱕ0 ∀ p, which ensures the decay of the exponential terms exp共⫿i␻qjx3兲 for large p 关see Eqs. 共A1兲 and 共A4兲兴.

The branch cuts depart from the branch points associated with the body-wave slownesses sj, as shown in Fig. 2. The qF-branch cut reduces to the imaginary axis and part of the real axis since the slowness of the fluid wave共sF兲 is real.

The current branch cuts are referred to as the “funda-mental” branch cuts.28 The corresponding Riemann sheet is referred to as the principal Riemann sheet17or the “physical” Riemann sheet.4

In Fig. 2 the closed contour is also displayed. It is formed by the entire real axis, the loops along the branch cuts and around the branch points, and an arc of infinite radius in the lower half-plane. For Re共p兲ⱕ0, the horizontal part of the contour lies just below the axis due to the pres-ence of a branch cut of the Hankel function at the negative real axis.29

Applying Cauchy’s residue theorem,30 we obtain gˆ+=

冕

−⬁ ⬁ f ˜+_{dp = − 2␲i}

兺

m Resp=s_m˜f+−

兺

j

冕

C_j f ˜+_dp, f ˜+₌ ␻ 2 4␲ n ˜+_共p,⳵ ␣,x3,␻兲 ⌬St H₀共2兲共␻pr兲p, 共21兲

where every sm denotes a first-order pole of the integrand inside the integration contour and every Cj denotes a loop along the specific branch cut. In Eq.共21兲, the contribution of the arc vanishes because of Jordan’s lemma.31The contribu-tions around the branch points are also zero.

The poles smresult from zeroes of the poroelastic Stone-ley dispersion equation 共⌬St= 0兲 on the principal Riemann sheet. The number of poles N present inside contour C is determined by applying the principle of the argument to the Stoneley equation30 N = 1 2␲i

.

C ⳵p⌬St ⌬St dp. 共22兲

The residue of the integrand at a first-order pole is given as Resp=s_m˜f+=

冋

␻2 4␲ n ˜+_共p,⳵ ␣,x3,␻兲 ⳵p⌬St H₀共2兲共␻pr兲p

册

p=sm . 共23兲

IV. NUMERICAL IMPLEMENTATION

To perform the integration along the hyperbolic branch cuts, we choose pi= Im共p兲 as the variable of integration ac-cording to

冕

C_j f ˜+⳵p ⳵pi dpi, 共24兲 where p =Re共sj兲Im共sj兲 pi + ipi, ⳵p ⳵pi = −Re共sj兲Im共sj兲 pi 2 + i. 共25兲

For the qF-branch cut, the integration path is the imaginary axis and part of the real axis, which follows from Eq.共25兲for vanishing imaginary part of the slowness Im共sF兲↑0 共Fig.2兲. Along the cut of qj, at the left side Re共qj兲⬎0 and at the right side Re共qj兲⬍0. At the specific cut Im共qj兲=0 and everywhere else Im共qj兲⬍0.

The numerical integration is performed using an adap-tive eight-point Legendre–Gauss algorithm,29,32 which can handle integrable singularities such as branch points.

For the numerical implementation of the principle of the argument, we apply Eq.共22兲separately for the areas between the various parts of the integration contour共branch cuts, real axis, and arc, see Fig.2兲 to find out where the poles can be expected. Subsequently, the pole locations are found numeri-cally by minimizing the left-hand side of equation 兩⌬St兩=0. Since it contains local minima and branch-cut discontinui-ties, it is important to choose a proper starting value. This requires some manual iteration. The accuracy, as expressed by兩⌬St共p=sm兲兩/兩⌬St共p=0兲兩, is typically O共10−10兲. Here sm de-notes the numerical value of the pole location.

V. NUMERICAL RESULTS AND DISCUSSION

In this section, we investigate the transient responses for four different fluid/poroelastic-medium configurations 共see Table I兲. In the first three configurations, water-saturated Bentheimer sandstone 共see Table II兲 is used as porous me-dium. The upper half-space is subsequently filled with water, air, or a light fluid. In the fourth configuration, which is the one of Feng and Johnson,7,8the porous medium is formed by water-saturated fused glass beads, while the upper half-space is filled with water.

For every configuration, we will show the vertical com-ponent of particle velocityv₃and the fluid pressure pf for an

1 P s S s 2 P s F s Re p Im p pR s * 2 P C * Pa s * Pb s 1 P s S s 2 P s F s Re p Im p pR s * 2 P C * Pa s * Pb s

FIG. 2. Complex p-plane with共– –兲 branch cuts, 共•兲 branch points sj, j =兵P1, P2, F , S其, and 共ⴱ兲 poles spR 共pseudo-Rayleigh兲, s¯ aP and sP¯ b 共addi-tional兲, for the calculation of the Green’s functions for Bentheimer/air con-figuration 2共see TableI兲. The branch points are formed by the body-wave slownesses specified in Appendix A共TableIII兲. The hyperbolic branch cuts are described by Im共qj兲=0. Poles are zeroes of the poroelastic Stoneley denominator, see Eq.共21兲. Only part of the closed integration contour共–兲 is displayed: real axis, arc in lower half-plane, and loop C_P2 along the qP2-branch cut. The direction of integration is indicated.

observation point at the interface x3= 0 at offset r = x1

= 0.1 m. Fluid pressure is related to dilatation only关see Eq. 共4兲兴 and hence, v3and pfcontain different information. Also, the comparison between the full response and a pole residue can be different inv3and pf, as will be shown.

The point force has Ricker signature,33

F共t兲 = Fmax

共

1 2␻0

2

t¯2− 1

兲

exp

共

−₄1␻₀2t¯2

兲

, 共26兲 where t¯= t − ts,␻0= 2␲f0, and center frequency f0= 500 kHz

共see Fig. 3兲. The magnitude F_max= 1 N and time shift ts = 5 ␮s. We perform the integration for the frequency range 0⬍ f ⱕ2 MHz. The full response is obtained by multiplica-tion of the spectra of the Green’s funcmultiplica-tions and the source 关see Eqs. 共17兲 and共18兲兴, and using a standard fast Fourier transform algorithm.

A. Residue contribution vs full response

First, we address the relation between a pole and a pseudo interface wave, as raised in point共1兲 in the Introduc-tion 共Sec. I兲. For configurations 1–3, the full transient re-sponses and separate pole residues 关see Eq. 共21兲兴 are dis-played in Figs.3–5. We identified the different arrivals in the full responses using the propagation velocities as obtained from the modal slownesses. Head waves are identified geo-metrically using the pertaining modal velocities and are in-dicated with double-mode symbols 关e.g., SP1: the shear共S兲

wave radiated by the fast共P₁兲 compressional wave兴. For the sake of clarity, a schematic snapshot of the full response with the different arrivals is shown in Fig.1共b兲.

We first note that the P₁-wave is present quite strongly inv₃although this component is perpendicular to the direc-tion of propagadirec-tion of this longitudinal wave共Figs.4and5兲. This is due to the contraction in vertical direction that can easily take place at the air/sandstone or light-fluid/sandstone interface. Remarkably, there is an arrival present in pf at the S-wave arrival time 共Figs. 4 and 5兲. This is not an

S-wavefront but radiated slow compressional 共P2兲 and fluid

共F兲 head waves, see Fig.1共b兲.

Now, we focus on the comparison of interface waves in the full responses and corresponding pole residues. The pole共s兲 that are present on the principal Riemann sheet con-tributing a residue are given in TableI, for each configuration separately. We found the pseudo-Stoneley共pSt兲, the pseudo-Rayleigh (pR), and two additional共P¯a, P¯b兲 poles. The latter ones are discussed in Sec. V C.

TABLE I. Various configurations for which the transient response is calcu-lated. The type of sandstone is Bentheimer. For fused glass beads, the bulk modulus of the drained matrix is chosen as Kb= 10 GPa and the permeabil-ity is chosen as k0= 10−11 m2. The upper half-space is filled with either

water共KF= Kf, ␳F=␳f兲, or air 共KF= 1.42⫻102 kPa, ␳F= 1.25 kg m−3兲, or a light fluid 共KF= Kf/10, ␳F=␳f/8兲. For every configuration, the poles present on the principal Riemann sheet are indicated: pseudo-Stoneley共pSt兲, pseudo-Rayleigh共pR兲, and two additional 共P¯a, P¯b兲 poles.

Porous solid Saturating fluid Upper half-space Poles

1 Sandstone Water Water pSt

2 Sandstone Water Air P¯a, P¯b,apR

3 Sandstone Water Light fluid pR, pSta

4 Fused glass beads Water Water pSta

a_{Its residue is not shown.}

TABLE II. Material parameters as used for water-saturated Bentheimer sandstone共Ref.34兲. The bulk modulus of the matrix Kbis found according to Kb= Kp−

4 3G.

Solid共frame兲 density␳s 2630 kg m−3

Fluid density␳f 1000 kg m−3

Tortuosity␣_⬁ 2.4

Porosity␾ 0.23

Permeability k0 3.7 ␮m2

Dynamic fluid viscosity␩ 0.001 Pa s

Shear modulus G 6.8 GPa

Constrained modulus Kp 14 GPa

Grain bulk modulus Ks 36.5 GPa

Fluid bulk modulus Kf 2.22 GPa

0.03 0.04 0.05 0.06 −4 −2 0 2 4x 10 −3 v3 [ms −1 ] 0.03 0.04 0.05 0.06 0.07 0.08 −1 0 1x 10 4 pf [Nm −2 ] t [ms] pSt pSt F(t) Full response pSt-pole

FIG. 3. Full response and pSt-pole residue at x2= x3= 0 and offset x1

= 0.10 m for Bentheimer/water configuration 1. The pSt-pole residue coin-cides with the pSt-waveform in the full response. Other wavemodes are too weak to be observed in this figure. The source signature F共t兲 is also dis-played. 0.03 0.04 0.05 0.06 −5 0 5 x 10−5 v3 [ms −1 ] 0.03 0.04 0.05 0.06 0.07 0.08 −0.02 −0.01 0 0.01 0.02 p f [Nm −2] t [ms] Full response pR-pole Pa-pole (x 1/200) P1 S pR P1 pR FS, P2S

FIG. 4. Full response and residues of pR-pole and P¯a-pole at x2= x3= 0 and

offset x1= 0.10 m for Bentheimer/air configuration 2. The P¯a-pole residue has been scaled down by a factor 200 to make it entirely visible. The

pR-pole residue coincides with the pR-waveform in the full response ofv3.

For configuration 1 共water as upper fluid兲, only the

pSt-pole is found on the principal Riemann sheet. From Fig.

3, we observe that its residue yields the entire pSt-waveform. For configuration 2共air as upper fluid兲, the pR-pole is found on the principal Riemann sheet. From Fig. 4, it is observed that its residue coincides with the pR-waveform in the full response of v3 共actually, the difference is nonzero but too

small to be observed兲. However, it does not coincide with that in the full response of pf. Its contribution is opposite, which means that the loop integrals along the branch cuts also contribute to the pseudo interface waveform. This was also found by Allard et al.11It implies that part of the perti-nent physical properties of the pseudo interface wave is cap-tured by the loop integrals.

This is more pronounced for configuration 3共light upper fluid兲, as shown in Fig.5, in which both the pR-pole and the

pSt-pole are found on the principal Riemann sheet. In both

components 共v3and pf兲, the pR-pole residue does not

coin-cide with the pR-waveform in the full response. The pSt-pole residue is not displayed separately because the pSt-wave strongly interferes with the F-wave.

To investigate how the residues and the interface wave-forms in the full responses compare for an observation point that lies off the interface, we calculated the responses for configurations 1 to 3 at x3= 0.01 m and offset r = x1= 0.1 m.

The corresponding results are displayed in Figs. 6–8. Com-pared to the previous responses at x3= 0, various head waves

can now be distinguished as separate arrivals, generated by the body wavefronts that propagate along the interface, cf. Fig.1共b兲. From Figs.6–8, we also observe that the residues now yield two waveforms in the full responses. The first one 共pR or pSt兲 is the waveform of the specific interface wave itself, while the second 共P2pR or P2pSt兲 corresponds to the P2-mode that is radiated by the propagating pseudo interface

wave. For configurations 1共water as upper fluid, Fig.6兲 and 2 共air as upper fluid, Fig.7兲, it is observed that both wave-forms are now captured entirely by the residue of the corre-sponding pole. For configuration 3共light upper fluid, Fig.8兲, this is not the case, as for x3= 0.

0.03 0.04 0.05 0.06 0.07 −2 −1 0 1 2 x 10−5 v3 [ms −1 ] 0.03 0.04 0.05 0.06 0.07 0.08 0.09 −50 0 50 p f [Nm −2] t [ms] Full response pSt-pole P1 _S P1 SP1 P2S P2P1 P2pSt P2pSt P2F P2P1 pSt

FIG. 6. Full response and pSt-pole residue at x2= 0, x3= 0.01 m, and offset

x1= 0.10 m for Bentheimer/water configuration 1. The pSt-pole residue

co-incides with the pSt- and P2pSt-waveforms in the full response. The

double-mode symbols are explained in Fig.1共b兲.

0.03 0.04 0.05 0.06 −5 0 5 x 10−6 v3 [ms −1] 0.03 0.04 0.05 0.06 0.07 0.08 0.09 −100 −50 0 50 100 p f [Nm −2] t [ms] P1 S P1 SP1 P2S P2P1 P2P1 pR P2pR P2S Full response pR-pole Pa-pole (x 1/525) P2pR

FIG. 7. Full response and residues of pR-pole and P¯a-pole at x2= 0, x3

= 0.01 m, and offset x1= 0.10 m for Bentheimer/air configuration 2. The

pR-pole residue coincides with the pR- and the P2pR-waveforms in the full

response. The P¯a-pole residue has been scaled down by a factor 525 to make it entirely visible. The double-mode symbols are explained in Fig.1共b兲.

0.03 0.04 0.05 0.06 −2 −1 0 1 2 x 10−4 v3 [ms −1 ] 0.03 0.04 0.05 0.06 0.07 0.08 −50 0 50 p f [Nm −2] t [ms] x 1/10 pR P1 S pR FS, P2S P1 F, pSt F, pSt Full response pR-pole

FIG. 5. Full response and pR-pole residue at x2= x3= 0 and offset x1

= 0.10 m for Bentheimer/light-fluid configuration 3. From the dashed verti-cal line onwards, the response pf has been scaled down by a factor 10 to make it entirely visible. The double-mode symbols are explained in Fig.

1共b兲. 0.03 0.04 0.05 0.06 0.07 −1 0 1 x 10−5 v3 [ms −1] 0.03 0.04 0.05 0.06 0.07 0.08 0.09 −200 −100 0 100 200 p f [Nm −2] t [ms] Full response pR-pole P1 S P2pR, pSt P1 SP1 P2S P2P1 P2P1 pR P2pR, pSt P2F, P2pSt P2F, P2pSt

FIG. 8. Full response and pR-pole residue at x2= 0, x3= 0.01 m, and offset

x1= 0.10 m for Bentheimer/light-fluid configuration 3. The double-mode

In addition to the observations on responses with entire waveforms, we give attention to the characteristics of a pseudo interface wave. With regard to the propagation veloc-ity, we observe that it is predicted properly by the residue of the corresponding pole for all presented numerical results. Concerning the attenuation, it was proposed by van der Hijden17 to quantify the true attenuation of a pseudo inter-face wave based on the full transient response. This is also done by Rosenbaum,5but he only showed the decay of the total waveform, which would result in one value for the at-tenuation. This is quite restrictive and therefore we use the following method to retrieve the frequency-dependent at-tenuation from a windowed pseudo interface waveform in the full response. Here, attenuation is defined by Im共sm

tr_兲,

where sm tr

represents the true wave slowness and m =兵pR,pSt其. As a starting point, we consider the pseudo in-terface wave in the far field where it does not interfere with other wavemodes, and we assume that it is described by

vˆm,3共r兲 ⬀ r−1/2exp共− i␻sm tr

r兲, 共27兲

which is found from the asymptotic behavior of the Hankel function.29The imaginary part of the wave slowness can be retrieved by comparing the amplitude spectra of the win-dowed waveform兩vˆm,3共r兲兩 at two different observation points r = raand r = rb, according to Im共sm tr_{共f兲兲 =} 1 2␲f共rb− ra兲 ln

冉

rb 1/2_兩vˆ m,3共rb兲兩 ra1/2兩vˆm,3共ra兲兩

冊

. 共28兲

For configurations 1共water as upper fluid兲 and 3 共light upper fluid兲, the attenuations are displayed in Figs. 9 and 10, re-spectively, together with the corresponding predictions ob-tained from the poles p = sm. The limited frequency range is due to the limited bandwidth of the retrieved spectra. For configuration 1 共water as upper fluid兲, we observe that the attenuation is described very well by the pSt-pole, except for the low frequencies where the far-field approximation of the Hankel function in Eq.共27兲is not valid. For configuration 3 共light upper fluid兲, however, the true attenuation of the

pR-wave is much greater than the value obtained from the

pR-pole residue. Obviously, the loop integrals along the

branch cuts cannot only affect the waveform but also the spatial decay of a pseudo interface wave.

Sometimes, a residue of a pole is共implicitly兲 considered to represent the corresponding interface-wave part of the spectrum of the Green’s function共see e.g., Refs. 7,12, and 14–16兲, while the loop integrals are considered to constitute the part related to body waves and head waves 共if present兲. This can be true but we emphasize that the choice of branch cuts is not unique. Therefore, the integration can be per-formed on another physically allowed Riemann sheet, i.e., a Riemann sheet that also meets the requirement of Im共qj兲 ⱕ0 for real p,2,4

which is the original path of integration关see Eq. 共20兲兴. This was done by Allard et al.11and Tsang,35 and clarified by Harris and Achenbach.36 Then, the construction of the 共x,␻兲-domain Green’s function is different as other poles have to be taken into account and different loop inte-grals are to be evaluated. Therefore, it might very well be that 共part of兲 the pertinent physical properties of a true or pseudo interface wave are captured by the integrals along the closed contour, rather than by the residue of a specific pole 共alone兲.

From the current observations, we conclude that a resi-due of a pole on the principal Riemann sheet does not nec-essarily yield all the pertinent physical properties of the cor-responding pseudo interface wave.

B. Presence of pR-pole and pSt-pole on Riemann sheets

Now we address the issue concerning the origin of a pseudo interface wave, as raised in point共2兲 of the Introduc-tion 共Sec. I兲. In the computations in Sec. V A, we already found that a pole related to a pseudo interface wave can be located on the principal Riemann sheet and, obviously, con-tribute a residue to the full response共see TableI兲. This con-tradicts the conventional explanation that a pseudo interface wave originates from a pole on a different Riemann sheet and is accounted for only by the loop integrals along branch cuts by causing a local maximum in the integrand. Allard et

0 0.2 0.4 0.6 0.8 1 1.2 1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 x 10−5 Im (spSt ) f [MHz] pSt-waveform pSt-pole

FIG. 9. True attenuation, defined as Im共spSt

tr _{兲, retrieved from windowed}

pSt-waveforms for Bentheimer/water configuration 1, at x2= x3= 0 and

fro-moffsets x1= 0.24– 0.26 m. The attenuation Im共spSt兲 obtained from the cor-responding pole residue is also displayed.

0 50 100 150 200 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0 x 10−5 Im (spR ) f [kHz] pR-waveform pR-pole

FIG. 10. True attenuation, defined as Im共spR

tr_{兲, retrieved from windowed}

pR-waveforms for Bentheimer/light-fluid configuration 3, at x2= x3= 0 and

from offsets x1= 0.24– 0.26 m. The attenuation Im共spR兲 obtained from the corresponding pole is also displayed.

al.11 already found this contradiction, but they did not refer to this as such because their concern was to determine whether or not a pole is related to a separate arrival in the full response.

Surprisingly, in case of fused glass beads saturated with water and covered with water共configuration 4, Fig.11兲 the

pSt-pole is present on the principal Riemann sheet only for a

limited frequency range. In Fig.12, the position of the pole in the complex plane is given, as expressed by Im共spSt兲. Also the position of the qP2-branch cut is displayed, as expressed by its imaginary part at Re共p兲=Re共spSt兲. As frequency in-creases, the pSt-pole moves toward the branch cut and as soon as it reaches the cut, it vanishes from the sheet. The pole is not present on the principal sheet for 310 kHz⬍ f ⱕ2 MHz. Therefore, the residue of the pSt-pole is not shown in Fig. 11. For the material properties used by Gubaidullin et al.,12 exactly the same situation occurs, al-though the transition takes place at a different frequency. Obviously, the presence of a pole on a certain Riemann sheet

is not only a matter of the contrast in material parameters of the half-spaces37but can also depend on frequency in case of viscous poroelastic media.

The behavior of the pSt-pole illustrates both the noncon-ventional and the connoncon-ventional explanation about the origin of a pseudo interface wave. The pole does contribute a resi-due over a certain frequency range and not outside that spe-cific range. For the pR-wave present in the full response of configuration 4共Fig.11兲, only the conventional explanation holds as the pR-pole is not found on the principal Riemann sheet and the entire waveform is captured by the loop inte-grals.

C. Physical significance of additional poles

Finally, we give attention to the physical significance of two additional zeroes of the poroelastic Stoneley dispersion equation 共⌬St= 0兲 as raised in point 共3兲 of the Introduction 共Sec. I兲. In configuration 2 共air as upper fluid兲, these zeroes show up as poles on the principal Riemann sheet at p = sP¯ a and p = sP¯ b. They are located to the left of the fast compressional-wave slowness关Re共sP¯ a,P¯b兲⬍Re共sP1兲兴 close to the qP1-branch cut共see Fig.2; p = sP¯ asignifies the pole that lies the closest to p = sP1兲. The additional 共P¯a, P¯b兲 poles are comparable with the so-called P¯ -poles that occur in nonpo-rous elastic solids with an interface, as described by Gilbert and Laster38 and Aki and Richards.4 The scaled real and imaginary parts of the poles are displayed in Fig. 13. The

P ¯

b-pole is only present on the principal Riemann sheet for limited frequency range 818.75 kHzⱕ f ⱕ2 MHz.

Allard et al.18also found one of the poles and refer to it as an improper surface mode. Feng and Johnson7stated that poles located to the left of shear-wave branch point 关Re共p兲 ⬍Re共sS兲兴 have lost all physical significance as pseudo inter-face modes. In the latter paper, the authors consider pseudo interface modes in the conventional way. In their configura-tion, the additional poles might indeed lie on a different Rie-mann sheet, but we find that they can also show up on the principal Riemann sheet. From Fig. 4 we observe that the

P ¯

a-pole has a substantial residue contribution to the full re-0.03 0.04 0.05 0.06 0.07 0.08 −5 0 5 x 10−5 v3 [ms −1 ] 0.03 0.04 0.05 0.06 0.07 0.08 −20 −10 0 10 20 p f [Nm −2] t [ms] pR P1 pR FS, P2S P1 F pSt x 1/10 pSt F S

FIG. 11. Full response at x2= x3= 0 and offset x1= 0.10 m for fused

glassbeads/water configuration 4. From the dashed vertical line onwards, the response pfhas been scaled down by a factor 10 to make it entirely visible. The double-mode symbols are explained in Fig.1共b兲.

0 50 100 150 200 250 300 350 −2 −1.8 −1.6 −1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2 x 10−4 Im (p ) f [kHz] qP2-branch cut pSt-pole

FIG. 12. The location of the pSt-pole in the complex p-plane, as expressed by Im共p兲=Im共spSt兲, with respect to the location of the qP2-branch cut at p = Re共spSt兲, as expressed by Im共p兲=Re共sP2兲Im共sP2兲/Re共spSt兲, for fused glass beads/water configuration 4. The 共–·兲 line indicates the frequency at the intersection. 0 0.5 1 1.5 2 0 0.5 1 1.5 Locat ion of add .p ol es f [MHz] Re

( )

sPa Re( )sP1 1 100Im

( )

sPa Im( )sP1 13Im

( )

sPb Im( )sP1 Re

( )

sPb Re( )sP1

FIG. 13. The location of the additional poles p = sP¯ aand p = sP¯ bin the com-plex p-plane with respect to the P1-wave slowness sP1. The P¯b-pole is only present on the principal Riemann sheet for limited frequency range 818.75 kHzⱕ f ⱕ2 MHz.

sponse, although it does not correspond to an interface wave-mode 共P¯b-pole similarly兲. Any pole that contributes to the full response should be considered as physically significant. Gilbert and Laster38 and Aki and Richards4 related the

P

¯ -poles in elastic solids to a separate arrival. Van der

Hijden,17however, stated that the concept of a separate pulse should be dismissed because it is just a peculiar tail to the compressional head-wave arrival. Harris and Achenbach36 confirmed this by stating that the poles yield features of the lateral waves. The observations in the current computations for poroelastic media also confirm this. From Fig.7, we ob-serve that the P¯a-pole contributes to the head waves gener-ated by the P1-wavefront. It also contributes to the P1-wavefront itself because it yields a strong pulse that

ar-rives even earlier 共Figs. 4 and 7兲, which is obviously ex-plained by the pole lying to the left of the compressional-wave slowness. The same is true for the P¯b-pole. The early-arriving parts are not present in the full responses and hence, the P1-waveform is constituted by both the residues of the P

¯ -poles and the loop integrals along the branch cuts. The fact

that a pole contributes to the P1-waveform illustrates that it

lies in the vicinity of the saddle point of the body wave, as used in asymptotic ray theory.35,39

There is one remarkable difference between the P¯ -poles in elastic and the ones in poroelastic media. In former, the poles never show up on the principal Riemann sheet4while this is possible for the latter. A similarity lies in the fact that in elastic solids共with rather small values of Poisson’s ratio兲, the poles lie also to the left of the compressional-wave slowness.35

VI. CONCLUSIONS

In this paper, we analyze the three-dimensional transient response of a fluid/poroelastic-medium configuration that is subjected to a vertical point force at the interface. For differ-ent materials, we quantitatively compare the full transidiffer-ent response with the residue contributions of pole singularities present on the so-called principal or physical Riemann sheet of integration. The poles are formed by zeroes of the po-roelastic Stoneley dispersion equation, i.e., the pseudo-Stoneley (pSt) and the pseudo-Rayleigh (pR) poles.

We find that the residues of these poles do not necessar-ily contain all pertinent physical properties of the corre-sponding pseudo interface waves. Part of them can be cap-tured by the loop integrals along the branch cuts. Therefore, it can be erroneous to use only the location of a zero of the Stoneley dispersion equation on the principal Riemann sheet, to predict the entire waveform, the propagation velocity, and attenuation of the corresponding pseudo interface wave.

According to the generally accepted explanation about the origin of a pseudo interface wave, it originates from a pole that lies on a nonprincipal Riemann sheet. The influence of the pole is only indirect in the sense that it causes a local maximum in the integrand of the Green’s function when its location is close to the real axis. We find, however, that this conventional explanation is not necessarily confirmed in the context of Biot’s theory for poroelasticity. The poles can

show up on the principal Riemann sheet. For the pSt-pole, we even show that its presence on the principal Riemann sheet is not only determined by the contrast in the material properties, but also by frequency.

Finally, we find that two additional zeroes of the po-roelastic Stoneley dispersion equation do have physical sig-nificance due to their residue contributions to the fast com-pressional wavefront and to the head waves that are radiated by this wavefront. In the literature the additional poles are, however, referred to as nonphysical because they are not related to pseudo interface waves. The poles are comparable with the P¯ -poles known in nonporous elastic solids, Refs.38 and 4. Depending on the specific material parameters and frequency, they can be present on the principal Riemann sheet or on another one.

ACKNOWLEDGMENTS

This research is supported by The Netherlands Research Centre for Integrated Solid Earth Sciences 共ISES兲. The au-thors are also grateful to E. C. Slob and A. V. Metrikine for the valuable discussions, and to B. Rossen for textual re-marks and suggestions.

APPENDIX A: TRANSFORM-DOMAIN RESPONSE In this Appendix, we derive the共p,x₃,␻兲-domain solu-tion to Eqs. 共5兲, 共6兲, and 共9兲–共14兲. Many of the involved symbols are explained in TableIII.

The general solution to the acoustic wave equation关Eq. 共9兲兴 in the 共p,x3,␻兲-domain can be readily found by

apply-ing the Fourier transform 关Eqs. 共15兲 and 共16兲兴 and solving the obtained ordinary differential equation. The result is

TABLE III. Symbols used in Appendix A The various indices are defined as j =兵P1, P2, F , S其, k=兵P1, P2其, and ᐉ=兵P1, P2, S其. Behind a number of

defini-tions, references are displayed where the specific expressions originate from.

␳ˆ₁₁ = ␳11− ibˆ/␻ Ref.23 ␳ˆ₂₂ = ␳₂₂− ibˆ/␻ Ref.23 ␳ˆ12 = ␳12+ ibˆ/␻ Ref.23 d0 = ␳ˆ11␳ˆ22−␳ˆ12 2 _Ref.23 d1 = −共R␳ˆ11+ P␳ˆ22− 2Q␳ˆ12兲 Ref.23 d2 = PR − Q2 Ref.23 sk2 = 共−d1⫿共d12− 4d0d2兲1/2兲/共2d2兲, Im共sk兲ⱕ0 Ref.23 sS2 = d0/共G␳ˆ22兲, Im共sS兲ⱕ0 Ref.23 sF2 = ␳F/KF p = 共p12+ p22兲1/2ⱖ0 qj = 共sj2− p2兲1/2, Im共qj兲ⱕ0 ␤ˆ k = −共␳ˆ11− Psk 2_兲/共␳ˆ 12− Qsk 2_{兲 Ref.23} ␤ˆ S = −␳ˆ12/␳ˆ22 Ref.23 A⬘ = A −共1−␾兲Q/␾ Ref.12 Q⬘ = Q −共1−␾兲R/␾ Ref.12 Hk = Q + R␤ˆk Ref.12 Kk = A⬘+ Q⬘␤ˆk+ 2G Ref.12 ␾ᐉ = 1 −␾+␾␤ˆᐉ Ref.12 ⌬1 = sP22HP2− sP12 HP1 ⌬2 = qP1sP22HP2− qP2sP12HP1 ⌬3 = −4p4␾␳ˆ22−1共qP1sP12 HP1− qP2s2P2HP2兲+4p2qSqP1qP2␾␳ˆ22−1⌬1 +2p2_s S 2_共q P1共␾P1+␾d2−1HP1KP2兲−qP2共␾P2+␾d2−1HP2KP1兲兲 − sS2G−1共qP1␾P1sP22 KP2− qP2␾P2s2P1KP1兲

p

˜F= i␻␳FA˜Fexp共+ i␻qFx3兲, x3⬍ 0, 共A1兲

where A˜F is the complex plain-wave amplitude of the fluid 共F兲 wave and qF=共sF2− p2兲1/2is the vertical slowness. It con-tains the wave slowness sF and the magnitude of the hori-zontal slowness p that are defined in TableIII.

The general solution to the Biot equations关Eqs.共5兲and 共6兲兴 can be derived by applying Helmholtz decomposition in the共x,␻兲-domain to these equations, according to12,23

uˆ =ⵜ␸ˆ_P1+ⵜ␸ˆ_P2+ⵜ ⫻ ␺ˆ, 共A2兲

Uˆ =␤ˆP1ⵜ␸ˆP1+␤ˆP2ⵜ␸ˆP2+␤ˆSⵜ ⫻ ␺ˆ, 共A3兲 where␸ˆ_P1 and␸ˆ_P2 denote the scalar potentials for the fast 共P1兲 and slow 共P2兲 compressional waves, respectively, and ␺ˆ

denotes the shear-wave共S兲 vector potential.␤ˆ_P1,␤ˆ_P2, and␤ˆS are the well-known fluid-solid共frame兲 amplitude ratios23for the separate body wavemodes共TableIII兲.

Applying the Helmholtz decomposition, the governing equations are decoupled and once the spatial Fourier trans-form关Eq.共16兲兴 is applied, the decoupled equations turn into ordinary differential equations for ␸˜_P1 and␸˜_P2, and ␺˜ that can be solved separately. The general solution for the dis-placements is obtained by adding the separate contributions according to Eqs.共A2兲and共A3兲. When the shear-wave term is split into a vertically polarized 共SV兲 and a horizontally polarized 共SH兲 part, the result for the wave vector w˜ =共u˜₁, u˜₂, u˜₃, −␾˜pf兲Tcan be written as

w˜ =

冤

p1 p1 qS p1 p sS 2p2 p2 p2 p2 qS p2 p − sS 2p1 p2 qP1 qP2 − p 0 − i␻s2_P1HP1 − i␻sP2 2 HP2 0 0

冥

⫻

冤

A ˜ P1exp共− i␻qP1x3兲 A ˜ P2exp共− i␻qP2x3兲 A ˜_{SVexp共− i␻qSx} 3兲 A ˜ SHexp共− i␻qSx3兲

冥

, x3⬎ 0. 共A4兲

Next to the solid displacements u˜ , the wave vector w˜ con-tains the fluid pressure p˜frather then the fluid displacements

U˜ because the four components of w˜ describe the wave field totally: there are only four independent variables, see Ref. 40. In Eq.共A4兲, A˜_P1, A˜_P2, A˜SV, and A˜SH denote the complex plain-wave amplitudes of the corresponding body wave-modes. In TableIII, the vertical slownesses qP1, qP2, and qS are defined共together with qF兲, as well as the fluid compress-ibility terms HP1and HP2.

The body-wave slownesses have Im共sj兲ⱕ0 and Som-merfeld’s radiation condition requires that Im共qj兲ⱕ0 for all body modes, j =兵P₁, P₂, F , S其.

The complex plane-wave amplitudes are determined by the boundary conditions at the interface x3= 0. Applying the

transforms 关Eqs. 共15兲 and 共16兲兴 to the boundary conditions 关Eqs. 共10兲–共14兲兴 and substituting the wave fields 关Eqs.共A1兲 and共A4兲兴, the following set of equations is obtained

冤

2Gp2_{− s} P1 2 KP1 2Gp2− sP2 2 KP2 0 2GpqS 0 sP1 2 HP1 sP2 2 HP2 −␾␳F 0 0 qP1␾P1 qP2␾P2 qF − p␾S 0 2pqP1 2pqP2 0 sS 2 − 2p2 ₊ p2 p1p qSsS 2 2pqP1 2pqP2 0 sS 2 − 2p2 ₋ p1 p2p qSsS 2

冥

⫻

冤

A ˜ P1 A ˜ P2 A ˜ F A ˜ SV A ˜ SH

冥

=

冤

Fˆ i␻ 0 0 0 0

冥

, 共A5兲

which is similar to that in Ref. 12, but extended to three dimensions. The constrained moduli KP1and KP2are defined in Table III. The solution is calculated analytically using

MAPLE©: A ˜ P1= − Fˆ i␻G⌬1⌬St 共␾␳FqP2共sS2␾P2− 2p2␾␳ˆ22−1sP22 HP2兲 + qF共sS2− 2p2兲sP22 HP2兲, A ˜ P2= Fˆ i␻G⌬1⌬St 共␾␳FqP1共sS 2_␾ P1− 2p2␾␳ˆ22 −1 s_P12 HP1兲 + qF共sS 2 − 2p2兲s_P12 HP1兲, A ˜ F= Fˆ i␻G⌬1⌬St 共共qP1␾P1sP2 2 HP2− qP2␾P2sP1 2 HP1兲 ⫻ 共sS 2 − 2p2兲 + 2p2⌬2␾S兲, A ˜ SV= 2pFˆ i␻G⌬₁⌬St 共␾2_␳ F␳ˆ22−1qP1qP2⌬1+ qF⌬2兲, 共A6兲

and A˜SH= 0. Here, the “poroelastic Stoneley-wave denomina-tor”共see Sec. III兲 is defined as

⌬St= qF⌬R+␾␳F⌬3/⌬1, 共A7兲

which is associated with interface waves along the fluid/ poroelastic-medium interface. It is very similar to the “Scholte-wave denominator” for a fluid/elastic-solid interface,27and equivalent to the one as given by Denneman

et al.41 It contains the “poroelastic Rayleigh-wave denomi-nator” that is associated with interface waves along a vacuum/poroelastic-medium interface

⌬R=共sS2− 2p2兲2+ 4p2qS⌬2/⌬1, 共A8兲

which is very similar to the one for a vacuum/elastic-solid interface.4,31

Now the plain-wave amplitudes are known, the 共p,x3,␻兲-domain solution to Eqs. 共5兲, 共6兲, and 共9兲–共14兲 is

determined and given by Eqs.共A1兲 and共A4兲. APPENDIX B: INVERSE FOURIER INTEGRAL

In this Appendix, we show how Eq.共19兲can be written in terms of a single integral according to Eq.共20兲, following Ref. 4. Transforming Eq.共19兲 to cylindrical coordinates ac-cording to p₁= p cos␸, p₂= p sin␸, and x₁= r cos␽, x₂ = r sin␽, where r =共x₁2+ x₂2兲1/2, it can be written as

gˆ+= ␻ 2 共2␲兲2

冕

0 ⬁

冕

0 2␲_˜n+_共p,␸,x 3,␻兲 ⌬St ⫻exp共− i␻pr cos共␸−␽兲兲pd␸dp. 共B1兲

The ␸-dependence of n˜+ _{can be replaced by 共horizontal兲}

partial-derivative operators⳵_␣,␣=兵1,2其, since the factors p_␣ that appear in n˜+ _{关see Eqs. 共A4兲 and 共17兲兴 correspond to}

horizontal derivatives 共−i␻p_␣↔⳵_␣兲 in the 共x,␻兲-domain.

Therefore, n˜+_{共p,␸, x}

3,␻兲 is defined such that it contains the

appropriate derivative operators according to gˆ+= ␻ 2 共2␲兲2

冕

0 ⬁_n˜+_共p,⳵ ␣,x3,␻兲 ⌬St ⫻

冕

0 2␲ exp共− i␻pr cos共␸−␽兲兲d␸pdp = ␻ 2 2␲

冕

0 ⬁_n˜+_共p,⳵ ␣,x3,␻兲 ⌬St J0共␻pr兲pdp, 共B2兲

where we used the integral representation of the zeroth-order Bessel function J0共¯兲, see Ref. 42. The Bessel function is

replaced by the sum of two zeroth-order Hankel functions of the first and second kinds,42 i.e., J0共z兲=

1 2共H0

共1兲_共z兲+H 0 共2兲_共z兲兲.

Using the equality H₀共1兲共z兲=−H₀共2兲共−z兲 and the evenness of the 共p,x3,␻兲-domain Green’s functions in p, Eq. 共B2兲 can be

written as gˆ+= ␻ 2 4␲

冕

−⬁ ⬁ _˜n+_共p,⳵ ␣,x3,␻兲 ⌬St H₀共2兲共␻pr兲pdp, 共B3兲

where the horizontal derivatives are applied to the Hankel function before the integration is performed, according to

⳵␣H0共2兲共␻pr兲 = −␻p x_␣

r H1

共2兲_{共␻pr兲. 共B4兲}

1_{D. R. Burns, “Acoustic waveform logs and the in-situ measurement of}

permeability—A review,” in Geophysical Applications for Geotechnical

Investigations, edited by F. L. Paillet and W. T. Saunders共ASTM,

Phila-delphia, 1990兲.

2_{I. A. Viktorov, Rayleigh and Lamb Waves: Physical Theory and}

Applica-tions共Plenum, New York, 1967兲.

3_{W. M. Ewing, W. S. Jardetzky, and F. Press, Elastic Waves in Layered}

Media共McGraw-Hill, New York, 1957兲.

4_{K. Aki and P. G. Richards, Quantitative Seismology共Freeman and Co.,}

New York, 1980兲.

5_{J. H. Rosenbaum, “Synthetic microseismograms: Logging in porous}

for-mations,” Geophysics 39, 14–32共1974兲.

6_{H. Deresiewicz, “The effect of boundaries on wave propagation in a}

liquid-filled porous solid: IV. Surface waves in a half-space,” Bull. Seis-mol. Soc. Am. 52, 627–638共1962兲.

7_{S. Feng and D. L. Johnson, “High-frequency acoustic properties of a fluid/}

porous solid interface. I. New surface mode,” J. Acoust. Soc. Am. 74, 906–914共1983兲.

8_{S. Feng and D. L. Johnson, “High-frequency acoustic properties of a fluid/}

porous solid interface. II. The 2D reflection Green’s function,” J. Acoust. Soc. Am. 74, 915–924共1983兲.

9_{M. J. Mayes, P. B. Nagy, L. Adler, B. P. Bonner, and R. Streit, “Excitation}

of surface waves of different modes at fluid-porous solid interface,” J. Acoust. Soc. Am. 79, 249–252共1986兲.

10_{L. Adler and P. B. Nagy, “Measurements of acoustic surface waves on}

fluid-filled porous rocks,” J. Geophys. Res. 99,共B9兲, 17863–17869 共1994兲.

11_{J. F. Allard, M. Henry, C. Glorieux, W. Lauriks, and S. Petillon,}

“Laser-induced surface modes at water-elastic and poroelastic interfaces,” J. Appl. Phys. 95, 528–535共2004兲.

12_{A. A. Gubaidullin, O. Y. Kuchugurina, D. M. J. Smeulders, and C. J.}

Wisse, “Frequency-dependent acoustic properties of a fluid/porous solid interface,” J. Acoust. Soc. Am. 116, 1474–1480共2004兲.

13_{D. L. Johnson, J. Koplik, and R. Dashen, “Theory of dynamic}

permeabil-ity and tortuospermeabil-ity in fluid-saturated porous-media,” J. Fluid Mech. 176, 379–402共1987兲.

14_{I. Edelman and K. Wilmanski, “Asymptotic analysis of surface waves at}

vacuum/porous medium and liquid/porous medium interfaces,” Con-tinuum Mech. Thermodyn. 14, 25–44共2002兲.

15_{B. Albers, “Monochromatic surface waves at the interface between}

po-roelastic and fluid half-spaces,” Proc. R. Soc. London, Ser. A 462, 701– 723共2006兲.

16_{M. G. Markov, “Low-frequency Stoneley wave propagation at the}

inter-face of two porous half-spaces,” Geophys. J. Int. 177, 603–608共2009兲.

17_{J. H. M. T. van der Hijden, “Quantitative analysis of the pseudo-Rayleigh}

phenomenon,” J. Acoust. Soc. Am. 75, 1041–1047共1984兲.

18_{J. F. Allard, G. Jansens, G. Vermeir, and W. Lauriks, “Frame-borne surface}

waves in air-saturated porous media,” J. Acoust. Soc. Am. 111, 690–696 共2002兲.

19_{J. F. Allard, M. Henry, C. Glorieux, S. Petillon, and W. Lauriks,}

“Laser-induced surface modes at an air-porous medium interface,” J. Appl. Phys.

93, 1298–1304共2003兲.

20_{M. A. Biot, “Theory of propagation of elastic waves in a fluid-saturated}

porous solid. I. Low-frequency range,” J. Acoust. Soc. Am. 28, 168–178 共1956兲.

21_{M. A. Biot, “Theory of propagation of elastic waves in a fluid-saturated}

porous solid. II. Higher frequency range,” J. Acoust. Soc. Am. 28, 179– 191共1956兲.

22_{M. A. Biot and D. G. Willis, “The elastic coefficients of the theory of}

consolidation,” J. Appl. Mech. 24, 594–601共1957兲.

23_{J. F. Allard, Propagation of Sound in Porous Media: Modelling of Sound}

Absorbing Materials共Elsevier Applied Science, London, 1993兲.

24_{D. M. J. Smeulders, R. L. G. M. Eggels, and M. E. H. van Dongen,}

“Dynamic permeability: Reformulation of theory and new experimental and numerical data,” J. Fluid Mech. 245, 211–227共1992兲.

25_{H. Deresiewicz and R. Skalak, “On uniqueness in dynamic poroelasticity,”}

Bull. Seismol. Soc. Am. 53, 783–788共1963兲.

26_{B. Gurevich and M. Schoenberg, “Interface conditions for Biot’s equations}

of poroelasticity,” J. Acoust. Soc. Am. 105, 2585–2589共1999兲.

27_{A. T. de Hoop and J. H. M. T. van der Hijden, “Generation of acoustic}

waves by an impulsive line source in a fluid-solid configuration with a plane boundary,” J. Acoust. Soc. Am. 74, 333–342共1983兲.

28_{K. Attenborough, S. I. Hayek, and J. M. Lawther, “Propagation of sound}

above a porous half-space,” J. Acoust. Soc. Am. 68, 1493–1501共1980兲.

29_{M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions}

共Dover, New York, 1972兲.

30_{B. A. Fuchs, B. V. Shabat, and J. Berry, Functions of a Complex Variable}

and Some of Their Applications共Pergamon, Oxford, 1964兲.

31_{J. D. Achenbach, Wave Propagation in Elastic Solids共North-Holland,}

Am-sterdam, 1973兲.

32_{P. J. Davies and P. Rabinowitz, Methods of Numerical Integration}

共Aca-demic, New York, 1975兲.

33_{N. Ricker, “Wavelet contraction, wavelet expansion, and the control of}

seismic resolution,” Geophysics 18, 769–792共1953兲.

34_{C. J. Wisse, “On frequency dependence of acoustic waves in porous}

cyl-inders,” Ph.D. thesis, Delft University of Technology, The Netherlands 共1999兲.

35_{L. Tsang, “Time-harmonic solution of the elastic head wave problem}

in-corporating the influence of Rayleigh poles,” J. Acoust. Soc. Am. 63, 1302–1309共1978兲.

36_{J. G. Harris and J. D. Achenbach, “Comment on ‘On the complex}

conju-gate roots of the Rayleigh equation: The leaky surface wave’,” J. Acoust. Soc. Am. 112, 1747–1748共2002兲.

37_{J. G. Scholte, “The range of existence of Rayleigh and Stoneley waves,”}

Mon. Not. R. Astron. Soc. 5, 120–126共1947兲.

38_{F. Gilbert and S. J. Laster, “Excitation and propagation of pulses on an}

interface,” Bull. Seismol. Soc. Am. 52, 299–319共1962兲.

39_{B. L. van der Waerden, “On the method of saddle points,” Appl. Sci. Res.,}

Sect. B 2, 33–45共1952兲.

40_{G. Bonnet, “Basic singular solutions for a poroelastic medium in the}

dy-namic range,” J. Acoust. Soc. Am. 82, 1758–1762共1987兲.

41_{A. I. M. Denneman, G. G. Drijkoningen, D. M. J. Smeulders, and K.}

Wapenaar, “Reflection and transmission of waves at a fluid/porous-medium interface,” Geophysics 67, 282–291共2002兲.

42_{I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series and Products}