Numerical simulation of ultrasonic waves in reservoir rocks with patchy saturation and fractal petrophysical properties

We simulate the propagation of ultrasonic waves in heterogeneous poroviscoelastic media saturated by immiscible ﬂuids. Our model takes into account capillary forces and viscous and mass coupling effects between the ﬂuid phases under variable saturation and pore ﬂuid pressure. The numerical problem is solved in the space–frequency domain using a ﬁnite element procedure and the time–domain solution is obtained by a numerical Fourier transform. Het-erogeneities due to ﬂuid distribution and rock porosity–permeability are modeled as stochastic fractals, whose spectral densities reproduce saturation an petrophysical variations similar to those observed in reservoir rocks. The numerical experiments are performed at a central frequency of 500 kHz, and show clearly the effects of the different heterogeneities on the amplitudes of shear and compressional waves and the importance of wave mode conversions at the different interfaces.


Introduction
During the past few years an important number of theoretical and experimental studies have demonstrated the significant influence that spatial heterogeneities within porous rocks have on elastic moduli, seismic wave velocities and other related quantities.In the context of hydrocarbon exploration, understanding the relation between the seismic response, the distribution pore fluids and the petrophysical properties may provide useful information about the reservoir.
As pointed out in [22], the heterogeneous nature of porous rocks often results in the heterogeneity of fluid distribution on scales greater than pore or grain size.In a rock saturated by immiscible fluids (such as water and gas) at macroscopic scales two simplified models of fluid distribution are generally considered in most published works in the subject: (i) the proportion of both fluids within the pores is the same everywhere, i.e., homogeneous saturation, and (ii) the fluids are arranged in patches, i.e., macroscopic regions (which may include thousands of grains), fully saturated with one of the two fluids.
However, these must be considered as limiting cases since saturation usually exhibits irregularities of different scale and in most cases there exists a residual saturation [21] as well as capillary forces.In many cases the existence of patchy saturation in reservoir sandstones is closely related with variations in lithology and clay content, which may cause small effects on elastic properties but have a very important influence on permeability and capillary pressure curves [22,47].
The computation of bulk elastic moduli and compressional wave velocities for this situation was studied by different authors by means of empirical relations [8] and effective medium theories.This latter approach is generally based in the validity of Gassmann's equations [24] within each patch under the assumption that their characteristic length is significantly smaller than a wavelength [22,31].It is also assumed that the fluids are distributed in patches of 100% saturation of either gas or water.
Using a different approach, White [45] analyzed the physics of wave propagation through patchy partially-saturated porous media using Biot's theory [5].This model considers spherical gas pockets embedded in a water-saturated porous medium.White found that one of the main attenuation/dispersion mechanisms is the conversion of fast P wave to slow modes, diffusive or wave-like, depending on the frequency range.Two other important mechanisms are the Biot loss, generally occurring beyond the sonic frequency range, and scattering attenuation, whose relaxation frequency depends on the size of the patches and frequency content of the source.The simulation of acoustic wave fields in this kind of media was recently treated in [11,27].
While the analysis of the acoustic response of reservoir rocks for heterogeneous saturation has been studied by different authors, to our knowledge, the influence of spatially variable petrophysical properties on the wave fields has not received much attention yet.This encouraged us to develop a numerical model to investigate the influence of all these type of heterogeneities on wave propagation in porous saturated reservoir rocks.
The use of stochastic fractals for the statistical description of variable scale heterogeneities arising in porous media has become widely accepted for flow studies and reservoir simulations [12,46].This approach will be adopted to model the spatial variability in saturation and porosity-permeability.
Our model is based on a generalization of Biot's theory [5][6][7] for porous rocks saturated by two immiscible fluids [36,39,40].It takes into account the existence of capillary forces at the pore scale, assuming that each fluid phase has a continuous distribution within the tortuous pore space and that both fluids can flow (funicular saturation regime).The formulation includes viscous and mass coupling interaction coefficients between the solid and fluid phases under variable saturation and pore fluid pressure conditions and frequency dependent dissipative effects associated with viscous forces and viscoelasticity.This theory predicts the existence of one shear wave and three compressional waves (one fast and two slow compressional waves).It has been found that for the low frequency range, the two slow modes are diffuse waves due to viscous effects and for very high frequencies, such as those used in laboratory testing, these two modes behave as propagating waves [39,40].
The numerical modeling of wave propagation in dissipative media can be efficiently done in the space-frequency domain because the solution at a given time can be obtained without the knowledge of the time history of the system.The generalized Biot model used in this work, besides introducing linear viscoelastic behavior in the rock frame (that may be stated and solved in the space-time domain using time-convolutions), employs frequency dependent mass and viscous coupling coefficients that are better described in the space-frequency domain, which is consequently the natural domain for setting our problem.Following the ideas given in [16,18,25,26] the space-frequency domain solution is computed for a finite number of temporal frequencies, and the spacetime domain solution is obtained via an approximate inverse Fourier transform.Numerical simulation of waves in this type of media is computationally very expensive due to the large number of degrees of freedom needed to represent the different types of waves accurately.This is due to the fact that there is at least one order of magnitude difference between the phase velocities of the fast and compressional waves, as can be seen in table 1.Thus employing parallel architectures combined with an iterative nonoverlapping domain decomposition algorithm in a hybridized formulation is an efficient way to approach the problem [14,15,17,26].To approximate the displacement vectors we use a nonconforming rectangular finite element [19] for the solid phase and the vector part of the Raviart-Thomas-Nedelec mixed finite element space of order zero for the two fluid phases, which are conforming spaces, see [33,37].There are two advantages when using this type of nonconforming elements, one is that considerably reduces the amount of information exchanged among processors, and the other is that almost halves the number of degrees of freedom needed to have a prescribed level of numerical dispersion [48], as compared when using standard conforming bilinear elements.At the artificial boundaries of the computational domain the first order absorbing boundary condition presented in [44] is employed.Future work includes the derivation of higher order absorbing boundary conditions for generalized Biot media.This algorithm is applied to perform numerical experiments in a real sandstone saturated by two fluids (gas-water).In previous works we dealt with homogeneous models [10,41,42], in which we verified numerically the propagation of fast and slow waves in this kind of media.
The structure of the paper is as follows.First, we give a brief description of the model and the equations of motion and formulate the iterative finite element procedure.Next we perform numerical experiments to analyze the partition of wave energy at the plane interfaces defined by an abrupt change in fluid saturation and by a jump in porosity-permeability.Then we consider the cases of patchy saturation distribution and fractal porosity-permeability, using different degrees of heterogeneity.The analysis of the snapshots obtained at ultrasonic frequencies indicate that important wave mode conversions and scattering effects take place at the plane interfaces and at the heterogeneities defined by the patchy fluid distribution and the fractal porosity-permeability fields.The attenuation-dispersion effects are function of the correlation length and the level of heterogeneity present in the rock.We also found a clear correlation between the location and size of gas patches and amplitude maxima of the wavefronts.

