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