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.)