Description of the model
In a porous solid saturated by two immiscible fluids, we consider a wetting phase and a non-wetting one, which will be indicated with the subscripts (or superscripts) "w" and "n", respectively.Let S w and S n denote the averaged wetting and non-wetting fluid saturations, respectively, with S rw and S rn being the corresponding residual saturations, whose physical significance is as follows.S rw is the amount of wetting fluid that will always remain in the pore space even at very high capillary pressures when the wetting fluid is being displaced by the nonwetting fluid (drainage regime).On the other hand, when the nonwetting fluid is being displaced by the wetting fluid (imbibition), at zero capillary pressure a certain amount S rn of nonwetting fluid remains; this is the residual saturation of the nonwetting fluid.At S n < S rn the nonwetting phase ceases to flow [2].
We assume that the two fluid phases completely saturate the porous part of the bulk material so that S n +S w = 1.We further assume that the flow is in the funicular saturation regime, in which each fluid phase occupies a continuous network of tortuous (funicular) paths where simultaneous flow of both fluids is possible, so that S rn < S n < 1 − S rw [2,35,43].Let u s , u n and u w be the Fourier transform of the averaged displacement vectors for the solid and fluid phases at the angular frequency ω.
The governing equations for wave propagation in the space-frequency domain as given in [42] are: Here S n and S w denote the reference fluid saturations while f s , f n and f w represent the external source in the solid and fluid phases, respectively.The coefficient ρ is the density of the bulk material given by ρ = (1 − φ)ρ s + φ(S n ρ n + S w ρ w ), φ is the matrix effective porosity, ρ s is the mass density of the solid grains and ρ n and ρ w are the mass densities of the two fluids.
The mass coupling coefficients g n , g w , g nw represent the inertial effects associated with dynamic interactions between the three different phases, while the coefficients d n , d w and d nw include the viscous coupling effects between the solid and fluid phases.Let K, K rn (S n ), K rw (S n ) and K rnw (S n ) denote the absolute and relative permeability functions, respectively and set Then we take these coefficients to be of the form: ) The functions K rn (S n ), K rw (S n ) and K rnw (S n ) and the parameter ² are chosen as indicated in section 4. The constants µ n , µ w are the fluid viscosities and G = 1 2 (1 + 1/ φ).The complex valued frequency dependent function F (θ l ) = F R (θ l ) + iF I (θ l ), l = n, w, nw is the frequency correction function defined by Biot [6] in the high-frequency range: , T(θ l ) = ber 0 (θ l ) + i bei 0 (θ l ) ber(θ l ) + i bei(θ l ) , with ber(θ l ) and bei(θ l ) being the Kelvin functions of the first kind and zero order and θ l = a l p √ ωρ l /µ l , a l p = 2 p KK rl A 0 / φ, where A 0 denotes the Kozeny-Carman constant [29,42].
Let P w and P n denote infinitesimal changes in the pressures of the wetting and non-wetting fluids, respectively, with respect to corresponding reference values P n and P w associated with the initial equilibrium state with corresponding nonwetting fluid saturation S n .Recall that P n and P w are related through the capillary relation [2,35,43] Based on experimental data and ignoring hysteresis, the function P ca is a positive and strictly increasing function of the non-wetting fluid saturation.The stress-strain relations are given by Ravazzoli et al. [36]: where τ , T n and T w are the generalized stresses of the system, β = P ca (S n )/P 0 ca (S n ), ζ = P w /P 0 ca (S n ), ε ij and e b = ε ii are the Fourier transforms of the strain tensor of the solid and its linear invariant, respectively and ξ θ = −∇ • u θ , for θ = n, w.
The coefficients N , λ c , B 1 , B 2 , M 1 , M 2 , M 3 are the elastic moduli of the medium, and can be determined as follows.The coefficient N is the shear modulus of the dry rock, while λ c = K c − 2 3 N in 3D space and λ c = K c − N in 2D space, with K c being the undrained bulk modulus.Following [39] K c is computed using the formulae where K m , K s , K n and K w are the bulk modulus of the empty matrix, the solid grains and the non-wetting and wetting fluid phases, respectively, with corresponding compressibilities C l = K −1 l , l = m, s, n, w.The remaining coefficients can be obtained by using the following relations: ) with To introduce viscoelasticity we use the correspondence principle stated by Biot [4,7], i.e., we replace the real poroelastic coefficients in the constitutive relations by complex frequency dependent poroviscoelastic moduli satisfying the same relations as in the elastic case, with some necessary thermodynamic restrictions.In this work we use the linear viscoelastic model presented in [32] to make the undrained bulk modulus K c and the shear modulus N complex and frequency dependent, while all other coefficients in (2.6) are real.Thus, we take (2.10) The coefficients K * c and N * are reference values of the closed bulk and shear moduli properly chosen to fit high frequency velocities usually measured in laboratory.The functions R l and T l , l = K c , N, associated with a continuous spectrum of relaxation times, characterize the viscoelastic behaviour and are given by [9,32] The model parameters b Q l , T 1,l and T 2,l are taken such that the quality factors Q l (ω) = T l /R l are approximately equal to b Q l in the range of frequencies where the equations are solved, which makes this model convenient for geophysical applications.
The plane wave analysis performed in [40] shows that in these type of media, three different compressional waves (P 1 , P 2 and P 3 ) and one shear wave (or S-wave) can propagate.The P 1 wave is the classical fast P-wave propagating in elastic or viscoelastic isotropic solids, while the P 2 and P 3 waves are slow waves strongly attenuated in the low frequency range, corresponding to motions out of phase of the solid and fluid phases.In the high-frequency range they become truly propagating modes.
For simplicity in notation, the spatial dependency of all coefficients is omitted throughout the text.However in the numerical applications we will consider spatial variations in saturation, porosity and permeability, which induce heterogeneity also in bulk density, elastic moduli, mass coupling and viscous drag coefficients.This, in turn, produces local variations in phase velocities and intrinsic attenuation, modifying the patterns of wave energy propagation.

