A Preconditioned Conjugate Gradient Based Algorithm for Coupling Geomechanical-Reservoir Simulations

In this article, we introduce a new coupling algorithm between the reservoir multiphase Darcy flow simulator and the geomechanical code accounting for the compaction of the porous medium. The coupling is defined on time periods in such a way that the reservoir unknowns are computed for time steps which are small enough subdivisions of the time period whereas the mechanical problem is solved at the end of the period. Our new approach is based on a nonlinear preconditioned conjugate gradient method which is applied to the mechanical displacement variable. This algorithm is compared, on a one dimensional example, to the staggered coupling algorithm with a porous medium compressibility as relaxation parameter. The main conclusion is that the preconditioned conjugate gradient method is much more robust and converge much faster than the staggered algorithm, while the additional cost per iteration should remain in practical situations very small.


INTRODUCTION
The production of hydrocarbon reservoirs involves pressure and temperature variations over long periods of time.These variations trigger off a modification of the reservoir geomechanical behavior during the production.A wellknown example of the stress equilibrium changes in and around the reservoir is the subsidence phenomenon that has been observed on different fields.The most famous case is the Ekofisk field in the North Sea, Norway, where a sea floor subsidence rate of 42 cm/year has been reached at the end of 1993 (Sylte et al., 1999).The cases of the Valhall field in Norway (Pattillo et al., 1998) and the Bachaquero (Merle et al., 1976) and Tia Juana (Mc Lendon and Sawyer, 1991) fields in Venezuela also illustrate the importance of the subsidence phenomenon in oil reservoir.
The subsidence phenomenon strongly depends on the reservoir compaction.In conventional reservoir simulators, the reservoir compaction is directly deduced from the reservoir pore pressure and temperature changes using porous medium compressibility with pressure and porous medium dilatability with temperature coefficients.However, several authors (Tortike and Farouq, 1993;Settari and Mourits, 1994;Ruisten et al., 1996;Gutierrez and Lewis, 1998) demonstrate that porous medium compressibility and dilatability parameters (even pressure-dependent) are not sufficient to reproduce pore volume changes induced by pressure and temperature changes.Therefore an accurate modeling of stress dependent reservoir compaction and other related geomechanical effects require the coupling of a geomechanical model accounting for the rock deformation with a reservoir model describing the fluid flow in the porous medium.The porosity change due to strain variation and the total stress change associated with the pressure variation are the reservoir-geomechanical couplings considered in this paper.
There are two different approaches for solving the geomechanical-reservoir problem.The first one consists in simultaneously solving all the equations in one simulator, which is referred to as the fully coupled algorithm (Gutierrez and Lewis, 1998;Chin et al., 1998), whereas the second one uses two different simulators in order to solve the two sets of equations (mechanical equilibrium and flow problem); each simulator solves its own system independently, and information is passed in both directions between the simulators.This second technique, which is referred to as the simulators coupling method, can be performed with two conventional reservoir and geomechanical simulators (Tortike and Farouq, 1993;Settari and Mourits, 1998;Settari and Walters, 1999).
The objective of this paper is to define and compare two algorithms coupling the reservoir simulations performed over time periods to the mechanical computations at the end of each period.The reservoir simulation over each period is computed for time steps which are small enough subdivisions of the time period.
The first coupling algorithm defined in Section 3 is the so called staggered algorithm.It is a fixed point method which performs a sequence of reservoir simulations with the porosity field fixed by the previous iteration, followed by an update of the mechanical displacement and porosity field.In the single phase linear case, this algorithm amounts to a Gauss-Seidel type iterative method.
In Section 4, a new coupling algorithm is introduced which is based on a preconditioned conjugate gradient approach to solve the coupled system.Each iteration of the conjugate gradient method involves one mechanical computation which plays the role of the preconditioning step, and two sets of reservoir computations on the whole time period to compute the conjugate gradient residual.
In Section 5, these algorithms are compared in terms of robustness and convergence rate on a one dimensional example.
For the sake of simplicity, throughout this article, the reservoir is assumed to be an isothermal dead-oil model whereas the geomechanics is considered as a linear poroelastic model with an uncompressible rock matrix.Both models are introduced in Section 1, and their discretizations, by a finite volume scheme for the reservoir equations and by a finite element method for the mechanical equilibrium, are briefly outlined in Section 2.

