# Coupled Geomechanics and Multiphase Flow

The main focus of this section is to describe the coupled system of geomechanics and multiphase flow equations.

## Governing Equations

The governing equation of coupled geomechanics and multiphase flow can be described in terms of solid and fluid equations [geomechanics].

**[geomechanics]**

Wan, Jing. *Stabilized finite element methods for coupled geomechanics and multiphase flow*. Diss. stanford university, 2003.

**Equations for the Solid**

We assume that the solid density is constant, and therefore, *the mass balance equation of the deforming porous medium* is

The solid velocity $\mathbf{v}_s$ is given by

and therefore we have

where the volumetric strain

It can be shown that Equation 1 leads to the *variation of porosity*:

**Equations for the Fluids**

The *mass balance equations of multiphase multicomponent fluid* are given by

The *linear momentum balance equation* is given by

Here, $K$ is the permeability tensor, but in our case we assume it is a space varying scalar value. $k_{ri}(S_i)$ is a function of $S_i$, and typically the higher the saturation, the easier the corresponding phase is to flow. $\tilde \mu_i$ is the viscosity. $Z$ is the depth cordinate, $\rho_i$ is the density, $\phi$ is the porosity, $q_i$ is the source, $P_i$ is the fluid pressure and $g$ is the velocity constant. We assume the movement of the solid is slow in this study, therefore Darcy’s law is still valid without acceleration terms.

is the relative velocity of the phase $i$ with respect to $\mathbf{v}_s$ (also called *interstitial velocity*).

**Fluid and Mechanics Coupling**

Here $\sigma'$ should be understood as the *effective stress*, which allows us to treat a multiphase porous medium as a mechanically equivalent single‐phase continuum. We use the concept of equivalent pressure in (5)

**Constitutive Relation**

The constitutive relation connects $\sigma'$ and the displacement $\mathbf{u}$. For example, the linear elastic relation is expressed as

Here, the strain is the Cauchy strain

Instead of assuming a linear elasticity model for the geomechanics, we can also model the subsurface solid material by a viscoelasticity model (see here for details). For example, the constitutive relation for the two dimensiona Maxwell material is as follows:

## Numerical Scheme

We use an iterative algorithm to solve the coupled equation; namely, we alternatively solve the mechanics equation and flow equation.

**Fluid Equation**

We define the fluid potential

and the capillary potential

Here the capillary pressure $P_c = P_1-P_2$ is assumed to be small. We define mobilities

We have the following formula from Equations 3-4:

We can solve for $\Psi_2$ using a Poisson solver.

Next, we have from Equations 3-4

Note we have an extra term $q_1 \frac{m_2}{m_1}$ to account for the assumption of zero capillary pressure.

Equation 8 is a nonlinear equation in $S_2$ ($m_2$ is defined in terms of $S_2=1-S_1$) and requires a Newton-Raphson solver.

**Solid Equation**

Upon solving the fluid equation, we obtain $S_1, S_2, \Psi_2$. We can use $\Psi_2$ to estimate the fluid pressure $P_1$ and $P_2$. Use Equations 5 and 6, we solve for $\mathbf{u}$ using

Here $\varepsilon_v = \varepsilon_{xx} + \varepsilon_{yy} = u_x + u_y$.

## Example

We simulate the coupled geomechanics (both the linear elastic and the viscoelastic material) and two phase flow model. The two models share the same Lamé constants $\lambda$ and $\mu$ (or equivalently, $E$ and $\nu$), except that the viscoelasticity model has one more viscosity parameter $\eta$. We show an animation of the von Mises stress, displacement (magnified by 50x), the fluid potential for the wet phase ($\Psi_2$), and the saturation. Particularly, we plot the $x$-direction displacements on the surface. The displacement data will be used as observation in the inverse problem.

The setting is as follows: we consider a layer model for the relative permeability $K$ (the values of the space varying $K$ are known). The blue triangle denotes the injection well, and the orange triangle denotes the production well. The blue triangles denote receivers, where we can collect horizontal displacement–-for example, we can obtain these data from satellite images.

The two flow phases are the oil and the water (the wet phase), and we use number 1 and 2 to denote them respectively. The computational domain is has 450 depth and 900 width, and it is discretized into $15\times 30$ uniform squares for both the finite element method (for the mechanics equation) and the finite volume method (for the fluid equation). The simulation has time horizon $[0,1000]$ and is temporally discretized into 50 intervals. The other parameters for simulations are: $\rho_1=1053$, $\rho_2=501.9$, $\tilde\mu_1=1$, $\tilde\mu_2=0.1$, $g =9.8$. In the linear elasticity case, we use $E=6\times 10^9$ and $\nu=0.35$. In the viscoelasticity case, in addition to the Lamé constants, the viscosity parameter $\eta=6.7\times 10^{11}$.

Model | Linear Elasticity | Viscoelasticity |
---|---|---|

Displacement $u$ | ||

Saturation $S_2$ | ||

Potential $\Psi_2$ | ||

Pressure $p$ | ||

Von Mises Stress $\sigma_{\mathrm{vm}}$ | ||

Observation $u_{\mathrm{obs}}$ |