Coupling Fluid Flow and Rock Mechanics: Formulations of the Partial Coupling Between Reservoir and Geomechanical Simulators

Résumé — Formulation du couplage partiel entre un simulateur réservoir et un simulateur de géomécanique — La compaction des roches-réservoirs à forte porosité est un phénomène complexe qui dépend du comportement constitutif de la roche-réservoir, du chemin de contrainte suivi par le réservoir, etc. La compaction des réservoirs pétroliers faiblement consolidés ne peut donc pas être modélisée à partir des simulateurs réservoirs traditionnels. En effet, ces simulateurs réduisent les effets mécaniques à un paramètre de compressibilité de pore qui ne reproduit pas parfaitement la complexité des phénomènes mécaniques. Une analyse plus précise de la compaction d’un réservoir pétrolier nécessite donc le recours à une modélisation thermo-hydro-mécanique du comportement du réservoir et de ses épontes. Les équations résultant d’un tel modèle peuvent être soit résolues simultanément dans un même simulateur (couplage total), soit résolues séparément dans un simulateur réservoir et un simulateur géomécanique avec échange de données entre les deux (couplage partiel). L’article présente trois formulations de couplage partiel obtenues sous les hypothèses d’un comportement linéaire élastique de la roche-réservoir et d’un écoulement monophasique. Ce cadre poroélastique permet Abstract — Coupling Fluid Flow and Rock Mechanics: Formulations of the Partial Coupling between Reservoir and Geomechanical Simulators — During high porosity reservoir production, the rock compaction is a complex phenomenon that depends on the rock constitutive behavior, the reservoir stress path, etc. Reservoir compaction can hardly be analyzed with conventional reservoir simulators as the pore factor, the only mechanical parameter used in such simulators, is not sufficient to represent the complex phenomena involved. In order to solve correctly this problem, the full thermo-hydro-mechanical equations must be addressed. The corresponding set of equations can be either solved simultaneously (fully coupled scheme) or using a conventional reservoir simulator in conjunction with a geomechanical simulator and information exchanges between the two simulators (partial coupling). The paper presents three formulations of the partial coupling, which are obtained in the framework of single-phase flow and a linear elastic isotropic rock behavior. This simple framework makes possible an easy and rigorous derivation of the porosity correction to be appended to the reservoir Lagrange’s porosity used in the reservoir simulator. The porosity correction depends on the pore compressibility factor and a mechanical contribution that can be expressed either in terms of volumetric strain, pore volume change, or the mean total stress change. One formulation is tested on a numerical test that depicts the water flood through a laboratory core sample initially saturated with oil and constrained to uniaxial strain. The numerical test illustrates the importance of the mechanical effects on the fluid flow problem and validates the partial coupling proposed. The example also highlights the role of the pore compressibility factor in the partially coupled reservoir simulation. Actually, in the partially (iteratively) coupled approach, the pore compressibility factor can be interpreted as a relaxation parameter controlling the convergence speed of the iterative process between reservoir simulation and geomechanical simulation.


