Chapter 06
Solving it
The coupled system has two properties that make it awkward: it is a saddle-point problem in space, and it is stiff in time. Both have well-understood consequences, and both punish the obvious implementation. This chapter is about the two traps and the standard ways out.
Discretisation
Almost every production THM code discretises displacement with continuous finite elements, because the momentum balance is naturally variational and the boundary conditions are tractions. For pressure and temperature the field splits: finite elements are common in geotechnics, finite volumes in reservoir simulation, because they conserve mass locally on each cell and handle strong heterogeneity in $k$ more gracefully.
After discretising in space and applying backward Euler in time, one step of the coupled system reads:
This matrix is symmetric but indefinite. That single fact drives most of what follows: it forbids naive iterative solvers such as conjugate gradients without a block preconditioner, and — more insidiously — it means that not every combination of displacement and pressure interpolation is admissible.
Trap 1 — spurious pressure oscillations
In the incompressible, undrained limit — small $\Delta t$, small permeability, small storage — the bottom-right block $\ten{S}+\Delta t\,\ten{H}$ nearly vanishes and the system becomes a pure saddle-point problem, structurally identical to Stokes flow. It then inherits the inf–sup (LBB) condition: equal-order interpolation of displacement and pressure is not stable, and the pressure field decorates itself with checkerboard oscillations near drained boundaries.
This is not a bug you can refine away. Halving the element size while holding $\Delta t$ makes it worse, because the offending dimensionless group is $\Delta t\,c/h^2$: the oscillations appear when the time step is too small for the mesh. The figure below assembles the actual 1-D finite-element system and solves the first time step, so the oscillation you see is the real discrete solution, not an illustration of one.
Equal-order elements, and what stabilisation does
Linear displacement and linear pressure on the same mesh, one backward-Euler step from a sudden surface load. The dashed line is the Terzaghi analytic solution at the same instant. Reduce the time step and the numerical pressure develops a saw-tooth near the drained surface; switch on the fluid-pressure-Laplacian stabilisation and it disappears without moving the physical answer.
- Satisfy inf–sup. Use Taylor–Hood elements: quadratic displacement, linear pressure. Correct, robust, and the most expensive per node.
- Stabilise. Add a consistent pressure-Laplacian or polynomial-pressure-projection term scaled by $h^2/(4G)$. Cheap, keeps equal order, and is what most codes actually do.
- Respect the time-step limit. Vermeer and Verruijt’s criterion $\Delta t \ge h^2/(6c)$ keeps an unstabilised equal-order scheme clean. Simple, but it forbids resolving the early transient — which is often exactly what you wanted.
Trap 2 — solving the blocks separately
Solving the full block system at once — monolithic coupling — is the unconditionally stable choice, but it means writing one solver that owns all the physics. In practice people very often want to reuse an existing flow simulator and an existing geomechanics code, exchanging fields between them. That is a sequential scheme, and whether it converges depends entirely on which quantity you hold fixed while solving each sub-problem.
The four classical choices, in a scalar model problem where the coupling strength is $\tau = \alpha^2 M/K_{od}$:
| Scheme | Flow solved at fixed… | Stability | Contraction factor |
|---|---|---|---|
| Drained split | displacement (mechanics first) | conditional | τ |
| Undrained split | fluid mass (mechanics uses ) | unconditional, slow | → 1 for large τ |
| Fixed-strain split | strain (flow first) | conditional | τ |
| Fixed-stress split | total stress (flow first, with a stress-corrected storage) | unconditional | τ /(1+τ) |
Convergence of the four splitting schemes
The iteration is run for real on the scalar model problem $K_{od}\varepsilon - \alpha p = \sigma$, $\;p/M + \alpha\varepsilon = \zeta$, starting from a unit error. Raise the coupling strength past $\tau = 1$ and the drained and fixed-strain splits diverge outright; the fixed-stress split keeps converging no matter how strongly coupled the problem is, which is why it has become the default.
Why fixed-stress works
The trick is to solve the flow problem with a storage coefficient that already anticipates the mechanical response. Instead of $S_\varepsilon$, use $S_\varepsilon + \alpha^2/K_{od}$ — the storage you would measure at constant total stress rather than constant strain — and add the corresponding stress-rate term to the right-hand side:
The extra storage is precisely the mechanical stiffness the flow solver would otherwise ignore, so the first iterate is already nearly right. It costs one line in the flow assembly and converts a conditionally-stable scheme into an unconditionally stable one.
Time stepping and nonlinearity
Start small, then grow
Every THM problem has an initial transient — undrained loading, a heater switch-on — over which the pressure changes by orders of magnitude. Geometric time-step growth from a small first step, with a target on the Newton iteration count, is standard and effective. But respect the equal-order lower bound if you are not stabilised.
Upwind the advection, or pay
The $(\rho c)_f\vec q\cdot\nabla T$ term is hyperbolic. Central differencing it on a grid with cell Péclet number above 2 produces oscillations that look exactly like a physical instability. Upwinding, SUPG, or a flux-limited scheme is not optional in advection-dominated problems.
Lag the weak couplings
Temperature feeds back on the others weakly enough that a T→H→M sequential order, with one pass per step, is usually accurate. Verify it by running one case monolithically — not by assuming.
Nonlinear coefficients need iteration
$k(\phi)$, $\mu(T)$, $k_r(S)$ and any plasticity make the system genuinely nonlinear. Newton with a consistent tangent converges quadratically; a fixed-point update of the coefficients converges linearly and can stall when the retention curve is steep.
Verifying a coupled code
A checklist that catches most implementation errors, in the order that isolates them fastest:
- Patch tests on each physics separately — uniform stress, uniform flux, uniform temperature gradient — must be exact to machine precision.
- Terzaghi. One-dimensional, analytic, exercises the $\ten{Q}$ blocks and the time integration. If the early-time pressure is not $\sigma_0$ everywhere, the undrained limit is wrong.
- Mandel. Two-dimensional, analytic, and the only common benchmark whose answer is non-monotonic. It catches a lagged or transposed coupling block that Terzaghi passes happily.
- Booker & Savvidou — a point or spherical heat source in a poroelastic half-space, analytic — brings the thermal blocks in.
- A convergence study in both $h$ and $\Delta t$, checking that the observed orders match the scheme. A coupled code that is first-order accurate where it should be second-order usually has a term evaluated at the wrong time level.
- Energy and mass audits over the whole domain, every step. Cheap, and they catch boundary-condition sign errors that convergence studies miss.