Chapter 02 · H
Hydraulics
The fluid half of the problem is a conservation statement and a constitutive law: water is neither created nor destroyed, and it moves down its energy gradient at a rate set by the pore geometry. Everything else — storage, diffusion, capillarity — follows from combining the two.
Darcy’s law
In 1856 Henry Darcy, tasked with supplying Dijon with clean water, measured the discharge through sand filters and found it proportional to the head loss and inversely proportional to the filter length. That empirical result survives essentially unchanged, because it is derivable: average the Stokes equations over a periodic pore geometry and the interfacial drag term collapses to a linear resistance.
Hover or tap a highlighted group.
- $\vec q$ Darcy flux m s⁻¹ Volume of fluid crossing unit bulk area per second — not the speed of a water molecule. Divide by porosity to get the pore (seepage) velocity that a tracer actually experiences.
- $\ten{k}/\mu$ Mobility m² Pa⁻¹ s⁻¹ Intrinsic permeability — a property of the porous medium alone — divided by the fluid viscosity. Splitting them this way is what lets a single material property serve water, brine, oil and CO₂ — and it is where temperature enters the hydraulics.
- $\nabla p$ Pressure gradient Pa m⁻¹ The driving force. In groundwater practice it is folded together with gravity into the hydraulic head $h = z + p/(\rho_f g)$.
- $\rho_f\vec g$ Buoyancy term Subtracting the hydrostatic gradient leaves only the part of the pressure field that drives flow. When density varies — with temperature or salinity — this term no longer cancels, and free convection becomes possible.
Two vocabularies coexist. Hydrogeology prefers hydraulic conductivity $K = k\rho_f g/\mu$ in metres per second, which bundles the fluid in; reservoir engineering and poromechanics prefer intrinsic permeability $k$ in square metres (or darcy, $1\ \mathrm{D} \approx 10^{-12}\ \mathrm{m^2}$). Any THM model must use the intrinsic form, because the fluid properties are unknowns.
Darcy’s experiment, and where his law stops working
The column is solved with the Forchheimer extension $i = \tfrac{\mu}{\rho g k}q + \tfrac{\beta}{g}q^2$ using $\beta = 0.55/\sqrt{k}$, so the departure from linearity is computed rather than sketched. Darcy’s law is the tangent at the origin; it holds while the pore Reynolds number stays below about 1–10. Coarse gravel near a pumping well is the usual place it fails.
What sets the permeability
Permeability has units of area, and the area in question is roughly the cross-section of a typical pore throat. That is why it scales with the square of the grain size and why it spans fifteen orders of magnitude between gravel and intact granite — a wider range than any other material property in engineering.
The Kozeny–Carman relation is a semi-empirical closure of exactly the interfacial integral that averaging left behind: $\phi^3/(1-\phi)^2$ stands in for connectivity and specific surface, and the factor 180 absorbs tortuosity. It works well for uniform sands and badly for clays, fractured rock, and anything with a bimodal pore-size distribution — but it captures the two dependences that matter, and it is the reason a small change in porosity produces a large change in permeability.
Permeability: grain size, porosity, and the range of real materials
Left: Kozeny–Carman at four porosities. Right: the observed range for common geological materials, with your current $(d_{50},\phi)$ marked. The cube dependence on porosity is the quiet villain of coupled modelling — a few per cent of compaction can halve the permeability, which is a M→H coupling hiding inside a coefficient.
Storage, and why pressure diffuses
Mass conservation for the fluid in a deforming skeleton says that the rate of change of fluid content equals minus the divergence of the flux. Fluid content changes for three reasons: the fluid compresses, the pore space compresses, and the skeleton strains. Collecting the first two into a storage coefficient and substituting Darcy’s law gives a diffusion equation.
Hold the skeleton still — $\varepsilon_v$ constant — and what remains is the classical groundwater equation with hydraulic diffusivity $c = k/(\mu S_\varepsilon)$. It is a genuine diffusion coefficient in m² s⁻¹, and it means pressure signals spread as $\sqrt{ct}$ rather than travelling at a speed. That single square root explains most of the timescales in the subject.
Pressure diffusion from a step change at the boundary
Analytic solution $p(x,t) = \Delta p\;\mathrm{erfc}\!\big(x/2\sqrt{ct}\big)$ for a half-space. The marker tracks the penetration depth $2\sqrt{ct}$, where the pressure has risen to about 5 % of the applied step. Quadrupling the time only doubles the depth — which is why injection pressure fronts stall, and why deep repositories stay dry for a very long time.
Pressure diffusion, heat conduction and Terzaghi consolidation are the same partial differential equation with different coefficients. Any solution you know for one — the $\mathrm{erfc}$ half-space, the Theis well function $W(u) = E_1(u)$, the Fourier series in a finite layer — transfers directly to the others. That is worth a great deal when checking a coupled code.
When the pores are only partly wet
Above the water table, air and water share the pore space and the water is held under suction: its pressure is below atmospheric because curved menisci pull it into the narrow throats. The capillary pressure $p_c = p_g - p_w = 2\sigma\cos\theta/r$ is larger in smaller pores, so as a soil dries the water retreats into progressively finer pores and the suction rises — steeply.
Two constitutive curves close the problem. The retention curve relates saturation to capillary pressure; the relative permeability says how much of the saturated conductivity survives when only part of the network is water-filled. Van Genuchten’s retention model and Mualem’s conductivity model, used together, are the standard pairing:
Retention and relative permeability
Parameters after Carsel & Parrish (1988). Note the vertical scale on the right: at 70 % saturation a sand has already lost more than 90 % of its conductivity. This nonlinearity — four orders of magnitude of $k_r$ over the wetting range — is what makes Richards’ equation hard, and it is why unsaturated solvers need adaptive time stepping.
Substituting the two curves into mass conservation gives Richards’ equation, a nonlinear diffusion–advection equation in which both the diffusivity and the storage coefficient depend on the unknown:
Infiltration into a dry column — Richards’ equation, solved live
Solved with the mass-conservative modified Picard scheme of Celia, Bouloutas & Zarba (1990): backward Euler in time, central differences in space, the storage term written in mixed form so that mass is conserved to machine precision, and a tridiagonal solve at each iteration. The front is sharp because $k_r$ collapses ahead of it; cumulative infiltration follows Philip’s $\sqrt{t}$ law early on and becomes linear once gravity takes over.
Retention curves are hysteretic: a soil holds more water while draining than while wetting at the same suction, because of ink-bottle pore geometry and contact-angle differences. Most engineering models ignore it, which is defensible for monotonic processes and indefensible for wetting–drying cycles. In THM problems suction is also temperature-dependent, since surface tension falls with heat — a T→H coupling that matters in near-field repository analysis.