On effective characteristics of wave propagation in a porous fluidsaturated medium containing entirely fluid inclusions

S U M M A R Y To investigate the behaviour of seismic wave propagation in natural rocks and soils, a model for the elastic properties of a heterogeneous porous medium containing uniform fluid inclusions is useful. The problem has already been addressed before, utilizing Biot’s theory of poroelasticity, and there have been some analytical and numerical results obtained for the different types of cavities. However, as it is shown in this paper, immediate passage within heterogeneous porous media from spherical porous heterogeneities to continuous fluid or solid inclusions cannot be performed on the basis of previously developed models. In this paper, we extend the results for wave scattering in porous media due to a random ensemble of spherical porous inclusions to allow for continuous fluid inclusions instead. To consider this case, we first obtain the single scattering series. Then we utilize multiple scattering theory to find the effective wavenumber due to an ensemble of inclusions. Our results make it possible to analyse wave attenuation and velocity dispersion at low frequencies due to the presence of entirely fluid inclusion heterogeneities within a porous material. Explicit analytical expressions for solution of single scattering problem for a fluid cavity, and effective propagation velocity and attenuation factor in a porous medium with random ensemble of these cavities constitute main results of this work.


I N T RO D U C T I O N
Recent publications (Gurevich & Lopatnikov 1995;Gelinsky et al. 1998;Pride & Berryman 2003;Ciz & Gurevich 2005;Muller & Gurevich 2005a,b) have shown that the presence of heterogeneities within a porous fluid-saturated medium can cause significant wave attenuation and velocity dispersion. This is a result of wave-induced fluid flow occurring between regions having different elastic compliance. There were number of different models proposed to describe wave propagation in media with mesoscopic heterogeneities (Toms et al. 2006). Of our particular interest is the so-called discrete random model of entirely fluid inclusions.
Numerically, the behaviour of longitudinal waves in media with random ensembles of fluid-filled cavities was investigated before for broad frequency range (Krutin et al. 1984;Markov et al. 2005). Moreover, the role of surface tension on the elastic waves scattering in an inhomogeneous poroelastic medium was analysed in the paper by Markov & Levin (2007).
However, there were no explicit analytical results obtained for this case. In this work, we derive the effective wavenumber for a heterogeneous porous material containing non-porous regions consisting entirely of a continuous fluid. The effective wavenumber yields the important properties: the inverse quality factor characterizing attenuation and the effective wave propagation velocity.
To obtain our results, we employ a two-step procedure. We first solve the single scattering problem for a plane elastic wave incident on a spherical inclusion (Berryman 1985;Ciz & Gurevich 2005). Then we follow Ciz et al. (2006) by applying multiple scattering theory (Foldy 1945;Waterman & Truell 1961) to obtain the effective wavenumber for a porous material containing an ensemble of randomly distributed spherical fluid volumes.
The paper is organized as follows. In Section 2, we introduce a simple model of a porous medium and define all notations subsequently utilized throughout the paper. Then, in Section 3, the single scattering problem is formulated. Here we also discuss the mathematical transition 1044 D. V. Ponomarev and O. V. Nagornov from fluid-saturated porous inclusion to non-porous fluid-filled cavity. We utilize the limit passage for the expressions published by Ciz & Gurevich (2005), and partly solve the problem analytically in a straightforward way to obtain the second piece of the result replacing formulas derived for the purely porous case (Berryman 1985), which do not admit passing to the entirely fluid inclusion limit. All obtained results are then verified numerically and demonstrated on the plots. In Section 4, we outline basic concepts of the multiple scattering theory as applicable to our problem, deriving the effective wavenumber for the compression wave of the first kind in a manner similar to other publications (Markov 2005;Ciz et al. 2006). This gives explicit expressions for inverse quality factor characterizing attenuation and effective propagation velocity in the media. Although the main results of this work given in Sections 3 and 4 are also formulated in alternative notation that is commonly used in the cited papers, in addition, in the appendix, we list expressions showing the equivalence between our poroelastic formulation and other approaches.

