Numerical simulation of two-phase fluid flow

We simulate two-phase fluid flow using a stress–strain relation based on Biot’s theory of poroelasticity for partial saturation combined with the mass conservation equations. To uncouple flow and elastic strain, we use a correction to the stiffness of the medium under conditions of uniaxial strain. The pressure and saturation differential equations are then solved with an explicit time stepping scheme and the Fourier pseudospectral method to compute the spatial derivatives. We assume an initial pressure state and at each time step compute the wetting- and non wetting-fluid pressures at a given saturation. Then, we solve Richards’s equation for the non wetting-fluid saturation and proceed to the next time step with the updated saturations values. The pressure and saturation equations are first solved separately and the results compared to known analytical solutions showing the accuracy of the algorithm. Then, the coupled system is solved. In all the cases, the non-wetting fluid is injected at a given point in space as a boundary condition and capillarity effects are taken into account. The examples consider oil injection in a water-saturated porous medium.


Introduction
Diffusion equations can be obtained in poroelasticity at low frequencies, when the inertial terms are neglected (e.g. Carcione 2007). Chandler and Johnson (1981) have shown the equivalence between quasi-static fluid flow and Biot's diffusive wave (see also Carcione 2007). Hence, fluid flow and pressure diffusion are phenomena described by the same differential equation. In hydrology and hydrocarbon exploration, diffusion equations are mainly used to model fluid flow in reservoir rocks (Peaceman 1977). Shapiro et al. (2002) describe the phenomenon of micro-seismicity caused by fluid injection in boreholes, while Müller (2006) provides a detailed physical analysis of the pore pressure induced by a fluid mass point source. His results support the hypothesis that the diffusive slow P-wave is mainly responsible for the triggering of microearthquakes. Carcione and Gei (2009) simulated fluid-pressure diffusion in inhomogeneous media by considering a single fluid. The purpose of the present work is to generalise that approach to more than one fluid and use a similar method to solve the equations governing the flow. Possible applications involve simulation of fluid depletion and injection of hydrocarbons, geomechanical analysis of reservoirs (Settari and Mourits 1994) and micro-seismicity, as mentioned above. A convenient equation to define the stress-strain relation involving two fluids can be obtained from Biot's theory of poroelasticity (Santos et al. 1990a, b;Ravazzoli et al. 2003). The diffusion can be fully uncoupled from the elastic deformations by neglecting the strain, but a less stringent condition can be obtained under uniaxial strain conditions, where the effect of the elastic deformations requires the modification of the rock stiffness involved in the diffusion equation. Uniaxial-strain conditions hold for the case of a fluid-injection source in a borehole, where the porous medium deforms uniaxially only in the vertical direction and the horizontal strains are zero. This approximation is valid for geological formations whose lateral dimensions are large compared to their thickness. In this case, the lateral deformations are small compared to the horizontal dimensions, and hence the lateral strain is approximately zero (e.g., Gutierrez and Lewis 2002).
In this work, we combine Biot's poroelastic equations and the fluid mass conservation equations averaged to obtain a system of differential equations for the pore pressures and saturations of the wetting and non-wetting phases. The theory used here to define the stress-strain relations is a generalization of Biot's theory and has been developed by Santos et al. (1990a, b) (see also Carcione et al. 2004). It is important to note here that the constitutive equations introduced by Santos et al. (1990a, b) have never been used to simulate fluid-flow phenomena, but only to describe wave propagation. Biot's wave propagation theory, neglecting the acceleration terms, is equivalent to the consolidation theory if the infinitesimal stresses and pressures are assumed to be their absolute counterparts, since the stress-strain relations can be interpreted as a relation between incremental fields, where stress and strain are increments with respect to a reference stress and strainthe case of wave propagation-or, as relations between the absolute fields (consolidation theory; Biot 1941). Other models accounting for the coupling between the mechanical behaviour of the matrix and the fluid dynamics include the Hassanizadeh model (Hassanizadeh et al. 2002) and the Barenblatt-Biot model (Barenblatt and Gilman 1987).
The poroelastic stress-strain relation for the two fluids can be recasted as a first-order differential equation in time and allows us to obtain the fluid pressures and total flow velocity. On the other hand, the mass conservation equation can be recasted to obtain the advection-diffusion Richards equation which is coupled to the fluid pressures through the flow velocity. The diffusion term is determined by the capillarity effects. We first solve the pressure and saturation equation separately and test the results with known analytical solutions for the advection and diffusion equation in 1D and 2D space. Then, we investigate the effects of partial saturation and capillary forces on the pressure and saturation profiles and finally solve the coupled system, which consists in computing the flow velocity and then Richards's equation for the fluid saturations. The algorithm is fully explicit and based on the Euler-Picard time stepping and on the Fourier pseudospectral method to compute the spatial derivatives (e.g., Carcione 2007). It has been shown that the Picard iteration method demonstrated perfect mass balance in fluid-flow equations when used with direct-grid methods such as the finite difference and finite element approximations in space (Zarba 1988). It is expected that the same performance can be achieved with the Fourier pseudospectral method. In this case, the model is discretized on a mesh the spatial derivatives are calculated with the fast Fourier transform. The spatial derivatives computed with the Fourier method have spectral accuracy and require coarser grids compared to finite-difference methods. The use of this spectral method overcomes two drawbacks: low accuracy and stringent stability conditions, since the error in time decays exponentially. The method has proved to be efficient when solving diffusion equations. (Carcione 2006(Carcione , 2010Carcione and Gei 2009).

