Understanding the Mathematics of the Richards as a limit paper

In Part 1 we built the toolkit on finite matrices: a Laplacian-like operator with \(\ker = \mathrm{span}\{\mathbf{1}\}\), the Fredholm alternative as the source of macroscopic equations, and the pseudo-inverse as the source of transport coefficients. Now we let the matrix indices become continuous and watch the same algebra, verbatim, derive Richards' equation. This post is a plain-language companion to the second of our two Physical Review E manuscripts, where Richards' equation is obtained as the hydrodynamic (Chapman–Enskog) limit of a kinetic theory of unsaturated soil water.

From nodes to pore classes

In the kinetic theory the state of the soil water at a macroscopic point \(x\) and time \(t\) is not a single number \(\theta(x,t)\) but a whole distribution: the pore-occupancy function \(g(r, x, t)\), telling us how water is apportioned among pore classes of radius \(r\). The water content is recovered as a moment,

$$ \theta(x,t) = \int g(r,x,t)\, \mu(dr), $$

with \(\mu\) the pore-size measure of the medium. The “nodes” of Part 1 have become the continuum of pore radii \(r\); a “vector” \(\mathbf{f}\) has become a function \(f(r)\); the dot product has become an integral, \(\langle f, h \rangle = \int f(r)\, h(r)\, \mu(dr)\). Nothing else changes.

Water is exchanged between pore classes — capillary rearrangement, film flow, local equilibration — and this exchange is encoded by a linear(ized) operator \(\mathcal{I}\) acting on functions of \(r\). Schematically, and up to the details spelled out in the papers,

