Six Notebooks to represent water budgets from compartment models (with EPNs inside)

Since the paper on GEOtop (Rigon et al., 2006), which claimed that a hydrological model should look at the water and energy budgets, rather than at single fluxes, in order to gain real knowledge of the system, I have been struggling to build models that respond to this effort. Estimating a discharge, an evaporation, a snow water equivalent one at a time is not the same thing as closing a budget: it is the budget — the simultaneous accounting of inputs, outputs and storage variations — that constrains the pieces to be mutually coherent, and it is from that coherence that understanding comes.

Detail of a Giacomo Balla artwork

The effort has grown in various directions since then. On the conceptual and strategic side, it evolved into the vision of Digital eARth Twin Hydrology systems, the DARTHs (Rigon et al., 2022), discussed in several posts on this blog (see the posts on DARTHs). On the software side, it required a complete rebuilding of the informatics of these systems, deployed in the GEOframe system: a component-based infrastructure built on the OMS3 framework (David et al., 2013; Formetta et al., 2014), whose NewAGE branch (Bancheri et al., 2020) is the one we routinely use for operational and research modelling (see the GEOframe posts). On the formal side, it produced the Extended Petri Nets (EPN) formalism (Bancheri, Serafin and Rigon, 2019), which gives a graphical and mathematical grammar for writing budget-based models: places for storages, transitions for fluxes, controllers for the variables that regulate rates without exchanging mass, splitters for the partition of fluxes.

In all of this effort, one persistent issue has been how to represent the water budget in a compact and coherent way (on the energy budget we are still working). Compact, because a budget involves many time series at once — inputs, outputs, several storages — and they must be readable on a single temporal abscissa. Coherent, because the representation should not merely juxtapose the series but make their mutual constraint visible: at every time step, inputs minus outputs must equal the sum of storage variations, and a figure of the budget should let the reader see whether this closure holds, column by column.

