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, https://arxiv.org/abs/2607.17358v2, 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: https://arxiv.org/abs/2607.09416
• 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.

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.

The Statistical physics of unsaturated soil water: kinetic theory and non commutative pore water dynamics

I am giving this talk at the EGU General Assembly 2026 in Vienna last week, in the Hydrological Sciences division. The argument, in a single sentence: Richards' equation is not wrong, but it is the equilibrium limit of a deeper kinetic theory — in the same sense that the Navier–Stokes equations are the hydrodynamic limit of the Boltzmann equation for a gas. Mario Putti twenty years ago once asked me, "if not Richards, what else?"; this is my attempt at an answer that arrives after year dedicated to properly solve Richards equation, before with GEOtop and later with WHETGEO. 

The core object is a filling distribution g(r, x, t) : ℝ⁺ → [0, 1] that gives the volume fraction of pores of radius r that are water-filled at position x and time t. Theta is recovered as θ[g] = φ ∫₀^∞ g(r) f(r) dr. Hysteresis becomes the non-commutativity [W, D] ≠ 0 of the wetting and drying operators — geometry, not memory. Richards' equation is recovered as the small Damköhler limit Da → 0, with K(ψ) emerging as a derived transport coefficient built from the connectivity kernel C(r, r') rather than being postulated.

Materials

• Slides (PDF)  the deck I'll use in the presentation.
• Storyboard (DOCX)  the slide-by-slide reading guide, in five columns: spoken text, visual content, speaker notes, mounting comments. Useful if you want to present the same material yourself, or if you just want to follow along with what I actually said.
• Extended version of the slides — give me a few days — an annotated version with the full speaker text, more references, and the bits I had to cut for time.

Notebooks

These are the Jupyter notebooks I used to generate some of the figures in the slides, plus a few that produce supporting evidence in the supplementary material of the upcoming PRE papers. All run on top of OpenPNM 3.x and a small custom Y–L percolation code.