The domain decomposition iteration
We consider the solution of equations (2.1)-(2.3) in a rectangular poro-viscoelastic domain Ä in the (x, z)-plane using a domain decomposition procedure.Let N h be a nonoverlapping partition of Ä into rectangles Ä j of diameter bounded by h such that Ä = S J j =1 Ä j .Set 0 j = ∂Ä ∩ ∂Ä j , 0 jk = ∂Ä j ∩ ∂Ä k , and denote by ξ j and ξ jk the midpoints of 0 j and 0 jk , respectively.Let us denote by ν jk the unit outer normal on 0 jk from Ä j to Ä k and by ν j the unit outer normal to 0 j .Let χ j and χ jk be two unit tangents on 0 j and 0 jk so that {ν j , χ j } and {ν jk , χ jk } are orthonormal systems on 0 j and 0 jk , respectively.
Let the positive definite mass matrix P ∈ R 6×6 and the nonnegative dissipation matrix C ∈ R 6×6 be defined by where I denotes the identity matrix in R 2×2 .Also set F = (f s , f n , f w ) and for u = (u s , u n , u w ) let L(u) be the second order differential operator defined by Then we seek the solution of our differential problem over each subdomain Ä j as follows: for j = 1, . . ., J , find u j (x, z, ω) such that with ω * an upper temporal frequency of interest.Problem (3.1) needs a set of boundary conditions chosen as indicated below.Set If Ä j has a part 0 j of its boundary contained in ∂Ä, we impose the absorbing boundary condition (see [44]) where the symmetric positive definite matrix B is given by [44] Furthermore, as in [17,26], at the interior interface 0 jk we use the Robin transmission boundary conditions: Here β jk is a positive definite matrix function defined on the interior boundaries 0 jk .The Robin transmission conditions (3.3), (3.4) impose the continuity of the solid displacement, the normal component of the fluid displacements and the generalized stresses at the interior interfaces 0 jk .The spatial discretization is performed as follows.To approximate each component of the solid displacement vector we employ the nonconforming finite element space as in [19], while to approximate the fluid displacement vectors we choose the vector part of the Raviart-Thomas-Nedelec space [33,37] of zero order.More specifically, set with the degrees of freedom being the values at the midpoint of each edge of b R. Also, if For each Ä j , let F Ä j : b R → Ä j be an invertible affine mapping such that F Ä j ( b R) = Ä j , and define The choice of θ( x) as above implies that for an orthogonality property that was used in [19] to derive optimal apriori error estimates for a global nonconforming procedure for second order elliptic problems.A similar requirement is needed in our case.
Next, following [17,19,26], we introduce a set of Lagrange multipliers η jk = (η s,ν jk , η s,χ jk , −η n jk , −η w jk ) associated with the values of the generalized forces at the mid points ξ jk of 0 jk in the sense that η jk ∼ G jk (U j )(ξ jk ).The Lagrange multipliers η jk belong to the following space of functions defined on the interior interfaces 0 jk : where P 0 (0 jk ) denotes the constant functions on 0 jk .Next, we state a domain decomposition iteration using a variational formulation.Let us denote by (•, •) j the usual complex inner product in L 2 (Ä j ).Moreover, for 0 = 0 j or 0 = 0 jk let h•, •i 0 denote the complex inner product in L 2 (0), and let hhu, vii 0 denote its approximation by the midpoint quadrature: hhu, vii 0 = (uv)(ξ jk )|0|, where |0| is the measure of 0.
For t = 0, 1, 2, . . ., let U t j = (U s,t j , U n,t j , U w,t j ) and η t jk be the discrete displacement vectors and the Lagrange multipliers at the t-iteration level.Then, multiply (3.1) by v ∈ [V h j ] 2 ×W h j ×W h j , and integrate over Ä j , using integration by parts in the (L(u), v)term.Next, apply the boundary conditions (3.2) and (3.3) and approximate the boundary integrals on 0 j and 0 jk using the midpoint quadrature rule.Then the domain decomposition iteration can be stated as follows: given Equation (3.6), used to update the Lagrange multipliers, is obtained directly from (3.3) evaluated at the mid point ξ jk .Equation (3.5) yields a 16 × 16 linear system of equations for the degrees of freedom associated with the vector displacements of the three phases on each subdomain Ä j at the t-iteration level.After solving these systems, the Lagrange multipliers are updated using (3.6).The iteration (3.5), (3.6) is a Jacobi type iteration.A twice as fast iteration may also be defined by using a red-black type iteration (see [17,26]).Besides, our Jacobi iteration was combined with relaxation in order to improve the convergence.The arguments given in [17,26] can be used here to show that the iteration (3.5), (3.6) converges and it is first order correct in the spatial discretization.The choice of the nonconforming element used to compute the solid displacement is based on the dispersion analysis performed in [48] showing that it almost halves the number of points per wavelength needed to have a desired accuracy as compared with the standard conforming bilinear element.The iteration parameter matrix β jk is chosen to have the same form of the matrix B in (3.2) and defined in terms of averages of the coefficients of the differential system on both sides of the interfaces 0 jk .The space-time solution is obtained by solving (3.5), (3.6) for a finite number of frequencies and an approximate inverse Fourier transform [18].The definition of the iteration (3.5), (3.6) can be extended to the case of larger subdomains Ä j , see [42].