COUPLING GEOMECHANICAL AND RESERVOIR SIMULATIONS
The reservoir is modeled as a porous nonlinear elastic medium with an instantaneous deformation response coupled with an unstationary compressible two phase Darcy flow.Assuming small deformations of the medium, both the fluid flow and the geomechanical equilibrium equations are written in Euler coordinates, and the computational geometrical domain is a fixed domain Ω.The coupling between the geomechanics and the Darcy flow is given by a simplified form of Biot's law connecting the variation of the porosity of the medium to its mechanical deformation, which amounts to assume that the solid matrix is incompressible.Conversely, we suppose that the medium is submitted to the fluid pore pressure.For the sake of simplicity we do not take into account the gravity forces.In addition, the reservoir is supposed to have zero permeability outside the subdomain ω ⊂ Ω.In other words, the subdomain ω represents the flow computational domain of the reservoir code, while the whole domain Ω is the geometrical domain of the geomechanical code (Fig. 1).Note that Ω is much larger than ω since the mechanical effects can be observed on a much larger scale than the scale of the reservoir domain; in particular the boundary conditions for the mechanical unknowns must be set sufficiently far away from ω.The coupled system consists in the following set of equations governing the time evolution of the fluid pore pressure p and the displacement field u of the porous medium.
The geomechanical equations are given by: (1) where σ and ε are the stress and the strain tensors for the porous medium with: λ and µ and are the Lamé constants in drained conditions, and the function b (x) is defined by: The reservoir equations state the conservation of each phase with fluxes defined by the two phase Darcy's law: (2) where we have neglected the capillary effects so that the water and oil pressures are equal.
The simplified form of Biot's law corresponding to the case of an incompressible rock matrix is given by: (3) In Equations ( 2) and (3) the following notations have been used: S denotes the saturation of the water phase, and 1 -S is the saturation of the oil phase; p is the pressure; is the porosity of the porous medium; κ is the intrinsic permeability; kr i,j = kr i,j (S) is the relative permeability of the phase i in presence of the phase j ≠ i; µ j is the dynamical viscosity of the phase i; ρ i = ρ i (p) is the density of the phase i.

THE DISCRETE PROBLEM
In order to approximate the fluxes and ensure the mass conservation, the reservoir set of equations is discretized using a finite volume method.As in most mechanical codes, the geomechanical problem is discretized using a finite element method.The displacement u is computed at times T 0 = 0, T 1 , ..., T k , ..., T p = T f .Also we suppose that T k+1 -T k = T, k = 0, 1, ..., p -1 and refer to T as the period.The unknown functions p, S, and φ are computed at times t k 0 , ..., t k n , ..., t k q = t 0 k+1 , where In fact, ∆t and T vary in the course of the computations; for the sake of simplicity, we shall consider them as fixed in Sections 2, 3, 4.
Furthermore, for the sake of conciseness, whenever it is not ambiguous, we shall use the notations p n K , S n K and φ n K rather than p K n,k , S K n,k and φ K n,k for the approximations of p, S, and φ at time t k n on the control volume K. Similarly let ρ n i, K , i = o or w, denote the cell values of the oil and water phase densities at time t k n .
The reservoir and mechanical geometrical domains.

Finite Volume Scheme for the Flow Problem
Let τ h be the set of control volumes, m (K) the measure of the cell K, and N (K) the set of its neighboring cells.The transmissibility between two neighboring cells K and L is defined by: with m (e K, L ) denoting the measure of the common face e K, L between the cells K and L, XK being a point associated to each cell K, and d (XK, XL) denoting the distance between the points XK and XL.The reservoir conservation Equations ( 2) are integrated on each space time element with a semi-implicit integration in time to obtain the following discrete system: with the "phase by phase upwinding" defined by the relations: The porosities , depend on the mechanical displacement through the Biot's law (3) and are defined in Subsection 2.3 below.

Finite Element Computation of the Geomechanical Equilibrium
Let us define the approximate pressure function by the following piecewise constant interpolation: The geomechanical equation is discretized by a finite element Galerkin approximation of Equation ( 1) with a continuous piecewise quadratic finite element space on a given mesh of the domain Ω.
Let u h, T (., T k ) denote this finite element approximation of u (., T k ), then it is defined by the following Galerkin variational approximation: (6) for all finite element test functionv h , where 〈f h ,.〉 is a right hand side accounting for the load at the boundary of the domain Ω.