A S I M P L E P O RO E L A S T I C M E D I U M M O D E L
The fluid-saturated porous medium consists of a framework of solid grains and a connected pore space filled entirely by fluid. The porosity of the porous medium is denoted by β (see Table 1 for the definition of β and all other notations used in this paper) and is equivalent to the ratio of the fluid volume divided by the total volume of the fluid-saturated porous medium.
The porous medium is a continuum having two displacement vectors at each point (details on such an averaging process are discussed in Molotkov (2001)): u refers to solid matrix displacement, while U refers to fluid displacement. Usually, w = β( u − U ) is defined which characterizes relative solid-fluid displacement.
The stress tensor for the porous material is given by where s = −βp is the force acting on the fluid component of the surface of the cubic volume, σ i j is the stress tensor acting on the framework (i.e. the solid component of the surface of the cubic volume), p is the fluid pressure, σ 0 i j is the stress tensor of the solid grains and δ i j is the Kronecker delta. We assume the porous media is isotropic; the stress and strain relations are a generalization of the Hooke's law given by where e i j is the solid strain tensor and e = ∇ · u, ε = ∇ · U are dilatations of the solid matrix and fluid, respectively. A, N, Q and R are coefficients listed in the appendix and defined in terms of elastic moduli (Biot & Willis 1957). The paper of Biot (1962) provides a rigorous derivation of the stress-strain relations listed above. However, Biot's derived results are self evident when considering that the fluid component can only make an isotropic contribution to the Hooke's law for a solid elastic material.
To obtain dynamic poroelastic equations, let us introduce the kinetic energy of unit volume, which we can write as a positively defined quadratic form of velocities T = ρ 11 2 u 2 x +u 2 y +u 2 z + ρ 12 u xUx +u yUy +u zUz + ρ 22 2 U 2 x +U 2 y +U 2 z and the simplest dissipation function where b = β 2 η/κ is a constant chosen to have the agreement with the Darcy law, η is fluid viscosity, and κ is permeability. Although we write the dissipation function here in the simplest form making it turn to zero when there is no relative fluid speed, it is eligible in the low-frequency range only. Now we can write Lagrange equations with dissipation ∂ ∂t where L = T − W is the Lagrange function and W is strain energy. When we substitute L, we should take into account equalities ∂ W and ∂ W ∂U i = ∂s ∂ x i , both sides of which give i-components of total forces per unit volume acting on solid and fluid, respectively. Finally, we obtain equations governing the motion of the porous medium Expressions for A, N, Q, R, ρ 11 , ρ 12 and ρ 22 in terms of measurable medium quantities and alternative coefficients are given in the appendix.