Numerical applications
We use the iterative procedure (3.5), (3.6) to simulate the propagation of waves in a sample of Nivelsteiner sandstone, a friable sandstone mainly composed of quartz with small percentages of rock fragments and potash-feldspar (see [30] for more details).Its material properties, taken from [1], are φ = 0.33, K = 5000 mD, ρ s = 2.65 g/cm 3 , grain bulk modulus K s = 36 GPa, frame bulk modulus K m = 6.21GPa and frame shear modulus, N = 4.55 GPa.The pore space is assumed to be filled by water (as the wetting phase) and a free hydrocarbon gas.Their properties are: ρ w = 1 g/cm 3 , µ w = 0.01 poise, K w = 2.223 GPa, ρ n = 0.1 g/cm 3 , µ n = 0.00015 poise, K n = 0.022 GPa.The reference fluid pressure P w is taken 30 MPa, corresponding to a typical hydrostatic pressure at a burial depth of about 3 km.
The relative permeability functions K rn (S n ) and K rw (S n ) and the capillary pressure function P ca (S n ) needed to describe our system are taken to be [13]: These relations are based on laboratory experiments performed on different porous rocks during imbibition and drainage processes (neglecting hysteresis effects).In the numerical experiment, we chose S rw = S rn = 0.05, and A = 30 kPa.The resulting capillary pressure at S n = 0.1 is about 3.4 kPa.In the absence of proper experimental data, the coupling permeability function K rnw (S n ) used in this work is assumed to be The parameter ² in (2.4) and the equation above is equal to 0.1, as in [40].The viscoelastic parameters describing the dissipative behaviour of the saturated sandstone are b Q l = 30, 20, for l = K c , N, respectively, T 1,l = 10 ms, T 1,l = 10 9 ms, for l = K c , N.
Values of the phase velocities c j and attenuation coefficients α j (in dB) for j = P 1 , P 2 , P 3 , S, computed as in [9,20], at the central frequency f 0 = 500 kHz are given in table 1 for S n = 0.1 and S n = 0.9.The value of the Kozeny-Carman constant A 0 used in the definition of the pore size parameters is equal to 5 [29].
In the following numerical experiments, we simulate a laboratory test of generation and propagation of body waves at ultrasonic frequencies.The domain for the numerical simulation is a square of side length L d = 6 cm with a uniform partition N h of Ä into squares of side length h = L d /N x , with N x = N z = 640.The source function (f s , f n , f w ) is a compressional point source located at (x s , z s ) = (3 cm, 2 cm) applied to the solid and fluid phases given by where δ x s ,z s denotes the Dirac distribution at (x s , z s ).The function g(ω) is the Fourier transform of the waveform g(t) = −2ξ(t − t 0 )e −ξ(t−t 0 ) 2 , with f 0 = 500 kHz denoting the source central (dominant) frequency and ξ = 8f 2 0 , t 0 = 1.25/f 0 .The spectrum of g(ω) is negligible for frequencies ω above ω * = 2π 1000 kHz.Thus, the iterative procedure (3.5), (3.6) is used to compute the particle velocities V l (x, z, ω) = iωU l (x, z, ω), l = s, n, w, at 500 equally spaced temporal frequencies in the interval (0, ω * ).The solution V l (x, z, t), l = s, n, w, at the desired discrete times is obtained by using a discrete time Fourier transform.The maximum simulation time in all the experiments was equal to 0.06 ms.See [18] for more details about the numerical procedure and error estimations.