Darcy's equations
We consider a porous medium (a rock for instance) saturated with two immiscible fluids and denote with the subscripts (and superscripts) w and n quantities related to the wetting and non-wetting phases, respectively. Let v s , v n ; and v w denote the interstitial particle-velocity vectors of the solid grains, non-wetting fluid and wetting fluid, respectively. First of all, we assume that the rock deforms slowly compared to multiphase flow. In this case, we have the relative particlevelocity components   Santos et al. 1990a), where S, q, g and j denote saturation, density, viscosity and permeability, respectively, g is the gravity constant, i indicates the spatial variables x(1), y(2) and z(3), d ij is Kronecker's delta and the subindex ''i'' denotes spatial differentiation. Darcy's velocities, as given in Gai (2004, Eq. 2.16) where N a are the so-called (dimensionless) concentrations.

Stress-strain relations
Following Santos et al. (1990a, b), we define where n is the variation of fluid content and a dot above a variable denotes time differentiation. Let ij denote the strain components of the solid grains let # ¼ 11 þ 22 þ 33 denote the dilatation field. Then, the total stress components and fluid stresses are where N is the dry-rock shear modulus and K c is the undrained (closed) bulk modulus. It is where K s is the bulk modulus of the solid grains, K m is the dry-rock bulk modulus, where K s is the bulk modulus of the grains, K n and K w are the bulk moduli of the non-wetting and wetting fluids, respectively, p c is the capillary pressure (i.e., the difference between the non-wet and wetting-fluid absolute pressures), p w is an average hydrostatic pressure (assumed to be zero in this work) and p 0 c the derivative of p c with respect to S n ; moreover, Uncoupling fluid flow and deformation The fluid pressure is coupled with the strain of the matrix in Eq.
(3). This fact makes the problem much more difficult to solve, but there are situations where these field variables can be uncoupled. They occur when the displacement field is irrotational or when the fluid is very compressible (e.g., Detournay and Cheng 1993). We may avoid such approximation by using a less stringent one. Let us assume the case of a fluid-injection source in a borehole, uniaxial strain conditions and vertical deformations only. In this case, the only non-zero differential strain is d# ¼ d 33 : Assuming no changes in the vertical stress, we obtain from Eq. (3): Using this result, the second and third Eq.
(3) become More general approaches involving the coupling of fluid flow and deformation based on Biot's equations are given in Gutierrez and Lewis (2002).

