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

Drying porous media

Evapotranspiration and soil water retention curves (SWRC) are intimately related. This was already apparent in the recent works of Dani Or and coworkers. Recently,  caused by studies to support the Master Thesis work by Niccolò Tubini on freezing soil, all these issues came back to my attention. As a result, I did a little of literary review (see below) which covered non traditional journals for a hydrologist, and which, however, is equally well illuminating of the aspects of the drying phenomenon.

Being as simple as possible, what I undeerstood in my readings is that the drying process (e.g. Metzger and Tsotsass, 2010) proceeds in a interplay between evaporation and capillary fluxes. Water evaporate from any meniscus, at different rates which are commanded locally by the Fick’s law of diffusion. The gradient moving water around is the difference between the saturation water vapor content (vapor tension) close to the meniscus and the vapor content already present in air gas. Since vapor tension in larger menisci (as described by Kelvin’s effect) is higher, larger pores evaporate faster than smaller pores. Because, water vapor above smaller menisci is less than over larger menisci, we must have an equilibrating vapor flow from larger pores to smaller pores locations (the same Fick’s law is acting). However, a second, and probably faster,  water transport is happening through the liquid phase. Internal pressures among pores of different size are imbalanced and a laminar fluid flows happens between them. This flow is controlled by the liquid viscosity. If viscosity is low, larger pores replenish smaller pores until they are emptied very fast. If viscosity is high smaller pores are also emptied.  At the extreme cases, low viscosity means (if large pores span the whole control volume) no mean gradient of water content, high  viscosity the formation of an evaporation front.
During process, besides temperature (which control both vapor tension, according to Clausius-Clapeyron relation and water viscosity - wi) what determines the evolution of the drying phenomenon are the topology and the geometry of the pores space, and the amount of liquid water itself.  The latter determines the continuity (or the discontinuity) of the water phase. When continuity of the water phase is broken (I imagine by the effect of heterogenous nucleation of air bubbles under high water tension values) smaller pores cannot be anymore supplied by water from larger pores, and macroscopic evaporation rate decays fast (I was tempted to say exponentially: but this has a quantitative interpretation that could not be true).

Overall, the outer water vapor content (which previously I called  "the vapor content already present in air gas") is determined by the turbulent transport of vapor. Therefore, the picture is pretty dynamic, heterogeneous, and quite in disequilibrium. The three equation of mass, energy and momentum conservation must be used to account properly for the phenomenon.

The poor hydrologist problem (php) is how to reduce this complexity at the macro (representative elementary volume - REV) scale. It is clear that playing around with all the above factors very different dynamics can be obtained which could appear weird or unlikely to happen. Pores in soil, at the end, do not have all the variety of structures and connection a mathematically oriented (or perverted) mind can image, and therefore, I hope, a macroscopic synthesis is possible.

Some would argue that this synthesis is already present and, being the actual theory that constitutes the core of the modern treatment of the Richards equation (starting from Mualem, 1976). Others remind the series of papers by Whitaker in late seventies. A must read, indeed.
However, I am not sure that this cannot be improved. The point were the latter theories are weaker is the description of the geometry and topology of porous media, which in both cases is simplified.

References

Defraeye, T. (2014). Advanced computational modelling for drying processes – A review. Applied Energy, 131, 323–344. http://doi.org/10.1016/j.apenergy.2014.06.027

Kemp, I. C. (2007). Drying software: past,present, and fuure. Drying Technology, 25, 1249–1263. http://doi.org/10.1080/07373930701438709

Metzger, T., & Tsotsas, E. (2010). Network models for capillary porous media: application to drying technology. Chemie Ingenieur Technik, 82(6), 869–879. http://doi.org/10.1002/cite.201000023

Metzger, T., Vu, T. H., Irawan, A., Surasami, V. K., & Tsotsas, E. (2008). Pore-Scale Modelling of Transport phenomena in drying, 187–204.

Mualem, Y. (1976). A new model fro predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research, 12(3), 513–522.

Neethirajan, S., Jayas, D. S., White, N. D. G., & Zhang, H. (2008). Investigation of 3D geometry of bulk wheat and pea pores using X-ray computed tomography images. Computers and Electronics in
Agriculture, 63(2), 104–111. http://doi.org/10.1016/j.compag.2008.01.019

Prat, M. (2011). Pore Network Models of Drying, Contact Angle, and Film Flows. Chem. Eng. Technol., 34(7), 1029–1038. http://doi.org/10.1002/ceat.201100056

Shokri, N., & Sahimi, M. (2012). Structure of drying fronts in three-dimensional porous media. Physical Review E, 85(6), 066312–1–8. http://doi.org/10.1103/PhysRevE.85.066312

Vorhauer, N., Metzger, T., & Tsotsas, E. (2010). Empirical Macroscopic Model for Drying of Porous Media Based on Pore Networks and Scaling Theory. Drying Technology, 28(8), 991–1000. http://doi.org/10.1080/07373937.2010.497088

Vu, T. H. (2006, June 1). Influence of pore size distribution on drying behaviour of porous media by a continuous model. Universität Magdeburg.

Whitaker, S. Simultaneous heat, mass and momentum transfer in porous media. Advances in Heat Transfer 1977, 13, 119–203.

  Yang, F., Griffa, M., Bonnin, A., Mokso, R., Di Bella, C., Munch, B., et al. (2016). Visualization of water drying in porous materials by X-ray phase contrast imaging. Journal of Microscopy, 261(1), 88–104. http://doi.org/10.1111/jmi.12319

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