Saturation interface
We first analyze the propagation of the wave fields in the presence of a plane horizontal interface defined by an abrupt variation in saturation.This interface divides the computational domain in two layers of equal thickness (3 cm).In the lower part of the domain (where the source is located) S n = 0.9 and in the upper part S n = 0.1.As expected, when the three compressional waves generated by the point source arrive at the interface, reflection, refraction and wave conversion phenomena can be observed.
In figure 1 we present a snapshot of V s z (x, z, t) (i.e., the normal component to the saturation interface) at t = 0.007 msec, showing the propagation of the fast P 1 wave and the development of a slow wave front still very close to the source.Figure 2 shows a snapshot of V s z (x, z, t) at t = 0.012 msec, where the effect of the saturation interface Figure 1.Saturation interface.Snapshot of V s z (x, z, t) at t = 0.007 msec, showing the fast compressional P 1 waves, direct (P 1d ), transmitted (P 1t ) and a slow direct wave (P 3d ).can be clearly observed.In this picture we can see the direct P 1 and P 3 waves, P 1 reflected and transmitted waves and shear S waves (reflected and transmitted), converted from the incident P 1 wave.No slow P 2 wave is observed due to its high attenuation coefficient, as can be seen in table 1 for S n = 0.9.The presence of these shear waves was also checked by computing the curl and divergence of the particle velocity field in the solid.The wave energy splitting at the saturation interface can be quantified by measuring for the solid phase the peak amplitudes of the different wave fronts incident (i), reflected (r) and transmitted (t) through the interface, which will be denoted as A m (r), m = i, r, t for r = P 1 , P 2 , P 3 , S. For this purpose we computed snapshots of the modulus of the two-dimensional particle velocity fields.The amplitude of the incident P 1 wave was measured just before arriving to the interface.For the reflected and transmitted P 1 waves we found the ratios A r (P 1 )/A i (P 1 ) ≈ 0.06 and A t (P 1 )/A i (P 1 ) ≈ 0.2, respectively.For the shear waves the corresponding ratios are A r (S)/A i (P 1 ) ≈ 0.09 and A t (S)/A i (P 1 ) ≈ 0.29.This shows the important wave energy partition that may take place due to saturation heterogeneities.The corresponding snapshots for V n z (x, z, t) and V w z (x, z, t) show a similar behaviour and are omitted.The wave field corresponding to V s z (x, z, t) at t = 0.045 msec is shown in figure 3, when the direct fast waves exited from the domain.The spherical wavefront corresponds to the direct slow P 3 wave arriving at the interface.Significant wave mode conversion from the incident P 3 wave to fast P 1 , S and slow reflected and transmitted waves can be observed.In particular, a P 2 wave is transmitted to the upper medium since for S n = 0.1 the attenuation coefficient shows a significant decrease; the corresponding transmitted P 3 wave is not observed due to its high attenuation for this saturation (see table 1).From the snapshots we obtained the following ratios for the reflected waves: A r (P 3 )/A i (P 3 ) ≈ 0.122, A r (P 1 )/A i (P 3 ) ≈ 0.1, A r (S)/A i (P 3 ) ≈ 0.044, and for the transmitted: A t (P 1 )/A i (P 3 ) ≈ 0.019, A t (P 2 )/A i (P 3 ) ≈ 0.055, A t (S)/A i (P 3 ) ≈ 0.006.Similar mode conversion phenomena in single phase Biot media were theoretically studied by some of the authors in [38].
Some numerical artifacts can be observed in the figures due to the low number of points per wavelengths used to represent the slow waves at high frequencies (about 3).No spurious reflections from the artificial boundaries are observed, showing that the absorbing boundary conditions is performing quite well in the simulation.