1045
To simplify these equations, let us express displacements in terms of scalar and vector potentials separating the time-harmonic multiplier u = e iωt (∇ϕ + ∇ × ψ), By applying divergence and curl operations on (3) and substituting the newly defined displacements, we obtain ⎧ ⎨ ⎩ χ = ψ, and respectively, with the following notations As one can see, the equations for scalar potentials (5) are coupled. To decouple them, we can perform simple linear algebra transformations introducing new independent variables If we choose γ 1 = D − a + (D − a) 2 + 4Bc/(2c), γ 2 = D − a − (D − a) 2 + 4Bc/(2c), then we will have two independent equations (written as one using commas in the subscript) where k 2 1,2 = −a − γ 1,2 c. It is easy to show that at low frequencies the wavenumbers behave as |k 1 | ∝ ω, |k 2 | ∝ √ ω, |Im k 1 | << |Re k 1 |, |Im k 2 | ≈ |Re k 2 | (to be more accurate we should write Im k 2 ≈ −Re k 2 for chosen time-harmonic dependence exp(iωt)). Thus, there are two types of longitudinal waves propagating in the fluid-saturated porous medium. Besides the fast wave with weak dispersion and attenuation, there exists the so-called Biot's slow wave, which is highly dissipative and has strong dispersion. The existence of this second type of compression wave predicted by Biot's poroelasticity theory (Biot 1956a(Biot ,b, 1962 was afterwards confirmed experimentally (Plona 1980). A slow wave demonstrates out-of-phase behaviour in solid-fluid motion and has a dissipative nature.
We should note one more thing: the equations obtained are eligible only under an assumption of low frequencies. However, the value of the critical frequency ω c = βη/(κρ f ) is high enough, so this limitation does not affect the most practical applications in seismology.

S C AT T E R I N G O N A S I N G L E S P H E R I C A L I N C L U S I O N
To utilize multiple scattering theory, it is necessary first to solve the single scattering problem. That is, the objective is to determine the amplitudes of reflected elastic waves, when a plane elastic wave is incident on a single spherical inclusion (Fig. 1). The single scattering problem in elastic media was considered by Yamakawa (1962). Then results were obtained for porous media (Berryman 1985) for small inclusion size and afterwards they were extended to a broader range of applicability by Ciz & Gurevich (2005).
In this section, we formulate the problem in terms introduced earlier, that is using potentials F 1 , F 2 and ψ, assuming harmonic time dependence giving by multiplier exp(iωt). Also here, as we utilize results obtained in another publication (Ciz & Gurevich 2005), we establish a link with the formulation in terms of dilatations and a time multiplier inverse to ours.
The single scattering problem is formulated as follows.
Let the incident fast plane wave amplitude be A 01 and the spherical inclusion radius R a , which we assume to be small enough to satisfy the condition |k 1 R a | 1. Choosing spherical coordinates that are more appropriate for the problem, the eqs (6) and (4) for external porous  Porosity, β 1.6 × 10 −1 Permeability, κ 9.5 × 10 −13 (m 2 ) Dynamic viscosity, η 1.0 × 10 −3 (Pa · s) Tortuosity factor, s 3 Density of solid matrix material, ρ s 2.76 × 10 3 (kg/m 3 ) Density of fluid in porous matrix, ρ f 1.0 × 10 3 (kg/m 3 ) Density of fluid filling cavity,ρ f 0.7 × 10 3 (kg/m 3 ) Shear modulus of solid matrix material, μ 3.4 × 10 9 (Pa) Bulk modulus of fluid in porous matrix, K f 2.2 × 10 9 (Pa) Bulk modulus of drained porous medium, K b 5.8 × 10 9 (Pa) Bulk modulus of solid matrix material, K s 4.0 × 10 10 (Pa) Bulk modulus of fluid filling cavity,K f 1.5 × 10 9 (Pa) Radius of cavity, R a 1.0 × 10 −3 (m) Here ψ is the azimuthal component of the vector potential ψ = {0, 0, ψ}, chosen in such form due to the axial symmetry of the problem. Following Ciz & Gurevich (2005), we first consider wave scattering from a single spherical porous inclusion embedded in a background porous medium having different physical properties. Later in the section we will derive results appropriate for a spherical inclusion, which is entirely fluid (i.e. non-porous).
We will denote all parameters inside this porous inhomogeneity with arc symbols. So the equations inside the porous inclusion (r < R a ) are absolutely similar to those written earlier We may write the general solutions for these wave equations as the spherical harmonics expansions (for instance, Landau & Lifshits 1989) where we simply write F instead of F 1 , F 2 , F 1 and F 2 , and instead of ψ and ψ; P l (cos θ) are the Legendre's polynomials. Solutions inside the inclusion must be regular at r = 0, thus two constants must be dependent so that a combination of exponents forms a sinc function term sin kr kr . Wavefields at r → ∞ must satisfy the attenuation conditions. Taking into account Im k 1 < 0 and Im k 2 < 0, it yields that only e −ikr kr terms must present. Also there is the Sommerfeld radiation condition ∂ψ ∂r + ik s ψ = o( 1 r ) at r → ∞ that leads us to a conclusion that there must be no e iks r ksr term in shear wavefield. Allowing for all these conditions and using expansion of the incident wave as follows: Tikhonov & Samarskiy 2004), we obtain total wavefields where we prolonged the equalities to write the expressions in a short form referring to the spherical Bessel and Hankel functions notations. We should note periodical far-field behaviour of spherical Hankel functions (Abramowitz & Stegun 1965) h Therefore, to truncate a series of scattered waves, we need to show smallness of higher harmonic amplitudes. All constants mentioned can be found from the interface conditions (r = R a ), which are the continuity of both normal and tangential total stresses, the equality of fluid pressures and the continuity of solid normal and tangential displacements and normal relative fluid displacements (Deresiewicz & Skalak 1963): 1048

D. V. Ponomarev and O. V. Nagornov
Using the expressions for stresses, pressures and displacements in terms of the potentials, we develop each of these interface conditions: Substituting here the wavefields given above and taking into account the independence property of the Legendre's polynomials, we come to the set of algebraic equations for the amplitudes of the harmonics.
It is not feasible to obtain the exact solutions because all the series are infinite, but it is possible to obtain the effective scattered wavefield by truncating the series assuming |k 1 R a | << 1.
The numerical solution of this set at low frequencies reveals that all terms in the scattered fields starting with l = 3 can be neglected because of higher degrees of amplitude-frequency dependence. Indeed, we can easily find out by drawing plots in logarithmic scale and calculating angular coefficients (that are multiples of 0.5) which amplitudes are significant and which ones can be neglected. Thus, the effective scattered wavefields that are of interests read as follows: Analytical expressions for scattering amplitudes were earlier obtained by Berryman (1985) and Ciz & Gurevich (2005) but they used dilatation terms e = ∇ · u, ξ = β∇ · ( u − U ) instead of potentials like we do to have all the equations in scalar form, writing the scattered wavefields in the form (the correspondences for H , C and further mentioned elastic constants K , M, μ, λ and σ can be found in the appendix).

1049
The other difference with the cited authors is that they initially implied time-dependence given by the multiplier exp(−iωt), which results in second-to-first-kind replacement of spherical Hankel functions in comparison to this paper where we assume exp(iωt) time-dependence.
At low frequencies, the Berryman's solution for scattered field amplitudes gives where ξ + = k + a = k 1 R a , ξ − = k − a = k 2 R a are dimensionless wavenumbers and A 0 is the amplitude of the incident wave; all quantities marked with a prime are related to the inclusion. As we do, Ciz and Gurevich also stressed that higher harmonics were negligible and that the l = 1 harmonic did not contribute to the scattered slow wavefield. Here are the analytical solutions they obtained: where ρ = βρ f + (1 − β)ρ s is the overall density of two-phase porous medium, and ρ f and ρ s are the fluid and matrix solid (grain) densities, respectively.

Taking into account
, we may establish correspondence with the notation we use through this paper. A where l = 0, 1, 2, ... . We should also allow for differences between the amplitudes of incident waves written in terms of potentials and dilatations; this gives A 01 = A 0 (γ 1 −γ 2 ) γ 2 k 2 1 . Thereby, we come to the expressions helping to verify the equality of results obtained using the different approaches for an lth harmonic. Next, we proceed with scattering on a continuous fluid obstacle. Expressing internal ideal fluid displacements in terms of potentials Ũ = e iωt ∇φ, we may apply the standard procedure to reduce the Euler equationρ f ∂ ṽ ∂t = −∇p and the linearized continuity equation ∂ρ ∂t = −ρ f ∇ · ṽ to the simple wave equation governing motion inside the inclusion (r < R a ) φ +k 2 lφ = 0, wherek 2 l = ω 2 c 2 l is the fluid wavenumber. Following the same steps as in the case of porous inclusion, we write the total wavefields as Interface conditions (r = R a ) are continuity of both normal and tangential total stresses, equality of fluid pressures (capillary effects in the present work are not taken into account), and continuity of normal relative fluid displacements (1 − β)u r + βU r =Ũ r .
The substitutions into each of these conditions yield a corresponding equation for potentials Thus, after plugging in all of the wavefields here, we end up with a set of algebraic equations for the harmonic amplitudes. Because there is no shear wave in an ideal fluid, we have only three equations for the zero harmonic amplitudes. Although calculations are quite tedious and too cumbersome to be given here, the following low-frequency results can be obtained in a straightforward way: l is the bulk modulus of a fluid filling the inclusion. These expressions can be developed further taking into account that for low frequencies we have γ 1 ≈ R+Q A+2N +Q and γ 2 ≈ −1. Thus, we arrive at Also we give here the same results expressed in terms of an alternative formulation of the problem These expressions are verified with the numerical solution of the set of eqs (25)-(28) in the low-frequency range (Figs 2 and 3). For higher harmonics it is not feasible to obtain analytical expressions directly. Therefore, the perturbation method was used by Berryman (1985) and Ciz & Gurevich (2005) to utilize the results of the scattering problem in a non-porous medium published by Yamakawa (1962). In this way the expressions for B + 1 , B − 2 , B + 2 given earlier were obtained. Mathematical passage to the limit of entirely fluid inclusion can be carried out as follows.  First, one should consider the non-porous solid limit of the expressions (19)-(21) by setting C → 0, M → 0, as was mentioned by Berryman (1985). This is implemented as β → 0, K f → 0, (K f /β ) → 0, K b → K s in the expressions for elastic medium moduli (A8)-(A12).
Then, to perform the transition to the fluid, we simply replace all material characteristics of the solid with those of the fluid filling the cavity, in particular setting μ = 0.
Berryman has been using complete (i.e. in both the inclusion and the background medium) transition to the non-porous medium limit to compare his results with those of Yamakawa (1962). Obviously, it cannot be done only for the medium of the inclusion leaving the background medium still porous. Indeed, one can immediately notice that the expression (17) has M in the denominator, and therefore is not defined after passing the limit. It follows that the expression (18) cannot be computed either. However, the expressions (19)-(21) are still eligible once the limit has been passed. We can directly perform it as follows: ρ =ρ f , Then, in terms of this work we have .
The validity of these formulas can be ensured by comparison with results obtained by means of the numerical solution of the set of linear algebraic equations for amplitudes (Figs 4-6). Small discrepancies for the amplitudes of the first and second harmonics are explained by the inexactness of the original expressions (19)-(21) that were derived in the cited papers using perturbation technique.
Numerical analysis also shows that expressions for the effective wavefields have the same form as those for scattering from porous media inclusion. That is to say, the infinite series can be truncated to the same terms giving the main contribution to the total field, only the harmonic amplitudes are different.
Thus, effective scattered fields of compression waves are where the amplitudes A (1)

S C AT T E R I N G O N M U LT I P L E I N C L U S I O N S
Having computed the effective wavefields, we can now proceed with considering a medium that contains a big ensemble of N (one should not confuse it with the notation of an elastic constant we introduced before) randomly distributed spherical inclusions that we will refer to as scatterers here. We assume scatterers' positions are independent and their density in the medium is low enough so we can use the obtained effective wavefields which are valid in the far field of each scatterer.  If we denote the scatterers' density as n( r ), the number of scatterers in an elementary volume d r will be d N( r ) = n( r )d r . The probability of finding a single scatterer in this volume is d P( r ) = d N( r ) N = n( r) N d r and the corresponding probability density at r is p( r ) = n( r) N . Then the overall probability density for N scatterers located at r 1 , ..., r N is p overall ( r 1 , ..., r N ) = p( r 1 ) · ... · p( r N ) = n( r 1 )·...·n( r N ) N N (as we noted earlier we consider the situation where all scatterers are independent so there is no correlation in their positions).
We showed that in a porous medium without scatterers' longitudinal wave motion was governed by Helmholtz equations of the form F 0 + k 2 0 F 0 = 0 with a constant k 0 .

D. V. Ponomarev and O. V. Nagornov
Now, involving N inclusions in the medium, we have the total wavefield as a superposition of fields scattered from each of them where E( r , r j ) = e −ik 0 | r − r j | | r − r j | . Each coefficient A j is defined by an incident field F j ( r ) acting on a jth scatterer. So denoting the reflection coefficient as T we have A j = T F j ( r j ).
According to the superposition principle we may write Combining (38) and (39), closed set of equations can be obtained for all j = 1, ..., N . However, because the scatterers' positions are random, we should resort to an averaging process of the total wavefield (38) and obtain where F j ( r j ) j is the field acting on a jth scatterer which is averaged over the positions of the rest. Due to a big number of scatterers N , we may assume that F j ( r j ) j ≈ F( r j ) . Also, because all scatterers are independent we are able to replace the sum in (41) with a multiplication by N and hereby derive an integral equation for the averaged total wavefield For this integral equation an exact solution can be found. Indeed, let us apply a ( + k 2 0 ) operator to the both sides of eq. (42) and notice that we will have the Dirac delta function under the integral sign.
Validity of the expression ( is the Green's function for the radiation problem, that is, G(r ) = − e −ik 0 | r − r | 4π| r − r | is the solution of ( + k 2 0 )G( r , r ) = δ( r − r ) which satisfies the Sommerfeld condition at infinity.
Thus, we get ( + k 2 0 ) F( r ) = ( + k 2 0 )F 0 ( r ) − 4π T n( r ) F( r ) or allowing for ( + k 2 0 )F 0 ( r ) = 0 we can finally write it as where is the effective wavenumber for the heterogeneous porous medium. The procedure that we have just followed is based on the paper by Foldy (1945), where the isotropic property of each scatterer was assumed. Though, as we showed earlier, in our case besides the zero harmonic there are also, generally speaking, the first and the second harmonics which cannot be neglected. However, according to a later publication (Waterman & Truell 1961), if the scatterers' density may be considered constant (i.e. n( r ) = n 0 ), a similar expression for the effective wavenumber is valid: where f (0) and f (π ) are total forward and backward scattered amplitudes.
Substitution of this result into (46) yields an effective wavenumber for the compression wave of the first kind: Because the amplitudes of the harmonics are found, we can explicitly write real and imaginary parts of (48) and thereby study attenuation and dispersion.
Attenuation is characterized by the inverse quality factor Effective velocity can be calculated as , that is, for low density of scatterers equivalent to where v 1 = ω k 1 is velocity of the fast compression wave. To get explicit expressions for these quantities, it remains to substitute the harmonic amplitudes given by (29), (33) and (34) into (49) and (50). This yields (52) The last expression can be also written in a slightly more compact form using notations utilized in the papers Ciz & Gurevich (2005) and Ciz et al. (2006)

C O N C L U S I O N S
In this paper, we have derived results that allow the study of elastic wave propagation through porous media containing entirely fluid non-porous spherical inclusions. We started by introducing general concepts of wave propagation in porous fluid-saturated media appropriate for low-frequency wave propagation. Then we solved the single scattering problem for a plane elastic wave propagating through a porous medium incident on a single spherical poroelastic inclusion having different physical properties. It was done by expressing the solutions in terms of spherical harmonic expansions, analysing the low-frequency behaviour of the harmonic amplitudes and then formulating the effective scattered wavefield at distances far from the inclusion.
Before, the analytical solutions for an effective scattered field were obtained only for the case of a porous medium spherical inclusion (Berryman 1985;Ciz & Gurevich 2005).
In this work, previous results (Ciz & Gurevich 2005) were partly used to obtain solutions for scattering on an entirely fluid spherical cavity by means of mathematical transition from a porous to a purely fluid inclusion case in expressions for higher harmonics. Also, the details and peculiarities of this passage to the limit of non-porous inclusion were pointed out. Namely, this transition cannot be performed for expressions for zero-harmonics amplitudes published by Berryman (1985), therefore these expressions were substituted with those newly derived and constitute another part of the results of this work.
All the results were verified numerically and are shown on the plots. Then, according to Foldy (1945) and Waterman & Truell (1961), we applied the single-scattering results to describe wave propagation in a porous medium with an ensemble of randomly distributed inclusions by introducing the effective wavenumber concept. Similar work has been done by Markov (2005) and Ciz et al. (2006), but for a medium containing only porous inclusions. For the case of fluid cavities only numerical results have been obtained before Markov et al. 2005;Markov & Levin 2007).
Obtained analytical results for the effective wavenumber yield explicit expressions for attenuation (characterized by inverse quality factor) and effective velocity of propagation of fast compression waves in the heterogeneous media containing a random ensemble of low-concentrated, entirely fluid, non-porous spherical inclusions.
where the new quantity ρ a was introduced. Additional mass ρ a can be expressed in terms of fluid density, porosity and the so-called tortuosity ρ a = βρ f (s − 1) where 1 ≤ s ≤ 5 (Zwikker & Kosten, 1949;Dutta & Ode, 1979). Below we give the correspondence between Biot's original elastic coefficients and those used in the cited papers: