# Algorithms for coupled mechanical deformations and fluid flow in a porous medium with different time scales

Fatima-Zahra Daïm, Danielle Hilhorst, Jacques Laminie, Robert Eymard

► To cite this version:

Fatima-Zahra Daïm, Danielle Hilhorst, Jacques Laminie, Robert Eymard. Algorithms for coupled mechanical deformations and fluid flow in a porous medium with different time scales. International Journal of Numerical Analysis and Modeling, 2009, 5 (4), pp.635-658. hal-01449039

HAL Id: hal-01449039

https://hal.science/hal-01449039

Submitted on 30 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ALGORITHMS FOR COUPLED MECHANICAL DEFORMATIONS
AND FLUID FLOW IN A POROUS MEDIUM
WITH DIFFERENT TIME SCALES

FATIMA-ZOHRA DAÏM, DANIELLE HILHORST, JACQUES LAMINIE, AND ROBERT EYMARD

**Abstract.** In this paper, we solve a problem describing the mechanical de-
formations of a porous medium in the presence of a monophasic linear flow
or a two phase nonlinear flow with the purpose of modelizing subsidence of
hydrocarbon reservoirs. An essential characteristics of this problem is that the
mechanical deformation and the flow have different time scales. In petroleum
industry, one uses different very efficient simulators for the flow problem and the
mechanical deformations, which enables to handle complex models. Therefore
it is necessary to be able to combine as efficiently as possible the exploitation
of these simulators. We propose two alternative splitting approaches. The first
one is the staggered algorithm used by engineers, which amounts to a Gauss-
Seidel method in the one phase linear case. The second approach is based upon
the preconditioned conjugate gradient method. We use a numerical multi-scale
method in both of these algorithms. We compare these two approaches and we
show that the preconditioned conjugate gradient algorithm is faster and more
robust than the staggered algorithm. Applying the preconditioned conjugate
gradient algorithm therefore seems to compensate for the fact that the inf-sup
condition for the mixed discretization method is not satisfied when combining
the simulators for the mechanical deformations and for the flow computations.

**Key Words.** Porous media, Darcy flow, Mechanical deformations, Gauss Seidel method, Nonlinear conjugate gradient method, Inf-Sup condition, Mixed formulation, Multiscale algorithm.

1. Introduction

The production of oil and gas in soft highly compacting reservoirs induces an
important reduction of the pore volume, which increases the oil productivity. This
compaction leads to undesirable effects such as surface subsidence or damage of well
equipments. A well- known example of subsidence is the Ekofisk field in the North
Sea in Norway, where a sea floor subsidence rate of 42 cm/year has been reached
at the end of 1993 (see [11]). The cases of the Valhall field in Norway (see [9]) and
the Bachaquero (see [8]) and Tia Juana (see [7]) fields in Venezuela also illustrate
the importance of the subsidence phenomenon in oil production. The purpose of
this paper is to simulate the mechanical deformations of the porous media in the
presence of a Darcy flow in two space dimensions, taking into account that they
have different time scales.

The behavior of the skeleton is described through a linear elastic equation, whereas we consider two alternative flow models, a monophasic linear case and a two phase nonlinear one. The coupling between these equations is given by the Biot's law which connects the variation of the porosity of the ground to the variation of the mechanical deformations. The main physical unknowns for the linear monophasic flow are the pressure and the porosity; as for the two phase nonlinear flow, the unknown functions are the pressure, saturation and porosity, and the geomechanical problem is written in terms of the displacement. A possible way is to write a complete simulator involving all the flow and mechanical unknowns (see [2]), but this seems to be too expensive. In practice the variation of the displacement is much smaller than the variation of the flow unknowns which means that we are dealing with a multi-scale phenomenon. In the petroleum industry, there already exist different simulators which permit to separately solve the flow problem and the geomechanical one, and since those simulators are very efficient and able to handle complex models, it seems natural to combine them. Typically the flow simulators are based upon a finite volume method whereas one applies a standard finite element method for the discretization of the geomechanical problem. In practice, one has to deal with the fact that geomechanical simulators are more expensive than flow simulators; further one should keep in mind the fact that the inf-sup condition is not satisfied by the mixed discretization method imposed by using the separate simulators. The purpose of this work is to compare two splitting approaches in the case of a two dimensional prototype, extending a previous study performed in a one-dimensional case (see [4]):

* the first one, which has been introduced by engineers, is based upon the computation of the flow unknowns and the displacement by means of a fixed point method; more precisely one makes use of an iterative procedure and we show that in the monophasic linear case this amounts to a Gauss-Seidel type iteration method;

* the new splitting approach which we propose is based upon a preconditioned conjugate gradient method.

In fact, we use a multi-scale coupling algorithm where we study the impact of ratio of the time steps of the mechanical and fluid flow computations on the precision of the numerical results.

The new algorithm turns out to be more robust than the first one; in particular it works fine even though the Inf-Sup condition is not satisfied in the mixed discretization method and it is faster and less expensive. Note that the precision of the computations is satisfactory as soon as the ratio of the two time steps is not too large. These approaches are applied in both the cases where the mechanical deformations of the ground are coupled to a monophasic linear flow and to a two phase nonlinear flow.

We present the two problems in Section 2. We will refer to the model with the monophasic linear flow as to the linear model and we will call nonlinear model the system involving the nonlinear two phase flow. In Section 3 we discretize the two models and present the complete discretized problems. In Section 4 we show that the Inf-Sup condition is not satisfied because of the mixed discretization which is imposed by using separate simulators. We introduce in Section 5 the first splitting approach which we refer to as the multi-scale staggered algorithm and we show that
in the linear case it amounts to a Gauss-Seidel method. In Section 6 we present our new approach based upon the preconditioned conjugate gradient method in both the linear and the nonlinear models. Finally, we present in Section 7 the numerical results: the two approaches are compared in terms of robustness and convergence rates.

## 2. Mathematical model

In this section, we assume small deformations of the medium and write the equations for the flow and the mechanical deformations of the ground.

We denote by $\Omega$ the computational domain of the geomechanics and by $\omega \subset \Omega$ the computational domain for the fluid flow so that the reservoir is completely impermeable in $\Omega \setminus \omega$ (see Figure 1). Furthermore we assume that $\Omega$ is much larger than $\omega$ 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 $\omega$. Moreover we suppose that the solid matrix is incompressible. The coupling between flow and mechanics is given by Biot's law which connects the variation of the porosity of the medium to its mechanical deformation. We also suppose that the medium is submitted to the fluid pore pressure. For the sake of simplicity we do not take into account the gravity force and the capillary effects in the case of the nonlinear two-phase flow.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. The reservoir and the geomechanical domains.</figcaption></figure>

We will study two models:

1. the first one describes the mechanical deformations of the porous media in the presence of a monophasic linear flow and is referred to as **the linear model**;

2. the second one, which describes the same mechanical deformations as the linear model, is coupled with a two-phase compressible flow; we refer to it as **the nonlinear model**.

In both cases the coupling between the flow and the geomechanics is given by Biot's law. The full system consists in the equations governing the time evolution of the fluid pore pressure $p = p(\mathbf{x}, t)$, the porosity of the medium $\phi = \phi(\mathbf{x}, t)$, and in the nonlinear model the fluid saturation $S = S(\mathbf{x}, t)$ and the displacement field $\mathbf{u} = \mathbf{u}(\mathbf{x}, t)$ of the porous medium.
**2.1. Monophasic linear flow.** A monophasic linear flow of compressible fluid is modelled by an evolution equation such that the unknown functions are the fluid pressure $p$ and the porosity of the medium $\phi$. The equation of mass conservation is uniform parabolic and is given by

$$ (1) \qquad \phi_0 \rho_o^0 c_o \frac{\partial p}{\partial t} - \frac{\rho_o^0}{\eta_o} \mathrm{div} \left( \kappa \nabla p \right) + \rho_o^0 \frac{\partial \phi}{\partial t} = Q \quad \text{in } \omega \times (0, T], $$

with the boundary and initial conditions