Immiscible patchy saturation
Here we abandon the hypothesis of homogeneous saturation and we consider a patchy distribution of the fluids (water and gas) throughout the domain.For this set of experiments porosity and permeability are taken constants, with the values given at the beginning of this section.We remark that, unlike other known models for patchy saturation, this one takes into account the capillary pressure and residual saturation effects that take place when the pore space is occupied by two (or more) immiscible fluids.
To generate the spatial distribution of the fluids we used a modified fractal field, based on the so-called von Karman self-similar correlation functions.These models are widely used in the statistical characterization of heterogeneities for different applications.Following [23,27], we consider a particular case for which the spectral density is given by: where k = p k 2 x + k 2 z is the radial wavenumber, a the correlation length, H is a selfsimilarity coefficient (0 < H < 1), S 0 is a normalization constant and E is the Euclidean dimension.Equation (4.3) corresponds to a fractal process of dimension D = E +1−H at scales smaller than a.For this application we take E = 2, a = 0.2 cm and D = 2.5.
The first step was to assign to each point of the mesh a random number using a two-dimensional random number generator with uniform distribution (white noise).This field was then Fourier transformed to the spatial wavenumber domain and its amplitude spectrum was then filtered using (4.3).The result was then transformed to the spatial domain, forced to have zero mean and normalized to a range appropriate to produce saturation fluctuations around a given value S * subject to the restriction S rn < S n < 1 − S rw .Finally, to construct the patches, we modified the saturation S n obtained at each point of the grid so that for cells having S n 6 S * we changed the saturation to a value close to the non-wetting residual (i.e., almost fully water saturated) and where S n > S * we assigned a saturation close to that corresponding to the wetting residual saturation.The result of this procedure is illustrated in figure 4. The black zones represent cells Figures 5(a) and (b) show snapshots of V s z (x, z, t) at t = 0.007 msec for the heterogeneous and the homogeneous saturation, respectively.As can be seen from figures 4 and 5(a), a remarkable correlation exists between the locations and size of the larger gas patches and the maxima and minima of the particle velocity fields, which can be associated with compressional and shear waves generated at the irregular saturation interfaces after the passage of the fast P 1 wave.The presence of shear waves was verified by computing the corresponding snapshot of the curl of the velocity field, which is not shown for brevity.On the right side of figure 5(a) a spherical wave front can be observed around the source, due to the propagation of the P 1 wave train through a zone of homogeneous saturation.For comparison, figure 5(b) shows the propagation of this wave mode at the same time under homogeneous saturation S w = 0.9.The peak amplitude relation between figures 5(b) and (a) is 1.372 which shows the effect of the presence of the gas bubbles on the amplitudes of the waves.
The next graphs show the corresponding snapshots after 0.02 msec, where we can observe the propagation of the slow wave P 3 through the heterogeneities of the domain (figure 6(a)) and its counterpart in the homogeneous case (figure 6(b)).Once again we can observe important wave-mode conversion from slow to fast compressional and shear waves at the irregular saturation interfaces, a fact that was also verified observing the snapshot of the curl for this time.The peak amplitude relation between figures 6(b) and (a) is 2.63; this amplitude increase is again associated to the presence of the gas bubbles.