INTRODUCTION
Hydrocarbon or gas reservoir production induces variations in time and space of reservoir pressure, saturation and temperature.In turn, changes in thermal and hydraulic reservoir properties may cause a modification of the stress state in and around the reservoir.The stress changes may then alter the reservoir fluid flow parameters and then the reservoir production scenario.Furthermore, they may also cause a threat to well integrity when the change in stress state leads to significant ground strains around wells.Geomechanical effects can be particularly pronounced in stress sensitive reservoirs as poorly compacted reservoirs and fractured or faulted consolidated reservoirs.For poorly compacted reservoirs, the stress changes may have benefit effects on fluid recovery due to the reservoir compaction.However, the reservoir compaction may also reduce the reservoir permeability, cause surface subsidence and create damage on well equipment.Fractured and faulted reservoirs are particularly affected by stress changes induced by reservoir thermal variations (cold water injection, steam injection, etc.).The resulting stress changes may enhance or reduce the fracture conductivity or create preferred flooding directions.Most of the time, in standard reservoir simulators, the description of the mechanical behavior of the reservoir and its possible adjacent formations are neglected or reduced to a weak contribution.This is for instance the case for reservoir compaction that is only accounted for in standard reservoir simulator with a pore compressibility factor.This simplified approach has been found to be insufficient to provide an exact description of the realistic phenomena that may occur in a stress dependent reservoir.The influence of geomechanics in reservoir simulations has been illustrated by Rhett and Teufel (1992), Ruistuen et al. (1996) and Gutierrez the pore compressibility factor, the only mechanical parameter used in such simulators, is not sufficient to represent the complex phenomena involved.In order to solve correctly this problem, the full thermohydro-mechanical equations must be addressed.The corresponding set of equations can be either solved simultaneously (fully coupled scheme) or using a conventional reservoir simulator in conjunction with a geomechanical simulator and information exchanges between the two simulators (partial coupling).The paper presents three formulations of the partial coupling, which are obtained in the framework of single-phase flow and a linear elastic isotropic rock behavior.This simple framework makes possible an easy and rigorous derivation of the porosity correction to be appended to the reservoir Lagrange's porosity used in the reservoir simulator.The porosity correction depends on the pore compressibility factor and a mechanical contribution that can be expressed either in terms of volumetric strain, pore volume change, or the mean total stress change.One formulation is tested on a numerical test that depicts the water flood through a laboratory core sample initially saturated with oil and constrained to uniaxial strain.The numerical test illustrates the importance of the mechanical effects on the fluid flow problem and validates the partial coupling proposed.The example also highlights the role of the pore compressibility factor in the partially coupled reservoir simulation.Actually, in the partially (iteratively) coupled approach, the pore compressibility factor can be interpreted as a relaxation parameter controlling the convergence speed of the iterative process between reservoir simulation and geomechanical simulation.
and Lewis (1998) who provided evidence that reservoir compressibility and permeability strongly depend on the stress path followed by the reservoir.This result demonstrates that the pore compressibility used in reservoir simulators is not sufficient enough to reproduce pore volume changes induced by pressure and temperature variations (see Tortike and Farouq, 1993;Settari and Mourits, 1998).Osorio et al. (1998Osorio et al. ( , 1999) ) emphasised the role of the reservoir surrounding domain (the reservoir extended stress-disturbed region induced by reservoir production) in their coupled fluid-flow/geomechanical modeling.They particularly investigated the effects of the mechanical boundary conditions (applied on the reservoir surrounding domain) and the mechanical parameters on the reservoir productivity.The water weakening effect described by Sylte et al. (1999) and Hermansen et al. (2000) (reduction of the reservoir mechanical strength with increasing water saturation) can't be accounted for with a conventional reservoir simulator.Chin and Thomas (1999) investigated the role of the rock constitutive behavior, stress-dependent permeability, overburden/reservoir interaction and water weakening effect on the oil recovery of a poorly compacted reservoir.In the case of fractured reservoirs, Gutierrez and Makurat (1997) showed that fracture permeability changes due to temperature and fluid pressure variations induced by reservoir production can significantly alter the water flood behaviour of the reservoir.The evidence of preferred flooding directions during water flood in both "naturally fractured" and "unfractured reservoirs" have been illustrated by Koutsabeloulis et al. (1994) and Heffer et al. (1994).
The previous references show that, for stress dependent reservoir, a conventional reservoir simulation may lead to significant inaccuracies in the oil production forecast given by the numerical simulations.These inaccuracies can be reduced when the problem is analysed in its whole from the thermo-hydro-mechanical point of view.There are two different approaches for solving the thermo-hydromechanical equations: the fully coupled approach simultaneously solves the all set of equations in one simulator, whereas the partially coupled approach uses a conventional reservoir simulator in conjunction with a conventional geomechanical simulator.This paper presents different formulations of the partially coupled approach, which describe the porosity changes resulting from the stress computation and that have to be implemented in a conventional reservoir simulator to account for compaction phenomenon.Partial coupling formulations are derived under isothermal conditions and assuming a linear poroelastic behavior of the rock saturated with one fluid.The linear poroelastic framework used in this paper makes possible an easy presentation of the partial coupling methodology and formulations.Finally the last section of the paper compares conventional reservoir, partially coupled and fully coupled modeling on a dead oil case with porosity changes induced by pressure changes.This simple onedimensional example illustrates the effects of the rock mechanical behavior on the fluid flow problem and the role of the pore compressibility factor on the partially coupled approach.