$$ \left\{ \begin{array}{ll} \frac{\kappa}{\eta_o} \nabla p \cdot \mathbf{n} & = 0 & \text{on } \gamma_N \times (0, T], \\ p & = p_D & \text{on } \partial\omega, \\ p(\cdot, 0) & = p_0 & \text{in } \omega. \end{array} \right. $$

**2.2. Two phase nonlinear flow.** The two phase compressible flow, where we neglect the capillary effects so that the water and oil pressures are equal, is given by the system of equations

$$ (2) \left\{ \begin{array}{ll} \frac{\partial}{\partial t} (\rho_w(p) S \phi) - \mathrm{div} \left( \frac{\rho_w(p) kr_{w,o}(S)}{\eta_w} \kappa \nabla p \right) & = 0 \quad \text{in } \omega \times (0, T], \\ \frac{\partial}{\partial t} (\rho_o(p)(1-S)\phi) - \mathrm{div} \left( \frac{\rho_o(p) kr_{o,w}(S)}{\eta_o} \kappa \nabla p \right) & = 0 \quad \text{in } \omega \times (0, T]. \end{array} \right. $$

with the boundary and initial conditions

$$ \left\{ \begin{array}{ll} p & = p_D \quad \text{on } \partial\omega \times (0, T], \\ \frac{\kappa}{\eta_i} \nabla p \cdot \mathbf{n} & = 0 \quad \text{on } \gamma_N \times (0, T] \text{ for } i = o \text{ and } w, \\ p(\cdot, 0) & = p_0 \quad \text{in } \omega, \\ S(\cdot, 0) & = S_0 \quad \text{in } \omega. \end{array} \right. $$

The indices $w$ and $o$ denote respectively the water and oil phases.

Next we define all the quantities introduced above:

* $\gamma_N$ a subset of the boundary of $\omega$;

* $S = S(\mathbf{x}, t)$ the saturation of the water phase, and $1-S$ is the saturation of the oil phase;

* $p = p(\mathbf{x}, t)$ the pressure of both phases;

* $\phi = \phi(\mathbf{u})$ the porosity of the porous medium and $\phi_0$ is the initial value of the porosity;

* $\kappa$ the intrinsic permeability;

* $c_i$ the fluid compressibility of the phase $i$ where $i \in \{o, w\}$;

* $\rho_i = \rho_i(p)$ the density of the phase $i$, which linearly depends on the pressure, such that it is given by $\rho_i(p) = \rho_i^0 (1 + c_i p)$, where $\rho_i^0$ denotes the initial density of the phase $i$;

* $kr_{i,j} = kr_{i,j}(S)$ the relative permeability of the phase $i$ in presence of the phase $j$;

* $\eta_i$ the dynamical viscosity of the phase $i$;

* $Q$ a term source.

**2.3. The geomechanical model.** We suppose that the rock has an elastic behavior so that the geomechanics is modeled by the elastic linear equation. We obtain the elliptic equation

$$ (3) \qquad \left\{ \begin{array}{ll} \mathrm{div}\sigma(\mathbf{u}) - b(\mathbf{x})\nabla p = 0 & \text{in } \Omega \times (0, T], \\ \sigma(\mathbf{u}) = 2\mu\epsilon(\mathbf{u}) + \lambda(\mathrm{div}\mathbf{u})\mathbf{I} & \text{in } \Omega \times (0, T], \end{array} \right. $$
with the boundary conditions

$$ \begin{cases} \mathbf{u} = \mathbf{0} & \text{on } \Gamma_D, \\ \sigma \cdot \mathbf{n} = \mathbf{q} & \text{on } \Gamma_N, \end{cases} $$

where $\Gamma_D \cap \Gamma_N = \emptyset$, $\Gamma_D \cup \Gamma_N = \partial\Omega$ and $\partial\Omega$ is the boundary of the domain $\Omega$. The functions $\sigma$ and $\epsilon$ are respectively the stress and the strain tensors of the porous medium with

$$ \epsilon(\mathbf{u}) = \frac{1}{2}(\nabla \mathbf{u} + \nabla^t \mathbf{u}). $$

$\lambda$ and $\mu$ are the Lamé constants in drained conditions, $\mathbf{q}$ is a given function, and the function $b$ is defined by

$$ b(\mathbf{x}) = \begin{cases} 0 & \text{in } \Omega \setminus \omega, \\ 1 & \text{in } \omega. \end{cases} $$

**2.4. Biot's law.** The geomechanics and the flow model are coupled by Biot's law which connects the variation of the porosity of the porous medium to its mechanical deformations. The corresponding equation is given by

$$ (4) \qquad \frac{\partial \phi}{\partial t} = b(\mathbf{x}) \operatorname{div} \frac{\partial \mathbf{u}}{\partial t} \quad \text{in } \Omega \times (0, T]. $$

**3. The discretization**

In order to approximate the fluxes and to insure the mass conservation, we apply the finite volume method, which we refer to as FV, for the discretization of the reservoir problem, and the finite element method, which we refer to as FE, to discretize the geomechanical problem. The choice of the mesh points, which is standard for such problems, is indicated in the Figure 2.

<figure><img src="image_3.png" /></figure>

○ The computation of the displacement on the vertexes of the triangles.

+ The computation of the reservoir unknowns in the centers of the control volumes.

<caption>FIGURE 2. The finite volume and the finite element mesh points.</caption>
In practice the variations of the displacement are slower than that of the unknown functions of the flow problem; in other words we have to deal with a multiscale phenomenon. The displacement **u** is computed at the times {$T^k$}$_{k=1,...,p}$ such that

$$T^0 = 0 < T^1, \dots < T^k < \dots < T^p = T,$$

where {$T^k$}$_{k=1,...,p}$ is a partition of the interval $[0, T]$. We also suppose that $T^{k+1} - T^k = \Delta T$ for $k = 0, ..., p-1$ and refer to $\Delta T$ as the period of the geomechanical computations. The unknown functions $p, S$ and $\phi$ are computed at the times $t^{1,k+1}, ..., t^{n,k+1}, ..., t^{q,k+1}$, where $t^{0,k+1} = T^k < t^{1,k+1} < ... < t^{n,k+1} < ... < t^{q,k+1} = T^{k+1}$, for all $k \in \{0, ..., p-1\}$. We denote by {$T^k$}$_{k=0,...,p}$ the meeting times of the geomechanical and the reservoir simulations. We set $t^{n+1,k+1} - t^{n,k+1} = \Delta t, n = 0, ..., q-1, k = 0, ..., p-1$. The parameters $\Delta t$ and $\Delta T$ can vary in the course of the computations; for the sake of simplicity we consider them as fixed in our presentation.

**3.1. Finite volume scheme for the flow problem.** Let $\tau_h$ be the set of the control volumes; we denote by $m(K)$ the measure of the control volume $K$ and by $N(K)$ the set of its neighbors. We define the transitivity between two neighboring volume elements $K$ and $L$ by $T_{K,L} = \frac{m(e_{K,L})}{d(\mathbf{x}_K, \mathbf{x}_L)}$ where $m(e_{K,L})$ is the measure of their common interface; we remark that there is a point $\mathbf{x}_K$ in each volume element $K$ and that the segment $[\mathbf{x}_K, \mathbf{x}_L]$ is orthogonal to the interface $e_{K,L}$. Further we denote by $d(\mathbf{x}_K, \mathbf{x}_L)$ the distance between the points $\mathbf{x}_K$ and $\mathbf{x}_L$. Finally we define the approximations of the pressure, the saturation and the porosity by

$$ (5) \qquad \begin{aligned} p_{h,\Delta t}(\mathbf{x}, t) &= p_K^{n+1,k+1} & \text{if } (\mathbf{x}, t) \in K \times (t^{n,k+1}, t^{n+1,k+1}], \\ S_{h,\Delta t}(\mathbf{x}, t) &= S_K^{n+1,k+1} & \text{if } (\mathbf{x}, t) \in K \times (t^{n,k+1}, t^{n+1,k+1}], \\ \phi_{h,\Delta t}(\mathbf{x}, t) &= \phi_K^{n+1,k+1} & \text{if } (\mathbf{x}, t) \in K \times (t^{n,k+1}, t^{n+1,k+1}], \end{aligned} $$

for all $n = 0, ..., q-1$ and for all $k = 0, ..., p-1$. In a similar way, we define the quantities $\rho_{i,K}^n$, ($i = o, w$).

**3.1.1. Discretization of the monophasic linear flow.** In order to obtain the finite volume scheme one formally integrates the equation (1) on the cell $K \times (t^{n,k+1}, t^{n+1,k+1}]$; this yields

$$ (6) \qquad \begin{aligned} & m(K) \phi_0 \rho_o^0 c_o \frac{p_K^{n+1,k+1} - p_K^{n,k+1}}{\Delta t} - \frac{\rho_o^0}{\eta_o} \kappa \sum_{L \in N(K)} T_{K,L} (p_L^{n+1,k+1} - p_K^{n+1,k+1}) + \\ & \qquad + m(K) \rho_o^0 \frac{\phi_K^{n+1,k+1} - \phi_K^{n,k+1}}{\Delta t} = Q_K^{n+1,k+1}, \end{aligned} $$

where $Q_K^{n+1,k+1}$ is the approximation of $\int_K \int_{t^{n,k+1}}^{t^{n+1,k+1}} Q \, dx \, dt$. We suppose that $\phi_K^{n+1,k+1}$ is obtained by linear interpolation, namely

$$ \phi_K^{n+1,k+1} = \phi_K^{0,k+1} + \Delta\phi_K^{k+1} \frac{(n+1)\Delta t}{\Delta T}, $$

for all $n = 0, ..., q-1, k = 0, ..., p-1$ and for all $K \in \tau_h$ where the variation of the porosity $\Delta\phi_K^{k+1}$ on $K$ over a period $T^{k+1} - T^k = \Delta T$, which is given by

$$ \Delta\phi_K^{k+1} = \phi_K^{q,k+1} - \phi_K^{0,k+1}, $$
is computed by a discrete approximation of Biot's law (see (14) below). Further we introduce the vector notations

$$ \mathbf{p}^{n,k+1} = (p_K^{n,k+1})_{K \in \tau_h}, \mathbf{S}^{n,k+1} = (S_K^{n,k+1})_{K \in \tau_h} \text{ and } \Delta\Phi^{k+1} = (\Delta\phi_K^{k+1})_{K \in \tau_h}. $$

Let $N_v$ denote the number of the volume elements. The previous notations permit to rewrite (6) as

$$ (7) \qquad \mathcal{R}_h \mathbf{p}^{n+1,k+1} + \mathbf{D}_h \Delta\Phi^{k+1} = \phi_0 c_o \mathbf{D}_h \mathbf{p}^{n,k+1} + \Delta t \mathbf{Q}_h^{n+1,k+1}, $$

where

* $\mathcal{R}_h = \phi_0 c_o \mathbf{D}_h + \Delta t \mathbf{R}_h$;

* The symmetric positive definite matrix $\mathbf{R}_h$ comes from the discretization of the gradient operator;

* $\mathbf{D}_h$ is a diagonal matrix such that its diagonal elements are given by $d_{KK} = m(K)$ for $K = 1, \dots, N_v$;

* The right hand side $\Delta t \mathbf{Q}_h^{n+1}$ of the reservoir model (6) is such that

$$ \mathbf{Q}_h^{n+1} = \left( \frac{1}{\rho_o^0} \mathbf{Q}_K^{n+1,k+1} \right)_K. $$

We successively substitute the pressures $\mathbf{p}^{n,k+1}$ for $n = 1, \dots, q-1$ in (7) to obtain the equation for $\mathbf{p}^{q,k+1}$

$$ \begin{aligned} & \mathcal{R}_h^q \mathbf{p}^{q,k+1} + \left( \sum_{i=0}^{q-1} \mathcal{R}_h^i (c_o \phi_0 \mathbf{D}_h)^{q-1-i} \right) \frac{\mathbf{D}_h}{q} \Delta\Phi^{k+1} = \\ & = \Delta t \sum_{i=0}^{q-1} (c_o \phi_0 \mathbf{D}_h)^{q-1-i} \mathcal{R}_h^i \mathbf{Q}_h^{i,k+1} + (c_o \phi_0 \mathbf{D}_h)^q \mathbf{p}^{0,k+1}. \end{aligned} $$

**3.1.2. Discretization of the nonlinear two-phase flow.** We integrate (2) on the space time element $K \times (t^{n,k+1}, t^{n+1,k+1}]$ to obtain

$$ (8) \qquad \left\{ \begin{array}{l} m(K) \frac{\phi_K^{n+1,k+1} \rho_{w,K}^{n+1,k+1} S_K^{n+1,k+1} - \phi_K^{n,k+1} \rho_{w,K}^{n,k+1} S_K^{n,k+1}}{\Delta t} - \sum_{L \in \mathcal{N}(K)} \Lambda_{w,K,L}^{n+1,k+1} = 0, \\ \\ m(K) \frac{\phi_K^{n+1,k+1} \rho_{o,K}^{n+1,k+1} (1 - S_K^{n+1,k+1}) - \phi_K^{n,k+1} \rho_{o,K}^{n,k+1} (1 - S_K^{n,k+1})}{\Delta t} - \sum_{L \in \mathcal{N}(K)} \Lambda_{o,K,L}^{n+1,k+1} = 0, \end{array} \right. $$

where the terms of the form $\Lambda_{w,K,L}^{n+1,k+1}$ and $\Lambda_{o,K,L}^{n+1,k+1}$ correspond to an “upstream weighting scheme” (see [5]):

$$ \Lambda_{w,K,L}^{n+1,k+1} = \begin{cases} \kappa T_{K,L} \frac{\rho_{w,K}^{n,k+1} kr_{w,o}(S_K^{n,k+1})}{\eta_w} & (p_L^{n+1,k+1} - p_K^{n+1,k+1}) \text{ if } p_L^{n+1,k+1} > p_K^{n+1,k+1}, \\ \kappa T_{K,L} \frac{\rho_{w,L}^{n,k+1} kr_{w,o}(S_L^{n,k+1})}{\eta_w} & (p_L^{n+1,k+1} - p_K^{n+1,k+1}) \text{ otherwise,} \end{cases} $$

$$ \Lambda_{o,K,L}^{n+1,k+1} = \begin{cases} \kappa T_{K,L} \frac{\rho_{o,K}^{n,k+1} kr_{o,w}(S_K^{n,k+1})}{\eta_o} & (p_L^{n+1,k+1} - p_K^{n+1,k+1}) \text{ if } p_L^{n+1,k+1} > p_K^{n+1,k+1}, \\ \kappa T_{K,L} \frac{\rho_{o,L}^{n,k+1} kr_{o,w}(S_L^{n,k+1})}{\eta_o} & (p_L^{n+1,k+1} - p_K^{n+1,k+1}) \text{ otherwise.} \end{cases} $$
Next we rewrite (8) in the more compact form

$$ (9) \qquad \mathcal{R}_h(\mathbf{D}_h \Delta\Phi^{k+1}, \mathbf{S}^{n,k+1}, \mathbf{p}^{n,k+1}; \mathbf{S}^{n+1,k+1}, \mathbf{p}^{n+1,k+1}) = 0, $$

where the operator $\mathcal{R}_h$ is nonlinear. Given a variation of the porosity $\Delta\Phi^{k+1}$ over the period $T^{k+1} - T^k = \Delta T$ and a pair $(\mathbf{S}^{n,k+1}, \mathbf{p}^{n,k+1})$, one computes the solution pair $(\mathbf{S}^{n+1,k+1}, \mathbf{p}^{n+1,k+1})$ of Problem (9) by means of the Newton method. We use the notation $\mathcal{R}_h^{-1}$ in the expression of the solution in order to have the same notations as in the linear model; for all $n = 0, \dots, q-1$, we write

$$ (10) \qquad (\mathbf{S}^{n+1,k+1}, \mathbf{p}^{n+1,k+1}) = \mathcal{R}_h^{-1}(\mathbf{D}_h \Delta\Phi^{k+1}, \mathbf{S}^{n,k+1}, \mathbf{p}^{n,k+1}). $$

In order to compute the last value of the pressure $\mathbf{p}^{q,k+1}$ which in turn will permit to calculate the displacement at this time with the geomechanical module, we need to perform $q$ computations by means of the Newton method. We use the notation $\mathbf{p}^{q,k+1} = \mathcal{R}_h^{-q}(\mathbf{D}_h \Delta\Phi^{k+1})$ for the successive results of the simulations:

$$ (11) \qquad \left\{ \begin{array}{l} \text{For } n = 0, \dots, q-1 \\ (\mathbf{S}^{n+1,k+1}, \mathbf{p}^{n+1,k+1}) = \mathcal{R}_h^{-1}(\mathbf{D}_h \Delta\Phi^{k+1}, \mathbf{S}^{n,k+1}, \mathbf{p}^{n,k+1}). \end{array} \right. $$

**3.2. Finite element computation of the geomechanical equilibrium.** We denote by $\mathbf{u}_{h,\Delta T}(\cdot, T^k)$, $k = 0, \dots, p$, the piece-wise linear approximation of the function $\mathbf{u}(\cdot, T^k)$, which one obtains by solving the discrete linear system

$$ (12) \int_{\Omega} \sigma(\mathbf{u}_{h,\Delta T}(\cdot, T^k)) : \epsilon(\mathbf{v}_h) \, d\mathbf{x} - \int_{\omega} (\operatorname{div} \mathbf{v}_h) p_{h,\Delta t}(\cdot, T^k) \, d\mathbf{x} = \langle \mathbf{f}_{h,\Delta T}, \mathbf{v}_h \rangle, $$

where

$$ \langle \mathbf{f}_{h,\Delta T}, \mathbf{v}_h \rangle = \int_{\Gamma_N} (\sigma(\mathbf{u}_{h,\Delta T}(\cdot, T^k))\mathbf{n}) \cdot \mathbf{v}_h \, ds - \int_{\partial\omega} p_D \mathbf{n} \cdot \mathbf{v}_h \, ds, $$

for all $k = 0, \dots, p$ and for suitable test functions $\mathbf{v}_h$; here

$$ \sigma(\mathbf{u}) : \epsilon(\mathbf{u}) = \sum_{i,j=1}^{3} \sigma_{ij}(\mathbf{u}) \epsilon_{ij}(\mathbf{u}), $$

with $\sigma(\mathbf{u}) = (\sigma_{ij})_{i,j}$ and $\epsilon(\mathbf{u}) = (\epsilon_{ij})_{i,j}$. Further the components of $\mathbf{u}^k$ correspond to the values of the vector function $\mathbf{u}_{h,\Delta T}(\cdot, T^k)$ at the vertices of the triangularization. Using these notations and in view of (12) we obtain the matrix expression

$$ (13) \qquad \mathbf{G}_h \mathbf{u}^{k+1} - \mathbf{B}_h \mathbf{p}^{q,n+1} = \mathbf{F}_h^{k+1}, $$

for all $k = 1, \dots, p-1$, where

* $\mathbf{G}_h$ is the symmetric positive definite stiffness matrix, and

* The sparse matrix $\mathbf{B}_h = (b_{mn})_{m=1,\dots,2N_s, n=1,\dots,N_v}$ corresponds to the coupling between flow and mechanical deformations, $N_s$ is the number of the vertices of the triangularization, and $\mathbf{B}_h^T$ is the transpose matrix of $\mathbf{B}_h$.

**3.3. Discretization of Biot's law.** Next let us express the relation between the discrete displacement and the discrete porosity, namely

$$ (14) \qquad m(K) \Delta\phi_K^{k+1} = \int_K \operatorname{div}\{\mathbf{u}_{h,\Delta T}(\mathbf{x}, T^{k+1}) - \mathbf{u}_{h,\Delta T}(\mathbf{x}, T^k)\} \, d\mathbf{x}, $$

for all $k = 0, \dots, p-1$, and for all $K \in \tau_h$. With the vector notations introduced above, (14) can be rewritten as

$$ (15) \qquad \mathbf{D}_h \Delta\Phi^{k+1} = \mathbf{B}_h^T (\mathbf{u}^{k+1} - \mathbf{u}^k). $$
**3.4. The complete discrete problem.** Finally we state the complete discrete problem for both the linear and nonlinear models:

(16)

(i) in the case of the linear flow we have that

$$
\begin{align*}
\mathcal{R}_h^q \mathbf{p}^{q,k+1} &= -\frac{1}{q} \left( \sum_{i=0}^{q-1} \mathcal{R}_h^i (c_o \phi_0 \mathbf{D}_h)^{q-1-i} \right) \mathbf{D}_h \Delta\Phi^{k+1} + \\
&\quad + (c_o \phi_0 \mathbf{D}_h)^q \mathbf{p}^{0,k+1} + \Delta t \sum_{i=0}^{q-1} (c_o \phi_0 \mathbf{D}_h)^{q-1-i} \mathcal{R}_h^i \mathbf{Q}_h^{i,k+1},
\end{align*}
$$

(ii) whereas in the case of the nonlinear flow we have

$$
\mathbf{p}^{q,k+1} = \mathcal{R}_h^q (\mathbf{D}_h \Delta\Phi^{k+1}),
$$

together with the system

$$
(17) \qquad \left\{
\begin{array}{l}
\mathbf{G}_h \mathbf{u}^{k+1} - \mathbf{B}_h \mathbf{p}^{q,k+1} = \mathbf{F}_h^{k+1}, \\
\mathbf{D}_h \Delta\Phi^{k+1} = \mathbf{B}_h^T (\mathbf{u}^{k+1} - \mathbf{u}^k),
\end{array}
\right.
$$

for all $k = 1, \dots, p-1$.

**4. Inf-Sup Condition**

In this section we show that the mixed method used in the discretization of
this coupled problem (Finite Elements-Finite Volumes) is not stable for the linear
model as well as in the nonlinear model. For the sake of simplicity, we consider
below the case that $q = 1$ which means that $\Delta t = \Delta T$. When the fluid has a small
$c_o$ compressibility and the time step is small too, the linear system (16-17) is close
to

$$
(18) \quad \begin{pmatrix} \mathbf{G}_h & -\mathbf{B}_h \\ \mathbf{B}_h^T & 0 \end{pmatrix} \begin{pmatrix} \mathbf{u}^{k+1} \\ \mathbf{p}^{k+1} \end{pmatrix} = \begin{pmatrix} \mathbf{F}_h^{k+1} \\ \mathbf{B}_h^T \mathbf{u}^k \end{pmatrix}.
$$

In view of the theorem 4.3, p. 127 of [1] we have that the mixed problem (18)
admits a solution if and only if the **LBB condition** (Ladyshenskaja-Babushka-
Brezzi condition) is satisfied. Here this condition reduces to where **V**<sub>h</sub> is the space
of the piecewise linear functions. The condition (19) implies that the application
associated to the matrix **B**<sub>h</sub> is surjective. The Schur complement S<sub>h</sub> associated to
the system (18) is given by S<sub>h</sub> = **B**<sub>h</sub><sup>T</sup> **G**<sub>h</sub><sup>-1</sup> **B**<sub>h</sub>. In this case, as we show in figure 3,
this operator is only positive semi-definite.

We remark that **the inf-sup condition** implies that the operator $S_h$ is definite.
Indeed, let $\{\lambda_i\}_{i=1,\dots,n}$ be the eigenvalues of the Schur complement such that $\lambda_1 \le$
$\lambda_2 \le \dots \le \lambda_n$ and let $\{\mathbf{v}_i\}_{i=1,\dots,n}$ be the associated eigenvectors. We denote
respectively by $\eta_{\min}$ and $\eta_{\max}$ the smallest and largest eigenvalues of $\mathbf{G}_h^{-1}$. We
have for all $i=1,\dots,n$,

$$
(19) \quad \exists \beta > 0 \text{ such that } \|\mathbf{B}_h \eta_h\| \ge \beta \|\eta_h\| \quad \text{for all } \eta_h \in \mathbf{V}_h,
$$

$$
(20) \quad \begin{array}{l} \eta_{\min} \|\mathbf{B}_h \mathbf{v}_i\|^2 \le (S_h \mathbf{v}_i, \mathbf{v}_i) \le \eta_{\max} \|\mathbf{B}_h \mathbf{v}_i\|^2, \\ \eta_{\min} \|\mathbf{B}_h \mathbf{v}_i\|^2 \le \lambda_i \|\mathbf{v}_i\|^2 \le \eta_{\max} \|\mathbf{B}_h \mathbf{v}_i\|^2. \end{array}
$$

Therefore if the condition (19) is satisfied, it follows that

$$
(21) \qquad \eta_{\min} \beta^2 \le \lambda_1,
$$

which implies that the operator $S_h$ is positive definite.
<figure><img src="image_4.png" /><figcaption>FIGURE 3. The dependence of the largest eigenvalue of $S_h^{-1}$ depending on the mesh size.</figcaption></figure>

Moreover if the fluid has a too small compressibility, the time step $\Delta t$ should also have a minimum value, which is rather constraining. However the constraint would disappear with a different choice of the finite element method in the mechanical simulator.

## 5. Multi-step staggered algorithm

We show in Figure 4 below a graphical representation of the algorithm which is valid in both the cases of the linear and the nonlinear models; it clearly expresses the fact that the mechanical deformations happen in a much slower time scale than the porous medium flow.

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Sketch of the iterations between reservoir and mechanical computations on the time period $\Delta T$.</figcaption></figure>
**5.1. The case of the linear model.** We solve the set of equations (16-17), where $c_r$ is a relaxation parameter, by means of a fixed point method. The corresponding algorithm, where the iteration parameter is denoted by $l$, is given by (22)

**Initialization of the period:**

$\mathbf{p}_0^{0,k+1} = \mathbf{p}_1^{0,k+1} = \mathbf{p}^{q,k}$, are given by the previous period,

$$ \Delta\Phi_0^{k+1} = 0, $$

**Iterations on the period $\Delta T = q\Delta t$:**

For $l = 0, \dots, \dots$, until convergence,

Reservoir simulations:

For $n = 0, \dots, q-1$,

$$ \mathcal{R}_h \mathbf{p}_{l+1}^{n+1,k+1} + c_r \phi_0 \mathbf{D}_h (\mathbf{p}_{l+1}^{n+1,k+1} - \mathbf{p}_{l+1}^{n,k+1}) + \frac{\Delta t}{\Delta T} \mathbf{D}_h \Delta\Phi_l^{k+1} = \phi_0 c_o \mathbf{D}_h \mathbf{p}_{l+1}^{n,k+1} + \Delta t \mathbf{Q}_h^{n+1,k+1}. $$

Geomechanics computation:

$$ \mathbf{G}_h \mathbf{u}_{l+1}^{k+1} + \mathbf{B}_h \mathbf{p}_{l+1}^{q,k+1} = \mathbf{F}_h^k, $$

$$ \mathbf{D}_h \Delta\Phi_{l+1}^{k+1} = \mathbf{B}_h^T (\mathbf{u}_{l+1}^{k+1} - \mathbf{u}^k) - c_r c_o \phi_0 \mathbf{D}_h (\mathbf{p}_{l+1}^{q,k+1} - \mathbf{p}^{0,k}). $$

The parameter $c_r$, which is often used by the engineers in petroleum industry, is refered to as the rock compressibility and the algorithm (22) can be seen as a Gauss-Seidel method with a relaxation parameter $\gamma = (1 + \frac{c_r}{c_o})^{-1}$.

**Remark 1.** The staggered algorithm (22) when $c_r = 0$ and $q = 1$, which is a classical splitting method, can be seen as a Block Gauss-Seidel method for solving the pressure/displacement coupled system. Indeed this method consists in solving system $\mathbf{M}\mathbf{x}_{l+1} = \mathbf{N}\mathbf{x}_l + \mathbf{b}$ for each iteration $l$, where $\mathbf{N}$ is upper triangular matrix of $\mathbf{A}$ and $\mathbf{M} = \mathbf{A} - \mathbf{N}$. In our case

$$ \mathbf{A} = \begin{pmatrix} \mathcal{R}_h & \mathbf{B}_h^T \\ -\mathbf{B}_h & \mathbf{G}_h \end{pmatrix}, \mathbf{M} = \begin{pmatrix} \mathcal{R}_h & 0 \\ -\mathbf{B}_h & \mathbf{G}_h \end{pmatrix}, \mathbf{N} = \begin{pmatrix} 0 & -\mathbf{B}_h^T \\ 0 & 0 \end{pmatrix}, $$

and

$$ \mathbf{b} = \begin{pmatrix} \mathbf{D}_h \mathbf{p}^k + \mathbf{B}_h^T \mathbf{u}^k + \Delta t \mathbf{Q}^{k+1} \\ \mathbf{F}_h^{k+1} \end{pmatrix}, \quad \mathbf{x}_{l+1} = \begin{pmatrix} \mathbf{p}_{l+1}^{k+1} \\ \mathbf{u}_{l+1}^{k+1} \end{pmatrix}, $$

This method converges if and only if the spectre $\rho(\mathbf{M}^{-1}\mathbf{N}) < 1$.

**5.2. The case of the nonlinear model.** As in the case of the linear model, we use a fixed point method; we obtain the following algorithm

**Period initialisation:**

$\mathbf{S}_0^{0,k+1} = \mathbf{S}^{0,k+1} = \mathbf{S}^{q,k}; \mathbf{p}_0^{0,k+1} = \mathbf{p}^{0,k+1} = \mathbf{p}^{q,k}$; are given by the previous period

$$ \Delta\Phi_0^{k+1} = 0. $$

**Iterations on the period $\Delta T = q\Delta t$:**

For $l = 0, \dots, \dots$, until convergence,

Reservoir simulation:

For $n = 0, \dots, q-1$,

Newton method :

$$ \begin{aligned} & \mathcal{R}_h(\Delta\Phi^{k+1} + c_r \frac{\Delta T}{(n+1)\Delta t}(\mathbf{p}_{l+1}^{n+1,k+1} - \mathbf{p}^{0,k+1}), \mathbf{S}^{n,k+1}, \mathbf{p}^{n,k+1}; \mathbf{S}^{n+1,k+1}, \mathbf{p}^{n+1,k+1}) \\ & = 0. \end{aligned} $$

Geomechanical computations:

$$ \mathbf{G}_h \mathbf{u}_{l+1}^{k+1} - \mathbf{B}_h \mathbf{p}_{l+1}^{q,k+1} = \mathbf{F}_h^{k+1}, $$

$$ \mathbf{D}_h \Delta\Phi_{l+1}^{k+1} = \mathbf{B}_h^T (\mathbf{u}_{l+1}^{k+1} - \mathbf{u}^k) - c_r c_o \phi_0 \mathbf{D}_h (\mathbf{p}_{l+1}^{q,k+1} - \mathbf{p}^{0,k}). $$
## 6. Preconditioned conjugate gradient method

**6.1. Linear model.** To apply the preconditioned conjugate gradient method to the linear model, we have two alternatives to construct the Schur complement.

1. Extract the expression of the displacement $\mathbf{u}^{k+1}$ from the first equation of (17) and substitute it into the second one, and then replace the result in the first one of (16).

2. Substitute the second equation of (17) in the first one of (16), extract than the expression of the pressure $\mathbf{p}^{q,k+1}$ from the result and substitute it into the first equation of (17).

It seems difficult to apply the preconditioned conjugate gradient method in the first approach because we can not define the matrix vector product, which is immediate in the second case. Hence, we obtain using the second approach

$$ (23) \left\{ \begin{array}{l} \left( \mathbf{G}_h + \mathbf{B}_h \mathcal{R}_h^{-q} \sum_{i=0}^{q-1} \mathcal{R}_h^i (c_o \phi_0 \mathbf{D}_h)^{q-1-i} \frac{\mathbf{B}_h^T}{q} \right) \mathbf{u}^{k+1} = \\ = \mathbf{F}_h^{k+1} + (c_o \phi_0 \mathbf{D}_h)^q \mathbf{p}^{0,k+1} + \mathbf{B}_h \mathcal{R}_h^{-q} \left( \sum_{i=0}^{q-1} (c_o \phi_0 \mathbf{D}_h)^{q-1-i} \mathcal{R}_h^i \mathbf{S}_h^{i+1,k+1} \right), \\ \\ \mathcal{R}_h^q \mathbf{p}^{q,k+1} = - \sum_{i=0}^{q-1} \mathcal{R}_h^i (c_o \phi_0 \mathbf{D}_h)^{q-1-i} \frac{\mathbf{B}_h^T}{q} \mathbf{u}^{k+1} + \\ \qquad \sum_{i=0}^{q-1} (c_o \phi_0 \mathbf{D}_h)^{q-1-i} \mathcal{R}_h^i \mathbf{S}_h^{i+1,k+1} + c_o \phi_0 \mathbf{D}_h)^q \mathbf{p}^{0,k+1}. \end{array} \right. $$

where $\mathbf{S}_h^{i+1,k+1} = \Delta t \mathbf{Q}_h^{i,k+1} + \frac{\mathbf{B}_h^T}{q} \mathbf{u}^k$ for all $i \in \{0, \dots, q-1\}$. We apply the preconditioned conjugate gradient algorithm to the first equation of (23) with $\mathbf{G}_h$ as preconditioning matrix and apply a fixed point algorithm to obtain the solution of the coupled system (23).

**6.2. Nonlinear model.** We now extend the preconditioned conjugate gradient method to the nonlinear case. From the formulas (16)-(17) we obtain

$$ (24) \left\{ \begin{array}{l} \mathbf{G}_h \mathbf{u}_{l+1}^{k+1} - \mathbf{B}_h \mathcal{R}_h^{-q} (\mathbf{B}_h^T (\mathbf{u}_{l+1}^{k+1} - \mathbf{u}^k)) = \mathbf{F}_h^{k+1}, \\ \mathbf{p}^{q,k+1} = \mathcal{R}_h^{-q} (\mathbf{B}_h^T (\mathbf{u}_{l+1}^{k+1} - \mathbf{u}^k)). \end{array} \right. $$

We apply the nonlinear preconditioned conjugate gradient method to (24). The residual at each iteration is given by

$$ (25) \quad \mathbf{G}_h \mathbf{u}_{l+1}^{k+1} - \mathbf{B}_h \mathcal{R}_h^{-q} (\mathbf{B}_h^T (\mathbf{u}_{l+1}^{k+1} - \mathbf{u}^k)) - \mathbf{F}_h^{k+1} = \mathbf{G}_h \mathbf{u}_{l+1}^{k+1} - \mathbf{B}_h \mathbf{p}_{l+1}^{q,k+1} - \mathbf{F}_h^{k+1}. $$

In order to obtain the descent coefficient $\alpha_l$, we have to compute the gradient of the left hand side of (25)

$$ (26) \quad \nabla_x \left( \mathbf{G}_h \mathbf{x} - \mathbf{B}_h \mathcal{R}_h^{-q} (\mathbf{B}_h^T (\mathbf{x} - \mathbf{u}^k)) \right) \cdot \mathbf{d}_l = \mathbf{G}_h \cdot \mathbf{d}_l - \nabla_x \left( \mathbf{B}_h \mathcal{R}_h^{-q} (\mathbf{B}_h^T (\mathbf{x} - \mathbf{u}^k)) \right) \cdot \mathbf{d}_l, $$

and use the following finite difference approximation with $\mathbf{x} = \mathbf{u}_{l+1}^{k+1}$

$$ \nabla_x \left( \mathbf{B}_h \mathcal{R}_h^{-q} (\mathbf{B}_h^T (\mathbf{u}_{l+1}^{k+1} - \mathbf{u}^k)) \right) \cdot \mathbf{d}_l \approx \frac{\mathcal{R}_h^{-q} (\mathbf{B}_h^T (\mathbf{u}_{l+1}^{k+1} + \epsilon \mathbf{d}_l - \mathbf{u}^k)) - \mathcal{R}_h^{-q} (\mathbf{B}_h^T (\mathbf{u}_{l+1}^{k+1} - \mathbf{u}^k))}{\epsilon}, $$
where $\epsilon$ is a small parameter. Let us note by $\tilde{\mathbf{p}}_l^{q,k+1}$ the intermediate pressure given by

$$ \tilde{\mathbf{p}}_l^{q,k+1} = \mathcal{R}_h^{-q}(\mathbf{B}_h^T (\mathbf{u}_{l+1}^{k+1} + \epsilon \mathbf{d}_l - \mathbf{u}^k)) $$

We will speak of one reservoir simulation when we perform in fact a set of $q$ reservoir computations; in order to obtain the descent parameter $\alpha_l$, we have to perform one reservoir simulation.

(27)

**Initialization of the period:**

$\mathbf{u}_0^{k+1}$ is given,

$$ \mathbf{p}_0^{q,k+1} = \mathcal{R}_h^{-q}(\mathbf{B}_h^T (\mathbf{u}_0^{k+1} - \mathbf{u}^k)), \quad \rightarrow \text{Reservoir simulation} $$

$$ \mathbf{r}_0 = \mathbf{F}_h^{k+1} - \mathbf{G}_h \mathbf{u}_0^{k+1} + \mathbf{B}_h \mathbf{p}_0^{q,k+1}, $$

$$ \mathbf{C} \mathbf{d}_0 = \mathbf{r}_0, $$

$$ \mathbf{z}_0 = \mathbf{d}_0. $$

**Iterations on the period $\Delta T = q\Delta t$:**

For $l = 0, \dots, l$, until convergence

$$ \tilde{\mathbf{p}}_{l+1}^{q,k+1} = \mathcal{R}_h^{-q}(\mathbf{B}_h^T (\mathbf{u}_l^{k+1} + \epsilon \mathbf{d}_l - \mathbf{u}^k)), \quad \rightarrow \text{Reservoir simulation} $$

$$ \mathbf{y}_l = \mathbf{G}_h \mathbf{d}_{l+1} - \mathbf{B}_h \frac{\tilde{\mathbf{p}}_{l+1}^{q,k+1} - \mathbf{p}_l^{q,k+1}}{\epsilon}, $$

$$ \alpha_l = \frac{(\mathbf{r}_l, \mathbf{z}_l)}{(\mathbf{y}_l, \mathbf{d}_l)}, $$

$$ \mathbf{u}_{l+1}^{k+1} = \mathbf{u}_l^{k+1} + \alpha_l \mathbf{d}_l, $$

$$ \mathbf{p}_{l+1}^{q,k+1} = \mathcal{R}_h^{-q}(\mathbf{B}_h^T (\mathbf{u}_{l+1}^{k+1} - \mathbf{u}^k)), \quad \rightarrow \text{Reservoir simulation} $$

$$ \mathbf{r}_{l+1} = \mathbf{F}_h^{k+1} - \mathbf{G}_h \mathbf{u}_l^{k+1} + \mathbf{B}_h \mathbf{p}_{l+1}^{q,k+1}, $$

$$ \mathbf{G}_h \mathbf{z}_{l+1} = \mathbf{r}_{l+1}, \quad \rightarrow \text{Geomechanical computation} $$

$$ \beta_{l+1} = \frac{(\mathbf{r}_{l+1}, \mathbf{z}_{l+1})}{(\mathbf{r}_l, \mathbf{z}_l)}, $$

$$ \mathbf{d}_{l+1} = \mathbf{z}_{l+1} + \beta_{l+1} \mathbf{d}_l. $$

**Remark 2.** We remark that the algorithm (27) coincides with the preconditioned conjugate algorithm in the linear case (23).

## 7. Numerical results

First, we show numerically that the schemes which we use are of second order under the stability condition. For the linear case, we construct analytical solutions by choosing a pressure $\mathbf{p}$ and a displacement $\mathbf{u}$, and computing the corresponding right hand sides of (3), (1) and (4). We show that the staggered algorithm and the preconditioned conjugate gradient algorithm both give similar results and we also show that the preconditioned conjugate gradient algorithm is faster than the staggered algorithm and more robust, in particular when the Inf-Sup condition is not satisfied (see Section 4). Then, we study the influence of the ratio of the time steps $q = \Delta T/\Delta t$. In order to try to better simulate the physics, we also use a second analytic solution as a test case, where the displacement $\mathbf{u}$ varies much slower than the pressure $\mathbf{p}$. In particular we check that the errors between the analytical solution and the computational one are similar for the staggered algorithm and the preconditioned conjugate gradient algorithm independantly of the ratio $q$. Moreover it is not surprising that when $q$ is large these errors are very important whereas we observe that for smaller values of $q$, the errors are close to those with the schemes where $q = 1$, that is where $\Delta T = \Delta t$. This validates the multiscale approach. The
problem of the precise choice of $q$ is still open and will be the subject of future work.

We then validate both algorithms in the nonlinear case with physical parameters. We show the influence of the ratio $q$ and of the viscosity rates $\eta_w, \eta_o$ on the behavior of the solution; we note that the results which we obtain agree with the physical previsions.

We conclude that the multiscale algorithm is well adapted and gives accurate results and that the preconditioned conjugate gradient method is less costly and more robust than the staggered algorithm; in particular it permits to perform numerical computations even in the case that the Inf-Sup condition is not fulfilled.

For most numerical results we refere [3].

**7.1. Linear model.** We define the following test case. We set $\Omega = [0, 1]^2$, $\omega \subset \Omega$ with $\omega = \left[\frac{1}{3}, \frac{2}{3}\right]^2$ and

$$p(x, z, t) = 10^4 (x - \frac{1}{3})^2 (x - \frac{2}{3})^2 (z - \frac{1}{3})^2 (z - \frac{2}{3})^2 \cos(t),$$

$$\mathbf{u}(x, z, t) = 100 (x^2 - x)^2 (z^2 - z)^2 \cos(t) \begin{pmatrix} 1 \\ 1 \end{pmatrix},$$

where the time interval of the simulations is (0, 10]. We inject the functions $p$ and $\mathbf{u}$ in (1), (3), and (4) to obtain the right-hand-sides [3] and compare the staggered and the preconditioned conjugate gradient algorithms in the case that $\Delta t = \Delta T$ (i.e. $q = 1$) and $c_o = 1$; in this case one can check numerically that the Inf-Sup condition is satisfied (see section 4). In what follows we denote by $(p_h, \mathbf{u}_h)$ the approximate solution.

We observe in this case that the two approaches have a similar behavior (see Figures 5 and 6). The order of convergence is equal to 2 for the pressure as Figure 5 shows. For the displacement field $u$ (see figure 6) when $\Delta x > 1/45$ we also have second order accuracy, but this is not the case for smaller values of $\Delta x$. These is due to the violation of the stability condition which has the form $\Delta t \le C\Delta x^2$.

In Section 4, we have seen that in the case of small fluid compressibility the Inf-Sup condition imposes a lower bound on $\Delta t$, which is not usual. The Tables 2 and 1 show that when the fluid has a small compressibility ($c_o = 0.5$), the staggered algorithm converges only in cases that the time step $\Delta t$ is not too small. Moreover Table 1 shows the robustness of the preconditioned conjugate gradient algorithm with respect to the staggered algorithm. We observe that the staggered algorithm converges for the first period $T^1 = 1$ and then breaks down whereas the preconditioned conjugate gradient algorithm converges rather fast. Table 1 shows how to make the staggered algorithm converge. The idea is to perform the computations either with a positive rock compressibility $c_r$ or with a larger time step. In any case the robustness of the preconditioned conjugate gradient method can be observed in terms of costs as well as of quality (see Tables 1, 2, and 3).

We now wish to validate both algorithms in the multiscale case where $q > 1$. We set $\Omega = \omega = [0, 1]^2$ and we choose test functions such that the displacement field $\mathbf{u}$
varies less fast than the pressure function p.

$$ 
\begin{aligned}
p(x, z, t) &= 10 (x^2 - x)^2 (z^2 - z)^2 \cos(2t) + 10, \\
\mathbf{u}(x, z, t) &= ((x^2 - x)^2 (z^2 - z)^2 + 10) \cos\left(\frac{t}{5}\right) \begin{pmatrix} 1 \\ 1 \end{pmatrix}.
\end{aligned}
$$

In Figure 7, we remark that the error $p-p_h$ on the pressure for $q$ equal to 1, 10 or 20 is similar for the two algorithms, which validates the multiscale method whereas for $q = 50$ to $q = 200$, the error grows with $q$. The error on the displacement $\mathbf{u}$ (see Figure 8) behaves in a similar way.

<figure><img src="image_6.png" /><figcaption>FIGURE 5. $L^2(\omega)$-norm of the error $p-p_h$ as a function of $\Delta x$, on the left-hand-side with the staggered algorithm and on the right-hand-side with the preconditioned conjugate gradient algorithm with $q = 1$, $c_o = 1$ and $\Delta t = 0.01$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIGURE 6. $L^2(\Omega)$-norm of the error $u - u_h$ as a function of $\Delta x$, on the left-hand-side with the staggered algorithm and on the right-hand-side with the preconditioned conjugate gradient algorithm with $q = 1$, $c_o = 1$ and $\Delta t = 0.01$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7. The relative error $p - p_h$ in $L^2(\omega)$-norm as a function of time, on the left-hand-side with the staggered algorithm and on the right-hand-side with the preconditioned conjugate gradient algorithm for different values of $q$ with $c_o = 1$, $\Delta t = 0.01$ and $\Delta x = 1/40$.</figcaption></figure>
<table><caption>TABLE 1. The number of iterations to compute the solution at the times $T_i$ and the corresponding CPU times for the staggered algorithm and the preconditioned conjugate gradient algorithm with $q = 1$, $c_o = 0.5$ and $\Delta t = 0.01$ ($\times$ indicates the absence of convergence).</caption><thead><tr><th rowspan="2">T<sup>i</sup></th><th colspan="2">Staggered algorithm</th><th colspan="2">P.C.G Algorithm</th></tr><tr><th>Iterations</th><th>CPU</th><th>Iterations</th><th>CPU</th></tr></thead><tbody><tr><td>T<sup>1</sup> = 1</td><td>117</td><td>3.523E+1</td><td>7</td><td>2.079</td></tr><tr><td>T<sup>3</sup> = 3</td><td>×</td><td>×</td><td>10</td><td>4.190</td></tr><tr><td>T<sup>5</sup> = 5</td><td>×</td><td>×</td><td>10</td><td>2.929</td></tr><tr><td>T<sup>7</sup> = 7</td><td>×</td><td>×</td><td>9</td><td>6.649</td></tr><tr><td>T<sup>9</sup> = 9</td><td>×</td><td>×</td><td>9</td><td>2.680</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>FIGURE 8. The relative error $\mathbf{u} - \mathbf{u}_h$ in $L^2(\Omega)$-norm as a function of time, on the left-hand-side with the staggered algorithm and on the right-hand-side with the preconditioned conjugate gradient algorithm for different values of $q$ with $c_o = 1$, $\Delta t = 0.01$ and $\Delta x = 1/40$.</figcaption></figure>

**7.2. Nonlinear model.** In this section, we consider the nonlinear model which describes the mechanical deformations of the porous media in the presence of the two phase flow for water and oil. We take $\Omega = \omega = [0, 1]^2$ as space domains and suppose that water is injected in order to extract the oil. We fix the pressure $p$ at the outcoming phase whereas we fix a value for the water flux $Q_w$ and we prescribe $Q_o = 0$ for the oil flux at the incoming phase. The phases water and oil are assumed to be compressible. The physical parameters are $\eta_o = 5.10^{-3}$ Pa.s, $\eta_w = 1.10^{-3}$ Pa.s, $\kappa = 2.10^{-14}$, $\phi_0 = 0.30$, $\rho_o^0 = 950$ kg.m$^{-3}$, $\rho_w^0 = 1000$ kg.m$^{-3}$,
<table><caption>TABLE 2. The number of iterations and the corresponding CPU times to compute the solution at the times $T_i$ for the staggered algorithm with $q = 1$ and $c_o = 0.5$ on the one hand with a relaxed algorithm, and on the other hand by taking a larger value of $\Delta t$.</caption><thead><tr><th rowspan="2">T<sup>i</sup></th><th colspan="2">Relaxed by c<sub>r</sub> = 0.2, Δt = 0.01</th><th colspan="2">c<sub>r</sub> = 0 and Δt = 0.1</th></tr><tr><th>Iterations</th><th>CPU</th><th>Iterations</th><th>CPU</th></tr></thead><tbody><tr><td>T<sup>1</sup> = 1</td><td>13</td><td>4.176</td><td>23</td><td>6.740</td></tr><tr><td>T<sup>3</sup> = 3</td><td>12</td><td>3.872</td><td>25</td><td>7.300</td></tr><tr><td>T<sup>5</sup> = 5</td><td>13</td><td>4.196</td><td>24</td><td>7.160</td></tr><tr><td>T<sup>7</sup> = 7</td><td>13</td><td>4.185</td><td>24</td><td>6.920</td></tr><tr><td>T<sup>9</sup> = 9</td><td>13</td><td>4.208</td><td>26</td><td>7.650</td></tr></tbody></table>

<table><caption>TABLE 3. Comparison of the efficiency of both algorithms at different times with $q = 1$ in a case where the Inf-Sup condition is satisfied.</caption><thead><tr><th rowspan="2">T<sup>i</sup></th><th colspan="2">Staggered algorithm q = 1</th><th colspan="2">P. C. G. algorithm q = 1</th></tr><tr><th>Iterations</th><th>CPU</th><th>Iterations</th><th>CPU</th></tr></thead><tbody><tr><td>T<sup>1</sup> = 1</td><td>11</td><td>3.538</td><td>5</td><td>1.993</td></tr><tr><td>T<sup>3</sup> = 3</td><td>12</td><td>3.871</td><td>9</td><td>3.318</td></tr><tr><td>T<sup>5</sup> = 5</td><td>12</td><td>3.895</td><td>9</td><td>3.317</td></tr><tr><td>T<sup>7</sup> = 7</td><td>13</td><td>4.228</td><td>7</td><td>2.660</td></tr><tr><td>T<sup>9</sup> = 9</td><td>13</td><td>4.212</td><td>8</td><td>2.971</td></tr></tbody></table>

<table><caption>TABLE 4. Comparison of the efficiency of both algorithms at different times with $q = 10$ in a case that the Inf-Sup condition is satisfied.</caption><thead><tr><th rowspan="2">T<sup>i</sup></th><th colspan="2">Staggered algorithm q = 10</th><th colspan="2">P. C. G. algorithm q = 10</th></tr><tr><th>Iterations</th><th>CPU</th><th>Iterations</th><th>CPU</th></tr></thead><tbody><tr><td>T<sup>1</sup> = 1</td><td>11</td><td>5.362</td><td>6</td><td>3.950</td></tr><tr><td>T<sup>3</sup> = 3</td><td>12</td><td>5.862</td><td>9</td><td>5.790</td></tr><tr><td>T<sup>5</sup> = 5</td><td>13</td><td>6.356</td><td>10</td><td>6.399</td></tr><tr><td>T<sup>7</sup> = 7</td><td>13</td><td>6.378</td><td>8</td><td>5.170</td></tr><tr><td>T<sup>9</sup> = 9</td><td>13</td><td>6.361</td><td>9</td><td>5.800</td></tr></tbody></table>

$c_w = 4.10^{-8} Pa^{-1}$, $c_o = 1.10^{-8} Pa^{-1}$, $\lambda = 5.10^8 Pa$, $\mu = 2.10^8 Pa$, atmospheric pressure = $10^5 Pa$. The viscosities $\eta_o$ and $\eta_w$ and the compressibilities $c_o$ and $c_w$ of the fluids are prescribed in each simulation. We denote by $R$ the ratio of the viscosities, namely $R = \eta_o/\eta_w$. We perform tests with different values of this ratio to check its influence on the flow. The initial pressure is the atmospheric pressure $10^5 Pa$ and the outcoming pressure is $0 Pa$.

Figure 9 shows the time evolution of the $L^2(\Omega)$-norm of the difference of the saturation computed with different values of $q$ and the reference saturation computed
with $q = 1$, both for the staggered scheme and for the preconditioned conjugate gradient algorithm. The figures 10 and 11 show similar results for the pressure and for the displacement. We observe that all three errors grow with $q$. This is not surprising and it shows that one has to be slightly careful with the choice of $q$.

Next we discuss the cost of the computations in terms of the number of iterations and the CPU time. The tables 5 and 6 illustrate the comparison between the two algorithms in the multiscale case where $q = 20$. In the case that $c_o = c_w = 4.10^{-9}$ (see Table 5) we observe that the preconditioned conjugate gradient algorithm is less expensive than the staggered algorithm. On the other hand when $c_o = 1.10^{-8} Pa^{-1}$ and $c_w = 4.10^{-8} Pa^{-1}$ (see Table 6) the two algorithms require a similar number of iterations.

Figure 12 shows the oil saturation and the pressure when the viscosity rate $R = 100$ or in other words when the oil phase is 100 times more viscous than the water phase. We observe the presence of a preferred outcoming direction exhibited by the finger at the time $t = 100 s$; on the contrary the oil comes out in a symmetric fashion in the case of Figure 13 where $R = 1$. The figures 14 and 15 show a similar phenomenon at the time $t = 580 s$. Let us remark that this agrees with the physical predictions.

The previous results have been obtained on a squared mesh without any preferred directions. The figures 16 and 18 exhibit two grids with different preferred directions and their influence on the numerical results at the times $t = 100 s$ and $t = 580 s$; it turns out that one obtains completely different saturation profiles. Here the numerical scheme cannot cope with the physical instability.

<figure><img src="image_10.png" /><figcaption>FIGURE 9. The time evolution of the $L^2(\Omega)$-norm of the error between the reference solution ($q = 1$) and the approximate saturations computed for different values of $q$ on the left-hand-side by the staggered algorithm and on the right-hand-side by the preconditioned conjugate gradient algorithm.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIGURE 10. The time evolution of the $L^2(\Omega)$-norm of the error between the reference solution and the approximate pressures computed for different values of $q$ on the left-hand-side by the staggered algorithm and on the right-hand-side by the preconditioned conjugate gradient algorithm.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 11. The time evolution of the $(L^2(\Omega))^2$-norm of the error between the reference solution and the approximate displacements computed for different values of $q$ on the left-hand-side by the staggered algorithm and on the right-hand-side by the preconditioned conjugate gradient algorithm.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 12. Saturation on the left-hand-side and pressure on the right-hand-side at the time $t = 100$ s with $R = 100$ on a regular mesh.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIGURE 13. Saturation on the left-hand-side and pressure on the right-hand-side at the time $t = 100$ s with $R = 1$ on a regular mesh.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 14. Saturation on the left-hand-side and pressure on the right-hand-side at the time $t = 500$ s with $R = 100$ on a regular mesh.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 15. Saturation on the left-hand-side and pressure on the right-hand-side at the time $t = 500$ s with $R = 1$ on a regular mesh.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIGURE 16. Saturation at the time $t = 100$ s with $R = 100$ and two different meshes: on the left-hand-side with vertically stretched cells and on the right-hand-side with horizontally stretched cells ones.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 17. Saturation at the time $t = 580$ s where $R = 100$ and two different meshes: on the left-hand-side with vertically stretched cells and on the right-hand-side with horizontally stretched cells ones.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIGURE 18. Pressure at time $t = 100$ s where $R = 100$ and two different meshes: on the left-hand-side with vertically stretched cells and on the right-hand-side with horizontally stretched cells ones.</figcaption></figure>
<table><caption>TABLE 5. Comparison in terms of iterations and CPU times between the multiscale staggered algorithm and the preconditioned conjugate gradient algorithm at the times $T^1, T^3, T^7, T^9, T^{10}$ and $T^{20}$ with $c_o = c_w = 4.10^{-9} Pa^{-1}$.</caption><thead><tr><th rowspan="2">T<sup>i</sup></th><th colspan="2">Multiscale staggered algorithm q = 20</th><th colspan="2">P. C. G. algorithm q = 20</th></tr><tr><th>Iterations</th><th>CPU</th><th>Iterations</th><th>CPU</th></tr></thead><tbody><tr><td>T<sup>1</sup> = 20</td><td>126</td><td>0.611E+3</td><td>19</td><td>0.172E+3</td></tr><tr><td>T<sup>3</sup> = 60</td><td>122</td><td>0.507E+3</td><td>17</td><td>0.145E+3</td></tr><tr><td>T<sup>7</sup> = 140</td><td>113</td><td>0.368E+3</td><td>16</td><td>0.107E+3</td></tr><tr><td>T<sup>9</sup> = 180</td><td>108</td><td>0.350E+3</td><td>15</td><td>0.997E+2</td></tr><tr><td>T<sup>10</sup> = 200</td><td>104</td><td>0.378E+3</td><td>15</td><td>0.110E+3</td></tr><tr><td>T<sup>20</sup> = 400</td><td>52</td><td>0.110E+2</td><td>12</td><td>0.524E+2</td></tr></tbody></table>

<table><caption>TABLE 6. Comparison in terms of iterations and CPU times between the multiscale staggered algorithm and the preconditioned conjugate gradient algorithm at the times $T^1, T^3, T^7, T^9, T^{10}$ and $T^{20}$ with $c_o = 1.10^{-8} Pa^{-1}$ and $c_w = 4.10^{-8} Pa^{-1}$.</caption><thead><tr><th rowspan="2">T<sup>i</sup></th><th colspan="2">Multiscale staggered algorithm q = 20</th><th colspan="2">P. C. G. algorithm q = 20</th></tr><tr><th>Iterations</th><th>CPU</th><th>Iterations</th><th>CPU</th></tr></thead><tbody><tr><td>T<sup>1</sup> = 20</td><td>11</td><td>0.488E+2</td><td>14</td><td>0.319E+3</td></tr><tr><td>T<sup>3</sup> = 60</td><td>11</td><td>0.464E+2</td><td>12</td><td>0.285E+3</td></tr><tr><td>T<sup>7</sup> = 140</td><td>10</td><td>0.217E+2</td><td>12</td><td>0.142E+3</td></tr><tr><td>T<sup>9</sup> = 180</td><td>10</td><td>0.216E+2</td><td>11</td><td>0.129E+3</td></tr><tr><td>T<sup>10</sup> = 200</td><td>9</td><td>0.196E+2</td><td>11</td><td>0.130E+3</td></tr><tr><td>T<sup>20</sup> = 400</td><td>6</td><td>0.130E+2</td><td>7</td><td>0.742E+2</td></tr></tbody></table>

## 8. Operation cost

We apply the preconditioned conjugate gradient algorithm even in the case of a slightly nonsymmetric problem. We have at each period the following computation costs:

* We compute the initial residual, and the first descent direction, which requires one reservoir simulation and one geomechanical computation.

* At each iteration we compute $\alpha_l$ and the preconditioned direction $z_{l+1}$, which requires one reservoir simulation and one geomechanical computation.

* At the end of each period, we perform one reservoir simulation.

The fact that we need to perform two reservoir simulations at each iteration is not constraining since one simulation requires 15 mn of CPU time whereas 6 hours are necessary for one geomechanical computation.

**Remark 3.** The convergence criterium for both algorithms is based upon the relative variation of the pressure from one iteration to the next; it is given by

$$ (28) \qquad \frac{\|p_{l+1} - p_l\|_{\infty}}{p_{atm}} < 10^{-j} $$

where $p_{atm}$ is the atmospheric pressure with $j = 6$ for the linear case and $j = 3$ for the nonlinear model.
References

[1] D. Braess, Finite elements, Cambridge University Press, 2001.

[2] L. Y. Chin, and L. K. Thomas, Fully coupled analysis of improved oil recovery by reservoir compaction, SPE Phillips Petroleum Compagny, SPE Journal, 1999.

[3] F. Z. Daïm, Etude théorique et numérique de couplages entre écoulements et déformations mécaniques dans l'extraction d'hydrocarbures, Thèse de Doctorat, Université de Paris XI, 2004.

[4] F. Z. Daïm, R. Eymard, D. Hilhorst, R. Masson and M. Mainguy, A preconditioned conjugate gradient based algorithm for coupling geomechanical-reservoir simulations, Oil & Gas Science and Technology, Vol. 57, 515-523, 2002.

[5] R. Eymard, T. Gallouët and R. Herbin, Finite volume methods, Handbook of numerical analysis, Vol. VII, 713-1020, North-Holland, 2000.

[6] R. Eymard, T. Gallouët, R. Herbin and A. Michel, Convergence of a finite volume scheme for nonlinear degenerate parabolic equations, Numer. Math., Numerische Mathematik, Vol.92, 41-82, 2002.

[7] T. H. Mac Lendon and D. N. Sawyer, Performance prediction for the M-6 area of the Tia Juana field using a rate dependent pore volume compressibility model and extended material balance, SPE Paper 22939, SPE Annual Technical Conference and Exhibition, Dallas, USA, 553-564, 1991.

[8] H. A. Merle, C. J. P. Kentie, G. H. C. Van Opstal and G. M. G. Schneider, The Bachaquero study- A composite analysis of the behavior of a compaction drive/solution gas drive reservoir, SPE/JPT, Eurock 98 Rock Mechanics in Petrolium Engineering, Vol. 1, 1107-1115, 1976.

[9] P. D. Pattillo, T. G. Kristiansen, G. V. Sund and R. M. Kjestadli, Reservoir compaction and sea floor subsidence at Valhall, SPE/ISRM Paper 47274, Proc, Eurock 98 Rock Mechanics in Petrolium Engineering, Vol. 1, 377-386, 1998.

[10] A. Settari and D. A. Walters, Advances in coupled geomechanical and reservoir modeling with applications to reservoir compaction, SPE Reservoir Simulation Symposium, SPE Journal, 1999.

[11] J. E. Sylte, L. K. Thomas, D. W. Rhett, D. D. Bruning and N. B. Nagel, Water induced compaction in the Ekofisk field, SPE Phillips Petroleum Compagny, SPE Journal, 1999.

Laboratoire de Mathématique (UMR 8628), Université de Paris-Sud, 91405 Orsay Cedex, France.

E-mail: fatima.daim@math.u-psud.fr and danielle.hilhorst@math.u-psud.fr

GRIMAAG, Université des Antilles et de la Guyane, 97157 Pointe à Pitre, France.

E-mail: jacques.laminie@univ-ag.fr

Laboratoire d'Analyse et de Mathématiques Appliquées, Université de Marne-La-Vallée, 77454 Marne La Vallée Cedex, France.

E-mail: Robert.Eymard@univ-mlv.fr