Porosity-permeability interface
Here we consider the propagation of the waves when the computational domain has homogeneous saturation S n = 0.1, but there exist an abrupt change in both porosity and permeability.This change in the petrophysical properties of the rock defines a plane horizontal interface at the middle of the domain.For the lower half we keep the reference porosity and permeability values given at the beginning of this section, that is K = 5000 mD and φ = 0.33, while for the upper half we consider a much lower permeability, equal to 20 mD.We assume that the porosity and permeability values at both sides of the interface are related through the Kozeny-Carman equation in the form [3] Thus, solving the resulting cubic equation for the porosity, we found for the upper layer a porosity φ = 0.0654.Figure 7 shows a snapshot of the vertical component of the particle velocity field after 0.012 msec.To analyze the wave energy splitting at this interface once again we measured the ratios of peak amplitudes of the different wave   fronts in the corresponding snapshot of the modulus of the particle velocity fields.We also measured the peak amplitude of the incident P 1 wave just before arriving at the interface.For the reflected P 1 wave we found that the ratio A r (P 1 )/A i (P 1 ) ≈ 0.06, for the transmitted A t (P 1 )/A i (P 1 ) ≈ 0.09, and for the shear reflected and transmitted waves the ratios are A r (S)/A i (P 1 ) ≈ 0.12 and A t (S)/A i (P 1 ) ≈ 0.25, respectively.

Fractal porosity-permeability
Now we analyze the propagation of waves when the saturation is uniform (S n = 0.1) but the permeability behaves like a stochastic fractal.The scattering of slow waves in poroelastic media with random and periodic heterogeneities saturated by single-phase fluids has been analyzed in [28].For the present analysis we describe permeability in the form [12] K(x, z) = K * e f K (x,z) , (4.5) where K * is a constant reference permeability and f K (x, z) is a fractal field.In the numerical tests we consider different permeability distributions to analyze the resulting wave fields for variable degree of spatial heterogeneity.Following [12] this effect is quantified by means of the dimensionless coefficient C ν defined as the ratio between the standard deviation of K(x, z) and its mean value.
To generate the fractal fields f K (x, z) we used equation (4.3) with Euclidean dimension E = 2 and fractal dimension D = 2.2 for different values of the correlation length and the reference permeability.The numerical procedure is analogous to that described in the previous set of experiments.The field f K (x, z) is used in (4.5) in order to obtain a fractal permeability distribution and the corresponding porosity distribution is obtained by solving the Kozeny-Carman equation (4.4) at each point of the domain.More general relationships between porosity and permeability for fractal pore geometries are presented in [34].
In figure 8 we show the fractal permeability distribution obtained using a correlation length a = 0.2 cm, K * = 1100 mD and a C ν coefficient of about 0.58.The resulting permeability field varies from 138 to 7284 mD.The darker areas in the figure correspond to the higher values.The corresponding porosity map, with values in the range from 0.12 to 0.36 is not included for brevity.
Figure 9(a) shows V s z (x, z, t) after 0.009 msec, where the combined effect of the P 1 and S waves can be observed.It is important to observe the irregular character of the wave fronts due to the scattering effects at the heterogeneities and also to the spatial variations in the intrinsic attenuation.Next, we analyze the acoustic response of the medium for the same correlation length but changing K * to 500 mD and C ν ' 0.93 (i.e., a more heterogeneous medium).In this case the permeability varies from 25.35 mD to 7578 mD and the porosity from 0.07 to 0.36.The corresponding snapshot is presented in figure 9(b).To quantify the overall attenuation of the waves we computed the ratio between the maximum amplitudes in figures 9(b) and (a) resulting 0.95.Now we reduce one order in magnitude the scale of the heterogeneity, taking a correlation length a = 0.025 cm, but keeping the same parameters of figure 9(a) (i.e., C ν ' 0.58).Here the permeability varies from 135 mD to 10903 mD and the porosity from 0.12 to 0.39.The resulting snapshot is shown in figure 9(c), where it can be noted that the form of the wavefronts are less affected by scattering than in the two previous figures.This is due to the small size of the heterogeneities compared to the wavelength of the fast waves.However, the relative amplitude decrease in this case (referred to figure 9(a)) is about 11% (amplitude relation 0.88), showing the effect of the intrinsic attenuation.
In the next set of graphs we analyze the effects of the previously described fractal porosity-permeability fields on the propagation of the slow waves.Figures 10(a)-(c) show the corresponding snapshots of V s z (x, z, t) at t = 0.06 ms, where the relative amplitude decrease (normalized to figure 10(a)) is about 60% for figure 10(b) and 59% for figure 10(c).These higher dissipation rates are due to the fact that the slow waves have amplitude components only in the high frequency range, with wavelengths on the order of the heterogeneities, consequently suffering important scattering effects.The  faster waves have frequency components also in the low frequency range (with larger wavelengths) being less affected by scattering at the heterogeneities.