LINEAR POROELASTICITY
We recall in this section the basic equations that govern the problem of consolidation of a porous medium as firstly introduced by Biot (1941).The porous medium is considered as the superposition of a moving fluid saturating the connected porous space and a deformable skeleton (the solid matrix and the connected pore space without the fluid).The following assumptions on the fluid and the skeleton are introduced: -the fluid flows through the porous medium according to Darcy's law; -the fluid is compressible; -the porous medium is isotropic; -the skeleton transformations are infinitesimal; -the behavior of porous medium is poroelastic (i.e.linear and reversible).
In what follows, we use the sign convention of continuum mechanics for which stress and strain are positive in tension.According to the previous assumptions, one can obtain in Euler coordinates the governing equations describing the deformation of the skeleton and the motion of the fluid in the porosity.Using the skeleton displacement vector and the fluid pressure as primary variables, the governing equations read (Biot, 1941;Boutéca, 1992;Coussy, 1995;Lewis and Schrefler, 1998): (1) (2) with ∇ 2 = the Laplacian operator, and ∇ = the gradient operator.The displacement vector defines the strain tensor ε with: (3) The volumetric strain ε v appearing in Equation ( 1) is given by the trace of the strain tensor: (4) The Biot's modulus M is related to the rock and fluid characteristic with (Coussy, 1995): where φ 0 is the reference porous medium porosity (Lagrangian or Eulerian definition).The Biot's coefficient b is related to the matrix and drained bulk moduli K s and K d with the Biot relation: (6) Equation ( 1) arises from the fluid mass balance equation whereas Equation (2) characterises the equilibrium equations in case of zero body forces.Equations ( 1) and ( 2) are fully coupled because the displacement vector and the pore pressure appear in each equation.On the one hand, the displacement vector acts on the pressure Equation (1) by means of the volumetric strains.The resulting term (i.e. the right hand side of Equation ( 1)) accounts for the effects of the rock bulk volume variations induced by the skeleton deformation on the pressure changes.On the other hand, the pore pressure gradient affects the stress equilibrium equation through the right-hand side of Equation (2).The full set of Equations (1-4) with boundary conditions and initial condition can be solved analytically for simplified geometry or numerically for more general cases.