$$ (\mathcal{I} f)(r) = \int \kappa_s(r, r')\, \big[ f(r') - f(r) \big]\, \mu(dr'), $$

with a symmetric pair conductance built as the harmonic mean of the single-class conductances,

$$ \kappa_s(r, r') := \frac{2\, \kappa(r)\, \kappa(r')}{\kappa(r) + \kappa(r')}. $$

Compare this with \((L\mathbf{f})_i = \sum_j A_{ij}(f_j - f_i)\) (sign flipped): \(\mathcal{I}\) is a weighted graph Laplacian on a continuum of nodes, with \(-\mathcal{I}\) playing the role of \(L\). The harmonic mean is not decoration — it is the series-resistor composition rule of Part 1, guaranteeing that exchange between two classes is throttled by the less conductive of the two, and it makes \(\kappa_s\) manifestly symmetric, hence \(\mathcal{I}\) self-adjoint.

The three properties, revisited

Every structural fact from Part 1 now reappears with physical meaning attached.

(i) \(\ker \mathcal{I} = \mathrm{span}\{\mathbf{1}\}\) — one collision invariant. The exchange operator annihilates constants because pairwise exchange conserves total water: \(\int (\mathcal{I}f)\, \mu(dr) = 0\) identically, by the antisymmetry of the integrand. In the Boltzmann theory of gases the collision operator has a five-dimensional kernel (mass, three momenta, energy), and the hydrodynamic limit correspondingly produces five balance equations — the compressible Euler/Navier–Stokes system. In soil water the exchange between pore classes conserves only mass: the kernel is one-dimensional, and the hydrodynamic limit produces exactly one balance equation. That equation is Richards'. The dimension of a kernel dictates the size of your macroscopic PDE system — I find this one of the cleanest structural insights the kinetic viewpoint offers.

(ii) \(\mathcal{I}\) is self-adjoint and negative semi-definite — an H-theorem. The continuum version of the sum-of-squares identity of Part 1 reads

$$ \langle f, \mathcal{I} f \rangle = -\tfrac{1}{2} \iint \kappa_s(r,r')\, \big[f(r') - f(r)\big]^2\, \mu(dr)\, \mu(dr') \ \le\ 0, $$

with equality iff \(f\) is constant (on a “connected” pore network, in the sense that \(\kappa_s\) does not decompose the pore space into non-communicating blocks). Exchange strictly dissipates any non-uniformity: this is the H-theorem of the model, and the reason equilibrium exists and is unique. The equilibrium itself is a packing state — pores fill in order of capillary strength, a Fermi-sea-like picture — but for the linear algebra all we need is that fluctuations around it relax under a self-adjoint, negative semi-definite \(\mathcal{I}\).

(iii) Spectral gap — the small parameter exists. Because the zero eigenvalue is simple and isolated, there is a gap \(\lambda_1 > 0\), hence a fastest-conserved and slowest-decaying separation of time scales. Its ratio to the macroscopic time defines the Damköhler number, \(\mathrm{Da} := \tau_{eq}/\tau_{mac}\), and \(\mathrm{Da} \ll 1\) is the regime in which a hydrodynamic description can be honest. When \(\mathrm{Da}\) is not small — coarse structured soils, preferential flow, rapid forcing — the fast modes never fully slave to \(\theta\) and Richards' equation degrades. The linear algebra even tells you how it degrades: through the modes just above the gap.

The Chapman–Enskog march, order by order

Write the kinetic equation schematically as

$$ \partial_t g + (\text{transport in } x) = \frac{1}{\mathrm{Da}}\, \mathcal{I} g, $$

and expand \(g = g^{(0)} + \mathrm{Da}\, g^{(1)} + \dots\). The algebra of Part 1 now executes itself.

Order \(\mathrm{Da}^{-1}\): \(\mathcal{I} g^{(0)} = 0\), so \(g^{(0)}\) lies in the kernel — it is the local equilibrium distribution, parametrized by the single conserved moment \(\theta(x,t)\). The population of pore classes is enslaved to the water content.

Order \(\mathrm{Da}^{0}\): an equation of the form \(\mathcal{I} g^{(1)} = \mathcal{S}[g^{(0)}]\), with \(\mathcal{S}\) collecting the transport terms. This is exactly the singular problem \(L\mathbf{u} = \mathbf{b}\) of Part 1. The Fredholm alternative demands \(\langle \mathbf{1}, \mathcal{S} \rangle = 0\): projecting the transport terms onto the conserved direction. That projection is the continuity equation,

$$ \partial_t \theta + \nabla \cdot \mathbf{q} = 0. $$

The macroscopic balance law is not assumed; it is the solvability condition of a singular linear problem.

The flux and the conductivity: granted solvability, the correction is \(g^{(1)} = \mathcal{I}^{+} \mathcal{S}\) — the pseudo-inverse at work — and inserting it into the flux moment yields a Buckingham–Darcy law, \(\mathbf{q} = -K(\theta)\, \nabla (\psi(\theta) + z)\)-type, in which the hydraulic conductivity emerges as a bracket of the form

$$ K \ \sim\ -\,\big\langle \Phi,\ \mathcal{I}^{+} \Phi \big\rangle, $$

with \(\Phi(r, r') = -\Phi(r', r)\) the antisymmetric driving potential of the exchange. Readers of Part 1 will recognize the structure immediately: it is the Green–Kubo / effective-resistance formula, the continuum sibling of \(R_{ij}\) built from \(L^+\). The conductivity of a soil is, quite literally, the inverse Kirchhoff-type resistance of its pore-class network, evaluated on the mode forced by gravity and capillarity. This is where the empirical shapes of \(K(\theta)\) — Mualem, van Genuchten and relatives — acquire the status of approximations to a spectral object.

The variational subtlety, honestly told. In an earlier version of the manuscript we characterized this bracket by a single-field quadratic functional, Cercignani-style. That is legitimate when the bracket is a genuine quadratic form \(\langle \chi, \mathcal{I}\chi \rangle\); it silently fails when the object of interest is a bilinear pairing between two different functions. The fix is a two-field, primal–dual functional \(\mathcal{I}[\tilde\chi, \tilde\chi^*]\), stationary in each argument separately, whose stationary value returns the bilinear bracket without smuggling in a false symmetry. On a 3×3 matrix this distinction is invisible; in function space it decides whether your variational bound on \(K(\theta)\) is a theorem or wishful thinking. I mention it because it is a good example of how the finite-dimensional intuition of Part 1, taken too casually, can bite.

Coda: why bother

One can use Richards' equation for a lifetime without this machinery. The point of the derivation is not to re-obtain a 1931 result; it is that every object in the equation now has an address. \(\theta\) is the kernel coordinate; the continuity equation is a Fredholm solvability condition; \(K(\theta)\) is a pseudo-inverse bracket over the pore-class network; and the validity of the whole enterprise is a statement about a spectral gap, quantified by \(\mathrm{Da}\). When the equation fails — and we all know soils where it does — the linear algebra tells you which assumption broke, and what the next term in the expansion looks like.

References and further reading

• S. Chapman and T. G. Cowling, The Mathematical Theory of Non-Uniform Gases, 3rd ed., Cambridge University Press, 1970. (The original Chapman–Enskog method.)
• C. Cercignani, The Boltzmann Equation and Its Applications, Springer, 1988. (Linearized collision operator, Fredholm alternative, variational principles for transport coefficients — the template we adapt from gases to soils.)
• H. Grad, “Asymptotic theory of the Boltzmann equation,” Physics of Fluids 6:147–181, 1963. (The role of the collision invariants and the hydrodynamic projection.)
• F. Golse, “The Boltzmann equation and its hydrodynamic limits,” in Handbook of Differential Equations: Evolutionary Equations, Vol. 2, Elsevier, 2005. (A modern, rigorous survey of hydrodynamic limits.)
• L. A. Richards, “Capillary conduction of liquids through porous mediums,” Physics 1:318–333, 1931.
• E. Buckingham, “Studies on the movement of soil moisture,” USDA Bureau of Soils Bulletin 38, 1907.
• Y. Mualem, “A new model for predicting the hydraulic conductivity of unsaturated porous media,” Water Resources Research 12:513–522, 1976; M. Th. van Genuchten, Soil Science Society of America Journal 44:892–898, 1980. (The empirical \(K(\theta)\) shapes reinterpreted here as spectral approximations.)
• [PRE-1] R. Rigon et al., A kinetic theory of unsaturated soil water — manuscript, Physical Review E (submitted). (insert final title/DOI)
• [PRE-2] R. Rigon et al., Richards' equation as a Chapman–Enskog hydrodynamic limit — manuscript, Physical Review E (submitted). (insert final title/DOI)

Some about the Linear Algebra used in the Richards as a limit paper

This is the first of two posts written as a gentle companion to our recent work on the kinetic theory of unsaturated soil water and, in particular, to the derivation of Richards' equation as a Chapman–Enskog hydrodynamic limit. Before touching any soil physics, I want to isolate the linear algebra that makes the derivation work. It turns out to be the linear algebra of one very special family of matrices — Laplacians — together with the two classical tools that let us handle their singularity: the Fredholm alternative and the pseudo-inverse. If you understand a 3×3 example, you understand the skeleton of the whole paper.

A matrix from a picture

Take the simplest possible network: three nodes in a line,

1 — 2 — 3.

Two matrices encode this picture. The adjacency matrix \(A\) records who is connected to whom (\(A_{ij} = 1\) if \(i\) and \(j\) share an edge, \(0\) otherwise), and the degree matrix \(D\) is diagonal, with \(D_{ii}\) counting the edges at node \(i\). Their difference

$$ L := D - A $$

is the graph Laplacian. For our three-node line:

$$ L = \begin{pmatrix} 1 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 1 \end{pmatrix} $$

Now assign a number \(f_i\) to each node — think of it as a water content, a temperature, a head. Then a one-line computation shows

$$ (L\mathbf{f})_i = \sum_{j \sim i} (f_i - f_j) $$

where \(j \sim i\) runs over the neighbors of \(i\). The Laplacian measures how much each node differs from its neighborhood. It is the discrete cousin of \(-\nabla^2\): a peak gives a positive value, a valley a negative one, and a node in balance with its surroundings gives zero. This is why the discrete heat (diffusion) equation reads

$$ \frac{d\mathbf{f}}{dt} = -L\mathbf{f} $$

and drives any initial condition toward uniformity.

The three properties that matter

Everything we will need in Part 2 follows from three elementary facts.

(i) The rows sum to zero. By construction, each diagonal entry exactly cancels its off-diagonal row. Consequently the constant vector \(\mathbf{1} = (1, 1, \dots, 1)^{\mathsf T}\) satisfies \(L\mathbf{1} = \mathbf{0}\): the Laplacian forgets constants. In dynamical language, uniform states are equilibria — once everything is level, diffusion stops. In conservation language: if you sum the components of \(L\mathbf{f}\), you get zero, so the total quantity \(\sum_i f_i\) is conserved by the dynamics. Keep this pairing in mind — kernel of the operator ↔ conserved quantity — because it is the leitmotiv of the whole series.

(ii) \(L\) is symmetric and positive semi-definite. For any \(\mathbf{f}\),

$$ \mathbf{f}^{\mathsf T} L \, \mathbf{f} = \sum_{\text{edges } (i,j)} (f_i - f_j)^2 \ \ge\ 0, $$

a sum of squares over the edges. This little identity is the discrete Dirichlet form; it says the Laplacian measures the total “roughness” of \(\mathbf{f}\) over the network. Symmetric + semi-definite means all eigenvalues are real and non-negative, \(0 = \lambda_0 \le \lambda_1 \le \dots \le \lambda_{n-1}\), and the eigenvectors form an orthogonal basis. (In the papers we work with the sign flipped: our exchange operator is negative semi-definite, \(\mathcal{I} = -L\) morally, so that it dissipates. The content is identical.)

(iii) The zero eigenvalue counts connected components. If the graph is connected, \(\lambda = 0\) is simple: the kernel is exactly \(\mathrm{span}\{\mathbf{1}\}\), nothing more. If the graph splits into two islands, the kernel is two-dimensional — spanned by the indicator of each island, e.g. \((1,1,0,0)^{\mathsf T}\) and \((0,0,1,1)^{\mathsf T}\) — because each island can sit at its own uniform level with no communication between them. Connectivity is thus an algebraic statement:

$$ \ker L = \mathrm{span}\{\mathbf{1}\} \iff \text{the network is one connected piece}. $$

For a hydrologist this should already ring a bell: it is the algebraic shadow of the percolation question — does the wet pore space form a single connected cluster through which pressure can equilibrate?

Eigenvectors as standing waves

Decompose any state on the eigenvector basis, \(\mathbf{f}(0) = \sum_i c_i \mathbf{v}_i\). Under diffusion each mode evolves independently:

$$ \mathbf{f}(t) = \sum_i c_i\, e^{-\lambda_i t}\, \mathbf{v}_i. $$

The eigenvectors are the standing waves of the network, and the eigenvalues are their decay rates. The constant mode (\(\lambda_0 = 0\)) never decays — it is the conserved total. The first non-trivial mode, associated with \(\lambda_1\) (the Fiedler eigenvalue), is the slowest transient: it is smooth, splits the network into a positive and a negative half, and encodes the network's most reluctant global imbalance. High eigenvalues correspond to jagged, rapidly oscillating modes that diffusion flattens almost instantly.

This separation of time scales — one frozen mode, and everything else decaying at rates bounded below by \(\lambda_1 > 0\) (the spectral gap) — is exactly the structure that a Chapman–Enskog expansion exploits. Fast modes slave themselves to the slow one; the slow one becomes the hydrodynamic field. In our soil-water papers the “slow mode” is the locally conserved water content, and the spectral gap is what allows the small parameter \(\mathrm{Da} := \tau_{eq}/\tau_{mac}\) to exist at all.

Inverting the non-invertible: the problem

A kernel encodes a conservation law — but a kernel also means the matrix is singular. Suppose we want to solve the steady-state problem

$$ L\mathbf{u} = \mathbf{b}, $$

where \(\mathbf{b}\) is some forcing (a source/sink pattern, an injection of water or heat at the nodes). Since \(L\mathbf{1} = \mathbf{0}\), the matrix has no inverse: you cannot undo the crushing of the constant direction. Two things can go wrong, and they are dual to each other. The two classical tools that resolve the impasse — the Fredholm alternative and the pseudo-inverse — are not mere technicalities: in a Chapman–Enskog derivation they are precisely the steps that produce the macroscopic balance equation and the transport coefficient.

The Fredholm alternative: when does a solution exist?

For a symmetric operator, the range is the orthogonal complement of the kernel. So \(L\mathbf{u} = \mathbf{b}\) has a solution if and only if \(\mathbf{b} \perp \ker L\), i.e.

$$ \mathbf{b} \perp \mathbf{1} \iff \sum_i b_i = 0. $$

Physically: a connected, closed system can reach a steady state only if the forcing is globally balanced — what is pumped in somewhere must be extracted somewhere else. You cannot pour net water into a sealed, connected pore network and expect a stationary pressure field. And when a solution exists, it is unique only up to an additive constant \(c\mathbf{1}\): steady states are defined modulo a uniform offset, exactly as potentials are.

This solvability condition is the unglamorous hero of the Chapman–Enskog method. At each order of the expansion one has to solve an equation of the form “singular operator applied to unknown correction = known inhomogeneity”, and the demand that the inhomogeneity be orthogonal to the kernel is what spits out the macroscopic equation. In gas kinetics the orthogonality to the collision invariants (mass, momentum, energy) yields the Euler and Navier–Stokes equations. In our soil-water setting the kernel is one-dimensional — water mass is the only invariant — and the solvability condition yields the continuity equation for water content. We will see this in Part 2.

The pseudo-inverse: how to write down the solution

Granted \(\mathbf{b} \perp \mathbf{1}\), how do we express \(\mathbf{u}\)? Restrict attention to the subspace of zero-mean vectors, \((\ker L)^\perp = \{\mathbf{x} : \sum_i x_i = 0\}\). On this subspace \(L\) is strictly positive definite, hence genuinely invertible. The Moore–Penrose pseudo-inverse \(L^+\) implements this restricted inverse and extends it by zero on the kernel. Spectrally, if

$$ L = \sum_{i \ge 1} \lambda_i\, \mathbf{v}_i \mathbf{v}_i^{\mathsf T} \qquad (\text{the } \lambda_0 = 0 \text{ term absent}), $$

then

$$ L^+ = \sum_{i \ge 1} \frac{1}{\lambda_i}\, \mathbf{v}_i \mathbf{v}_i^{\mathsf T}, $$

and the zero-mean solution is \(\mathbf{u} = L^+ \mathbf{b}\). Notice the weighting: the slow, smooth, low-\(\lambda\) modes dominate \(L^+\), because \(1/\lambda\) is largest for them. Inversion amplifies exactly the structures that diffusion is most reluctant to erase. In the continuum, \(L^+\) becomes the Green's operator: the solution operator whose integral kernel is the Green's function of the Laplacian, i.e. the potential generated by a balanced distribution of sources.

Variational detour: the pseudo-inverse as a minimization

There is a second, equivalent way to characterize \(\mathbf{u} = L^+\mathbf{b}\) that deserves its own paragraph, because it is the route we ultimately follow in the second PRE paper. The solution of \(L\mathbf{u} = \mathbf{b}\) on \((\ker L)^\perp\) is the stationary point of the functional

$$ \mathcal{J}[\mathbf{u}] = \tfrac{1}{2}\,\mathbf{u}^{\mathsf T} L\, \mathbf{u} - \mathbf{b}^{\mathsf T}\mathbf{u}, $$

and the stationary value \(-\tfrac12 \mathbf{b}^{\mathsf T} L^+ \mathbf{b}\) directly evaluates the quadratic form of the pseudo-inverse. One caveat, learned the hard way: when the operator is applied between two different vectors — a bilinear bracket \(\mathbf{a}^{\mathsf T} L^+ \mathbf{b}\) with \(\mathbf{a} \ne \mathbf{b}\) — a single-field quadratic functional no longer suffices, and one must resort to a two-field (primal–dual) functional whose stationarity in each argument reproduces the bracket. This distinction between a quadratic form and a genuine bilinear form looks pedantic on a 3×3 matrix; in function space it is the difference between a correct and an incorrect variational bound on the hydraulic conductivity. (More on this in Part 2.)

What the entries of \(L^+\) mean

The pseudo-inverse is not just a formal device; on a network its entries are measurable quantities.

Treat every edge as a unit resistor. Then the effective resistance between nodes \(i\) and \(j\) is

$$ R_{ij} = (L^+)_{ii} + (L^+)_{jj} - 2 (L^+)_{ij}, $$

a genuine distance on the graph (the resistance distance of Klein and Randić), proportional to the mean commute time of a random walker between \(i\) and \(j\). The diagonal entry \((L^+)_{ii}\) measures how peripheral node \(i\) is with respect to the whole network — large for dead-ends, small for well-embedded hubs — and, if one reads \(L^+\) as the covariance of a Gaussian free field, it is the variance of the field at node \(i\). The trace,

$$ \mathrm{Tr}(L^+) = \sum_{i \ge 1} \frac{1}{\lambda_i}, $$

is proportional to the Kirchhoff index, a single number summarizing the global transport capacity of the network. For a hydrologist the dictionary is irresistible: pore networks as resistor networks, hydraulic conductance in place of electrical conductance, and \(L^+\) as the object that converts local pore-scale conductances into global, geometry-aware transport coefficients. That conversion is exactly what a Chapman–Enskog closure does.

The take-home diagram

$$ \text{kernel} \;\Rightarrow\; \text{conservation law} \;\Rightarrow\; \text{solvability condition (macroscopic equation)} $$ $$ \text{pseudo-inverse on } (\ker)^\perp \;\Rightarrow\; \text{Green's operator} \;\Rightarrow\; \text{transport coefficient}. $$

In Part 2 we let the nodes become a continuum of pore classes and watch this diagram become, line by line, the derivation of the Buckingham–Darcy flux and of Richards' equation.

References and further reading

• F. R. K. Chung, Spectral Graph Theory, CBMS Regional Conference Series in Mathematics 92, AMS, 1997. (The standard reference; the normalized-Laplacian point of view.)
• B. Mohar, “The Laplacian spectrum of graphs,” in Graph Theory, Combinatorics, and Applications, Wiley, 1991. (A very readable survey of the properties used here.)
• M. Fiedler, “Algebraic connectivity of graphs,” Czechoslovak Mathematical Journal 23:298–305, 1973. (Where \(\lambda_1\) got its name.)
• G. Strang, Introduction to Linear Algebra, Wellesley–Cambridge Press. (For the quadratic-form and spectral-theorem background; Strang is fond of exactly our 3×3 example.)
• R. Courant and D. Hilbert, Methods of Mathematical Physics, Vol. I, Interscience, 1953. (The continuous ancestor of everything above, including the variational characterization of eigenvalues we will meet again in Part 2.)
• E. B. Davies, Linear Operators and their Spectra, Cambridge University Press, 2007. (Self-adjointness, spectral theorem, Fredholm theory in one place.)
• D. J. Klein and M. Randić, “Resistance distance,” Journal of Mathematical Chemistry 12:81–95, 1993.
• A. Ghosh, S. Boyd, A. Saberi, “Minimizing effective resistance of a graph,” SIAM Review 50(1):37–66, 2008. (Kirchhoff index, variational characterizations, optimization view.)
• P. G. Doyle and J. L. Snell, Random Walks and Electric Networks, MAA, 1984 (freely available on arXiv). (The most enjoyable introduction to the resistor-network picture.)
• C. Cercignani, The Boltzmann Equation and Its Applications, Springer, 1988. (Chapter on the linearized collision operator: Fredholm alternative used exactly as above.)

Richards as a limit: a derivation of Richards' equation from the Continuum Kinetic soil Equation

Two quantities that have no business being related: the spectral gap of a pore network — a purely structural number, computed from the connectivity, with no flow solved anywhere — and the Stokes permeability of the same network, computed by actually solving viscous flow under a pressure drop. Drain the network step by step and they vanish at the same water content. Not approximately. Identically, both zero, at the same θ. That coincidence is the subject of this post: it is the point where a macroscopic constitutive law stops existing, and the theory says so by itself.

In July I posted the first of the two papers — the kinetic theory of the pore-occupancy g(r) — and promised the companion "shortly." Here it is:

Richards' equation as a hydrodynamic limit: Chapman–Enskog reduction of the continuum kinetic equation for unsaturated soil water
R. Rigon, arXiv:XXXX.XXXXX [physics.flu-dyn], 24 pages, 1 figure, 8 appendices, with numerical Supplemental Material.

It goes to Physical Review E, like its companion. The first paper said: θ is not enough, the state is g(r). This one says the other half, and it is the half that makes the first one respectable.

The argument in one sentence

Richards' equation is not an assumption of soil physics. It is a theorem — the solvability condition of a kinetic equation whose redistribution operator has exactly one invariant.

That is the whole paper. Everything below is a consequence.

Two limits that hydrology has always taken together

The reason Richards' equation has never been derived, only motivated, is that the passage from pores to fields conflates two entirely different operations. The paper's first move is simply to take them apart.

• The spatial limit, ε = L/Λ → 0, shrinks the representative volume to a point. It is pure kinematics: it says nothing whatsoever about time scales. What comes out is a closed continuum kinetic equation, ∂g/∂t + ∇·F = 𝒞[g], with F a pore-resolved flux that is still entangled — a transport coefficient and a driving gradient multiplied together inside one kernel, inseparable.
• The temporal limit is where the physics is. Redistribution among pores is fast; the forcing is slow; their ratio is the Damköhler number Da. When Da ≪ 1 the soil is pinned near local equilibrium, and one can expand — a Chapman–Enskog reduction, exactly as one passes from Boltzmann to Navier–Stokes.

Keeping them apart is what makes the structure visible, and it is why the paper can be honest about where each classical assumption enters.

The one operator everything depends on

Linearise the redistribution operator about equilibrium and you get 𝒥, and then the whole reduction is a statement about 𝒥. Three properties, and nothing else, do all the work:

1. 𝒥 is self-adjoint — in the mass inner product, ⟨u,v⟩ = ∫ u v f dr. And here is the small piece of algebra I find most satisfying in the paper: self-adjointness and conservation of water are literally the same statement. Detailed balance K(r,r′) f(r) = K(r′,r) f(r′) holds identically because the mobility is symmetric, and that single fact gives you both.
2. Its kernel is one-dimensional. There is exactly one thing redistribution cannot change: water. The Boltzmann gas has five invariants and therefore five macroscopic equations. Viscous pore flow has one invariant, and therefore one macroscopic equation.
3. It has a spectral gap — a slowest relaxation rate λ₁ > 0 — which is what makes the expansion asymptotic at all.

Then the derivation is almost mechanical. Take the f-moment of the first-order equation: redistribution drops out (that is the null space), and what remains is a condition on the source. That condition — the Fredholm alternative, the unglamorous requirement that the first-order correction should exist at all — is mass conservation. It is Richards' equation.

Five things that fall out, which I did not put in

1. Hydraulic conductivity is a transport coefficient. K is not a constitutive input. It is the first-order Chapman–Enskog coefficient, K = φ⟨κ|𝒥⁻¹|S⟩ — the exact structural counterpart of viscosity in the kinetic theory of gases. Which means, among other things, that K is a property of the medium's relaxation spectrum and is independent of the forcing, in the same sense that viscosity does not depend on the shear rate you apply.

2. The classical formulas are approximations of that coefficient. Take the mean field of it and you get back the standard Mualem–Burdine integral. Resum the serial paths — the fact that large pores must push through small-pore bottlenecks — and out comes Mualem's heterogeneity penalty exp(−4σ²), not as an empirical factor but as a Neumann series that collapses to a harmonic mean. I did not expect that to work as cleanly as it did.

3. Dual-permeability models are derived, not posited. This is the result with the widest practical reach. Give the operator a bimodal pore-size distribution and its relaxation spectrum splits into two bands separated by a gap. Apply Chapman–Enskog within each band, keep the slow cross-band relaxation, and two coupled Richards equations fall out — Gerke–van Genuchten, with Weiler's IN3M as the three-band case and mobile–immobile as the limit where one band is conductively dead. And the exchange coefficient Γ_w, which everybody fits, is computed from the cross-band connectivity. What was a modelling choice becomes a theorem with a formula.

4. The theory predicts its own breakdown — by two different routes. Raise the forcing and Da → 1: sharp fronts, the expansion fails, and you are in the preferential-flow regime the first paper described. But lower the water content and something else happens: at the percolation threshold the spectral gap closes, ‖𝒥⁻¹‖ diverges, and the closed conductivity ceases to exist. These are genuinely two different exits from the Richards regime — one by fast forcing, one by loss of connectivity — and the theory locates both. That is the figure at the top: the gap and the permeability going to zero together.

5. Hysteresis and dynamic capillarity have an operator address. The non-commutativity [𝒲,𝒟] ≠ 0 from the first paper turns out, at the operator level, to be the curvature of the projection onto local equilibrium — the same object that carries dynamic capillarity. I am not claiming to have solved hysteresis. I am claiming to know where in the mathematics it lives.

The part I am least able to hide behind

Point 4 above is the paper's most exposed claim, so it is the one I made numerical. The Supplement does three parameter-free computations, on soils I can name:

• A loam (unimodal, median 8 μm, draining around −1.9 m). Diagonalise the operator: one exact invariant (λ₀ ≈ 10⁻¹⁷), and the modes turn out to be localized by pore radius, each relaxing at the local rate of its own radius, to within 1%. They are not standing waves on the band — I had assumed they were, wrote it into an early draft, and the numerics said no. Ordering the modes by rate is ordering the pores by size: slow modes on small pores.
• The same loam, drained on a 24³ pore network. Gap and permeability vanish together at θ_c ≈ 3.5×10⁻³. Below it, no spanning cluster: the water is there, it simply cannot go anywhere.
• A structured soil (matrix at 4 μm plus macropores at 40 μm). The spectrum splits into two bands with a slow, sign-changing mode across the split — the exchange mode, appearing on its own. And projecting the operator onto the two bands gives an exchange rate that converges to the true one as the bands decouple. The split radius the operator produces falls at ψ ≈ −1 m ≈ −10 kPa, which is where soil physicists have been drawing the macropore boundary by hand for decades.

Where this connects

Nothing here replaces anything. Capillary-bundle models are the diagonal limit; Mualem and Burdine are the mean field; critical-path analysis is the spectral limit near θ_c; Gerke–van Genuchten is the two-band projection. The classical results are not overturned — they are located, each one identified as a particular approximation to a single operator inversion. That is the most useful thing a derivation can do for a field: not to declare the old results wrong, but to say precisely what they are approximations of, and therefore when they will fail.

What I am not claiming

Again, better said by me than to me.

This is a formal Chapman–Enskog reduction, in the sense the phrase carries in kinetic theory: the expansion is organised in powers of Da and closed order by order, but I do not prove convergence, and the higher-order remainders are not bounded rigorously. The closures come from the companion paper and are physical, not derived from molecular dynamics. The numerics are on synthetic networks, not on imaged soils. And the linearisation that makes 𝒥 an operator at all is exactly that — a linearisation, with the nonlinearity pushed into successive sources.

There is also one honest limitation inside the numerics that I have written into the Supplement rather than left for a referee to find: away from the threshold, the spectral gap of a growing cluster is increasingly dominated by its size rather than its connectivity, so the correlation between gap and permeability is only meaningful near θ_c. What is unambiguous — and all the argument needs — is that they vanish together, and that is exact rather than statistical.

Materials

• Paper: arXiv:XXXX.XXXXX (link when it lands) (The provisional pdf here until acceptance on arXiv)
• Companion paper: arXiv:2607.09416 — the kinetic theory this reduces
• The first post: If not Richards, what else ?
• Numerical Supplement (included in the main file) + Jupyter notebook: every number and figure above is reproducible; the notebook runs top to bottom in a few seconds.
• "The Real Book": a companion document I wrote for myself and then decided to keep — the entire derivation worked step by step, blackboard style, nothing skipped, with boxes reminding the reader of the linear algebra (Fredholm alternative, pseudo-inverse, graph Laplacians, Rayleigh quotients) and a glossary. If the paper looks forbidding, start there. It is the gentlest way in.
• On the same line of the Real Book, I also provide a little rehearsal on Linear Algebra and Linear operators, whose knowledge is necessary to the understanding of the paper calculations. 

As always: comments, objections and counterexamples are welcome — especially the counterexamples. This paper makes a falsifiable structural claim (Γ_w computable from connectivity, gap closure at θ_c) and I would rather find out early.

Water In Soil - A MOOC

I'm happy to share that my new MOOC, SOIL – The Hydrology of Soil, is now live on the University of Trento's MOOC platform. It's free, open, and self-paced, and it's aimed at anyone who wants to properly understand how water moves through unsaturated soil — not just as a set of formulas to memorize, but as a coherent chain of physical reasoning.

What the course is about

The course covers the hydrology of unsaturated soils and builds up, step by step, the mathematical tools needed to describe water flow through them — culminating in the Richards equation, the cornerstone of unsaturated flow theory.

It's organized into five chapters, each combining short video episodes, further readings, hands-on activities, and self-assessment quizzes:

1. What is Soil — the basic quantities used to describe water content and structure in soil.
2. The Energy of Water in Soil — how water's energy is distributed, capillary pressure, and the construction of soil-water retention curves.
3. Darcy and Buckingham's Laws — from Darcy's law in saturated soil to Buckingham's extension for unsaturated conditions, and how hydraulic conductivity is derived.
4. Soil Water Budget — the mass budget, conservation laws, and the derivation of Richards' equation itself (plus its groundwater counterpart).
5. Solving Richards' Equation and Its Limits — pedotransfer functions, numerical solution strategies, macropores, and where the classical theory breaks down.

What you'll be able to do by the end

By the end of the course, you should be able to explain the physical meaning of the water retention curve and the Richards equation and the reasoning that connects them, solve simple unsaturated flow problems by choosing the right constitutive relationships, and critically evaluate when the Richards-equation framework is — and isn't — a valid description of what's happening in real soil.

Who it's for

If you work in hydrology, agronomy, environmental engineering, or soil science — or you're a student who wants to go beyond a black-box use of Richards' equation — this course is built for you. The course can be seen as an introduction to the theory of the WHETGEO model.  WHETGEO provides an open source tool for solving 1D and 2D Richards equation. 

Try it

The course is free to enroll: mooc.unitn.it/course/view.php?id=32

If not Richards, what else ?

Two soils with the same water content θ are not in the same hydraulic state. The pore-occupancy g(r) — the fraction of pores of radius r that are water-filled — distinguishes them, while θ, being an integral of g, cannot. The navy step is the reference equilibrium g_eq = H(r* − r): water fills the small pores first. Everything else on the plot is a state that Richards' equation is blind to.

In May I posted the talk I gave at EGU 2026 in Vienna, and promised the two papers “in a couple of weeks after EGU.” It took a little longer than that — it always does — but the first one is now public:

The Statistical Physics of Unsaturated Soil Water: kinetic theory and non-commutative pore-water dynamics
R. Rigon, arXiv:2607.09416 [cond-mat.stat-mech], 22 pages, 9 figures, 2 appendices.

It is going to be submitted to Physical Review E. The companion paper — the Chapman–Enskog derivation that recovers Richards' equation as a hydrodynamic limit — follows shortly, and I will post it here when it lands.

The argument in one sentence

Unchanged from the talk, and worth repeating because everything else is a consequence of it:

Richards' equation is not wrong; it is the equilibrium limit of a deeper kinetic theory — in the same sense that Navier–Stokes is the hydrodynamic limit of Boltzmann's equation for a gas.

Mario Putti asked me, twenty years ago, “if not Richards, what else?” This is my attempt at an answer, and it arrives only after many years spent trying to solve Richards' equation properly — first with GEOtop, later with WHETGEO. You have to take an equation seriously for a long time before you earn the right to say what it is missing.

What the theory actually says

The state variable is not θ. It is the pore-occupancy g(r, x, t): the fraction of pores of radius r that are water-filled at position x and time t. Water content is recovered as a moment of it, θ[g] = φ ∫ g(r) f(r) dr — which is precisely the point: θ is an integral of g, so it throws information away. Two soils with the same θ, one wetted by rain (which fills pores by areal exposure, favouring the large ones) and one drained to the same θ (which empties the large ones first), are in genuinely different hydraulic states. They will conduct water differently, and they will respond to the next rainfall differently. Richards' equation cannot see the difference. That is the figure above, and that is the whole motivation.

The theory is built by passing through three scales, and I think this is the part hydrologists will find easiest to trust, because each step is ordinary physics:

• Microscale. A single water transfer between two pores is set by a Hagen–Poiseuille rate and driven by the difference of pore chemical potentials Φ(r, r′) — capillary and gravitational here, but open to adsorptive, osmotic, or thermal refinement without touching the structure of the theory.
• Mesoscale. Averaging over a representative volume gives a master equation — a gain–loss (Boltzmann-type) kinetic equation whose terms relax the occupancy toward its equilibrium, with a connectivity kernel C(r, r′) that encodes which pores can actually talk to which.
• Macroscale. A Chapman–Enskog reduction gives back Richards' equation in the quasi-static limit Da → 0.

Everything the theory needs as input is a geometric property of the pore network — measurable from micro-CT. Nothing is calibrated against macroscopic hydrological data. I want to be blunt about how unusual that is, and how exposed it leaves me: the theory makes parameter-free predictions, and parameter-free predictions can be wrong in public.

Four things that fall out, which I did not put in

This is the part I care about. These were not assumptions; they are consequences.

1. Matric potential and hydraulic conductivity exist only in the limit. ψ and K are not primitive quantities of the theory. They emerge at Da → 0, and K is derived from the connectivity kernel rather than postulated. Below the percolation threshold, K vanishes — not as a fitting choice, but because the water phase stops spanning the medium. Field capacity gets a geometric meaning: θ_FC ≈ θ_c.

2. Hysteresis is geometry, not memory. It is the holonomy of a forcing bundle — a geometric phase, arising from the non-commutativity [W, D] ≠ 0 of the wetting and drying operators. Wetting fills by areal exposure; drying empties by capillary ordering; the two operations do not commute, so a closed loop in the forcing does not return you to where you started. Independent-domain and Preisach models posit bistable pores and reproduce the loop. Here the loop is derived, and it comes with a falsifiable prediction: the loop area scales as H ∼ I² with the forcing intensity. Domain models are rate-independent and predict no such thing. That is a clean experimental discriminant, and I would very much like someone to go and measure it.

3. Preferential flow is not a separate process. It is what the same equation does when Da > 1. The molecular-chaos (Stosszahlansatz) closure that underlies the kinetic equation fails exactly when pore occupancies become correlated near the percolation threshold — and that correlated, channelized regime is fingering and preferential flow. So the Richards / preferential-flow dichotomy dissolves into a continuous, Da-controlled crossover. We do not need two domains and a phenomenological exchange term; we need one equation and an honest look at its Damköhler number.

4. Out of the quasi-static limit, g(r) is irreducible. No single scalar — not θ, not ψ — is a complete description. And, as I discovered while revising the manuscript (a lesson in the value of being asked a hard question at the right moment): this is true even at equilibrium. With gravity present in a finite volume, the equilibrium occupancy is not the sharp step H(r* − r) at all; it is a smeared step, because a large pore low in the profile can stay filled while a smaller pore higher up has already drained. Each pore holds water within its own Jurin rise. The retention curve — the last place where the classical scalar picture was supposed to be exact — is not exact either.

Where this connects

The framework absorbs rather than replaces. Capillary-bundle models are its diagonal limit; critical-path models are its spectral limit; Hassanizadeh–Gray is a thermodynamically consistent extension, here resolved pore-class by pore-class; phase-field methods are gradient flow on a free energy, here with explicit network connectivity; dual-permeability models are the Da > 1 regime, without the phenomenology. The same machinery, with capillary pressure replaced by freezing-point depression, is the freezing-soil problem I have worked on with Niccolò Tubini and John Mohd Wani.

This is not a parallel universe to Richards. It contains it.

What I am not claiming

I would rather say this myself than have it said to me. The paper is a construction, not a rigorous reduction from molecular dynamics: the closures are posited on physical grounds and judged by their consequences. The full kinetic equation has not yet been solved on a real soil — the numerics live in the companion paper and in the supplementary demonstrations. And the single most obvious next step is also the hardest and the most interesting one:

directly observing g(r).

Micro-CT can see it. Nobody, as far as I know, has yet used it to test a kinetic theory of soil water. If you work with imaging of pore-scale water and this sounds like a collaboration, write to me.

Materials

• Paper: arXiv:2607.09416
• The EGU 2026 talk (slides, storyboard, and the notebooks behind the figures): the May post — still the gentlest way in, if the paper looks forbidding.
• Code: will be added on GitHub; the OpenPNM notebooks that generate the supporting figures are already in the OSF repository linked from the talk.
• The companion paper: "Richards as a Limit"  where Richards equation is derived from the main general equation. Here.
• My MOOC (Massive Open Online Course) on Water in Soil. This, especially the part related to the energy of water in soil, can be though as an introduction to these more advanced topics. 

Comments, objections, and counterexamples are all welcome — especially the counterexamples. A theory that cannot be attacked is not saying anything.

Extended PETRI Net examples from the MARRMoT models collection

Five years ago, when we wrote that any lumped-parameter hydrological model can be represented as an Extended Petri Net, the statement had the flavor of a theorem asserted with a couple of worked examples. Now it has the flavor of a proof by exhaustion. The presentation below, prepared with Marialaura Bancheri and Anna De Nardi, contains the EPN translation of all forty-six conceptual models of the MARRMoT collection (Knoben et al., 2019), from the single-bucket Collie River Basin 1 up to SACRAMENTO, PRMS and CLASSIC. Anna carried out the bulk of this work as her graduation exercise, completing the task I had assigned to my class back in 2021, and I think the result deserves to be seen in its entirety.

The presentation can be found by clicking on the above image. 

I will not explain here what an EPN is: the definitive references remain Bancheri, Serafin and Rigon (2019), which introduces the formalism and its exact correspondence with the ordinary differential equations of the water budget, and Rigon and Bancheri (2021), which shows how the same topology carries, almost for free, the travel time, response time and tracer dynamics. A gentler entry point, with slides and a video tutorial, is this older post, and the conceptual background on the equivalences among the various hydrological dynamical systems is discussed here.

What the collection adds is something the papers could not give: the experience of seeing forty-six models side by side under a single graphical grammar. Some things become obvious that the original equations, or the traditional bucket sketches, keep hidden. Family resemblances jump out — the MOPEX series, the Flex variants, the Tank models reveal themselves as small mutations of a shared skeleton, and one starts to suspect that the space of conceptual models is much smaller than the number of their names suggests. Complexity becomes measurable at a glance: you can literally count places, transitions, splitters and see where a model concentrates its assumptions. And the pathologies show up too — when the wiring of a model resists a planar, readable drawing, as it happens with SACRAMENTO, that is telling you something about the model, not about the drawing. In this sense the EPN works as a diagnostic instrument, not merely an illustration.

There is also a forward-looking reason to care. Once a model is a graph with typed nodes, it is data: it can be stored, compared, composed with other graphs, translated automatically into code — which is what we pursue in the GEOframe/OMS3 world — and it connects naturally with the compositional, category-theoretic view of open systems I discussed apropos of stock and flow diagrams. The forty-six drawings below are therefore not an endpoint but a dataset.

All the previous material on the topic is collected under the EPN label of this blog, starting from the original announcement of the WRR paper.

References

Bancheri, Marialaura, Francesco Serafin, and Riccardo Rigon. 2019. "The Representation of Hydrological Dynamical Systems Using Extended Petri Nets (EPN)." Water Resources Research 55 (11): 8895–8921. https://doi.org/10.1029/2019WR025099

Rigon, Riccardo, and Marialaura Bancheri. 2021. "On the Relations between the Hydrological Dynamical Systems of Water Budget, Travel Time, Response Time and Tracer Concentrations." Hydrological Processes 35 (1). https://doi.org/10.1002/hyp.14007

Knoben, W. J. M., J. Freer, K. J. A. Fowler, M. C. Peel, and R. A. Woods. 2019. "Modular Assessment of Rainfall-Runoff Models Toolbox (MARRMoT) v1.0: An Open Source, Extendable Framework Providing Implementations of 46 Conceptual Hydrologic Models as Continuous Space-State Formulations." Geoscientific Model Development. https://gmd.copernicus.org/articles/12/2463/2019/

What is a DARTH, and which is its way to implement a digital twin of Hydrology ?

So what, exactly, is a DARTH? Not the software, not the cloud — the thing itself. This post is my attempt to answer, in thirteen tenets.

Digital twins are everywhere now, at least we see several contributions that name themselves such. They appear in the great international programmes, in the roadmaps of agencies, in the slides of almost every keynote about the future of Earth science. Many of them are magnificent feats of engineering: high-resolution solvers, elastic cloud, fast emulators. And yet, watching talk after talk, I kept feeling that something was missing. Technology is not the same thing as a vision. A high-resolution engine and a cloud are not, by themselves, a scientific instrument.

The precedent I keep returning to is FAIR. Years ago the FAIR principles changed the way our community treats data — not by inventing any new technology, but by writing down, compactly and citably, what good stewardship requires: that data be findable, accessible, interoperable, reusable. Their power was in the act of definition. We need something similar for twins. But a twin is neither only data nor only software. It is a model bound to a real basin, kept in correspondence with it, and used to make claims about its past, present, and possible futures. So I sat down and tried to write what makes such a twin trustworthy.

I came out with seven points. Talking them over with colleagues, the seven became thirteen, and they arranged themselves naturally into five movements — from what a twin is, through the record it keeps and the machinery that runs it, to the commons that builds it and the ends it serves.

I · Representation — what a twin is

1. A multi-resolution, multi-process, representation of the hydrological cycle, built from interchangeable modelling solutions that close mass and energy budgets and quantify their own error — at least in a subjective, Bayesian sense — propagated across scales.

2. Bound to a real system: calibrated, validated, and — where data allow — kept synchronised with the basin through data assimilation. This binding to reality is what distinguishes a twin from a mere simulator.

3. Transparent — no black boxes. Every component is inspectable, and even the machine-learning parts carry a stated structure and an interpretable role.

II · Data — the record it keeps

4. All data, consumed and produced, are geo-registered, versioned, and FAIR.

5. Every result is reproducible end to end — code, parameters, environment, and data lineage recorded as provenance. FAIR governs the data; this governs the computation.

III · Software — the machinery that runs it

6. Model parts are interchangeable by construction, so alternative descriptions can compete as falsifiable hypotheses.

7. Deployable over the web and, through shared standards and interfaces, interoperable across infrastructures — components and twins from different groups composing together, running from laptop to HPC.

8. Open source: built on, released as, and developed in the open, and free of lock-in.

IV · Community — the commons that builds it

9. A shared infrastructure, not isolated codes, that serves the widest community and grows with knowledge.

10. Organised, by construction, for participation and cooperation — for the collective action of scientists.

11. Stewarded: it carries its own documentation, training, and transparent governance, so that it outlives any single project or generation of its makers.

V · Purpose — the ends it serves

12. General-purpose by design — prediction, scenarios, decisions, education, or pure inquiry alike.

13. Bounded by an ethical purpose: to serve the public good and the stewardship of water and land, and to do no harm.

Please observe that Number 1, excludes practically all exiting model, which are certainly not multi-resolution, neither multiprocess. Therefore, for now, we have to consider DARTHs as an "organizing metaphor", i.e. a working program.  If you make me choose the tenets I care most about, they are these. Number 2, the binding to observation: a model that never meets a measurement is a beautiful animation, not a twin. Number 3, no black boxes: I am not against machine learning — I use it — but a learned component still has to say what it is and what role it plays. Number 7, interoperability among infrastructures: no platform is a twin of the planet on its own; the digital Earth appears only when many of them interoperate. One Earth, many infrastructures. And number 13, the ethical bound, which I added last and on purpose — openness and transparency are not only engineering choices, they are how we stay accountable.

Here is the part I most want to insist on: these are not wishes. The means to honour most of them already exist, in the open. The Basic Model Interface (CSDMS) makes components interchangeable and couplable; the OGC API – Processes family lets us deliver and chain models as web services across infrastructures; and the Object Modeling System with its cloud companion CSIP turns model components into open services at scale. GEOframe and GEOtop are simply our own way of walking that path. What no standard hands you for free is the last rung — transparency, honest uncertainty, components that behave as hypotheses rather than as fitted curves. That rung is the science, and it stays our job.

I am posting this as a draft, deliberately, because arguing about it in the open is rather the point of tenet 10. We are also writing it up — in the spirit of the FAIR paper — as a short article, anchored to those existing standards. If you think I have a tenet wrong, or that I have missed one, please tell me. That is how a manifesto earns the name.

---

The one-page manifesto (PDF): here.  ·  The short paper, in preparation: draft.  ·  The original blueprint: Rigon et al., HESS Opinions: Participatory Digital eARth Twin Hydrology systems (DARTHs), HESS 26, 4773–4800 (2022).