One step in this direction is the production of six Jupyter Notebooks, obtained with the help of Claude (Anthropic's AI assistant) and certainly improvable, which you find on OSF. As usual, this is a seed for improvements.

What the notebooks contain

The six parts share a single Python library (hydrobudget.py, accompanied by a test suite) — with a second module, epn_registry.py, for the model dictionary of Part 6 — and exchange data through feather files and JSON declarations, so that each notebook can be read — and modified — on its own.

Part 1 — Synthetic data. A minimal three-bucket model (snowpack, root zone, groundwater, with air temperature as forcing) generates the datasets used throughout: a one-year demonstration run, a four-year analysis run with an imposed dry year, and — this is a methodological point I care about — an independent thirty-year baseline from which all climatologies and drought thresholds are computed, so that the reference statistics are never contaminated by the period being analyzed. Everything is saved to feather, and every generated dataset is closure-checked at birth — including, now, a small hillslope-and-riparian system, which is the net whose drawing in Part 4 has a crossing: it could previously be drawn but never verified, for want of a table. This is shortcut since I was not having real data. 

Part 2 — Representing the budget. The central figure of the series: inputs drawn as bars hanging from the top of the panel (the hyetograph convention), outputs rising from the bottom, storage variations as signed stacked bars below, optionally the total storage in a third panel — all on one shared temporal abscissa, at any time step from hourly to yearly, with fluxes summed and states sampled at period end under aggregation. The budget check is built into the figure: a line marks in − out at every step, and each column of storage variations must reach it exactly. On the synthetic data the maximum column-wise closure error is of the order of 10⁻¹⁴ mm; on real data, the columns where the bars miss the line become a diagnostic of where measurements or model do not close. The bookkeeping is declared, not hard-coded — one learns, for instance, that snowfall is not a flux input to the soil but an input to the snowpack storage, and the figure enforces the distinction. What makes the residual informative, rather than merely small, is how it behaves when the control volume is declared wrongly: the notebook now tabulates four candidate volumes — the soil view as drawn, the same view with Δ(snowpack) added, the whole system, and the whole system with Δ(snowpack) removed. The two consistent pairings close to about 10⁻¹³ mm; both inconsistent ones fail by the same 116.32 mm, which is the largest monthly change in snowpack storage. Equal failure at the scale of exactly one term is the signature of a term placed on the wrong side of a boundary, not of an accumulating error.

Part 3 — Droughts. Once the budget is an object, droughts become excursions of that object: persistent negative anomalies of storages and accumulated fluxes with respect to their (baseline) climatology, identified by the classical threshold-level method of run theory (Yevjevich, 1967), with pooling, severities and intensities in explicit units. Because each drought type lives in a different variable of the budget, the notebook separates them — meteorological, snow, soil moisture, streamflow, groundwater — and displays their propagation chain on a timeline (in the spirit of Van Loon, 2015). Snow droughts are further classified as dry, warm, warm-and-dry or other, using precipitation ratios and temperature anomalies (Harpold et al., 2017): this is why the temperature is generated and stored in Part 1. The reference is now computed by day of year rather than by calendar month, and the notebook shows both: the step discontinuities of a monthly reference are visible in the anomaly it produces, and therefore in the spells the threshold method extracts from it — which is a way of seeing how much the choice of reference, and not only the choice of threshold, is doing. There is certainly more to work out, in this part.

Part 4 — Representing the EPNs. The same systems drawn as Extended Petri Nets, following the graphical conventions of Bancheri, Serafin and Rigon (2019): colored circles for places, squares for transitions inheriting the color of the place they exit, framed triangles for controllers with dashed information arcs, diamonds for splitters with their partition fractions, everything laid on an integer grid with mathematical symbols as labels. Two details I find pleasing: arc routing can circumnavigate nodes through waypoints, and when a crossing between arcs is genuinely unavoidable — the notebook distinguishes by computation the avoidable from the necessary ones — it is denoted by a small open circle, the old circuit-diagram convention for "crossing without exchange". The rules that place the nodes are worth stating, because they are computable and not aesthetic preferences: places and transitions stay on grid nodes, arrows are straight lines or compositions of straight lines except for the controllers, a multi-part arrow is permitted precisely in order to avoid a crossing, and among the admissible configurations the more compact is chosen. The justification for not trading one against the other is that a net which obeys the no-crossing rule in some layout can always be deformed onto the grid, so compactness and crossing-freedom are not in competition. Part 6 gives these rules a test they cannot have on a hand-made example: applied blind to 47 published nets, none of which declares a single coordinate, 34 come out with no crossing at all. The basic idea is that giving as input the EPN representation of the model, the program knows which data to expect and consequently which graphs to produce, see Part 5 below. There are various improvements that can be envisioned also here. 

Part 5 — The EPN as budget analysis. The point of the whole exercise: the EPN declaration is not only a drawing, it is the bookkeeping. From the same JSON that draws the net, the code derives the whole-system budget and one budget per place, verifies the closure node by node, and plots the two-panel budget figure of every place of the net. The same JSON then acts as a data contract for simulation outputs: we propose a format that GEOframe-NewAGE results should obey — one tidy table per hydrological response unit, one column per EPN transition (internal fluxes included, which is an argument for making them first-class outputs of any simulator), one column per place state — so that if simulator and analysis share the same EPN, every node must close to numerical precision, and any residual is either a format violation or a bug. Coherence between the GEOframe outputs and the EPN scheme of the simulator, in other words, becomes something one can check mechanically. One lesson from making the contract strict: the closure check used to skip the missing steps, so a table that was 95 % gaps reported a better residual than an honest one — precisely backwards. It now reports the sample size beside the error (on a deliberately broken demo table, n_steps = 4 of_possible = 364), so a good number obtained from four steps cannot pass for a good number, and a declared column that is present but mostly empty now fails the contract outright.

Part 6 — The EPN dictionary. If an EPN is a grammar, then it should be possible to write a vocabulary in it, and the vocabulary is worth more than any single sentence. Part 6 transcribes the models of the MARRMoT collection (Knoben et al., 2019) as EPNs — 47 nets over 46 model identifiers, since one model is drawn twice in the source with two different structures — and stores each as a JSON entry carrying the net, the governing equations as they are printed beside it, the constitutive closures, the parameters, and the page it was read from. The collection then becomes something one can query rather than leaf through: find the models with a snowpack, the ones with an IUH transition, the ones with at least four storages. Sizes run from one storage to nine, with a median of four; 43 of the nets have a junction, 33 a splitter, one an IUH.

The reason for storing the equations next to the drawing is that it makes the transcription falsifiable. The deck prints the ODE of every store beside its net, so dS/dt = Ps − Es − Perc is an independent statement of which fluxes touch that store and with which sign, and the code compares the two readings: 173 of the 187 places carry a printed equation, and of those 161 agree with the arcs drawn around them. Twelve places in ten models do not. I have not guessed at those: each is recorded in a table with the arc missing on one side or the other, because deciding which of the two readings is right means going back to the original paper, not to my transcription. Fourteen places print no equation at all in the deck, so they cannot be checked — a gap in the source, and it seemed better to count it than to hide it.

The second half of the notebook is about adding an entry, since a dictionary nobody can extend is a catalogue. A small two-bucket model is walked through the three gates an entry must pass — the schema, the equation cross-check, and a redraw — and each gate is broken first, so that the reader sees the error message rather than a description of it.

Two things I did not expect to learn from doing this. The first is that transcribing a drawing is the weak link in the whole exercise: the errors the cross-check caught in my own transcriptions were exactly the kind a careful reader makes — a flux attached to the wrong store, a splitter branch that should have been an evaporative loss, a melt term turned into a self-loop. The second is that the disagreements which survive are interesting in themselves. A published net and a published equation that do not say the same thing are a fact about the literature, not about my code.

A seed

The notebooks run top to bottom without errors, the figures are checked geometrically before they ship, and the library carries 112 tests — the dictionary's own invariants among them: that an entry's stored counts match its net, that its status flag is the computed verdict of the cross-check rather than an opinion, that the search index is recomputable from the entries it summarises — but none of this makes them finished. The synthetic model is deliberately minimal; the drought thresholds are heuristic; the energy budget is absent; the connection to real GEOframe runs is, for now, a contract waiting for its first signatures. As usual, this is a seed for improvements: take them, break them, and tell me where.

One practical warning for anyone who does. The .ipynb files are generated, by the scripts build_01.py … build_06.py: an edit made in Jupyter is discarded by the next build. This is deliberate — it is what makes the series reproducible, and the build is now byte-identical, so two builds of unchanged sources give identical files and a real change shows up as a real difference — but it means that the natural thing to do, which is to open a notebook and fix something in place, is the one thing that will quietly lose the fix. Edit the builder.

Files are here: https://osf.io/jt7z8/files/osfstorage

References

• Bancheri, M., Serafin, F., & Rigon, R. (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
• Bancheri, M., Rigon, R., & Manfreda, S. (2020). The GEOframe-NewAge modelling system applied in a data-scarce environment. Water, 12(1), 86. https://doi.org/10.3390/w12010086
• David, O., Ascough II, J. C., Lloyd, W., Green, T. R., Rojas, K. W., Leavesley, G. H., & Ahuja, L. R. (2013). A software engineering perspective on environmental modeling framework design: The Object Modeling System. Environmental Modelling & Software, 39, 201–213. https://doi.org/10.1016/j.envsoft.2012.03.006
• Formetta, G., Antonello, A., Franceschi, S., David, O., & Rigon, R. (2014). Hydrological modelling with components: A GIS-based open-source framework. Environmental Modelling & Software, 55, 190–200. https://doi.org/10.1016/j.envsoft.2014.01.019
• Harpold, A. A., Dettinger, M., & Rajagopal, S. (2017). Defining snow drought and why it matters. Eos, 98. https://doi.org/10.1029/2017EO068775
• Knoben, W. J. M., Freer, J. E., Fowler, K. J. A., Peel, M. C., & Woods, R. A. (2019). Modular Assessment of Rainfall–Runoff Models Toolbox (MARRMoT) v1.2: an open-source, extendable framework providing implementations of 46 conceptual hydrologic models as continuous state-space formulations. Geoscientific Model Development, 12(6), 2463–2480. https://doi.org/10.5194/gmd-12-2463-2019
• Rigon, R., Bertoldi, G., & Over, T. M. (2006). GEOtop: A distributed hydrological model with coupled water and energy budgets. Journal of Hydrometeorology, 7(3), 371–388. https://doi.org/10.1175/JHM497.1
• Rigon, R., Formetta, G., Bancheri, M., Tubini, N., D'Amato, C., David, O., & Massari, C. (2022). HESS Opinions: Participatory Digital eARth Twin Hydrology systems (DARTHs) for everyone. Hydrology and Earth System Sciences, 26, 4773–4800. https://doi.org/10.5194/hess-26-4773-2022
• Van Loon, A. F. (2015). Hydrological drought explained. WIREs Water, 2(4), 359–392. https://doi.org/10.1002/wat2.1085
• Yevjevich, V. (1967). An objective approach to definitions and investigations of continental hydrologic droughts. Hydrology Papers 23, Colorado State University.

The mathematician of slow manifolds, or: on work that waits for its readers (or paths for emergent properties)

A few weeks after we posted the second paper of our kinetic-theory series on arXiv — Richards' equation as a hydrodynamic limit,  I received a letter from A. J. (Tony) Roberts, of the University of Adelaide. He had read the preprint, recognized in it a structure he has been building for forty years, and wrote,  generously, precisely,  to offer "an alternative framework, one that provides complementary illumination." Attached was his recent paper in the Transactions of Mathematics and Its Applications (Roberts, 2025). Reading it, and then following the thread backwards through his earlier work, I had two reactions in quick succession. The first: this is exactly the rigorous scaffolding our derivation needed. The second, more uncomfortable: why had I, why had, as far as I can tell, essentially the whole hydrological and homogenization literature, never engaged with it?

This post is about both reactions. The sociology first, because it carries a lesson beyond this particular case; the substance after, because the substance is what hydrologists should actually take home.

How good work gets stranded

Roberts' program — using the modern theory of invariant manifolds to derive macroscale models from microscale dynamics, with proofs, at the finite scale separations of real physics — should have landed squarely in the homogenization mainstream: the community that computes effective properties of heterogeneous media, the RVE world that every pore-scale modeler implicitly inhabits. It did not, and the reasons are worth naming because none of them concerns the quality of the mathematics.

He arrived from the wrong direction: dynamical systems rather than the calculus of variations, publishing in journals (ANZIAM J., IMA J. Appl. Math., SIAM monographs) that the mechanics community does not routinely scan. His computer algebra runs on Reduce, a system few researchers under sixty have installed. And his 2025 paper confronts the mainstream head-on — by my count it contains thirty-one explicitly flagged points of contrast with standard homogenization practice, nearly all of them, as far as I can judge, technically warranted. But communities metabolize challenges more slowly than contributions, and an outsider's justified critique reads, sociologically, as an outsider's critique first and as justified much later.

There is a subtler reason too: his most distinctive results resist sloganization. "Macroscale models valid down to scale separations of two" contradicts folklore so entrenched — one or two orders of magnitude between micro and macro, says every textbook — that readers assume a special case. "An exact remainder term for the gradient expansion" sounds like bookkeeping until the day you need an error bar. And his "backwards theory" — your reduced model is exact, but for a system provably close to the one you specified (Hochs & Roberts, 2019) — is philosophically the correct validity statement and rhetorically a hard sell, because it sounds weaker than the false statement people prefer to make. The closest parallel I know is Gorban and Karlin's work on exact hydrodynamic manifolds for kinetic equations, which had the same semi-overlooked trajectory until the framing "Hilbert's sixth problem" finally gave it a banner. Roberts never found his banner. Perhaps hydrology, of all fields, can lend him one; we are, after all, professional users of the equation his theory certifies.

The substance, for hydrologists

When I described the framework to a colleague, the reaction was a version of a question I had asked myself: isn't this obvious? A linear operator has a null space; the null space becomes the model; the rest follows. It is worth saying carefully why the rest does not follow, because everything a practicing hydrologist would pay for lives precisely in the part that doesn't.

The null space tells you what the fast processes cannot erase — for soil water, exactly one thing, mass, hence the water content θ. That is a direction, not a model. A model requires that a curved surface exist in the space of all possible pore-filling configurations — one point per value of θ — onto which every soil state slides and along which it then travels. That this surface exists, attracts, and can be computed is a theorem with a hypothesis, and the hypothesis is not the null space: it is the spectral gap, the clean separation between the slowest internal redistribution rate and the rate of the forcing. When the gap holds, "local equilibrium" stops being an assumption and becomes a state the soil demonstrably reaches, at a computable rate — in principle a measurable spin-up time after every irrigation pulse. When the gap closes — at the percolation threshold, when the water phase fragments — the surface does not become inaccurate; it ceases to exist. Richards' equation fails there the way a rating curve fails when the river leaves its banks: the object being parametrized is gone.

Around this central fact, Roberts' theorems deliver things I have not seen stated anywhere in the hydrological literature. That the validity of a macroscale model is local and checkable: the gradient expansion carries an exact remainder (Bunder & Roberts, 2021), so the model is quantitatively fine in the drained profile and quantitatively suspect at the wetting front, with a number attached, instead of a global incantation about ε → 0. That the required scale separation is startlingly small — his worked examples hold down to about twice the microscale (Roberts, 2015) — which should give pause to every campaign that agonizes over REV support volumes. That a state variable of a reduced model need not correspond to any conservation law: the second variable of a dual-permeability model, seen clearly, is not the budget of a second continuum but a wetness contrast between pore populations, which does not balance but relaxes, like the overtone of a struck string dying away under the fundamental. And — the result that genuinely surprised me — that the admissible nonlinearity of a multi-domain model is capped by a ratio of two relaxation rates: if that ratio is modest, an elaborately nonlinear macropore–matrix exchange function is fitting structure the reduced description cannot resolve. One number, two eigenvalues, and a ceiling on how fancy your two-domain model is allowed to be.

There is even an answer to a question we never ask: when pore-network modelers impose periodic boundary conditions on a unit cell, who authorized them? Roberts' phase-shift construction — consider the ensemble of all shifted copies of the medium, in which periodicity becomes a theorem rather than an assumption — is the receipt. (Our kinetic theory, as it happens, never needed the trick: the pore-size axis is separate from space from the start, which is one of the small structural blessings of that formulation.)

What it did to our papers

The test of a framework is whether it changes what you write. Our series — the statistical physics of unsaturated soil water (the kinetic theory itself) and Richards' equation as its hydrodynamic limit (the Chapman–Enskog reduction) — derives the hydraulic conductivity as a Green–Kubo bracket over the relaxation spectrum of the pore network, and the dual-permeability models as a band projection of the same kinetic equation. Roberts' letter, and his theorems, sharpened both in ways that are now in revision. Where we wrote that multiple spectral gaps yield "several Richards equations," the correct statement — his correction, and he is right — is a nested family: each member of the hierarchy rests on a single gap, and reduces to the next by adiabatic elimination, with an exact bookkeeping identity (a sum rule) tying every level to the one conductivity K(θ). Where we invoked the limit Da → 0, the theorems permit the honest, stronger, finite-Da statement: existence of the slow manifold in a finite neighbourhood, attraction at a computable rate, error of the order of the residual. And his two-zone/two-mode equivalence (Roberts & Strunin, 2004) told us something we had not seen: our band description and the spectral description of dual permeability are the same object in two coordinate systems — Gerke–van Genuchten and the eigenmodes of the redistribution operator stop being rivals.

A closing thought on reading, and on how this post came to be

I will confess the obvious: metabolizing forty years of another person's mathematics is slow, and I have not done it alone. This post, and the revisions to our papers that preceded it, were worked out in sustained dialogue with Claude, Anthropic's AI assistant — not as an oracle, but as an interlocutor with whom I could transcribe Roberts' constructions into our own notation, step by step, asking at every turn "what, exactly, does this theorem consume, and what does it deliver?", testing objections, and letting the papers themselves remain the court of appeal. It is a different mode of study than the one I was trained in. It is dramatically faster at one specific thing: locating which of a framework's results are load-bearing for one's own problem, as opposed to true but idle. It does not replace reading the papers — nothing does, and the reading is slower and still under way — but it changed the order of operations: understand first, then read to verify and deepen, rather than read for months hoping understanding arrives.

And it is worth being precise about the causal chain, because none of its links was dispensable. Without Tony Roberts' email, none of this happens: I would not have found his work by searching, since — as the first half of this post argues — the literature's own structure had hidden it from where I was looking. Without the dialogue, his email would have produced a polite acknowledgment and a citation, not a restructured argument: the nested-family correction, the finite-Da certificate, the band–mode equivalence, the sum rule — each of these took days of back-and-forth to extract, verify, and fold into the manuscripts, work that by traditional means would have taken me months, if I had attempted it at all. And without the papers — his and ours — there would have been nothing to connect. A generous correspondent, a tireless interlocutor, and the primary literature: the triangle is the method, and I suspect it is quietly becoming the method of many of us. Better to say so openly than to let the acknowledgments pretend otherwise.

For those who want the patient version: start with the 2025 TMA paper for the panorama, the 2015 IMA paper for the spatial theory, and his SIAM book (Model Emergent Dynamics in Complex Systems, 2015). Hydrology runs, and has always run, on reduced models. It is a strange comfort to learn that there exists a body of theorems about when we are allowed to.

My thanks to Tony Roberts for writing, and for reading us first.

References

Roberts, A. J. (2025). Accurate families of multi-continuum micromorphic homogenisations in multi-D space-time via dynamical systems theory. Trans. Math. Appl. 9, tnaf001. doi:10.1093/imatrm/tnaf001
Roberts, A. J. (2015). Macroscale, slowly varying, models emerge from the microscale dynamics in long thin domains. IMA J. Appl. Math. 80, 1492–1518. doi:10.1093/imamat/hxv004
Roberts, A. J. (2015). Model Emergent Dynamics in Complex Systems. SIAM, Philadelphia. doi:10.1137/1.9781611973563
Bunder, J. E., Roberts, A. J. (2021). Nonlinear emergent macroscale PDEs, with error bound, for nonlinear microscale systems. SN Appl. Sci. 3, 1–28. doi:10.1007/s42452-021-04229-9
Hochs, P., Roberts, A. J. (2019). Normal forms and invariant manifolds for nonlinear, non-autonomous PDEs, viewed as ODEs in infinite dimensions. J. Differ. Equ. 267, 7263–7312. doi:10.1016/j.jde.2019.07.021
Roberts, A. J., Strunin, D. V. (2004). Two-zone model of shear dispersion in a channel using centre manifolds. Q. J. Mech. Appl. Math. 57, 363–378. doi:10.1093/qjmam/57.3.363
Gorban, A. N., Karlin, I. V. (2014). Hilbert's 6th problem: exact and approximate hydrodynamic manifolds for kinetic equations. Bull. Amer. Math. Soc. 51, 187–246. doi:10.1090/S0273-0979-2013-01439-3
Rigon, R. (2026). The statistical physics of unsaturated soil water. arXiv:2607.09416
Rigon, R. (2026). Richards' equation as a hydrodynamic limit: Chapman–Enskog reduction of the continuum kinetic equation for unsaturated soil water. arXiv:2607.17358

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.