Staggered Scheme

Project file on GitLab

Injection and Production in 1D Linear Poroelastic Medium

This benchmark simulates a soil column with fluid injection at the bottom and a production well at the top. It is taken from from Kim [1], in detail it coincides with one of his examples (case 2, coupling strength $\tau=1.21$). A brief description of the used staggered scheme follows at the end.

Simulation model with fluid source, sink, observation point and boundary conditions

The fluid enters and leaves only via the source and sink in the domain, there is no flow across the boundaries. The displacements at the bottom are fixed, whereas there is a vertical traction applied on top. Originally the problem is one-dimensional, for simulation with OpenGeoSys it is created in two dimensions with corresponding boundary conditions. All parameters are concluded in the following tables.

Material Properties
Property Value Unit
Fluid density $10^3$ kg/m$^3$
Viscosity $10^{-3}$ Pa$\cdot$s
Fluid compressibility $27.5\cdot 10^{-9}$ Pa$^{-1}$
Porosity $0.3$ -
Permeability $493.5\cdot 10^{-16}$ m$^2$
Young’s modulus (bulk) $300\cdot 10^6$ Pa
Poisson ratio (bulk) $0$ -
Biot coefficient $1.0$ -
Solid density $3\cdot 10^3$ kg/m$^3$
Solid compressibility $0$ Pa$^{-1}$
Dimensions and Discretization
Property Value Unit
Height ($y$) $150$ m
Width ($x$) $10$ m
Finite Elements $15$ Taylor-Hood quadrilateral elements 10 m $\times$ 10 m
Time step $86.4\cdot10^3$ s
Sources/Sinks, Initial and Boundary Conditions
Hydraulic Mechanical
injection over area $10$m$\times 10$m $+1.16\cdot 10^{-4}$ kg/(m$^3$s) -
production over area $10$m$\times 10$m $-1.16\cdot 10^{-4}$ kg/(m$^3$s) -
top no flow $\sigma_{yy}=-2.125\cdot 10^6$ Pa
left no flow $u_x=0$
right no flow $u_x=0$
bottom no flow $u_y=0$
initial state $p(x,y)=2.125\cdot 10^6$ Pa $u_x(x,y)=u_y(x,y)=0$

The gravity related terms are neglected in both: the Darcy velocity and the momentum balance equation.

Note that 100 time steps were used for the following results, whereas the provided input file is set to 1 time step (1 day = 86400 s). Kim plots his results over nondimensional time, referring to the time at which the produced fluid volume equals the pore volume of the domain (450 days).

Pressure at observation point (marked by circle) versus time (t=0…100 days) and spatial pressure distribution at t=100 days

Staggered Scheme: Fixed-stress splitting

For each time step run alternating simulations of the hydraulic (H) problem and the mechanical (M) problem until a convergence criterium is met. The fixed-stress split starts with the mass balance (H) followed by the momentum balance (M). These coupling iterations (H,M,H,M,…) add another iteration level compared to the monolithic formulation (HM). However, due to splitting into smaller problems this may result in a speedup.

The weak form of mass (scalar test function $v$) and momentum balance (vectorial test function $\mathbf{v}$) after discretization in time (Euler-backward) reads at the $n^\mathrm{th}$ time step for the field equations (boundary terms omitted) $$ \eqalign{ \int_V \left( \alpha \frac{ \varepsilon_V(\mathbf{u}^{n})- \varepsilon_V(\mathbf{u}^{n-1})}{\Delta t}v +S\frac{p^n-p^{n-1}}{\Delta t}v +\frac{k}{\mu}\nabla p^n \cdot \nabla v \right) \mathrm{d}V &=&\int_V \left(s v-\frac{k}{\mu}\varrho_\mathrm{f} \nabla \cdot \mathbf{g} v \right) \mathrm{d}V, \cr \int_V \left( \boldsymbol{\sigma}^\mathrm{eff}(\mathbf{u}^n):\boldsymbol{\varepsilon}(\mathbf{v}) -\alpha p^n \varepsilon_V(\mathbf{v}) \right) \mathrm{d}V &=& \int_V \ \varrho_\mathrm{b} \mathbf{g}\cdot \mathbf{v} \ \mathrm{d}V, } $$ where $\alpha$ denotes Biot coefficient, $S$ storage coefficient, $k$ permeability, $\mu$ viscosity, $\varrho_\mathrm{f}$ fluid density and $\varrho_\mathrm{b}$ bulk density. The fixed-stress split makes the assumption of constant volumetric stress to eliminate the displacement from the mass balance. At the $i^\mathrm{th}$ coupling iteration this assumption leads to $$ K_\mathrm{b} \varepsilon_V(\mathbf{u}^{n,i})-\alpha p^{n,i} = K_\mathrm{b} \varepsilon_V(\mathbf{u}^{n,i-1})-\alpha p^{n,i-1}, $$ where $K_\mathrm{b}$ denotes the drained bulk modulus. Thus the weak form of the mass balance (H) depends only on pressure $p^{n,i}$ as unknown $$ \int_V \left( \beta\frac{p^{n,i}-p^{n,i-1}}{\Delta t}v +S\frac{p^{n,i}-p^{n-1}}{\Delta t}v +\frac{k}{\mu}\nabla p^{n,i} \cdot \nabla v \right) \mathrm{d}V =\int_V \left( s v-\frac{k}{\mu}\varrho_\mathrm{f} \nabla \cdot \mathbf{g} v -\alpha\frac{ \varepsilon_V(\mathbf{u}^{n,i-1})- \varepsilon_V(\mathbf{u}^{n-1})}{\Delta t} v \right) \mathrm{d}V, $$ where $\beta=\beta^\mathrm{ph}=\frac{\alpha^2}{K_\mathrm{b}}$ would correspond to the physically motivated volumetric stress. However, choosing another artificial volumetric stress, i.e. using a different bulk modulus than $K_\mathrm{b}$, may accelerate the convergence of the coupling iterations. The analysis by Mikelic and Wheeler [2] revealed that $\beta^\mathrm{MW}=\frac{\alpha^2}{2K_\mathrm{b}}$ is generally close to the a-priori unknown optimal value $\beta^\mathrm{opt}$.

The obtained pressure $p^{n,i}$ is then inserted into the momentum equation (M) with $\mathbf{u}^{n,i}$ as unknown $$ \int_V \ \boldsymbol{\sigma}^\mathrm{eff}(\mathbf{u}^{n,i}):\boldsymbol{\varepsilon}(\mathbf{v}) \ \mathrm{d}V = \int_V \left( \varrho_\mathrm{b} \mathbf{g}\cdot \mathbf{v} + \alpha p^{n,i} \varepsilon_V(\mathbf{v}) \right) \mathrm{d}V. $$ These two steps (H,M) are repeated until convergence is reached and thus the solution of the $n^\mathrm{th}$ time step is found.



Kim, J. and Tchelepi, H.A. and Juanes, R (2009): Stability, Accuracy and Efficiency of Sequential Methods for Coupled Flow and Geomechanics. SPE International, vol. 16, p. 249--262, DOI:


Mikelic, A. and Wheeler, M.F. (2013): Convergence of iterative coupling for coupled flow and geomechanics. Computational Geosciences, vol. 17, p. 455--461, DOI:

This article was written by Dominik Kern and Wenqing Wang. If you are missing something or you find an error please let us know. Generated with Hugo 0.79.0. Last revision: March 25, 2021
Commit: use unit mD precisely and correct nondim time a8a2144fa  | Edit this page on