• Physics of groundwater flow can be genuinely tricky.

Everything else today serves the above statements.

A bit of global hydrology

for geophysicists

• Movement of water between atmosphere, ocean, ice, surface and subsurface as gas, liquid and solid.

Discount the ocean and the ice.
groundwater is ~99% of the unfrozen fresh water on the planet.

• It supplies about half of domestic water use worldwide.
• It supplies over 40% of the water used for irrigation.

Given to us
P - precipitation.
Comes from the climate system. We take it as forcing. We are essentially powerless over it.

Mostly forced
E - evaporation.
Set by temperature, water availability and aerodynamics (wind). Largely dictated to us.

Ours to model
Q + \Delta S - what runs off and what is stored underground.
This is the subsurface. This is where the physics is, and where we come in.

The balance funnels us straight to the term we can actually compute: the water below the surface.

Two stores, very different behaviour:

• Surface water - rivers, dams, lakes. Highly visible, but set largely by topography and by this year’s rainfall. It rises and falls fast.
• Groundwater - the slow, buffered store. It integrates decades of recharge and keeps rivers alive through drought (baseflow).

Surface and subsurface flow are inseparable.

Groundwater flow cannot be observed directly.

• The subsurface is sampled only at scattered boreholes.
• Conductivity and storativity are known, if at all, at a handful of points.
• Through surrogate methods, variations can be observed. satellite geodesy, seismology.
• You cannot interpolate/extrapolate between boreholes: the water between them moves by physical law, not by an assumption of smoothness.

A forecast needs two components. And an efficient inverse needs the forward solver plus its sensitivities. So the forward solver is the foundation – build that first.

The forward physics

From Darcy to Richards

AI Translation: ”… the experiments I carried out at Dijon together with the engineer M. Charles Ritter, to determine the laws of the flow of water through sands. The experiments were repeated by the chief engineer M. Baumgarten.”

”… there is a linear relationship between the head and the volume you are getting out of the filter”

Darcy (1856), studying sand columns for Dijon’s water supply, gave the momentum balance for slow flow through a porous medium:

Combine it with conservation of mass and you get the classical groundwater-flow equation:

In the saturated zone K is constant, the equation is essentially linear, and mature tools like MODFLOW are built on it.

• Agriculture: soil moisture sets crop yield and irrigation demand.
• Recharge: it decides how much rain actually reaches the aquifer.
• Bushfire: dry soil means dry fuel; a leading indicator of fire danger.
• Contamination: nitrate, pesticides, PFAS, spills all enter groundwater through it.

The soil turns the constants into curves: moisture \theta(h) and conductivity K(h).

Lorenzo Richards

Above the water table lies the unsaturated (vadose) zone. Pores hold both air and water, and the water is held under negative (capillary) pressure.

Darcy’s law still holds. What changes is that the coefficients now depend on the state:

• Conductivity collapses as the soil dries: K(h) = K_s \, k_r(h), falling many orders of magnitude.
• Storage depends on the moisture the pressure sets up.

Darcy in the unsaturated zone, \vec q = -K(h)\,\nabla(h+z), combined with mass conservation, gives Richards’ equation. We solve the mixed form, because it is the only one that conserves mass under discretisation:

• \theta(h): volumetric moisture content \theta = \theta(h)
• K(h): hydraulic conductivity; K = K_s \, k_r(h)
• z: elevation head – the gravity term that drives percolation
• s: sources and sinks (recharge, extraction)

Reliability
Are the solutions accurate?

Scalability
Can we solve a very big system? (high # of DoF)

Flexibility
Terrain-following meshes and a choice of discretisation?

Connection to data
Real elevation, stratigraphy and borehole observations?

Inversion
Can we do gradient-based (adjoint) optimisation?

Built in G-ADOPT (the Geoscientific ADjoint Optimisation Platform) on Firedrake, PETSc and pyadjoint - the same stack we use for mantle and lithosphere dynamics, now pointed at the near surface.

• Discontinuous Galerkin discretisation: locally mass-conservative, tolerant of sharp fronts and material jumps.
• Newton iteration on an exact, symbolically derived Jacobian.
• Implicit Runge-Kutta time stepping (Irksome); Backward Euler here.
• Anisotropy-aware multigrid preconditioning that coarsens to match the thin-sheet aspect ratio.
• Terrain-following mesh assembled from a national elevation model and borehole stratigraphy.

Verification of the solvers

Or the part I convince you it does!

We check the solver against a suite of benchmarks:

• Tracy (2002) - an analytic solution. The strictest test: no reference-solution ambiguity.
• Vauclin (1979) - a measured water-table recharge experiment, compared against four independent codes.
• Mass conservation - local and global mass balance held to round-off.
• A heterogeneous 3D case (Cockett 2018) - a realistic test before we extend to continental scale.

Vauclin: localised recharge builds a wetting front and a groundwater mound – unsaturated percolation and saturated flow in one equation.

Head evolving to steady state at t = 0,\; 5\times10^4,\; 2.5\times10^6 s.

Why Tracy? An exact solution to Richards’ equation – so we measure the true error, not code-versus-code.

• Exponential soil, \theta,\,K \propto e^{\alpha h}, linearises the equation and admits an analytic solution.
• Square domain, 15.24 \times 15.24 m. A sinusoidal head on the top boundary drives infiltration into initially dry soil.
• \alpha = 0.25 /m, \theta_r = 0.15, \theta_s = 0.45, K_s = 10^{-5} m/s.
• Driven to steady state with Backward Euler.

Testing mass conservation

Most fundamental constraint in hydrology: No drop left behind

Discontinuous Galerkin conserves water locally, cell by cell, to the tolerance of the nonlinear solve.

Google Maps

Map of MDB

The point of the preconditioning is that the linear solves stay efficient as the problem grows.

• Weak scaling: hold work per processor fixed, grow the problem – runtime stays flat.
• Strong scaling: fix the problem, add processors – it keeps speeding up.

Reliability
Accurate, convergent, mass-conservative.

Scalability
Weak- and strong-scaling at basin scale.

Flexibility
Terrain-following meshes, choice of discretisation.

Connection to data
Real elevation, stratigraphy and borehole data.

Inversion
Adjoint is set up – this is the next piece.

We now have a reliable solver for the full Richards’ equation at basin scale.

Forward modelling: done. Inversion and surface coupling: next.