Conclusions
We have developed a model to describe the acoustic response of porous dissipative rocks saturated with two-phase immiscible fluids.The constitutive relations and the equations of motion include the effects of capillary forces, mass and viscous coupling coefficients and viscoelastic dissipation.Frequency dependent correction factors for the high-frequency range are taken into account, so that the formulation is valid for the full frequency range.
This model is applied for the numerical simulation of ultrasonic waves in a real clean sandstone, saturated with gas and water, with different kinds of heterogeneities.The numerical solution is obtained by means of an iterative domain decomposition finite element procedure formulated in the space-frequency domain, allowing for the simultaneous and independent solution of the equations for a finite number of frequencies.The space-time solution is obtained by using an approximate inverse Fourier transform.
We analyzed the response of the medium in the presence of plane horizontal interfaces defined by changes either in saturation or porosity-permeability and also in the case of more realistic distributions in these parameters given by stochastic fractals.Our main results can be summarized as follows: • At the plane interfaces, important energy splitting is observed due to mode conversion from fast compressional and shear waves to slow waves and vice versa.• In the case of patchy saturation, the presence of gas pockets produces local increases in the amplitudes of the wavefronts, both for fast and slow waves.• For the fast waves there exist a noticeable correlation between the location and size of the maximum amplitudes of the wavefronts in the snapshots and the location of the larger gas pockets in the patchy saturation map.• Important mode conversion from slow to fast waves occurs at the gas pockets.• In the case of fractal porosity-permebility, there exist attenuation and dispersion effects in all the wavefronts, associated to scattering at the irregular interfaces and also to variable intrinsic attenuation.However, the energy losses are more pronounced in the case of slow waves, because their wavelengths are in the order of the size of the heterogeneities.• When the correlation length becomes smaller, the wavefronts associated with fast waves become more regular, as in the presence of an equivalent effective medium.
All these results indicate that a careful analysis of the characteristics of ultrasonic waves in porous saturated reservoir rocks may provide useful information about the spatial variability of its petrophysical parameters at different scales, about the nature of the reservoir fluids and to discriminate the saturation state.

Figure 2 .
Figure 2. Saturation interface.Snapshot of V sz (x, z, t) at t = 0.012 msec, displaying direct, transmitted and reflected fast compressional waves (P 1d , P 1t , P 1r ), reflected and transmitted shear waves (S r , S t ) and the direct slow wave P 3d .

Figure 3 .
Figure 3. Saturation interface.Snapshot of V s z (x, z, t) at t = 0.045 msec, showing the direct and reflected slow waves (P 3d , P 3r ), a transmitted P 2t wave and transmitted and reflected fast compressional and shear waves (P 1t , P 1r , S r , S t ).

Figure 4 .
Figure 4. Patchy saturation distribution.Black zones represent regions of almost full water saturation.The macroscopic gas saturation is 10%.

Figure 5 .
Figure 5. Patchy saturation.Snapshots of V s z (x, z, t) at t = 0.007 msec for (a) patchy and (b) homogeneous 10% gas saturation.The peak amplitude relation between snapshots (b) and (a) is 1.372.

Figure 6 .
Figure 6.Patchy saturation.Snapshots of V s z (x, z, t) at t = 0.02 msec for (a) patchy and (b) homogeneous 10% gas saturation.The peak amplitude relation between snapshots (b) and (a) is 2.63.

Figure 7 .
Figure 7. Porosity-permeability interface.Snapshot of V s z (x, z, t) at t = 0.012 msec, displaying direct, transmitted and reflected fast compressional waves (P 1d , P 1t , P 1r ), reflected and transmitted shear waves (S r , S t ) and the direct slow wave P 3d .

Table 1
Phase velocities and attenuation factors at 500 kHz.