Finite Volume Discretization of Biot's Law
The relation between the discrete displacement and the discrete porosity φ Κ n+1 is given by: ( (8) For c r = 0, ∆φ k K , denotes the variation of porosity during the period [T k , T k+1 ] which is linearly distributed over the subdivision t k n , n = 1, ..., q.Furthermore, the Definitions (7) and (8) of φ Κ n+1 and ∆φ k K include the usual porous medium compressibility correction c r of reservoir simulations.For q = 1, this correction cancels out between both Equations ( 7) and ( 8) and only plays the role of a relaxation parameter which stabilizes and speeds up the convergence of iterative coupling schemes (see Section 3).For q > 1, it is expected in addition to provide a better interpolation of the porosity within the period.

STAGGERED COUPLING ALGORITHM
The staggered algorithm performs iteratively until convergence a sequence of one reservoir simulation over the period T = T k+1 -T k on the subdivision t k 1 , ..., t k q with porosity variations ∆φ k K , K ∈τ h , given by the previous iteration, followed by an update of the geomechanical displacement and the porosity variations ∆φ k K , K ∈τ h , using the computed pressure at time t k q = T k+1 (Fig. 2).
Sketch of the staggered coupling iterative algorithm between the reservoir simulation over one period and the mechanical computation.
Let us denote by: the vector of the pressure, saturation, and porosity unknowns at time t k n , by u k the vector of the displacement node values at time T k , and by the vector of the variation of porosity unknowns.The subscript l will be used to denote the staggered algorithm iterations.
Initialization of the period: .Iterations on the period Do l = 1, … until convergence Do n = 1, … , q Reservoir simulation over the period: p l n,k and S l n,k are computed from (4) and (5) using End Do Geomecanical simulation: u 1 k+1 computed from (7) using p l q,k ∆φ k l computed from (8) End Do End of the period.The convergence criterion is based on the maximum relative variation of the pressure between two successive iterations.In the case of single phase linear flow and linear elasticity, this algorithm can be seen as a Block Gauss-Seidel iterative method for solving the pressure displacement coupled system.In that case a convergence analysis can be carried out and the convergence rate is shown to depend strongly on the compressibility ratio between the fluid and the rock.Furthermore, the porous medium compressibility c r can be interpreted as a relaxation parameter which is required to ensure the convergence of the algorithm when the compressibility of the porous medium is greater than the compressibility of the fluid (Bévillon and Masson, 2000;Bévillon, 2000).
It is well known that Gauss-Seidel type algorithms converge slowly and that their convergence rate can be considerably improved by using acceleration techniques like conjugate gradient methods.The definition of such an algorithm that can apply to nonlinear fluid flow models as well as to the subdivision of the geomechanical periods (q > 1) is the purpose of the next section.

PRECONDITIONED CONJUGATE GRADIENT METHOD
We consider throughout this section that the porous medium compressibility parameter c r is set to zero.The geomechanical Equation ( 6) are rewritten in matrix form as: (9) where G h denotes the stiffness matrix, B h the discrete gradient matrix, and F h the vector of the right hand side node values.With these notations, Equations ( 7) and ( 8) can be rewritten as follows: (10) where t B h denotes the transpose matrix of B h .From Equations (4), ( 5), (10), it is clear that the reservoir simulation over the period [T k , T k+1 ] relates the pressure p q,k to t B h u k+1 and to other quantities depending only on the previous period (like p q,k-1 , S q,k-1 , φ q,k-1 , and u k ).Let us denote by: (11) this nonlinear operator, where the superscript k of R k h stands for all variables computed at the previous period Conjugate gradient algorithms apply to symmetric positive definite linear equations.They can be extended to nonlinear equations with symmetric or slightly non symmetric positive definite Jacobian matrices.In order to be as close as possible to these properties, we have chosen to apply the conjugate gradient algorithm to the Schur complement system defined on the displacement variable obtained by computation of the pressure from Equation ( 11) and substitution in equation ( 9): (12) The nonlinear system ( 12) is clearly symmetric positive definite for single phase flow with q = 1 and is expected to Reservoir computations at each time Mechanical computations during the period y T = period ∆t = time step remain slightly non symmetric and positive in more general practical simulations.A natural symmetric positive definite preconditioner for the system ( 12) is the mechanical inverse operator: The residual gradient of ( 12) along the descent direction will be computed by finite difference of step denoted by ε.All these choices lead to the following preconditioned conjugate gradient algorithm, where l denotes the conjugate gradient iterative subscript, and (X, Y), X, Y ∈ IR N , the canonical scalar product of IR N .
Initialization: Compared with the previous staggered algorithm, each iteration involves one additional reservoir simulation in order to compute the gradient along the descent direction, and two mechanical residual computations.The evaluation of the mechanical residual is always negligible compared to the inversion of the mechanical operator.Furthermore, in practical situations for which the geomechanical model is a nonlinear plasticity model, the cost of a reservoir simulation over one period is much smaller than the cost of one mechanical computation.It results that the additional cost of each conjugate gradient iteration will be very small.Note also that this algorithm readily extends to nonlinear mechanical models with no significant additional cost.
The conjugate gradient algorithm exhibits a lack of robustness to non symmetry and non positivity of the system.In particular, although it has not been observed in the following numerical experiments, the algorithm could blow up if (y l , d l ) = 0 which is not excluded when the Jacobian of the system is not positive definite.More robust related algorithms will be investigated in a future work to improve the robustness of our method especially for highly nonlinear operators.

NUMERICAL EXAMPLE
In this section, the staggered algorithm presented in Section 3 is compared with the conjugate gradient algorithm of Section 4. The comparison is performed for a onedimensional test case dealing with an isotropic porous cylinder of radius and length (Fig. 3).
Water is injected at the bottom (incoming phase) and oil is extracted at the top (outcoming phase).There is no lateral movement in the x-and y-directions and we assume that the strains follow the z-direction (uniaxial strains).The total deformation u is fixed to zero, and the pressure p is imposed at the outcoming phase.At the incoming phase, the water flux Q w is prescribed to a given value and the oil flux Q o is set to zero.
The relative permeabilities are given by the following quadratic laws: The water and oil phases are assumed to be compressible, and we set: w (1 + c w p) where c w and c o denote the water and oil compressibility, and ρ 0 w and ρ 0 o the reference water and oil densities.The mechanical and fluid flow physical parameters used for the simulations are given in Table 1.
The meeting times between the reservoir and the mechanical simulators (the T k 's which are not equally spaced here) are given in minutes by: T 1 = 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 45, 60, T 20 = 120 The saturation equation is a Buckley-Leverett type problem where capillary effects have been neglected.It satisfies a first order nonlinear conservation law whose exact solution may develop chocks.The saturation plots displayed Figure 4 exhibit a steep front which approximates the discontinuity of the exact solution.The extraction of oil, which is five times more viscous than water, requires a very high pressure gradient decreasing in time as it is illustrated on Figure 4.
The reference solution is obtained with a fully coupled code and a computation of the geomechanics at each reservoir time step (q = 1).
For both the staggered and the conjugate gradient algorithms, the errors plotted Figures 5 and 6 between the computed solutions and the reference solutions are small, showing that the choice of the successive periods provides enough accuracy.
Note that the main differences on the saturation profiles are located at the discontinuities of the solution of the continuous problem, which could be expected.The order of magnitude of the saturation error for the staggered algorithm is lower than that of the conjugate gradient algorithm (Fig. 5).This is explained by the use of the exact porous medium compressibility for the staggered algorithm, which helps to decrease the effects of the pressure increase as water is injected.On the contrary, the magnitude of the pressure error for the conjugate gradient method is lower than that of the staggered algorithm (Fig. 6).This is due to a poor convergence of the Gauss-Seidel like method for such slightly compressible fluids.This lack of robustness appears more clearly in Tables 2 and 3, where j denotes a convergence criterion of 10 -j on the maximum relative variation of pressure between two successive iterations.Note in particular on Table 2, the nonconvergence of the staggered algorithm at the 18th period over a total of 20 periods and for j = 3. TABLE 2 Total number of mechanical computations and reservoir period simulations for the staggered algorithm over the 20 periods for a given accuracy 10 -j on the maximum relative variation of pressure between two successive iterations Total number of mechanical computations and reservoir period simulations for the conjugate gradient algorithm over the 20 periods for a given accuracy 10 -j on the maximum relative variation of pressure between two successive iterations The pressure error between the previous reference solution and, on the left, the staggered algorithm solution, and, on the right, the conjugate gradient algorithm solution.

CONCLUSION
In this paper we have introduced a new coupling algorithm between the geomechanical and the reservoir simulators.The coupling is defined on time periods in such a way that the reservoir unknowns are computed for time steps which are small enough subdivisions of the time period whereas the mechanical problem is solved at the end of the period.
Our new approach is based on a nonlinear preconditioned conjugate gradient method which is applied to the mechanical displacement variable.For each period and each iteration of the conjugate gradient algorithm, two reservoir simulations over the period are performed to compute the residual of the conjugate gradient algorithm and the mechanical computation plays the role of the preconditioning step.
This algorithm is compared to the staggered coupling algorithm with a porous medium compressibility as relaxation parameter.The main conclusion is that the preconditioned conjugate gradient method is much more robust and converge much faster than the staggered algorithm, while the additional cost per iteration should remain in practical situations very small.

Figure 3
Figure 3Cylinder with no lateral movement.