Two-Phase Flow Effects on Seismic Wave Anelasticity in Anisotropic Poroelastic Media

: We study the wave anelasticity (attenuation and velocity dispersion) of a periodic set of three ﬂat porous layers saturated by two immiscible ﬂuids. The ﬂuids are very dissimilar in properties, namely gas, oil, and water, and, at most, three layers are required to study the problem from a general point of view. The sequence behaves as viscoelastic and transversely isotropic (VTI) at wavelengths much longer than the spatial period. Wave propagation causes ﬂuid ﬂow and slow P modes, inducing anelasticity. The ﬂuids are characterized by capillary forces and relative permeabilities, which allow for the existence of two slow modes and the presence of dissipation, respectively. The methodology to study the physics is based on a ﬁnite-element uspcaling approach to compute the complex and frequency-dependent stiffnesses of the effective VTI medium. The results of the experiments indicate that there is higher dissipation and anisotropy compared to the widely used model based on an effective ﬂuid that ignores the effects of surface tension and viscous ﬂow interference between the two ﬂuid phases.


Introduction
Wave anelasticity in porous media is a topic with applications in many areas of geophysics and material science [1]. Waves generate fluid flow in the pores and energy losses that can be observed in field and laboratory experiments [2][3][4]. Recent laboratory experiments performed in the seismic range have shown the frequency dependence of anelasticity in sandstones with partial gas or oil saturations [5][6][7], while experiments conducted in Reference [8] show a significant attenuation in the extensional and bulk deformation modes, as well as numerical simulations in close agreement with laboratory data.
The theoretical pioneering work of Biot [9][10][11] describes waves propagation in porous media saturated by single-phase fluids. The theory predicts a shear wave and two compressional waves, a fast one, where the solid and fluid move together, and a slow one, where the displacement is out of phase. At low frequencies, the slow wave is diffusive and becomes a propagating wave at high frequencies. Significant anelasticity of the fast wave is due to mode conversion at mesoscopic scale heterogeneities (mesoscopic loss). White et al. [12] were the first to introduce this loss mechanism based on Biot's theory for a periodic sequence of two flat and thin porous layers saturated with a single-phase fluid. This mechanism is also termed WIFF or wave-induced fluid flow loss.
The study of propagation in partially saturated porous media has been presented in Reference [13][14][15]. Lo et al. [16] present an Eulerian model for waves propagating in porous media saturated by two fluids. In this model, the stress-strain relations are obtained by taking into account capillary pressure assumed as a unique function of saturation and ignoring hysteresis, while the dynamic equations are formulated from a mass balance equation for each solid and fluid phase. Three compressional waves are predicted to exist, one wavelike and two of them of diffusive type. The model is applied to analyze waves in soils saturated by water with air or oil. A macroscopic model based on a two-component Biot model and a mixture theory is presented in Reference [17]. Other models to analyze the wave behavior in porous rocks saturated my immiscible fluids appeared in [1,15,[18][19][20][21][22].
Concerning studies on porous media by using numerical simulations, the work by Thovert et al. [23] analyzes wave propagation using an homogenization approach. The authors consider the existence of connected pores that may or may not percolate, obtain the macroscopic coefficients and apply the results to a synthetic porous medium. Hamzehpour et al. [24] study acoustic waves in two-dimensional fractured porous media using finite differences to discretize the differential equations. This analysis concludes that, near the source, waves amplitude decay exponentially. A generalization of the Biot theory, appearing in Reference [25], takes into account more than one fluid phase saturating the pores, such as brine-gas and oil-gas mixtures. This theory predicts three compressional waves and a shear wave. The stress-strain relations are derived from the complementary virtual work principle in order to include the capillary pressure relation, while the dissipation function is defined in terms of Darcy's law based on two fluids [26].
The presence of an additional slow wave is expected to induce more attenuation and velocity dispersion, as shown in Reference [27], for the case of vertical propagation across a periodic set of two thin layers saturated with two immiscible fluids. This work considers the more general case of a sequence of three layers, also leading to a VTI behavior at long wavelengths, but to a different frequency dependence of the velocity and attenuation. Qi et al. [28] consider capillary effects in random media of patchy saturation by using a membrane stiffness, but this approach leads to an increase in phase velocities and lower attenuation. An experimental study on how capillarity affects the compressional P-wave velocity for patchy saturation during an imbibition process is presented by Liu et al. [29], where the patch size depends on the fluid saturation. The effect of capillarity is related to the injection rate, changing with imbibition. The results exhibit an increase of the experimentally measured P-wave static velocity as capillary pressure increases.
This study is based on the theory of Cavallini et al. [30] valid for periodic N layers saturated for single-phase fluids and a FE procedure to obtain an equivalent VTI medium to an horizontally layered medium consisting of a three periodic poroelastic solid saturated by immiscible fluids. An FE procedure is used to determine a VTI medium equivalent at long wavelengths to a periodic sequence of thin poroelastic layers. It consists of defining five independent boundary-value problems, each one associated with either a compressibility or shear time-harmonic test, whose solutions are obtained by the FE method. The results are validated against the theory by Cavallini et al. [30] for effective single-phase fluids. Then, several examples are given, where the velocities and attenuation are obtained by using effective-single phase and two-phase fluids.

The Differential Model
The poroelastic medium is saturated by non-wetting and wetting fluid phases, whose properties and variables are indicated by subscripts or superscripts "n" and "w", respectively. The corresponding saturations are denoted by S w = S w (x) and S n = S n (x), with x = (x, y, z), respectively, with associated residual saturations S rw and S rn . It is assumed that both fluid phases occupy the whole pore space and that a continuous network of these phases exists (funicular regime) [31]. Thus, The Fourier transforms of the particle displacement of the three phases composing the material, i.e., the solid, non-wettting and phases phases, are denoted as u s = (u s i ), u n = (ũ n i ), andũ w = (ũ w i ), i = 1, 2, 3, respectively. Set u = (u s , u n , u w ).
as the matrix effective porosity, with the relative displacement and variation in fluid content of each fluid phase defined as Define P w and P n as the Fourier transforms of infinitesimal changes in the wetting and non-wetting fluid pressures, respectively, with respect to an initial equilibrium state of pressuresP w ,P n , porosityφ, and non-wetting saturationS n . The capillary relation is [26,31,32] P ca = P ca (S n +S n ) =P n + P n − (P w + P w ) = P ca (S n ) + P n − P w ≥ 0. (1) Ignoring hysteresis, P ca is an increasing function of S n and positive. Let ε ij (u s ) and e s = ε ii (u s ) be the Fourier transforms of the solid strain tensor and corresponding linear invariant, respectively. Denoting with τ ij (u) the stress-tensor components of the bulk material, the constitutive relations are: where In (2), λ u = K u − 2 3 µ, with K u and µ being the wet-rock bulk and dry-rock shear moduli, respectively. Thus, is the wet-rock P-wave modulus. The coefficients in (2)-(4) can be determined as indicated in Appendix A. The model describing the response of a poroelastic medium saturated with two fluids in the diffusive range of frequencies consists of a partial differential equations imposing the stress equilibrium of the bulk material (Equation (6)), together with a generalization of the two-phase Darcy law [26,31,32] (Equations (7) and (8)). Thus, if ω denotes the angular frequency, these equations are The coefficients in (7) and (8) depend on the absolute permeability κ, fluid viscosities η l and relative permeabilities K rl (S l ), l = n, w. They are defined as The coefficient d nw (S n ) in (10) is a dissipative function, where is small. It describes the viscous drag between the immiscible fluids. In the absence of experimental data, it has the form given in (10).

The Equivalent Viscoelastic Transversely-Isotropic Medium
For long wavelengths, compared to the layer thicknesses, a fluid-saturated porous medium is seen as an effective or equivalent VTI medium characterized by five frequency dependent and complex coefficients that can be determined by means of harmonic experiments performed on a representative sample of the medium. These experiments, given next for the 2D case, are defined as boundary value problems (BVPs), whose solutions are obtained by using an FE method. For single-phase fluids, the analytical solutions given in Appendix B are used to validate the numerical simulations.
Let us consider a representative squared sample Ω = (0, L) 2 with boundary Γ in the (x, z)-plane, where x and z denote the coordinates along the horizontal and vertical directions. Let Γ R , Γ L , Γ T , and Γ B denote the right, left, top, and bottom boundaries of Ω. In addition, let χ be a unit tangent on Γ oriented counterclockwise, and ν the unit outer normal on Γ.
Let E ( u s ) and T ( u s ) be the Fourier transforms of the strain and stress tensors at the macroscale, with u s denoting the macroscopic displacement vector of the solid. The stress-strain equations for a VTI medium can be stated as follows: where p 12 = p 11 − 2p 66 . To compute the five independent frequency dependent and complex stiffness coefficients in (11)-(16), we solve Equations (6)- (8) in Ω imposing no change in fluid content of both fluid phases, i.e., with the boundary conditions and additional boundary conditions stated below for each coefficient p I J .
The numerical experiment to determine p 13 is defined by applying the boundary conditions: In this experiment, Then, ε 11 and ε 33 do not vanish, and E 11 and E 33 can be obtained as Next, from Equations (11) and (13), it follows that Note that it follows from (27) that T 11 = T 33 = −∆P, so that using (31) p 13 is determined as The following boundary conditions are used to determine p 55 : where This BVP satisfies the conditions ε ij = 0 only for i = 1, j = 3, so that E ij = 0 only for i = 1, j = 3, and (15) reduces to Next, by computing the average of the local strain ε 13 over the sample the stiffness p 55 can be determined from (35). Finally, p 66 is computed by using the boundary conditions: where Then, p 66 is determined as indicated for p 55 . The five BVPs are solved with the FE method. The components of the displacement vector of the solid are represented with locally bilinear polynomials, while global continuity is imposed on Ω. Furthermore, only the normal components of the vector displacements of the two fluids need to be continuous across the common edges of adjacent computational cells. This condition is imposed by using polynomials linear in the x-direction and constant in the z-direction to represent the first component of each fluid phase, while polynomials constant in x and linear in z are used for the second component.

Results. Numerical Experiments
The solution of the five BVPs to determine the complex stiffnesses p I J stated in Section 3 is obtained as a function of the propagation direction and frequency with the FE method. Appendix A gives analytical expressions for the phase and energy velocities and dissipation factors of the qP, qSV, and SH waves.
In all the experiments, the numerical samples are discretized with a 90 × 90 uniform mesh representing six periods, each one with three layers of equal 20 cm thickness, saturated with two fluids. The material properties of the solid and fluid phases are listed in Tables 1 and 2. The behavior of the two fluids is described with relative permeabilities, K rn (S n ) and K rw (S n ), and capillary pressure function, P ca (S n ), defined by [33][34][35]: with A in (41) being the capillary pressure amplitude. The wetting phase is brine, and the non-wetting phase is oil or gas. In all the experiments, the residual saturations are S rw = 0.01, S rn = 0, the capillary-pressure amplitude is 30 kPa, and = 0.01 in the definition of the coefficient d nw in ( [26,31,32,36] for additional information on flow of two fluids in porous media). Capillary pressure and relative permeability functions may be obtained from log-well data, as indicated in Reference [37].
The experiments show dissipation factors and phase and energy velocities of waves computed by using two and one fluid phases. The properties of the effective fluid (viscosity η ( * ) , density ρ ( * ) , and bulk modulus K ( * ) ) are obtained with Reuss averages for the bulk moduli and arithmetic averages for densities and viscosities: The following cases are analyzed in the experiments • Case 1: Six periods of three layers (L 1 , L 2 , L 3 ) of layer thickness 20 cm, where the non-wetting fluid saturations in each layer are L 1 : 1% gas, L 2 : 5% gas, and L 3 : 5% oil. The curves labeled "single-phase" are obtained by using the classical Biot theory and effective single-phase fluids. • Case 2: Six periods of three layers (L 1 , L 2 , L 3 ) of layer thickness 20 cm, where the non-wetting fluid saturations in each layer are L 1 : 1% oil, L 2 : 5% oil, and L 3 : 5% gas.
The following notation is used in the figures: cpII and Q II = Re(p II )/Im(p II ) are the phase velocity and dissipation factor of the waves traveling parallel and perpendicular to the layering, "11-waves" and "33-waves", respectively. Figures 1a,b and 2a,b show the phase velocities and corresponding dissipation factors of "33-waves" and "11-waves" as a function of frequency for Cases 1 and 2 computed with the FE method. The plots compare the results for two fluids and one (effective) fluid, where the latter is defined in Equations (42).
At low frequencies, the more realistic case of two-phase fluids has lower velocities than for single-phase fluids, with similar asymptotic values at higher frequencies. As expected, phase velocities for "33-waves" are lower than for "11-waves". The dissipation factors corresponding to two-phase fluids ( Figure 2) exhibit a very different behavior as compared with the effective single-phase fluids, since there is an additional attenuation peak at lower frequencies, much stronger for Case 2. This new phenomenon can be explained as follows. Ignoring the cross-dissipative coefficient d nw , the particle velocity of the l-phase iωu l , l = o, g, b is equal to the gradient of the generalized fluid pressure T l times the absolute permeability κ and the l-phase mobility γ l (S l ), defined as γ l (S l ) = [K rl (S l )]/η l [26]. For the values of brine, oil and gas saturation is used in the examples, at 1% (5%) oil saturation, and γ o is five (three) orders of magnitude lower than γ b and γ g . Thus, fluid flows differently within the pore space, with the fluids interfering each other and inducing additional energy losses, which are not taken into account by the single-phase effective-fluid approach. Furthermore, Case 2 has more oil content than Case 1, thus inducing higher attenuation (see Figure 2). Moreover, the attenuation peaks for "33-waves" have higher amplitudes than for "11-waves", a known effect in finely layered porous media. An additional result observed in the experiments is that, at a given saturation, varying the amplitude of the capillary pressure does not affect the attenuation.  Figures 4 and 6 show a perfect agreement between the single-phase analytical and FE dissipation factors of the qP and qSV waves (SH waves are lossless), but the attenuation corresponding to the two-phase fluids, including capillarity effects, is higher, in agreement with the previous remarks regarding Figure 2.        (7) and (8). These vertical components, corresponding to Case 2, are denoted by v n,z and v w,z in the plots. The higher mobility γ b , as compared with the gas and oil mobilities γ g and γ o , explains the higher values of the wetting particle velocity in Figure 9 when compared with those of the nonwetting phase in Figure 8. As stated above, the interference between the two fluid phases due to their different particle velocities induces energy losses and velocity dispersion. Finally, Figure 10 exhibits the modulus of the gradient of the total fluid pressure T = T n + T w at 10 Hz for Case 2, which is clearly observed at the layer interfaces. This pressure gradient is another useful illustration of the WIFF mechanism for the case of two-phase fluid saturation.

Discussion
Wave-induced fluid flow in porous media is known to be a major cause of attenuation and velocity dispersion. Concerning the existence of additional slow waves when multiphase fluids saturate a porous medium, Albers [17] developed a macroscopic linear model to study wave propagation in unsaturated soils. This model, which assumes that capillary pressure and relative permeabilities depend on saturation, predicts the existence of three compressional waves, one of them fast and two slow, and a shear wave, with phase velocities and attenuation depending on saturation, and the coefficients in the stress-strain relations derived as in Reference [25].
Furthermore, Lo and Sposito [16] present differential equations to model dilatational waves in porous media saturated with two immiscible fluids. Their analysis predicts three compressional waves (P1, P2, P3), with the attenuation of the P1 wave highly affected by the pore fluids. On the other hand, the attenuation of the P2 and P3 waves is found to be related to the inverse of the sum of the relative mobilities of the two fluids, thus dominated by the fluid with larger relative mobility. This model also assumes that capillary pressure is a unique function of saturation (ignoring hysteresis), with the dissipation coefficients are defined in terms of the absolute permeability and the saturation-dependent relative permeability functions using the form given in Reference [25]. Quoting these authors "the model equations by Brutsaert [13], Garg and Nayfeh [14], Berryman et al. [15] and Tuncay and Corapcioglu [19] can be considered as special cases (of Reference [16]) and Santos et al. [25] (model) as a version based on a Lagrangian framework".

Conclusions
Finite-element time-harmonic quasistatic experiments were used to study the wave anelasticity of a periodic set of three thin and flat layers saturated with two fluids, focusing on the effects of surface tension (capillarity forces). The results are compared with the case in which the two fluids are replaced by an effective single-phase one. This later case is validated against the analytical solution, whereas there is no analytical solution for the more general and realistic scenario involving capillary forces. The results show that, in this case, the loss is higher for both waves traveling normal and parallel to the layers, and there are two attenuation peaks, as opposed to one peak for the effective single-phase fluid.
The energy velocities and dissipation factors of the three waves, shown as polar plots at 10 Hz versus the propagation angle, exhibit a perfect agreement between the analytical and FE single-phase fluid solutions. While the energy velocities coincide in all cases, the dissipation factors are always higher for the case of two-phase fluids, as mentioned above. The physics can be explained by combined effects of capillary pressure and different mobilities of each fluid phase in the two-phase Darcy law.
The existence of an additional slow wave compared to the classic slow wave of the Biot theory is due to surface tension effects represented by a saturation-dependent capillary pressure. On the other hand, relative mobilities induce relative motions between the fluids, causing losses absent when an effective fluid is considered, as it is mostly the case in the literature regarding wave propagation in porous media. Future work includes a more detailed analysis of the influence of wettability, phase mobilities, and saturation on wave anelasticity of rocks saturated with immiscible fluids.  matrix is symmetric and the flow is normal to the layering. Each layer is defined by the porosity: φ, grain and fluid bulk moduli: K s and K f , dry-rock bulk and shear moduli: K m and µ, grain and fluid densities: ρ s and ρ f , fluid viscosity: η and permeability: κ.
The White model was generalized by Krzikalla and Müller [38] to anisotropic media, to obtain the five stiffnesses, p I J , of the symmetric 6 × 6 matrix of the effective transversely isotropic medium, where c r I J and c I J are the relaxed and unrelaxed values, i.e., the low-and high-frequency values. Gelinsky and Shapiro [39] (their Equations (14) and (15)) give the expressions of these stiffnesses, where λ m = K m − (2/3)µ, E m = K m + (4/3)µ, and the brackets indicate the average ζ = d −1 ∑ i d i ζ i . In the unrelaxed regime, there is no interlayer flow, and the stiffnesses are c 66 = c r 66 , is the Gassmann P-wave modulus. Because the relaxed and unrelaxed shear moduli coincide, there is no shear dissipation along and normal to layering. The qSV wave suffers attenuation because it is coupled with the qP wave, and the SH wave is not dispersive, since c 55 = c r 55 and c 66 = c r 66 imply p 66 = c 66 and p 55 = c 55 . A set of thin layers saturated with different fluids but with the same shear modulus is isotropic. The average density of the medium is simply ρ = ρ .