Mass conservation equations
The system of equations are complemented with the balance of mass (Gai 2004, Eqs. 2.17-2.19;Peaceman 1977;Santos et al. 1990a), where v m = (v 1 m , v 2 m , v 3 m ) T , s m is the source or sink term and we have assumed constant volume factors.

The pressure and Richards equations
Replacing the flow velocities from Eq. (1) into the time derivative of Eq. (18) and into Eq. (19), we obtain where g is the gravity vector. These equations are subject to the constraints From Eq. (18) we have and the first constraint Eq. (21) becomes In particular, we set From Eq. (20), we obtain the system to be solved: where the third (Richards's) equation is obtained in ''Appendix A'' c t is the total Darcy velocity. This system is complemented with Eqs. (18) and (22), subject to the initial condition (p w0 , S w0 ). In the limit S w ! 1 (assuming the residual saturations S rw = S rn = 0), a 2 ! ÀM; where M is given in Eq. (10) of Carcione and Gei (2009) and the first Eq. (25) becomes Eq. (9) of that paper. The algorithm to solve these equations is explicit in the time variable, t = j dt, where dt is the time step and j ¼ 0; 1; . . .: It uses the forward Euler method to approximate the time derivatives and the Fourier pseudospectral method to compute the spatial derivatives (e.g., Carcione 2006Carcione , 2007Carcione and Gei 2009). We assume an initial state (p w0 , S w0 ), compute s w at j ? 1 (first equation) assuming the saturations obtained at j on the r.h.s. of the equations. Then, we compute s n , the pressures with Eq. (22) and solve the last two equation for S n at j ? 1 proceed to the next time step. To improve the Euler method we solve the saturation Eq. (25) by using the Picard method (Pang 2006) (see ''Appendix B'').

Fourier analysis
The characteristics of the pressure diffusion can be analysed with a Fourier analysis. Let us assume a kernel of the form p w ¼ exp½iðxt À k Á xÞ; where x is the angular frequency, k is the complex wavenumber vector and x is the position vector. Assuming homogeneous properties, we have k = k (l 1 , l 2 , l 3 ) T , where k ¼ ReðkÞ À i a; a is the attenuation factor l i are the direction cosines defining the propagation direction, where Re takes the real part. From Eqs. (18) and (22) we have s w ¼ ÀS w p w þ fp c and _ s w ¼ ÀS w _ p w at constant saturation. Substituting the above kernel into the first Eq. (25), in the absence of source, gives the dispersion equation We define the complex velocity as (the quantity inside the square root is positive). The same kinematic concepts used in wave propagation (acoustics and electromagnetism) are useful in this analysis (see Carcione 2007, Chap. 8). For S w ! 1; we obtain the complex velocity of Carcione and Gei (2009). The phase velocity of the pressure front and respective attenuation factor can be obtained from the complex velocity as respectively, where Im takes imaginary part. The skin depth is the distance d for which expðÀ adÞ = 1/e, where e is Euler's number. It is usually taken as the effective distance of penetration of the signal. Then, d ¼ 1= a:

Numerical simulations
Tests are performed in the cases where analytical solutions exist. First, we consider the first Eq. (25) in the homogeneous case, without gravity, S w = 1 and p w0 = 0. In this case, s w =p w , a 3 = 0 and a 2 b w = -Mj /g w , where M is given in Eq. (10) of Carcione and Gei (2009). We obtain the equation where D is the Laplacian, j is the absolute permeability and here d is the Dirac delta function, representing the source. The analytical solution is the Green's function given by Eq. (35) in Carcione and Gei (2009). We consider the 1D case, the material properties shown in Table 1, a number of grid points n x = 165 and a grid spacing dx = 4 m/ n x . Figure 1 compares numerical and analytical time histories (normalized) at 24 cm from the source computed with a time step dt = 10 ls. The agreement is excellent.
In the absence of sources and gravity effects and D ¼ h n =/ (constant), an analytical solution is known for the 1D version of the third Eq. (25), i.e., where D is a diffusivity, for the initial condition S n (x,0) = 0 and boundary conditions S n (x 0 , t) = 1 and S n ð1; tÞ ¼ 0: It is given by where erfc is the complementary error function (Peaceman 1977, Eq. 4.35). In 1D space c t is a constant since qc t / qx = 0. The velocity at which the saturation front travels is Equation (30)  x C x 0 , where S n (x 0 = 0 , t) = 1. An important parameter is the Courant number C, giving the distance traveled due to advection in a time step dt with respect to the grid spacing, i.e., C ¼ cf 0 n dt=dx: Being explicit, the algorithm is stable for C \ 1, which means that one cannot advect the saturation more than one grid cell in a single time step. Accuracy requires a time step smaller than that imposed by stability. Figure 2b shows the comparison between the numerical and analytical solutions at 1 s, where n x = 429, dx = 4 m/n x and dt = 0.2 ms. The analytical solution is shown for x [ x 0 . In 2D space, assuming diffusivity anisotropy ðD x ; D z Þ a flow velocity (v x , v z ) and a point-source injection, Eq. (30) can be generalized as where s 0 is the rate of injection per unit area and the solution to this equation yields the Green function. It is well known that this equation can be recast as a pure diffusion equation. Let us use an analogy and consider the x-direction. If we are traveling in a train moving at velocity f 0 n v x ; we do not see the advection. If the train starts to move at t = 0 and x = 0, its location is x ¼ f 0 n v x t and we may define a moving coordinate x 0 ¼ x À f 0 n v x t; so that x 0 ¼ 0 always. With respect to those moving coordinates, we have the following diffusion equation: whose solution, at the original coordinates, is (Carslaw and Jaeger 1984;Carcione 2007). If the injection occurs at (x 0 , z 0 ), the solution is obtained by simply replacing t by tt 0 , x by xx 0 and z by zz 0 . In 1D space the solution is S n ðxÞ ¼ ðs 0 = ffiffiffiffiffiffiffiffiffiffiffiffi 4ptD x p Þ exp Àðx À f 0 n v x tÞ 2 =ð4D x tÞ h i : We consider n x = n z = 429, dx = dz = 4 m/n x and dt = 0.2 ms. Recall that the discrete spatial delta is represented by 1/dx in 1D space and 1/(dx dz) in 2D space.
with v 0 = 1.2 m/s at all the grid points.
Next, we consider the boundary condition S n (x 0 , t) = 1, a radial velocity flow as shown in Fig. 4a and (v ¼ v 0r ), with v 0 =1.2, f 0 n = 1.2 and the same diffusivity of the previous examples. Figure 4b shows the snapshot of the saturation, where we have taken n x = n z = 1,287, dx = dz = 4 m/n x and dt = 0.2 ms. Now we consider the mobilities and the capillary pressure effects. The permeabilities are given by (van Genuchten 1980;Lo et al. 2007), where j is the absolute permeability, 0 \ p \ 1 S rn and S rw are residual saturations. These relations are based on laboratory experiments performed on various porous rocks during imbibition and drainage processes (neglecting hysteresis effects).
The capillary pressure as a function of saturation is the van Genuchten function where P 0 is a reference pressure (e.g., Goumiri et al. 2011) and S rwv = 0.1 S rw . Here, p 0 c ¼ dp c =dS n : Relative permeabilities and capillary pressure as a function of the wettingfluid saturation are shown in Fig. 5, where j = 0.8 D, S rw = 0.25, S rn = 0.05, P 0 = 10 kPa and p = 0.75. Figure 6 display the advective and diffusive terms f 0 n and h n appearing in the Richards Eq. (25).
We now solve the first Eq. (25) for partial saturation with an initial condition in the wetting-fluid pressure field, without gravity effects and considering the capillary pressure. The equations are _ s w ¼ a 3 r Á b n rp n þ a 2 r Á b w rp w ; s n ¼ 1 S w ½ð1 À S w Þs w À Qp c ; p w ¼ ½ðS w þ fÞs n À ðS n þ b þ fÞs w =Q p n ¼ ½fs n À ðb þ fÞs w =Q: The initial condition is where p 0 = 1 MPa, p 1 = 0.1 MPa k 0 = 3/m. We consider the properties of Table 1, S rw = 0.05, S rn = 0.05 and p = 0.5. We first test the case where there is an analytical solution, i.e. only the presence of the wetting phase (S w = 1) and an initial condition f ðxÞ ¼ p 1 exp½Àk 2 0 ðx À x 0 Þ 2 : The 1D Green's function is well known (Carslaw and Jaeger 1984;Carcione and Gei 2009) and is given by GðxÞ ¼ ð4pdt 0 Þ À1=2 exp½Àðx À x 0 Þ 2 =ð4dt 0 Þ; where d = jM/g w and t 0 is a given time. The numerical solution fits the analytical solution (not shown here), which is given by p w ðxÞ ¼ G Ã f ; where ''Ã'' denotes spatial convolution. The next simulation considers three saturations with and without the presence of capillary pressure. We take n x = n z = 165, dx = dz = 4 m/n x and dt = 50 ls. The results at 60 ms are shown in Fig. 7, where the reference initial condition is also shown. Increasing wetting-fluid saturation implies attenuation and dispersion with higher levels at higher wetting-fluid saturations. The presence of capillary pressure produces a shift of the curves downwards. Figure 8 shows the non-wetting fluid pressure, which is not affected by capillary effects. Now, we compute the numerical solution of Richards's Eq. (25) at different saturations, with and without capillary effects. We ignore sources and gravity. We have where We then apply the boundary condition S n (x 0 , t) = 1 -S rw . We consider a 1D medium with c =c 0 for x B x 0 and c = ?c 0 for x C x 0 , where c 0 = 0.4 m/s. The grid size is n x = 315, with dx = 4 m/n x . Figure 9 shows a snapshot of the saturation at 15 s as a function of distance, where we have used P 0 = 10 kPa, p = 0.75, a time step dt = 0.05 ms and K = 10 iterations of the Picard algorithm (see ''Appendix B''). The capillarity pressure has little effect on the shape of the saturation front. In order to appreciate the capillarity effects, we consider P 0 = 1 GPa, p = 0.5, c 0 = 0.4 m/s, S n (x 0 , t) = 1 -S rw -0.01, the previous mesh and a time step dt = 0.1 ms.  Finally, we solve the coupled system Eq. (25) in 1D space, where the flow velocity (last equation) is obtained by solving for the pressures (first and second equations). We set the initial condition Eq. (39), with p 0 = 1 MPa, p 1 = 0.1 MPa, and k 0 = 3/m, S rw = 0.25, S rn = 0.05 and 95 % water saturation at t = 0. Then, we apply the boundary condition S n (x 0 , t) = 0.95 -S rw , with P 0 = 1 GPa and p = 0.5. We take n x = n z = 165, dx = dz = 4 m/n x , K = 10 and dt = 10 ls. Figure 11 shows the saturations at various propagation times. Diffusion effects due to the capillarity forces dominate in this case.

Conclusions
We have developed a fully explicit algorithm to simulate the injection of a non-wetting fluid into a porous medium saturated with a wetting fluid. The fluid pressures are obtained from the stress-strain relation of the theory of poroelasticity, from which the total (Darcy) flow velocity is computed. This velocity determines the advection term of Richards's saturation equation. Then, the coupled pressure-saturation system of equations is solved with the Euler-Picard time stepping algorithm and the Fourier pseudospectral method. The example assumes oil injection in a water-saturated porous medium. Analytical solutions allow us to verify the accuracy of the results in various situations, where the two unknown variables (pressure and saturation) are uncoupled, with and without capillarity effects based on the van Genuchten functions. Strong capillary forces have to be assumed to see the effects as shown in the simulations of the saturation front. Here, we have considered two fluids. The extension of the theory to  include gases, which requires to include pressure-volumetemperature and chemical effects, will be developed in a future work.