Chapter 05

The coupled system

Three balance laws, three primary unknowns, and a set of cross terms that refuse to be ignored. This chapter assembles the system and then tests it against three problems whose exact solutions are known — because a coupled code that cannot reproduce them is not coupled, it is merely complicated.

Assembling the system

Take the momentum balance from chapter 04, the fluid mass balance from chapter 02 and the energy balance from chapter 03, and write them in terms of the primary unknowns $(\vec u, p, T)$. Nothing new is introduced; the couplings were already there, hiding in the constitutive laws.

Hover or tap a highlighted group. Colours mark which physics the term belongs to.

  • $\alpha\nabla p$ H → M Effective stress. Symmetric with term d below — the pair forms the classic Biot coupling and, in a discretised system, an exactly symmetric off-diagonal block.
  • $3K\alpha_T\nabla T$ T → M Thermal stress. Note there is no matching term in the energy balance: the reverse coupling (thermoelastic heating) is real but negligible, so this coupling is essentially one-way.
  • $S_\varepsilon$ Storage $\phi/K_f + (\alpha-\phi)/K_s$. Small in stiff rock — which is exactly why pressure responds so sharply to a small volumetric strain.
  • $\alpha\,\dot\varepsilon_v$ M → H Pore-volume change. Deleting this term (and only this term) turns a poroelastic model into an uncoupled flow model — and destroys the Mandel–Cryer effect entirely.
  • $\Gamma\,\dot T$ T → H Thermal pressurisation. $\Gamma = \phi(\beta_f-\beta_s)$ to leading order: the fluid expands about ten times more than the pore space, and if it cannot escape the pressure climbs.
  • $\vec q\!\cdot\!\nabla T$ H → T Advection. The only nonlinear coupling in the set that is not hidden in a coefficient, and the one that forces an upwinded or stabilised discretisation.

Two further couplings live inside the coefficients rather than as separate terms, and they are easy to forget: permeability depends on porosity, which depends on volumetric strain ($k(\phi(\varepsilon_v))$), and viscosity depends on temperature ($\mu(T)$). Both are nonlinear couplings, and both can matter more than the explicit terms — a factor of five drop in viscosity between 20 and 150 °C changes flow rates more than most pressure changes do.

Block structure

Discretised, the system has a recognisable shape. Writing the unknown vector as $(\vec U, \vec P, \vec\Theta)$ for displacement, pressure and temperature degrees of freedom:

The $\ten{Q}$ / $\ten{Q}^{\mathsf T}$ pair is symmetric, so the mechanical–hydraulic block is symmetric-indefinite — a saddle-point system, with all the discretisation and solver consequences that implies. The thermal row is not symmetric with the others: the $\ten{L}$ and $\ten{\Gamma}$ blocks have no transposes in the bottom row, because the mechanical and hydraulic feedback on temperature is negligible. That asymmetry is not sloppy modelling; it is a physical fact about the magnitudes, and it is what makes sequential T→H→M solution strategies viable. More on that in chapter 06.

Benchmark 1 — Terzaghi consolidation

A saturated layer of thickness $H$, drained at the top, loaded instantaneously with a surface pressure $\sigma_0$. At $t=0^+$ the water carries everything: $p = \sigma_0$ throughout. Then it drains, the effective stress rises, and the layer settles. Since the strain is uniaxial and the total stress constant, the excess pressure obeys pure diffusion, and the solution is a Fourier series:

Figure 5.1

Terzaghi isochrones and the settlement curve

Series summed to 200 terms, so the $t\to0$ limit is resolved. Left: the excess-pressure isochrones; each curve is a snapshot in dimensionless time. Right: the degree of consolidation $U = 1-\sum(2/M^2)e^{-M^2T_v}$, the fraction of the final settlement achieved. The useful rules of thumb — $U=50\%$ at $T_v=0.197$, $U=90\%$ at $T_v=0.848$ — are read straight off it.

Benchmark 2 — Mandel’s problem

Terzaghi’s solution is one-dimensional, and one-dimensionality hides the most interesting thing poroelasticity does. Take a specimen between two rigid, impermeable, frictionless plates, drained on its vertical sides, and load it vertically at $t=0$. The pressure jumps everywhere to the Skempton value — and then, at the centre, it rises further before it falls.

The reason is a genuine two-way coupling. Water escapes at the edges first, so the edges lose their pore pressure and soften. The rigid plate cannot tilt, so the load transfers inward onto the still-stiff, still-undrained core. That extra load squeezes the core, and by the $\alpha\dot\varepsilon_v$ term it raises the pressure there. Only when drainage catches up does the pressure finally decay.

Why every poroelastic code is tested on this

The Mandel–Cryer overshoot cannot be produced by a diffusion equation for pressure alone. It requires the mechanical response to feed back into the fluid mass balance. If a code shows a monotonic decay here, its $\ten{Q}^{\mathsf T}$ block is wrong, missing, or lagged — and no amount of mesh refinement will fix it.

Figure 5.2

The Mandel–Cryer effect

Analytic solution of Mandel (1953) in the form given by Abousleiman et al. (1996). The eigenvalues $\alpha_n$ solve $\tan\alpha_n = \tfrac{1-\nu}{\nu_u-\nu}\,\alpha_n$ and are found here by bisection; 60 terms are summed. The overshoot grows with the contrast between drained and undrained Poisson’s ratios — set $\nu_u=\nu$ and it disappears completely, which is the uncoupled limit.

Benchmark 3 — heating a saturated column

Now add temperature. A saturated column, drained at the top, is heated at its base — the canonical near-field problem for a waste repository. Three things happen in sequence, on three different timescales, and the full THM system is needed to get any of them right.

  1. Pressurisation. Heat arrives faster than the water can escape, and because the water expands roughly ten times more than the pore space, pressure climbs at $\Gamma/S_\varepsilon \approx 0.1$–$1$ MPa per kelvin.
  2. Dissipation. The pressure drives flow and decays diffusively, on the hydraulic timescale $H^2/c$.
  3. Expansion. The skeleton itself expands thermally, on the much longer thermal timescale $H^2/\kappa$ — so the sample first heaves from pressure, then settles as pressure dissipates, then heaves again as the skeleton itself warms through.
Figure 5.3

Thermal pressurisation and thermal consolidation, solved live

A fully coupled 1-D THM solver under uniaxial strain: equilibrium is enforced exactly ($\sigma_{zz}$ uniform), so the mechanics reduces to $\varepsilon_v = (\sigma_0 + \alpha p + 3K\alpha_T\Delta T)/K_{od}$, and pressure and temperature are advanced with backward-Euler tridiagonal solves. Watch the non-monotonic surface displacement — heave, settle, heave — which no single-physics model reproduces.

Permeability decides the whole story

Slide the permeability. Above about $10^{-16}\,\mathrm{m^2}$ the pressure escapes as fast as it is generated and thermal pressurisation is invisible. Below about $10^{-19}\,\mathrm{m^2}$ the response is effectively undrained and the pressure can approach — or exceed — the minimum principal stress, hydraulically fracturing the rock. Repository host rocks are chosen from that lower band, which is precisely why this coupling is a licensing issue.