CONVENTIONAL RESERVOIR APPROACH BASED ON RESERVOIR COMPRESSIBILITY
In reservoir engineering, the effects of the stress variations on the fluid behavior are neglected or simplified so that the reservoir engineer focuses on the fluid flow problem.Incorporating Darcy's law in the fluid mass balance equation, one obtains: (7) The porosity appearing in Equation ( 7) must be considered as an Euler's porosity, i.e. the ratio of the pore volume and the bulk volume defined in the deformed configuration: (8) where V p is the pore volume in the deformed configuration and V b is the bulk volume in the deformed configuration.However, the porosity used in reservoir simulators refers to the initial configuration (Settari and Mourits, 1994) and therefore must be considered as a Lagrange's porosity, i.e. the ratio of the pore volume V p defined in the deformed configuration and the bulk volume V b 0 in the initial configuration: (9) Within the assumption of infinitesimal deformation, Euler's and Lagrange's porosities are related with the volumetric strain by: (10) Note that the initial Euler's and Lagrange's porosities are equal as long as no deformation has occurred, and they are noted φ 0 in this paper.Furthermore, Equation (10) shows that the assumption that the Lagrange's porosity may be used in reservoir simulators only holds for non deformable porous media.In reservoir simulators, it is a common practise to introduce a pressure dependency to the reservoir porosity (i.e. the Lagrange's porosity).A typical equation for the Lagrange's porosity is (we recall that the temperature is supposed constant here): (11) According to Equation ( 11) and assuming a linear constitutive behavior, the reservoir porosity linearly depends on the pore pressure and the proportionality coefficient, the pore compressibility c p , can be computed from: where V p 0 is the pore volume in the reference configuration.This coefficient accounts for the in situ pore volume strains that follow from changes in the reservoir pore pressure.It is measured by reservoir engineers under a reservoir specific and constant stress path.Note that exponential dependency of the porosity-pressure relationship are also used in reservoir simulators.Equalling Expression (11) of the Lagrange's porosity used in reservoir simulator and the expression of the Euler's porosity deduced from linear poroelasticity equations for isotropic material, Settari and Mourits (1994) show that the pore compressibility should be given by (at the first order): (13) where c s and c b are material properties: c s = 1/K s is the matrix compressibility and c b = 1/K d is the drained compressibility.In addition to the pore compressibility term, Settari and Mourits (1994) propose to add to the reservoir porosity given by ( 11) with (13) a porosity correction deduced from a geomechanical simulation and given by: ( 14) where σ m is the mean total stress = σ ii /3, resulting from the geomechanical simulation.
The following developments show how Equation ( 7) used with a Lagrange's porosity can be transformed in a pressure diffusion equation representative of the problem solved in a reservoir simulator.Conventional reservoir simulators consider compressible fluid with a density depending on the pore pressure.The linearized form of the fluid state law reads: (15) where ρ f 0 is the initial fluid density and c f is the fluid compressibility accounting for fluid density changes with pore pressure and is defined by: ( 16) According to the reservoir approach, we use the Lagrange's porosity in Equation ( 7).Expanding the first term of the left-hand side of Equation ( 7), we get: Introducing expressions ( 11) and ( 16) of the Lagrange's porosity and the fluid density, and linearizing the resulting equation yields the following pressure diffusion equation (see also Boutéca, 1992;Cossé, 1993, in the case of a nondeformable porous media for which c p = 0).( 18) Equation ( 18) can be used with boundary and initial conditions to predict pressure variations in the reservoir.Note that according to this equation, the pore compressibility remains the only mechanical parameter involved in a reservoir simulator.

FROM POROELASTICITY TO A RESERVOIR SIMPLIFIED DIFFUSIVITY EQUATION
This section shows how the diffusion Equation ( 18) can be recovered from simplifying assumptions in the geomechanical model (i.e.Eqs. ( 1) and ( 2)).As for the reservoir engineering approach, the whole attention is paid to the pressure diffusion Equation (1) whereas the mechanical problem given by Equation ( 2) is simplified.First, we assume that the rock matrix is uncompressible (i.e.K s = ∞) so that Equations ( 5) and ( 6) give the following approximations of the Biot's modulus and coefficient: Next, the right-hand side of Equation (1) has to be approximated as a function of the pore pressure.Using the last approximation of the Biot's modulus and the fact that the pore and bulk volume variations are equalled for an uncompressible matrix, the right-hand side of Equation ( 1) is approximated on the form: (20) with the pore compressibility given by Equation ( 12).When introduced in Equation (1), approximations ( 19) and ( 20) provide the pressure diffusion Equation ( 18) usually encountered in reservoir engineering.
Despite the "realistic" hypothesis of an uncompressible rock matrix, the diffusion Equation ( 18) and more generally the standard reservoir approach disregarding geomechanical coupling has been found to be inappropriate in numerous studies concern with stress dependent reservoirs (e.g.Chin et al., 1998;Gutierrez and Lewis, 1998;Osorio et al., 1998;Settari and Mourits, 1998).This comes from the lack of consideration for the geomechanical problem that can significantly influence reservoir engineer predictions.

STRESS DEPENDENT RESERVOIR SIMULATORS
A stress dependent reservoir simulator is a reservoir simulator that integrates the geomechanical problem.Therefore, a stress dependent reservoir simulator solves the fluid pressure problem (i.e.Eq. ( 1)) with the stress problem (i.e.Eq. ( 2)) and possibly the thermal problem.The hydrothermo-mechanical problem can then be solved using different degrees of coupling as firstly described by Settari and Mourits (1994) and Settari and Walters (1999).
-The fully coupled approach consists in simultaneously solving the all set of equations that govern the hydrothermo-mechanical problem (i.e.Eqs. ( 1) and (2) in our case).However, in the fully coupled approach, the hydraulic or geomechanical mechanisms are often simplified by comparison with conventional uncoupled geomechanical and reservoir approaches.-In the partially coupled approach, the stress and flow equations are solved separately for each time step but information is passed between the reservoir and geomechanical simulators.Therefore the reservoir and geomechanical problems have to be reformulated according to the original fully coupled problem.The partial coupled approach is termed explicit if the information exchange between both simulators is only performed once for each time step, or iterative if it is repeated until convergence of the stress and fluid flow unknowns.Contrarily to the fully coupled approach, the partial coupling looks more flexible and benefits from the high developments in physics and numerical techniques of both the reservoir simulator and the mechanical software.
The fully coupled approach does not produced particular difficulties because the all equations are solved in the same simulator.This is not the same for the partial coupling, for which the coupling methodology is detailed in the next subsections.The key idea of the partially coupled approach is the reformulation of the stress-flow coupling such that a conventional stress analysis code can be used in conjunction with a standard reservoir simulator.We present in this section how this reformulation can be realised for Equations ( 1) and (2).

Formulation of the Stress-Flow Coupling in the Stress Simulator
This is the easy part of the partial coupling because the righthand side of Equation ( 2) is the only term that distinguishes Equation ( 2) from the stress equilibrium equation encountered in classical mechanics.As previously mentioned, the partial coupling is based on an exchange process between the reservoir and the stress simulators.Thus, assuming that the pressure is known from the last reservoir simulation, one can evaluate the right-hand side of Equation ( 2) with reservoir pressure values.Consequently, pore pressure variations (and more generally temperature changes) computed by the reservoir simulator are converted in distributed load in the stress simulator.
When applying the distributed loads, the stress analysis performed on the reservoir and its adjacent formations provides the rock stress and strain.The rock strain and thus the volumetric strain is in turn used in the standard reservoir simulator with the methodology described in the next section.

Formulation of the Stress-Flow Coupling in the Reservoir Simulator
According to the methodology of partial coupling, Equation (1) arising from the poroelasticity theory has to be reformulated with a left-hand side similar to the left-hand side of the diffusion Equation ( 18) that is representative of the problem solved in reservoir simulators.This section details different possible expressions of the right-hand side of the equation resulting from this process.These expressions use the values of the pressure computed at the previous iteration and, depending on the stress simulator outputs, the volumetric strain, the pore volume variation or the average total stress.Following the previous methodology, Expression (5) of the Biot's modulus M is introduced in Equation (1).One obtains: (21) Equation ( 21) constitutes an appropriate formulation for the partial coupling.According to this formulation, a source term given by the right-hand side of Equation ( 21) and depending on the pore compressibility factor must be considered in the standard reservoir simulator.When discretizing Equation ( 21), the pore pressure time derivative appearing in the source term has to be approximated with the pore pressure time derivative computed at the previous reservoir iteration (iteratively coupled) or time step (explicitly coupled).The volumetric strain derivative of Equation ( 21) is evaluated with the difference between the last volumetric strain computed with the stress simulator and the same quantity given at the beginning of the time step.Hence the partial coupling can be simply achieved with a modification of the second member of the linear system classically solved in conventional reservoir simulators.Furthermore, the second member modification due to the partial coupling can be simply carried out when adding the right-hand side of Equation ( 21) to the Lagrange's reservoir porosity at the current time step (Bévillon and Masson, 2000).Thus the porosity correction expressed in terms of pore pressure and volumetric strain variations to be appended at each iteration (iteratively coupled) or time step (explicitly coupled) reads: The pore volume compressibility c p appears in Equations ( 21) and ( 22).In conventional reservoir simulation this parameter is determined by the reservoir engineer.In partial coupling, it can be considered as a numerical parameter because whatever the value used by the reservoir engineer is, the geomechanical simulator gives the exact porosity from which the correction to be applied to the reservoir simulator is calculated.Hence the main role of the reservoir pore compressibility appears to act for the stabilisation of the partially coupled approach.Indeed, if no compressibility is considered, Equation (21) takes the form: (23) Therefore, for uncompressible fluid, the stability of the coupled approach is not always ensured (see Bévillon and Masson (2000) for a stability criterion).Despite the numerical nature of the pore compressibility parameter, the convergence rate of the method depends on the magnitude of the pore volume correction to be applied in the reservoir simulation.Most of the time, the stress path followed by the reservoir is estimated from ideal cases where either the total stresses are blocked, or lateral displacements are blocked and the vertical load is constant (i.e.oedometric deformation).
In the case where the total stresses are blocked (i.e.∂σ kk = 0), the pore compressibility is given by: (24) In case of oedometric deformation (i.e.∂σ z = 0 et ∂ε r = 0), the pore compressibility is given by: (25) Equation ( 21) can be reformulated introducing the pore volume correction or the mean total stress instead of the volumetric strain.The linear theory of poroelasticity (for an isotropic porous medium) gives the following relation between the Lagrange's porosity, the total volumetric strain and the pore pressure (Coussy, 1995;Dangla, 1999): Introducing Equation (26) in Equation ( 21), it comes: (27) Equation ( 27) provides a second possible formulation that can be used for the coupling between the reservoir and the geomechanical simulators.Contrarily to the first formulation, the modification of the second member of the conventional reservoir simulator now depends on the pore pressure derivative (taken at the previous iteration or time step) and the Lagrange's porosity variation (i.e. the pore volume change) given by the last stress simulation.With this formulation, the pore compressibility always appears as the numerical parameter and the right-hand side of Equation ( 27) can still be interpreted as a reservoir porosity correction of the form: (28) Note that for an uncompressible rock matrix for which b ≈ 1, c s ≈ 0 and dV p ≈ dV b , both Formulations ( 21) and ( 27) are the same.Finally, let us introduce a last formulation of the partial coupling approach using the mean reservoir stress.For an isotropic porous medium, the first law of behavior of linear poroelasticity leads to: (29) with σ m the mean total stress = σ ii /3.Subtracting Equation (29) to the same equation multiplied by c s /c b , we get with Equation ( 6): (30) This last equation can be introduced in (21) to give: (31) Equation ( 31) provides a third possible formulation for the coupling between the reservoir and the geomechanical simulators.This formulation defines a porosity correction to be appended to the reservoir porosity and depending on the pore pressure derivative (taken at the previous iteration or time step) and the mean total stress computed by the stress simulator: (32) Note that when using the pore compressibility factor (13) recommended by Settari and Mourits (1994), the term in factor of the pressure time derivative vanishes and the same reservoir porosity correction as the one proposes by the previous authors is obtained (see Eq. ( 14)).Contrarily to Formulations ( 21) and ( 27), the last formulation (32) using the geomechanical mean total stress variable does not easily generalise to complex rock constitutive behavior.The different formulations of the partial coupling approach for linear poroelasticity are gathered in Table 1.The main advantage of these formulations compared to the formulation proposed by Settari and Mourits (1994) is that the pore compressibility factor is not fixed but can be chosen by the reservoir engineer according to either a known physical value or either stability considerations.

EXAMPLE FOR A DEAD OIL CASE
The aim of this section is two folds: firstly, to highlight the differences between reservoir and hydro-thermo-mechanical simulators as mentioned in the Introduction, and secondly to illustrate the methodology of partial iterative coupling as introduced in Section 4. The section also analyses the influence of the pore compressibility factor on the iteration numbers required in the iteratively coupled simulation.
Numerical tests presented here are realised on a onedimensional example that makes possible an easy development and comparison between reservoir, fully coupled and iteratively coupled reservoir simulators.

Example Description
We consider an isotropic porous cylinder of radius R = 0.25 m and length L = 1.50 m (Fig. 1).For this given geometry, we successively describe the mechanical and fluid flow problems addressed.
Figure 1 Cylinder with no lateral movement.
For the mechanical problem, we consider that there is no lateral movement in the x-and y-directions and we assume that the strains only occurs in the z-direction (uniaxial strain hypothesis).The boundary condition consists in a total displacement ∆u z that is applied to the whole cylinder in the z-direction.According to the geometry of the problem the strain and stress tensors read in cylindrical coordinates: (33) Using the law of linear poroelasticity for an initial isotropic state given with a null fluid pressure and a total stress σ 0 , the stress and strain are related to the fluid pressure with: (34) (35) where K uni is the bulk modulus in uniaxial condition and λ d is the Lame's constant in drained conditions.K uni and λ d are related to the rock shear modulus G and the drained bulk modulus K d with: (36) (37) The stress equilibrium equation written with no body forces shows that the uniaxial total stress σ z is constant (i.e.does not depends on the cylindrical coordinates r, θ and z).Using this result and the boundary condition (i.e. we integrate Equation (34) in the zdirection.One obtains: (38) Equation ( 38) provides a useful relation that relates the uniaxial stress variation (σ z -σ 0 ) to the total pressure change in the cylinder.When the uniaxial stress variation is known, it is possible to compute the strain ε z as a function of the pressure with: (39) Note that if we assume that the vertical load is constant (i.e.σ z = σ 0 ), the volumetric strain is proportional to the pore pressure variation and the proportionality coefficient can be used to compute the pore compressibility factor in oedometric deformation (Eq.( 25)).As a consequence, in oedometric deformation, the change in Lagrange's porosity is proportional to the pore pressure change.This is not the case for the mechanical problem considered here where we assume that the total vertical displacement is blocked instead of a constant vertical load.Thus, for the present case, the Lagrange's porosity change depends not only on the pore pressure change but also on the vertical stress variation (Eq.( 39)).Furthermore, as the vertical stress variation depends on the total pressure variation (Eq.( 38)), the local porosity change can't be proportional to the local pressure change.This result explains the discrepancy between conventional reservoir simulation and fully coupled reservoir simulation presented in Section 5.2.
For the fluid flow problem, we consider that the cylinder is initially saturated with oil and that a constant water flux is imposed at the top of the cylinder, a null pressure being imposed at the bottom boundary.The fluid flow problem can be interpreted as an isothermal water flood through a laboratory core sample.The mass balance equations for the oil and water read (gravity effects are neglected): where S w is the water saturation, for fluid i (i = o for oil and i = w for water), ρ i is the fluid mass density, µ i is the fluid dynamic viscosity, k ri is the fluid relative permeability, and p i is the fluid pressure.For the simulations presented here, the following assumptions are introduced: -capillary effects are neglected, so that p o = p w = p.This assumption allows us to define a single pressure for the fluid mixture, which can be considered as the equivalent pressure of the single-phase flow model supposed in linear poroelasticity; -relative permeabilities follow the quadratic laws of the form: (42) -the oil is noncompressible; -the water is assumed compressible, with its mass density determined by (c w is the water compressibility and ρ w 0 is the initial water mass density): (43) Finally, Table 2 gives parameter values used in the simulations presented in the following sections.Fluid flow and mechanical parameters correspond to an elastic highly deformable rock with an uncompressible rock matrix and high porosity and permeability.The fluid flow problem is solved with a finite volume method on a 150 grid size mesh whereas mechanical unknowns are computed from the direct discretization of the mechanical Equations ( 38) and (39).

Comparison between Fully Coupled and Conventional Reservoir Simulations
This section compares the solutions obtained with a reservoir simulator and a fully coupled simulator.
-The reservoir simulator solves Equations ( 40) and (41) (with ( 42) and ( 43)) with the finite volume method (Eymard et al., 2000).The scheme derived for the resolution of Equations ( 40) and ( 41) uses a pore compressibility factor c p as described in Section 2 of this paper.Its value is fixed to the pore compressibility factor in oedometric condition.Taken into account the uncompressible rock matrix assumption, the pore compressibility factor equals 1/(φ 0 K uni ).The time step size of the reservoir simulator is not constant but computed according to a predictor-corrector criterion.-The fully coupled approach simultaneously solves the two mass balance Equations ( 40) and ( 41) with the finite volume method and Equation (38) discretized with the unknown values of the pore pressure.This means that, assuming that the number of finite volume cells is N, the fully coupled approach turns to simultaneously solve 2N+1 equations (i.e.Equations (40) (41) for each finite volume cell + Equation ( 38)) with 2N+1 unknowns (i.e.water saturation and pore pressure at each finite volume cell + the uniaxial stress).Due to the uncompressible rock matrix assumption, the porosity of each finite volume cell appearing in the mass balance Equations (40) ( 41) is computed as the initial porosity plus the volumetric strain (i.e.ε z ) of each finite volume cell.These last values are deduced from the discretized form of Equation ( 39). Figure 2 compares the oil production (i.e. the oil volume recovered from the sample) obtained with the two simulations.The oil productions obtained with each simulator significantly differ during the two hours of the "numerical" experiment.As previously mentioned, the discrepancy between the two simulations arises from the fact that the reservoir simulation assumes oedometric deformation whereas in the fully coupled simulation the total vertical displacement is supposed to be zero.Note that the gap in oil production slowly decreases as time flows since the final volume of oil produced must converge at infinity towards the volume of oil initially present in the core sample (i.e.φπR 2 L ≈ 0.088 m 3 ).Figure 2   1 , Figure 3 displays the pore pressure obtained with both simulators at different values of time.One can notice that the pressure plots strongly differ since the short time values.This difference arises from the pore volume reduction that is predicted in the bottom of the cylinder (Fig. 4) by the coupled simulation due to the boundary condition of zero displacement.The pore volume reduction leads to an instantaneous increase of the pressure for the fully coupled approach and that cannot be accounted for with the reservoir approach (Gutierrez and Lewis, 1998).The pore pressure increase together with the pore pressure boundary condition at the bottom of the cylinder initiate an early oil production in the fully coupled simulation than for the conventional reservoir simulation (Fig. 2).  Figure 5 displays the water saturation obtained with both approaches.The water saturation profiles are very closed at the beginning of the experiment.However, the water saturation "front" computed with the fully coupled simulator progresses more rapidly in the cylinder.This leads to an early water breakthrough for the fully coupled model than for the reservoir model (see also Fig. 2).Finally, Figure 4 presents the porosity obtained with the reservoir simulator (initial porosity + porosity change due to pressure change according to the pore volume compressibility factor) and the fully coupled simulator (initial porosity + volumetric strain computed from Equation ( 39)).On the one hand, the porosity change computed with the reservoir simulator is necessarily a porosity increase due to the pore pressure increase.On the other hand, the fully coupled simulator predicts a porosity reduction in the bottom zone of the cylinder.This is due to the null strain boundary condition imposed to the whole cylinder.Due to the porosity increase in the top zone of the cylinder, the null boundary displacement of the whole cylinder necessarily implies a porosity decrease (and then a fluid compression) in the bottom of the cylinder.This effect cannot be accounted for in the reservoir simulator where the mechanical problem is not addressed.

Comparison between Fully Coupled and Partially (Iteratively) Coupled Simulations
This section compares the solutions obtained with the fully coupled simulator and the iteratively coupled reservoir simulator.The fully coupled simulator is described above, so we just describe the iteratively coupled simulator as firstly introduced in Section 4. The iteratively coupled simulator is based on an iterative scheme between reservoir simulation and mechanical simulation (illustrated here by the resolution of Equations ( 38) and ( 39)).The reservoir simulator uses the oedometric pore compressibility c p = 1/(φ 0 K uni ) and a porosity correction given by Equation ( 22) with c s = 0 and b = 1 due to the incompressible rock matrix assumption.The mechanical simulator consists here in the computation of the uniaxial stress σ z (see Equation ( 38)) using the finite volume pore pressures given by the reservoir simulator.Using these last values and Equation (39), it is possible to compute the new volumetric strain that will be used to compute the porosity correction needed for the next reservoir simulation.
Reservoir simulator and mechanical simulator are run at each time step until a convergence criterion in porosity is reached: the difference between the "reservoir porosity" and the "mechanical porosity" must be less than 10 -6 times the initial porosity.This criterion is used for all the iteratively coupled simulations because it ensures an acceptable convergence of the pressure, and therefore, a high convergence of the porosity, saturation and oil production.Figure 6 compares the oil production obtained with the fully coupled simulator and the iteratively coupled simulator.As expected, the two approaches exactly give the same results.This can also be observed on Figure 7 were the pore pressures obtained with both simulators are superimposed.These comparisons validate the partially coupled approach introduced in this paper.Iteratively coupled simulations presented on Figures 6 and 7 are realised with a pore compressibility value fixed to the oedometric compressibility: c p = 1/(φ 0 K uni ).In the partially coupled approach, this parameter of the reservoir simulator is in fact a "numerical parameter" that can be interpreted as a relaxation factor.To illustrate the role of this parameter, Figure 8 plots, for different values of the pore compressibility factor, the iteration numbers required to reach the same convergence level of the partial coupling scheme at each time  Iteration numbers obtained with the iteratively coupled simulator for different values of the reservoir compressibility: ( ) c p = 1.5/(φ 0 K uni ), (s) c p = 1/(φ 0 K uni ), (O) c p = 0.6/(φ 0 K uni ).
step.The oil production history obtained for the different iteratively coupled simulations are superimposed and thus not compared on an additional figure.As shown by Figure 8, the beginning of the numerical experiment involves small time steps that necessitate high iteration levels due to the low stability of the system for small time step sizes.Figure 8 demonstrates that 1/(φ 0 K uni ) is not the optimum parameter for the iteratively coupling scheme because the iteration numbers can be reduced when using a pore compressibility equal to 0.6/(φ 0 K uni ).Note that it was not possible to reduce the iteration numbers when reducing the pore compressibility to 0.5/(φ 0 K uni ) for which divergence problem starts to occur.
also shows that the water breakthrough (characterised by the change in slope of the oil production curve) occurs first for the fully coupled simulation after approximately half an hour and next for the conventional reservoir simulation after 40 min of water injection.

Figure 2
Figure 2Comparison of the oil production for the fully coupled model (solid lines) and the reservoir model (solid line and ).

Figure 3
Figure 3 Comparison of the pressure in the cylinder for the fully coupled model (solid lines) and the reservoir model (dashed line) at different values of time: ( ) 1 min, (s) 15 min, (O) 30 min.

Figure 4
Figure 4 Comparison of the porosity in the cylinder for the fully coupled model (solid lines) and the reservoir model (dashed line) at different values of time: ( ) 1 min, (s) 20 min, (O) 30 min.

Figure 5
Figure 5Comparison of the saturation in the cylinder for the fully coupled model (solid lines) and the reservoir model (dashed line) at different values of time: ( ) 1 min, (s) 15 min, (O) 30 min.

Figure 6
Figure 6Comparison of the oil production obtained with the fully coupled simulator (solid lines) and the iteratively coupled simulator ( ).

Figure 7
Figure 7Comparison of the pressure in the cylinder for the fully coupled simulation (solid lines) and the iteratively coupled simulation at different values of time: (×) 30 s, ( ) 2 min, (s) 5 min, (O) 10 min.

TABLE 1
Porosity correction as a function of the geomechanical variable