• Hysteresis_SWRC.ipynb — drainage and wetting branches in the (ψ, S_e) plane on a 3D pore network, with internal scanning curves. The figure on slide 9 of the talk comes from here. The notebook also documents an algorithmic artifact near the air-entry value (the missing air-trapping term during imbibition) — which is honest enough that I left it in.
• OpenPNM_Da_overshoot.ipynb — non-equilibrium overshoot in (θ, ⟨r⟩) and the universality crossover when the pore-size distribution becomes bimodal, governed by the Bhattacharyya overlap of the two modes.
• Percolation_K_threshold.ipynb — the percolation scaling K ∝ (θ − θ_c)^t with t ≈ 2, with finite-size scaling on three lattice sizes.
• subsection_pnm_mapping.tex — a short LaTeX subsection on how a two-tier pore-network maps onto the kinetic theory through a bimodal f(r) and a block-structured C(r, r'). Background reading for the OpenPNM notebooks.

Please find them zipped at this link.

Two upcoming papers

The full theoretical development is in two manuscripts, going to arXiv soon and submitted thereafter to Physical Review E --- give me a couple of weeks after EGU26:

1. The Statistical Physics of Unsaturated Soil Water: kinetic theory and non-commutative pore-water dynamics — the long paper. Builds the kinetic equation from the network thermodynamics, identifies the Onsager–Rayleigh gradient-flow structure, and proves that hysteresis is a geometric property of the configuration bundle (not a memory effect).
2. Richards' equation as a hydrodynamic limit: Chapman–Enskog derivation from the kinetic equation for unsaturated soil water — the short companion. Walks through the Chapman–Enskog expansion that recovers Richards' equation in the Da → 0 limit, with K(ψ) derived from the connectivity kernel.

Where this connects

The framework absorbs and extends a number of existing approaches that have been circling the same physics from different angles:

• Mixed-form Richards as the Da → 0 limit, with K(ψ) derived rather than postulated.
• Hassanizadeh–Gray as a thermodynamically consistent extension — pore-class-resolved here.
• Phase-field methods (Cahn–Hilliard) as gradient flow on a free energy — with explicit pore-network connectivity through C(r, r').
• Lucas–Washburn and its fractal variants as the single-capillary kinetic building block of C(r, r').
• Percolation-based hillslope frameworks with Damköhler and Péclet, where macropore activation is the Da > 1 transition.
• Compressible statistical soil mechanics (Einav–Liu 2023) — same occupancy dynamics governs the (ψ, σ') coupling.
• Freezing soil thermodynamics (Rempel et al. 2023, and our own work with Wani and D'Amato) — same kinetic framework with capillary pressure replaced by freezing-point depression.

This kinetic theory is not a parallel universe to Richards. It absorbs the existing physics, and it opens new measurements — directly observing g(r) is the obvious next experimental challenge

Theoretical and Numerical Tools for Studying the Critical Zone from Plots to Catchments

 What a valuable work is the thesis by Niccolò Tubini, here presented in its draft.  It covers works in hydrology of the critical zone, numerics, programming, software engineering, open science methods. Having a so wide horizon of interest it could  not be easy to grasp in all of these details, but it is well written and, we hope inspiring. As the Author, Ph.D. candidate says: "In the following we suggest that studying the CZ requires tools that are not yet readily available to researchers; then we propose one of our own. These tools should be flexible enough to allow the quick embedding of advancements in science"

Who wants to access the draft, can click on the figure of the Thesis first page below.

The Thesis included the work present in the two submitted paper by Niccolò,  on The Cryosphere and a second one presented in GMDD regarding WHETGEO-1D.  Whilst a thesis being considered kind of a definitive work, this one remains very much a work in progress with the extensions of the codes foreseen to arrive soon and whose informatics has already been implemented. All the tools developed during the Thesis are open source and freely available both as executable and source codes on Github.  Any comment or suggestion to the  Thesis as well as to the papers are welcomed. 

The Video of the defense is here.

WHETGEO 1D is out

WHETGEO-1D (Water HEat and Transport in GEOframe) is a physically based model simulating the water and energy budgets in a soil column. WHETGEO-1D is developed as an open-source code, adopting the Object-Oriented paradigm and a generic programming approach to improve its usability and expandability. WHETGEO-1D is fully integrated in the GEOframe/OMS3 system allowing the use of the many ancillary tools it provides. It comes on top of several years of work on engineering software, discussing and debating about Richards equation, taking care of getting appropriate integration methods in a travel that crossed hydrology, mathematics, numerics and software engineering.  Click on the figure below to access the paper.

The code is really solid and was throughly tested over the last three year by my students of the Hydrology class and in various applications. Its inputs and outputs can be analyzed by using Python and some standard Notebooks prepared to help the user to do it.  The paper has been submitted to the Geoscientific model development discussions (GMDD) and it is available for discussion to everyone and can be obtained by clicking on the Figure above. Complementary material is present in the GEOframe blog. The paper comes with its software, documentations and test and the GEOframe blog post explains where all of it is. You did a nice job Niccolò Tubini!

Here the paper on GMD.

Video Lectures on Hydrology

I am collecting here my video lectures on Hydrology (in my broken English). These are mostly part of the two GEOframe Winter School held in 2019 and 2020 and from a Summer School on Landslides made a few years ago. Since video lectures on these topics are uncommon, I think it is useful to index them.  I also invite anyone who has similar contributions to share them. I will be happy to add them to my list here.

Here they are below subdivided by arguments with their companion slides:

Catchments Delineation and Geomorphometry

• The basic theory (YouTube video 2019,YouTube2020)
• The derived quantities (YouTube video 2019,YouTube2020)
• Extracting the hillslope (YouTube Video 2019,YouTube2020)
• A little on some geomorphic laws (YouTube video 2019, YouTube2020)

Data Interpolation with Kriging 

• The Kriging's equations (YouTube2019 - Slides were a little modified for 2020, YouTube2020)
• Variography (YouTube video 2019, YouTubeVideo2020)
• The GEOframe Spatial interpolation Package (SIK) (YouTube2019, YouTube2020)
• Error estimation (leave one out) within SIK (YouTube video 2019, YouTube2020)

Richards equation

• Part I
• Part II
• Part III

Radiation for Hydrologists

• The Sun
• Stefan-Boltzmann law
• From Sun to Earth
• Copying with Earth Surface
• Shortwave absorptions by the atmosphere (YouTube video)
• Considering clouds (YouTube video)
• Shortwave-Terrain geometries (YouTube video)
• Longwave radiation (YouTube video)
• Table of Symbols

Evaporation and Transpiration

• Momentum and water vapor transport in atmosphere (YouTube2020)
• Evaporation as energy flow (YouTube2020)
• Derivation of the simplified energy budget (a là Penman-Monteith after Schymanski and Or) (YouTube2020)
• Simplified relatives of PMSO (YouTube2020)
• Leaf Temperature feedbacks
• Transpiration (YouTube2020)
• Further discussions on resistances and the physiology of plants in this video.
• Evaporation from soil (YouTube2020)
• From leaves to canopies (YouTube2019 Video, YouTube2020)

Hydrological Dynamical Systems  (a.k.a. lumped hydrological rainfall-runoff models) 

• Models in Hydrology
• The GEOframe type of models (YouTube2020)
• Representing the Hydrological Dynamical Systems (YouTube2020)
• The Embedded Reservoir Model (YouTube2020)

Other resources

• Other Videos, that I am providing for my Hydrological Modelling Class are here.
• Video collected by Kevin McGuire (GS) are here.

If you do not want to be just a tourist, you can go deeper and exercise with  Jupyter lab and GEOframe. For the latter, please see the material of the GWS2020. To anyone requesting, I can provide the original slides.

Advances in Richards 2D presentation at the Italian Hydrological Society meeting

The work of Niccolò Tubini is going who already developed a very solid Richards1D (with ponding),  code coupled with the energy budget is proceeding towards a 2D version on unstructured grids coupled, at present with a 1D de Saint-Venant equation. This is the summary of the work done so far given at the Italian Hydrological society meeting in Bologna.

Actually the de Saint-Venant coupling is not yet ready but it will be very soon. Stay tuned. Click on the figure for getting the presentation.

EGU Wien 2019: Two numerical models to solve Richards and energy equations

This contribution discusses two methods of integration for Richards equation and the heat equation (for T > 0 centigrades).  Results are encouraging and show that temperature could be important to get the right runoff production.

The model uses new numerics based on work by Vincenzo Casulli and Paola Zanolli, called nested Newton. The original poster is obtained by clicking on the figure.

Peter Germann's open book on preferential flow

Peter German (GS) first raised the issue of how preferential flow in soils due to preferential flow in a famous paper with Keith Beven (GS) in 1982, Macropores flow in soils. A benchmark paper from which uncountable studies followed. Now retired, he wrote a book entitled "Preferential flow: Stokes approach to infiltration and drainage".

Since 1982 there was a great debate on preferential flow modelling which is not yet concluded. Germann's book is evidently  a contribution from a protagonist of this story and he had the kindness to make it open. I am not sure about the content. Stokesian flows are those that happens when all acceleration is dissipated, so the subtitle is not very informative on the content. More detailed comments will arrive when I will have read it. For now you can get it by clicking on the Figure.

A little on soil physics

This shows the lectures I gave this week on soil and soil water to my class of hydrology.

Soils

Texture and structure of soils

Definitions

Darcy and Buckingham laws

Soil Water retention curves

Hydraulic conductivity

Hydraulic conductivity at saturation

Richards equation (first part)

Solving Richards equation

Richards equation 1d

Information about solving Richards equation with GEOframe tools

Towards the simultaneous account of Water and Energy budgets (a.k.a. GEOtop 3.0)

There are several ways to categorize this work. One way is to see it inserted in the GEOtop 3.0 project that aims to rebuild GEOtop physics and informatics. Another way is to see it inserted in the work to get better cryospheric parameterisations. A third way is to see it as a part of the process delineate in my professorship talk.  More pragmatically, it represents a new implementation of the coupled Richards equation and energy budget by using sound numerics and OMS3.
No easy to read. It requires application.

You can see the presentation at the OSF site, here.

A new topic for a Ph.D. in Hydrology at University of Trento. Modelling water flows under phase transitions

This study starts from a pore scale view of flow in soil and aggregate it at the representative elementary volume, (REV), scale according to statistical assumptions, to obtain new forms of the Richards equation. Flows are assumed to happen under normal and/or freezing conditions and under evapotranspiration demand. Transitions from unsaturated to saturated conditions will be properly accounted in all types of flow. The theoretical work at the basis of this proposal is contained in Dall’Amico et al. 2011 and Tubini, 2017. At the beginning the system will be modeled by coupling the water budget equation and the energy budget equation, neglecting vapor mass budget, as usually done. The candidate should take care of integrating the equations with appropriate and sound numerical methods that guarantee mass and energy conservation, following the footsteps of the work by Casulli and Zanolli (2010) and work for possible extensions.

There are various possible further development of this research. One is to couple the water and energy budget with surface waters simultaneously solved, another is to deal with water vapor explicitly. Others developments could come ongoing.

The informatics behind the code will follow (and, in case co-develops) the developments pursued by dott. Serafin, Ph.D. work inside the Object Modelling System, version 3 or subsequent (OMS3, David et al., 2013), that will take care implicitly of execution of parallel processes and will provide various services to computation (e.g. Serafin, 2016).

All the code developed will be done in Github (or similar platform), inside the GEOframe community and will be Open Source according to the GPL v3 license.
The candidate will take care of implementing, besides the code, the appropriate procedures for continuous integration of the evolving source code, and s/he will be also asked to maintain a regular rate of commits to the common open platform. Despite these conditions, and being free and open source, the code will be intellectual property by the coder. This will be guaranteed also by the components-based infrastructure offered by OMS3, which allows to better define the contributions of anyone. (See also: For incoming students, The tales of open source codes).

The implementation part will be followed, accompanied by testing activities, either for mathematical consistency, than for physical consistency with experiments and field measurements. These will be made especially by Dr. Stephan Gruber (GS) group at Carleton University, where the candidate will be asked to spend some periods od his/her doctorate. Participation to experimental activities will not be intended to be purely passive, the candidate will be asked to actively participate as much as feasible and reasonable to any part of the research.

The Ph.D. student is intended to produce, besides working and tested codes, also at least three papers in major journals (VQR Class A), of which, at least one as first Author.  Duration of the doctoral studies could be three or four years.

This project can enter either the curriculum C (Environmental Engineering) or the curriculum A (Modelling and Simulation) of our doctoral school.

For information please refers to riccardo.rigon <at> unitn.it

Essential References

Casulli, V., & Zanolli (2010). A nested newton-type algorithm for finite volume methods solving Richards' equation in mixed form. SIAM J. SCI. Comput., 32(4), 2225–2273.

M. Dall’Amico, S. Endrizzi, S. Gruber, and R. Rigon, An energy-conserving model of freezing variably-saturated soil, The Cryosphere, 5, 469-484, 2011, doi:10.5194/tc-5-469-2011.

David, O., Ascough, J. C., II, Lloyd, W., Green, T. R., Rojas, K. W., Leavesley, G. H., & Ahuja, L. R. (2012). A software engineering perspective on environmental modeling framework design: The Object Modeling System. Environmental Modelling and Software, 39, 1–13. http://doi.org/10.1016/j.envsoft.2012.03.006

Serafin, F., About graphs, DSL and replicable research, 2016, http://abouthydrology.blogspot.co.at/2016/11/about-graphs-dsl-and-replicable.html

Tubini, N. (2017, March 31). Theoretical Progress in freezing-thawing process studies. (R. Rigon, F. Serafin, & S. Gruber, Advisors.).

Reservoirology #4: the case of Richards 1d

This is the follow up of Reservoirology #3 post, where (there) I used  Petri nets to represent ordinary differential equations (ODEs). Here, instead, I try to use the same graphical formalism to represent Partial Differential Equations (PDEs), by extending a little the graphics.

Let's say that the number of specifications needed is larger here and, therefore, a lot of ancillary information has to be conveyed through Tables (and some interpretation). However, click on the figure above if you are curious to know more.

The Soil Water Retention Curves

When dealing with soils you are forced to implement mass conservation dependent on two variables, the dimensionless water content, usually named $\theta$ and suction, $\psi$, i.e. The energy contained in a volume of soil per unit mass. Therefore, to solve the budget, you need (at least) to get a new relationship which connects them.  This relation is called soil water retention curve. The plural in the title means that there are many. At least one for any soil type. 

In fact,  the  relationship, and precisely $\theta(\psi)$, is dependent on soil types and structure (and some other factor probably, like temperature, organic content etc). It is a statistical quantity, which averages the behavior of many pores, and an ensamble of water injecting/extracting possibilities.  

The figure below from Lu (GS) and Godt (GS) book (2013) is a clear visualisation of the problem.

The same Ning Lu, in a recent paper (2015) tried to disentangle the various forces acting on water when in pores, and obtained what is shown below.

As expected, the forces acting are not all of the same type, at varying suction values. At very high suction, adsorption forces act in which single water molecules adhere to soils. When more layers of water molecule add, water constitute  thermodynamics compounds, whose equilibrium is globally determined in between adhesion forces, bulk water weights, surface of water and air gas interactions, and which is usually known as capillarity.  

Laws governing capillarity are described by Young-Laplace and  Kelvin laws.  Some insight of the therodynamics of these phenomena (an excellent explanation, indeed) can be found in the first pages of Steudle (2001) review about plant-root suction. 

At this stage liquid water seems to, constitute a disconnected phase, while air gas is continuous inside the medium pores.  

Increasing the water content water becomes a continuos medium and usual hydrodynamics laws become valid.  A recent review of parameterizations of the soil water retention curves (not particularly deep or brilliant though) is given by Too et al. (2014) that cites other older reviews.

When pressure increase, however, we can have two effect which partially depends on how wetting happens. If wetting happens through some sort of flooding then air can stay trapped in pores and decrease the space available for water. The net effect is  associable to a decrease of porosity. However, when water fills all the space (i.e. the soil is saturated) the soil matrix cannot be considered anymore rigid. 

Assume it would be rigid. Then water content could not increase, any pressure applied to the saturated soil would transmit instantaneously through the water volume and water would be expelled where pressure is not applied or there is less pressure in a sort of piston flow.  

Instead, because the medium is not rigid, any pressure is transmitted with a certain speed, and pressure waves can be measured. This fact implies that after saturation, the system behaves as porosity increases and, at the same time pressure varies.  

From a practical point of view, soil water retention curves can be extended to positive pressure (negative suctions) adding a term which is well known in groundwater literature and is called specific specific storage. 

These qualitative descriptions do not end the complex phenomenology of water retention curves.  

As Nunzio Romano (GS) and coworkers noticed, and Kosugi (1994) before them, soil water retention curves shape depend directly on the pores' distribution. This, however, is not necessarily a unimodal distribution but can be multimodal because of soil structure and soil "disturbances" in form of macropores due to animal or roots decay. In this case soil water retention curves (their integral) can be more complex than expected, as shown in Figure below.

This opens to a series of generalisation, but it would be the topic of some other post (and actually was already the topic of several posts on soil freezing).

References

Kosugi, K. 1994. Three-parameter log-normal distribution model for soil water retention. Water Resour. Res. 30:891–901. 

Lu N, Godt JW. Hillslope Hydrology and Stability. Cambridge: Cambridge University Press; 2013. 

Lu, N. (2016). Generalized Soil Water Retention Equation for Adsorption and Capillarity. Journal of Geotechnical and Geoenvironmental Engineering, 142(10), 04016051–15. http://doi.org/10.1061/(ASCE)GT.1943-5606.0001524

Romano, N., Nasta, P., Severino, G., & Hopmans, J. W. (2011). Using Bimodal Lognormal Functions to Describe Soil Hydraulic Properties. Soil Science Society of America Journal, 75(2), 468. http://doi.org/10.2136/sssaj2010.0084

Steudle, E. (2001). The Cohesion-Tension Mechanism and the Acquisition of Water by Plant Roots. Annual Review of Plant Physiology-Plant Molecular Biology, 847–877.

Too, V. K., Omuto, C. T., Biamah, E. K., & Obiero, J. P. (2014). Review of Soil Water Retention Characteristic (SWRC) Models between Saturation and Oven Dryness. Open Journal of Modern Hydrology, 04(04), 173–182. http://doi.org/10.4236/ojmh.2014.44017

Integration of a Three-Dimensional Process-Based Hydrological Model into the Object Modeling System

This paper represents a first step of the unavoidable integration of GEOtop into OMS. I spent a lot of words in favour of this, and I cannot repeat it.  the paper abstract says:

The integration of a spatial process model into an environmental modeling framework can enhance the model’s capabilities. This paper describes a general methodology for integrating environmental models into the Object Modeling System (OMS) regardless of the model’s complexity, the programming language, and the operating system used. We present the integration of the GEOtop model into the OMS version 3.0 and illustrate its application in a small watershed. OMS is an environmental modeling framework that facilitates model development, calibration, evaluation, and maintenance. It provides innovative techniques in software design such as multithreading, implicit parallelism, calibration and sensitivity analysis algorithms, and cloud-services. GEOtop is a physically based, spatially distributed rainfall-runoff model that performs three-dimensional finite volume calculations of water and energy budgets. Executing GEOtop as an OMS model component allows it to: (1) interact directly with the open-source geographical information system (GIS) uDig-JGrass to access geo-processing, visualization, and other modeling components; and (2) use OMS components for automatic calibration, sensitivity analysis, or meteorological data interpolation. A case study of the model in a semi-arid agricultural catchment is presented for illustration and proof-of-concept. Simulated soil water content and soil temperature results are compared with measured data, and model performance is evaluated using goodness-of-fit indices. This study serves as a template for future integration of process models into OMS.

How is sweet to begin a year with a new publication. Have a nice reading (by clicking on the figure or here). 

A few topics for a challenging Master thesis

It is quite a long time that I am thinking to assign some master thesis around GEOtop 3.0.

A Master Thesis could be done even on the "simple" Richards equation. In this case the idea would be implementing the nested-Newton Casulli-Zanolli's method (on unstructured grids).

Time ago I assigned a little grant on this topic, but unfortunately, the work was not completed/ The material produced, in any case, is here. (I have also some other material, in FORTRAN, anyway).  In this case, the idea would be to use Java and develop further what already made by Francesco Serafin in his thesis. 

My outstanding colleague Michael Dumbser already promised to help me to complete the precedent work, and, at that point the main work would be to translate the procedural concepts into a object-oriented framework. In any case, who does it, would place one of the first stones of GEOtop 3.0.

Just thinking loudly, once started the Richards' work (which constitutes, however,  well defined and challenging enterprise) one could think how to implement coherently different flavours of the equation, for using bimodal or other water retention curves; for extending Richards analysis to integrate also the groundwater 3D equation, or studying the coupling with the evaporation/transpiration sink.  All alternatives that are interesting either from the numerical and the physical point of view.

The original problem I had in mind when I started this post, was the more ambitious one connected with the numerics and the physics of soil freezing (see Matteo Dall'Amico Ph.D. Thesis).  Our reference paper in the topic is the 2011 Dall'Amico et al. In the thesis and in the paper we wrote the equations in  3D but solved them in 1D with a not so particularly elegant method, which the nested-Newton algorithm could surpass by far. 

Working on the cryospheric side of Richards equation open a series of opportunity and especially the collaboration with  Stephan Gruber (with Carleton we have an exchange agreement, and the candidate could also stay for a few months there). 

Actually all these topics suggest that a very basic trial could be made to envision a scheme and an infrastructure that can accomodate all of these Richards variants by minimising code rewriting. But this would be probably a theme that could be completely developed in a Ph.D. In fact all these topics' task can clearly produce journal papers, if completed, and certainly open the road to some Ph.D. carrier.

Someone can think that everything is too much challenging, and actually, it is. However, all the topics are pretty mature in my mind, and the path to the solution is pretty well designed.