Daily Forcings Datasets (Rainfall and Temperature) for hydrological Purposes over the Po river district

If you want to run a hydrological model in the Po basin, or anywhere in northern Italy, the first question is which meteorological forcing to use. A few years ago the choice was short. Today there are at least a dozen gridded products, regional and European, built with different networks, methods and purposes. This post is a review of what is available, written from the point of view of someone who needs daily forcing for water budgets with GEOframe or a similar model.

A disclosure: with Hossein Salehi and colleagues I have a dataset of my own in this list (Salehi et al., preprint, currently in opendiscussion at ESSD). I have tried to judge it with the same yardstick as the others.

What a hydrologist needs from a forcing dataset

Not every good climatological product is a good forcing. For water budgets, these are the requirements I use:

Requirement	Why it matters
Daily time step, with the day aligned to discharge records (00–24)	Snow, soil moisture and routing run daily; misaligned days shift peaks
Precipitation and temperature from the same network and method	Rain–snow partition, melt and ET depend on both; mixing products mixes biases
Tmin and Tmax, or radiation	Needed by Hargreaves or Priestley–Taylor ET and by melt and frost thresholds
Resolution close to the model units, about 1 km in Alpine sub-basins	Aggregating is exact; disaggregating needs assumptions
Coverage of the whole hydrological district, Swiss and French headwaters included	Budgets must close over the basin, not over administrative borders
Uncertainty for each cell and day	So that forcing error can be carried into budgets and calibration
Station data and metadata	So that users can check, re-interpolate or bias-correct
Documented behaviour with elevation and undercatch	High-elevation input errors dominate Alpine budgets
Open licence, persistent identifier, standard format	Reuse, citation, reproducibility (the FAIR principles)
Updates after 2020	Operational water-balance and drought monitoring

No product meets all of them.

The datasets

Dataset	Variables	Period, step	Grid	Method
ARCIS precipitation (Pavan et al., 2019)	P	1961–present, daily (9-to-9 local day)	~5 km	Modified Shepard with topographic distance, on 1,048 homogenized series
ARCIS temperature (Pavan et al., 2026)	Tmin, Tmax	1991–present, daily, monthly updates	~5 km	Piecewise lapse rate, urban and water fraction, in 28 sub-areas; Shepard on residuals
IH-GAR (Manara et al., 2026)	P	1951–2023, daily and monthly	~800 m	Anomaly method: elevation-based normals, monthly anomalies, daily fractions
APGD (Isotta et al., 2014)	P	1971–2008; extended version to 2019	5 km	Anomaly method with PRISM-type normals; Alpine part only
Crespi et al. (2018)	P normals	1961–1990 climatology	~800 m	Local weighted regression on elevation
Crespi et al. (2021)	P, mean T	1980–2018, daily	250 m	Anomaly method; Trentino–South Tyrol only
Salehi et al. (preprint)	P, mean T	1991–2020, daily	1 km	Kriging with daily variograms; daily lapse rate for temperature
EMO-1 (Salamon et al., preprint)	P, Tmin, Tmax, T, wind, radiation, vapour pressure	1990–2024, daily and 6-hourly	~1.5 km	Angular distance weighting
E-OBS	P, T, Tmin, Tmax and more	1950–present, daily	~11 km	Ensemble interpolation
CERRA-Land	P and land variables	1984–present	5.5 km	Regional reanalysis
EEAR-Clim (Bongiovanni et al., 2025)	Station P, T, Tmin, Tmax	to 2020, daily	stations	Quality-controlled, homogenized series, about 9,000 stations
BIGBANG (ISPRA)	Water-balance outputs	1951–2024, monthly	1 km	National monthly water-balance model

A few notes on each family.

ARCIS is the work of the regional meteorological services of north-central Italy. Its strengths are care for the data (homogeneity tests, synchronization of 9-to-9 readings) and the fact that it is updated operationally; its temperature companion (Pavan et al., 2026) handles Po Valley inversions with sub-regional, two-slope lapse rates. Its limits for Alpine hydrology are the 5 km grid and the absence of an elevation model for precipitation: Manara et al. (2026) find ARCIS lower than other products above about 750 m.

IH-GAR (Manara et al., 2026) uses the densest network (7,417 stations in its calibration domain, about 2,000 active each year) and the anomaly method, so its normals capture the increase of precipitation with elevation where stations are missing. Its daily values are built from monthly anomalies and interpolated daily fractions; the paper validates normals and monthly anomalies, but I could not find a validation of the daily fields themselves. Temperature has to come from another product.

APGD (Isotta et al., 2014) remains the Alpine reference for daily precipitation and shows the sharpest valley gradients. It covers only the Alpine part of the basin and is not updated.

EMO (Thiemig et al., 2022, and EMO-1) is the most complete in variables, open and updated. It is worth knowing how it is built in the Alps: EMO-5 imports about 5,000 APGD grid points as virtual stations, and its direct Italian station input is small. In the Alpine part of the Po basin it largely re-grids APGD.

E-OBS and CERRA-Land are useful for large-scale consistency, but too coarse, or too model-dependent, for Alpine sub-basins.

Our dataset was designed as hydrological forcing: precipitation and temperature from the same network (1,583 precipitation and 1,555 temperature stations, from the Po River District Authority and EEAR-Clim), on a 1 km grid, over the whole district plus a 20 km buffer. The review it received was useful and sometimes hard. The precipitation field is smoother than ARCIS, IH-GAR and APGD in some Alpine valleys, Valtellina above all. The lapse rate is one per day for the whole domain, so it misses local inversions. The record stops in 2020. What it adds is the daily variograms, the daily lapse rates and the per-cell uncertainty, which we will release with the fields.

How many stations, in the same area?

Comparing station totals is misleading, because each product covers a different domain. What matters is how many stations support the field inside the area you model. For the Po River District with a 20 km buffer (about 90,000 km²):

Dataset	Stations in the district	Basis
Salehi et al., precipitation	472 (1991–2000), 954 (2001–2010), 1,108 (2011–2020) active per day	Exact
ARCIS precipitation	roughly 500–580	Estimate, scaled by area from about 1,000 series per year
ARCIS temperature	roughly 570 (1991) to 930 (after 2013)	Estimate, scaled by area from their Fig. 1a
IH-GAR	roughly 1,100–1,400 active per year	Estimate, from their station density
EMO	few direct Italian stations; APGD points in the Alps	Thiemig et al. (2022)

The estimates assume uniform station density and should be replaced by real counts. They still show a point that is easy to miss: station support changes a lot over time. Our network more than doubles after 2000, and any network that grows this much makes trends from the gridded series suspect unless the series are homogenized and the network is controlled, as ARCIS and IH-GAR do.

Fitness for water-budget modelling

Dataset	P and T, same network	Tmin/Tmax	≤ 1.5 km	Whole district	Per-cell uncertainty	P with elevation	Homogenized	After 2020
ARCIS P + T	same agencies, different methods	yes	no	Italian part	no	no	yes	yes
IH-GAR	P only	no	yes	yes	no	normals	yes	to 2023
APGD	P only	no	no	Alps only	ensemble version on request	normals	no	to 2019
EMO-1	yes	yes	yes	yes	yes (EMO-5)	limited	no	yes
E-OBS	yes	yes	no	yes	ensemble spread	no	no	yes
Salehi et al.	yes	mean T only	yes	yes	yes, kriging variance and cross-validation	tested and rejected	no	no

My reading, by use:

• Long-term climate and trends: ARCIS and IH-GAR, which are homogenized and designed for it.
• Precipitation totals in Alpine basins: IH-GAR, or APGD up to 2019; check ARCIS above 750 m.
• Operational monitoring: ARCIS precipitation and temperature, updated monthly.
• Daily forcing for a distributed model of a Po sub-basin: EMO-1 if you need many variables and recent years; our dataset if you want precipitation and temperature from one regional network with their uncertainty, with the caveats above; and in any case a check of the totals against IH-GAR at altitude.

Data availability and FAIR principles

Dataset	Identifier and access	Licence	Station data
ARCIS	arcis.it, no DOI	check the site	no
IH-GAR	UNIMI Dataverse, doi:10.13130/RD_UNIMI/CSYYAP	check the repository	no
APGD	MeteoSwiss, doi:10.18751/Climate/Griddata/APGD/1.0, registration	non-commercial use	no
Crespi et al. (2018)	ISAC-CNR	CC BY-NC-ND 4.0	no
Crespi et al. (2021)	PANGAEA, doi:10.1594/PANGAEA.924502	CC BY 4.0	partly
EMO-1	JRC Data Catalogue	CC BY 4.0	no
E-OBS, CERRA-Land	Copernicus Climate Data Store	open	partly (E-OBS)
EEAR-Clim	Zenodo, doi:10.5281/zenodo.10951609	CC BY 4.0	yes
Salehi et al.	Zenodo, doi:10.5281/zenodo.19207256	CC BY 4.0	no (sources: Po District Authority, EEAR-Clim)

The weakest point is common to all regional grids: none releases the input station series, because they belong to the regional services. EEAR-Clim is the exception and a real step forward. A second gap is code: among these products, only EMO and ours (GEOframe Krigings, GPL-3) point to open interpolation code.

What is still missing

• A daily product that combines the elevation-aware normals of the anomaly method with an explicit uncertainty for every day and cell.
• Validation of daily fields at the daily scale, including wet-day frequency and extremes, for products built from monthly components.
• Independent checks of precipitation at altitude, through water budgets closed on gauged basins with storage and evapotranspiration estimated independently. Snow undercatch cannot be seen by comparing gauges with gauges.
• Open station data. Without it, every dataset is a black box that can only be compared with other black boxes.

References

• Bongiovanni, G., et al. (2025). EEAR-Clim: a high-density observational dataset of daily precipitation and air temperature for the extended European Alpine region. Earth Syst. Sci. Data, 17, 1367–1391. doi:10.5194/essd-17-1367-2025
• Crespi, A., Brunetti, M., Lentini, G., Maugeri, M. (2018). 1961–1990 high-resolution monthly precipitation climatologies for Italy. Int. J. Climatol., 38, 878–895. doi:10.1002/joc.5217
• Crespi, A., Matiu, M., Bertoldi, G., Petitta, M., Zebisch, M. (2021). A high-resolution gridded dataset of daily temperature and precipitation records (1980–2018) for Trentino-South Tyrol. Earth Syst. Sci. Data, 13, 2801–2818. doi:10.5194/essd-13-2801-2021
• Isotta, F. A., et al. (2014). The climate of daily precipitation in the Alps. Int. J. Climatol., 34, 1657–1675. doi:10.1002/joc.3794
• Manara, V., et al. (2026). A new daily high-resolution gridded precipitation dataset for the Italian Greater Alpine Region (1951–2023). J. Hydrol.: Reg. Stud., 67, 103776. doi:10.1016/j.ejrh.2026.103776
• Pavan, V., et al. (2019). High resolution climate precipitation analysis for north-central Italy, 1961–2015. Clim. Dyn., 52, 3435–3453. doi:10.1007/s00382-018-4337-6
• Pavan, V., et al. (2026). A new operational dataset of gridded minimum and maximum temperature over north-central Italy 1991 to present. Climate Services, 43, 100689. doi:10.1016/j.cliser.2026.100689
• Salamon, P., et al. EMO-1: an improved version of the high-resolution multi-variable gridded meteorological dataset for Europe. Earth Syst. Sci. Data Discuss., essd-2025-723
• Salehi, H., et al. (2026). A 30-year 1-km daily precipitation and air temperature dataset for the Po River District (Italy). Earth Syst. Sci. Data Discuss. doi:10.5194/essd-2026-467
• Thiemig, V., et al. (2022). EMO-5: a high-resolution multi-variable gridded meteorological dataset for Europe. Earth Syst. Sci. Data, 14, 3249–3272. doi:10.5194/essd-14-3249-2022
• Copernicus Climate Data Store: E-OBS and CERRA-Land
• ISPRA, BIGBANG national water balance: workshop presentation, March 2026

Introduction to the Random Forest algorithm (for Hydrologists)

Random forests are among the most used machine-learning tools in hydrology, and among the least understood by the people who use them. This post is a guided reading of how they work, written for hydrologists who want the logic before the library call.

I follow one source closely: Chapter 8, Tree-Based Methods, of An Introduction to Statistical Learning (ISL) by Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani, together with the Stanford video lectures in which Hastie and Tibshirani teach that chapter. The book exists in an R edition (2nd ed., 2021) and a Python edition (ISLP, 2023, with Jonathan Taylor); both are free as PDFs from statlearning.com. The chapter numbering and notation below are theirs. Where I depart from them, it is to add a hydrological reading, and I say so. The relevant lectures are embedded section by section; the whole Chapter 8 series is in this YouTube playlist.

The chapter builds in four steps, and so does this post: a single tree, then many trees averaged (bagging), then many decorrelated trees averaged (random forests), then trees grown sequentially (boosting). Each step fixes a weakness of the previous one.

1. A single regression tree (ISL §8.1)

A regression tree cuts the predictor space into \(J\) non-overlapping boxes \(R_1, \dots, R_J\) and predicts, in each box, the mean of the training responses that fall in it. Think of predicting mean annual runoff from basin area, mean elevation, aridity index and forest cover: the tree asks is aridity < 0.8?, then is elevation > 1500 m?, and so on, until each leaf holds a group of similar basins.

The boxes are chosen to minimise the residual sum of squares:

\[ \mathrm{RSS} = \sum_{j=1}^{J} \sum_{i \in R_j} \left( y_i - \hat{y}_{R_j} \right)^2 \]

Finding the best partition is computationally infeasible, so the tree is grown by recursive binary splitting: at each step, take the single predictor \(X_j\) and cutpoint \(s\) that most reduce RSS, split, and repeat inside each half. It is greedy and top-down; it never looks back to see whether an earlier split was a poor choice.

Hastie and Tibshirani, Decision Trees (Stanford Statistical Learning, Ch. 8, 14:37).

A tree grown until each leaf holds a handful of points fits the training data almost perfectly and predicts new data badly. ISL's remedy is cost-complexity pruning: grow a large tree \(T_0\), then for each value of a penalty \(\alpha\) find the subtree \(T\) that minimises

\[ \sum_{m=1}^{|T|} \sum_{i:\, x_i \in R_m} \left( y_i - \hat{y}_{R_m} \right)^2 + \alpha |T| \]

where \(|T|\) is the number of leaves. \(\alpha\) plays the same role as \(\lambda\) in the lasso, and it is chosen by cross-validation.

Pruning a Decision Tree (Ch. 8, 11:45).

For classification the same machinery applies, with the Gini index or cross-entropy replacing RSS; the lecture Classification Trees and Comparison with Linear Models covers it.

Hastie and Tibshirani are frank about the balance sheet. Trees are easy to explain, mirror human decision-making, handle qualitative predictors without dummy variables, and can be drawn. But a single tree is usually less accurate than other methods, and it is non-robust: change the data a little and the tree can change a lot. In statistical language, trees have high variance. Everything that follows in the chapter is a way of killing that variance.

2. Bagging (ISL §8.2.1)

The idea rests on an elementary fact. If \(Z_1, \dots, Z_n\) are independent, each with variance \(\sigma^2\), their mean has variance \(\sigma^2/n\). Averaging reduces variance. If we had \(B\) independent training sets, we could grow \(B\) trees and average their predictions.

We do not have \(B\) training sets, so we fake them with the bootstrap (ISL Chapter 5): draw \(B\) samples of size \(n\), with replacement, from the one training set we have; grow a deep, unpruned tree on each; average:

\[ \hat{f}_{\mathrm{bag}}(x) = \frac{1}{B} \sum_{b=1}^{B} \hat{f}^{*b}(x) \]

For classification, average becomes majority vote. Each tree is deliberately overfitted (low bias, high variance); the averaging removes the variance. \(B\) is not a tuning parameter in the usual sense: a large \(B\) does not overfit, it simply costs time. You stop when the error has stabilised, often a few hundred trees.

Bagging also gives an error estimate for free. Each bootstrap sample leaves out, on average, about one third of the observations (the fraction tends to \(1/e \approx 0.368\)). For each observation, predict it using only the roughly \(B/3\) trees that did not see it. The resulting out-of-bag (OOB) error is a valid test-error estimate, and for large \(B\) it is essentially equivalent to leave-one-out cross-validation, as the authors note in the lecture below.

A caution the book does not need to make, but hydrologists do: the bootstrap assumes the observations are exchangeable. Daily streamflow is not. If neighbouring days land in both the bootstrap sample and the OOB set, the OOB error is optimistic. With time series, a block or year-based hold-out is the honest test.

3. Random forests (ISL §8.2.2)

Random forests improve bagging with one small change: at each split, the tree may only choose among a random subset of \(m\) of the \(p\) predictors, drawn afresh at every split. Typical choices are \(m \approx \sqrt{p}\) for classification and \(m \approx p/3\) for regression; \(m = p\) gives back bagging.

Bootstrap Aggregation (Bagging) and Random Forests (Ch. 8, 13:45): the core lecture for this post.

Why would forbidding a tree from seeing the best predictor help? Because bagged trees are correlated. Suppose one predictor is very strong (in a runoff problem, precipitation almost always is). Every bagged tree will then split on it first, and the trees will look alike. Averaging similar things does not reduce variance much. The arithmetic, given in The Elements of Statistical Learning (§15.2) and implicit in ISL, is the variance of an average of \(B\) identically distributed variables with pairwise correlation \(\rho\):

\[ \mathrm{Var}\!\left( \frac{1}{B}\sum_{b=1}^{B} \hat{f}^{*b} \right) = \rho\,\sigma^2 + \frac{1-\rho}{B}\,\sigma^2 \]

As \(B\) grows the second term vanishes, but the first does not. Bagging attacks the second term; random forests attack \(\rho\). By hiding the strong predictor from most splits, the other predictors get a chance, the trees diversify, and the average becomes more stable. In the lecture Hastie illustrates this on a high-dimensional gene-expression example, where \(m = \sqrt{p}\) clearly beats bagging; the gain is largest when many predictors are correlated.

The practical consequences:

• Like bagging, a random forest does not overfit by adding trees; \(B\) only needs to be large enough.
• \(m\) is the one parameter worth tuning, using the OOB error.
• Trees are grown deep and unpruned. Minimum leaf size (typically 5 for regression) is a secondary control.
• The price is interpretability: a forest of 500 trees cannot be drawn.

4. Variable importance, and how to misread it

To recover some interpretability, ISL records, for each predictor, the total decrease in RSS (or Gini index) produced by splits on it, averaged over all \(B\) trees. Ranked and plotted, this is the familiar variable-importance chart that ends half the hydrological ML papers I review. The lecture in §5 below discusses it alongside boosting.

Three warnings, the first from ISL's spirit and the others from the later literature:

1. Importance is predictive, not causal. A high score says the forest found the variable useful for splitting this dataset. It does not say the variable drives the process.
2. Correlated predictors share or steal credit. Elevation, temperature and snow fraction carry overlapping information; the forest spreads importance among them in ways that depend on \(m\) and on chance.
3. Impurity-based importance is biased toward continuous predictors and categorical ones with many levels (Strobl et al., 2007). Permutation importance (Breiman, 2001), which measures how much the OOB error grows when one predictor's values are shuffled, is the safer default.

If the question is physical, the variable-importance plot is the beginning of an argument, not its conclusion.

5. The other branch: boosting and BART (ISL §8.2.3–8.2.4)

Chapter 8 does not stop at forests, and a reader should know why. Bagging and random forests grow trees independently and average them. Boosting grows them sequentially: each small tree is fitted to the residuals of the current model, and only a shrunken fraction of it is added:

\[ \hat{f}(x) \leftarrow \hat{f}(x) + \lambda\, \hat{f}^{b}(x), \qquad r_i \leftarrow r_i - \lambda\, \hat{f}^{b}(x_i) \]

The model learns slowly. It has three tuning parameters: the number of trees \(B\) (which, unlike in a forest, can overfit), the shrinkage \(\lambda\) (typically 0.01 or 0.001), and the tree depth \(d\) (often 1 or 2; with \(d = 1\), stumps, the model is additive). Boosting usually edges out random forests in accuracy but needs more care. XGBoost and LightGBM are its modern industrial descendants.

Boosting and Variable Importance (Ch. 8, 12:03).

The 2nd edition of the book adds Bayesian Additive Regression Trees (BART), which mixes both ideas: many trees, each repeatedly perturbed to fit the partial residual left by the others, with the output being a posterior distribution, hence uncertainty bands, at little tuning cost. BART is not in the 2014 video series.

6. A hydrologist's reading

What follows is mine, not the book's. Random forests are excellent at some hydrological jobs and structurally unable to do others, and the reasons follow directly from the construction above.

Where they shine. Regionalisation (predicting signatures or parameters in ungauged basins from catchment attributes); digital soil mapping; land-cover and snow-cover classification; gap-filling of station records; and post-processing, that is, learning the residuals of a process-based model such as one built in GEOframe. In all these cases the predictors are tabular, the relations are non-linear with interactions, and the forest needs almost no tuning.

Where they fail by design.

• No extrapolation. A forest predicts averages of training responses, so its output is bounded by what it has seen. A flood larger than any in the record, or a climate warmer than the calibration period, is outside its reach. This is the most important single fact about trees for climate-impact work.
• No conservation laws. Nothing ties predicted runoff, evapotranspiration and storage change to a closed water budget. Each output is fitted separately.
• No memory. A forest maps inputs to output; dynamics must be put in by hand through lagged predictors, antecedent indices and the like.
• Dependent data. As noted for OOB error, spatial and temporal autocorrelation make random cross-validation optimistic. Use spatially blocked or leave-basins-out schemes.

Uncertainty. A plain forest gives a point estimate. Quantile regression forests (Meinshausen, 2006) keep the full distribution of responses in each leaf and return predictive quantiles at no extra cost, which is often what a hydrological application actually needs.

My own conclusion is that forests are best used with physics, not instead of it: as a fast, honest benchmark that a process model must beat, and as a tool for learning what the process model gets wrong.

7. Try it: a synthetic Budyko test

The example below generates 500 synthetic basins whose runoff follows Fu's form of the Budyko curve, plus a small forest effect and noise. Elevation is a decoy with no effect. It fits a random forest with \(m = p/3\), reads the OOB score, computes permutation importance, and then asks the forest about a basin wetter than any it was trained on. ISLP's own lab (§8.3) does the same steps on standard datasets, and the R lab is shown in the video below.

import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.inspection import permutation_importance

rng = np.random.default_rng(42)
n = 500
P   = rng.uniform(400, 2000, n)    # annual precipitation (mm)
PET = rng.uniform(400, 1200, n)    # potential ET (mm)
F   = rng.uniform(0, 1, n)         # forest fraction (-)
Z   = rng.uniform(200, 3000, n)    # mean elevation (m): a decoy
w = 2.6                            # Fu parameter
ET = P * (1 + PET/P - (1 + (PET/P)**w)**(1/w))
Q  = P - ET - 50*F + rng.normal(0, 20, n)
X  = np.column_stack([P, PET, F, Z])
names = ["P", "PET", "forest", "elevation"]

rf = RandomForestRegressor(n_estimators=500, max_features=1/3,
                           min_samples_leaf=5, oob_score=True, random_state=0)
rf.fit(X, Q)
print(f"OOB R^2: {rf.oob_score_:.3f}")

imp = permutation_importance(rf, X, Q, n_repeats=10, random_state=0)
for name, m in sorted(zip(names, imp.importances_mean), key=lambda t: -t[1]):
    print(f"{name:>10s}: {m:.3f}")

# Extrapolation: a basin wetter than any in the training set
x_new = np.array([[3000, 800, 0.5, 1000]])
et_new = 3000 * (1 + 800/3000 - (1 + (800/3000)**w)**(1/w))
print(f"True Q: {3000 - et_new - 25:.0f} mm, RF Q: {rf.predict(x_new)[0]:.0f} mm")

The output (scikit-learn, random_state fixed):

OOB R^2: 0.896
         P: 1.219
       PET: 0.185
    forest: 0.020
 elevation: 0.020
True Q: 2212 mm, RF Q: 918 mm   (max Q in training: 1559 mm)

Three lessons in a few lines. The OOB \(R^2\) of 0.90 is a decent fit obtained with no tuning. The importance ranking recovers precipitation and PET, but the real forest effect is indistinguishable from the decoy: small physical effects drown in noise. And for \(P = 3000\) mm the forest predicts 918 mm against a true 2212 mm, below even the largest runoff in its training set. That is §6's extrapolation warning, in numbers.

The same experiment in Java, with Smile

For those of us whose models live on the JVM, as GEOframe does, Smile (Statistical Machine Intelligence and Learning Engine, by Haifeng Li) offers a random forest that can sit inside the same code base. The defaults match the book: \(m = p/3\) for regression, minimum leaf size 5, OOB metrics computed during training. Add the dependency (Smile 4 and later need Java 21 or newer):

<dependency>
  <groupId>com.github.haifengl</groupId>
  <artifactId>smile-core</artifactId>
  <version>6.3.0</version>
</dependency>

import java.util.Random;
import smile.data.DataFrame;
import smile.data.formula.Formula;
import smile.regression.RandomForest;

public class BudykoForest {

    // Fu's form of the Budyko curve: actual ET from P and PET
    static double fuET(double P, double PET, double w) {
        return P * (1 + PET / P - Math.pow(1 + Math.pow(PET / P, w), 1 / w));
    }

    public static void main(String[] args) {
        int n = 500;
        double w = 2.6;
        Random rng = new Random(42);

        double[][] data = new double[n][];
        double qMax = Double.NEGATIVE_INFINITY;
        for (int i = 0; i < n; i++) {
            double P   = 400 + 1600 * rng.nextDouble();   // annual precipitation (mm)
            double PET = 400 +  800 * rng.nextDouble();   // potential ET (mm)
            double F   = rng.nextDouble();                // forest fraction (-)
            double Z   = 200 + 2800 * rng.nextDouble();   // mean elevation (m): a decoy
            double Q   = P - fuET(P, PET, w) - 50 * F + 20 * rng.nextGaussian();
            data[i] = new double[] {P, PET, F, Z, Q};
            qMax = Math.max(qMax, Q);
        }
        String[] names = {"P", "PET", "forest", "elevation", "Q"};
        DataFrame df = DataFrame.of(data, names);

        // 500 trees; mtry defaults to p/3, nodeSize to 5
        RandomForest rf = RandomForest.fit(Formula.lhs("Q"), df,
                                           new RandomForest.Options(500));

        System.out.printf("OOB R^2: %.3f%n", rf.metrics().r2());

        // Impurity-based importance: total RSS decrease, summed over trees
        double[] imp = rf.importance();
        for (int j = 0; j < imp.length; j++) {
            System.out.printf("%10s: %.3e%n", names[j], imp[j]);
        }

        // Extrapolation: a basin wetter than any in the training set
        DataFrame wet = DataFrame.of(new double[][] {{3000, 800, 0.5, 1000, 0}}, names);
        double qTrue = 3000 - fuET(3000, 800, w) - 25;
        System.out.printf("True Q: %.0f mm, RF Q: %.0f mm (max Q in training: %.0f mm)%n",
                          qTrue, rf.predict(wet.get(0)), qMax);
    }
}

Java's random generator is not NumPy's, so the numbers will differ slightly from the Python run, but the story does not: a high OOB \(R^2\), precipitation dominating, and a wet-basin prediction capped near the training maximum. One difference is instructive. Smile's importance() is the impurity-based measure of §4, not permutation importance, so expect the decoy elevation to receive a visibly non-zero score. That is warning 3 of §4 at work.

Lab: Random Forests and Boosting (Ch. 8, 15:35), in R. A companion Lab: Decision Trees precedes it.

References and lectures

The basis of this post

• James, G., Witten, D., Hastie, T., Tibshirani, R. (2021). An Introduction to Statistical Learning, with Applications in R, 2nd ed. Springer. Chapter 8, Tree-Based Methods. Free PDF at statlearning.com.
• James, G., Witten, D., Hastie, T., Tibshirani, R., Taylor, J. (2023). An Introduction to Statistical Learning, with Applications in Python. Springer. Chapter 8 and its lab. Same site.
• Hastie, T., Tibshirani, R. Statistical Learning, Stanford online course, Chapter 8 lectures: Decision Trees, Pruning a Decision Tree, Classification Trees and Comparison with Linear Models, Bagging and Random Forests, Boosting and Variable Importance, Lab: Decision Trees, Lab: Random Forests and Boosting. Full Chapter 8 playlist; an index of all chapters. The current edition of the course runs on Stanford Online and edX.

For going further

• Hastie, T., Tibshirani, R., Friedman, J. (2009). The Elements of Statistical Learning, 2nd ed. Springer. Chapter 15, Random Forests.
• Breiman, L. (2001). Random forests. Machine Learning, 45, 5–32.
• Strobl, C., Boulesteix, A.-L., Zeileis, A., Hothorn, T. (2007). Bias in random forest variable importance measures. BMC Bioinformatics, 8, 25.
• Meinshausen, N. (2006). Quantile regression forests. Journal of Machine Learning Research, 7, 983–999.

Software

• Li, H. Smile: Statistical Machine Intelligence and Learning Engine. haifengl.github.io; artifact smile-core on Maven Central.
• Pedregosa, F. et al. (2011). Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12, 2825–2830.

Interpolating for hydrological purposes

Every hydrological model starts from a map of precipitation and temperature that nobody measured. We measured points, at rain gauges and thermometers, and we filled the space between them with an interpolation. How well we did that filling sets the quality of everything downstream: snow, evapotranspiration, discharge, the water budget itself.

This post collects what I have learned, from the literature and from building a 1 km daily dataset for the Po River District (Salehi et al., preprint), about doing this filling well. It is written as a set of practical guidelines, with the papers behind each one.

1. A fine grid is not fine information

The starting point is still Daly (2006): do not equate resolution with realism. A 1 km grid is a rendering choice. The information it contains is set by the stations and by the spatial structure of the field, not by the grid spacing.

Daly also tells us where the difficulty lies. Terrain and coastal effects matter least at scales above 100 km and most below 10 km. With a typical station spacing of 100 km they cannot be resolved, “except in densely populated regions of developed countries”. Northern Italy is one of those regions, so the question is not whether to interpolate below 10 km, but how to say where the result is informative and where it is not.

Two things are often read into Daly’s paper that it does not say:

• it does not reject interpolation below 10 km; its conclusions depend on the setting (flat land, mountains without rain shadows, complex terrain) and on station density;
• it does not reject cross-validation; it warns against comparing cross-validation errors computed on different networks, which is a different and correct point.

2. Know your network before choosing a method

Most interpolation problems are network problems. Before any method, describe the network:

• Stations active each day, not in total. A network of 1,500 stations with 50% missing data can mean 500 stations active in the 1990s and 1,100 after 2000. In the Po dataset, the active precipitation network grew from 472 stations (1991–2000) to 1,108 (2011–2020).
• Distances, not densities. For each cell and day, compute the distance to the nearest station and to the k-th nearest, where k is the number of stations your interpolation uses. These maps show where the field is supported and where it is extrapolated. In the Po dataset, the median distance to the seventh-nearest station fell from about 28 km to 15 km over the three decades.
• Hypsometry. Compare the elevation distribution of stations with that of the DEM, as Manara et al. (2026) and Pavan et al. (2026) do. High elevations are almost always under-sampled.
• Time of observation. Many manual gauges in Italy recorded 9-to-9 local time; ARCIS corrects for it (Pavan et al., 2019). Mixing 9-to-9 and 00–24 days shifts daily values and misaligns them with discharge.
• Record length and homogeneity. Short series are useful for spatial detail. Long homogenized series are needed for trends. You rarely get both from the same selection.

3. Temperature and precipitation need different treatment

An elevation trend works at the daily scale for temperature, and only at the climatological scale for precipitation. Every dataset I reviewed is consistent with this.

Temperature. The lapse rate is physically constrained and changes smoothly from day to day, so a daily regression on elevation, with kriging or weighting of the residuals, works well. The main difficulty is inversions, frequent in the Po Valley and in Alpine valleys in winter. A single lapse rate for a large domain misses them. Pavan et al. (2026) show that fitting a piecewise (two-slope) lapse rate in each of 28 sub-areas gives lower cross-validation errors, every year, than one regression for the whole domain. In our own dataset the lapse rate is domain-wide, re-estimated daily, and temperature errors grow from 1.4 °C below 500 m to 2.1 °C above 1,500 m. That is the price.

Precipitation. On a single day the relation with elevation is weak and changes with the weather. Kirshbaum et al. (2018) explain why:

• in unblocked flow, convection starts on windward slopes;
• in blocked flow, it starts upstream or downstream of the barrier;
• in summer, thermally driven flows make crests hot spots of convective initiation.

The same mountain can be wetter or drier than the plain depending on the day. When we tried a daily elevation trend for precipitation, the slopes were unstable and, multiplied by the elevation of 1 km cells, produced an obviously artificial topographic imprint. We dropped it.

At the climatological scale the relation is robust, mostly because it rains more often at higher elevations. This is why the best precipitation products use elevation only for the normals:

• the anomaly method: interpolate monthly normals with an elevation model, then interpolate the anomalies (or daily fractions) separately. APGD (Isotta et al., 2014), Crespi et al. (2021) and IH-GAR (Manara et al., 2026) work this way;
• PRISM-like local regression for the normals (Daly et al., 2008; Crespi et al., 2018), with weights that account for distance, elevation, slope, aspect and distance from the sea.

The anomaly method has a further advantage that a reviewer of our paper rightly stressed: anomalies are spatially smoother than absolute values, so they are much less sensitive to missing stations.

4. The main families of methods

Method	Weights from	Elevation	Uncertainty per cell	Typical use
Inverse distance, Shepard	Prescribed distance functions	Only if added as a trend	No	ARCIS precipitation (Pavan et al., 2019)
Angular distance weighting (SPHEREMAP)	Prescribed, with angular correction	No	Yes, with the Yamamoto (2000) method	EMO-5 and EMO-1 (Thiemig et al., 2022)
Local weighted regression (PRISM-like)	Prescribed kernels on distance and terrain similarity	Explicit	No	Precipitation normals (Crespi et al., 2018)
Anomaly method	Regression for normals, distance weighting for anomalies	Through the normals	No	APGD, IH-GAR
Ordinary kriging	Variogram fitted to the data	No	Kriging variance	Precipitation, when no drift is reliable
Detrended kriging or kriging with external drift	Variogram of residuals	Through the trend	Kriging variance	Daily temperature

Two remarks on kriging, since I am a long-time user of it in GEOframe (Bancheri et al., 2018):

• It is often lumped together with inverse distance weighting, as in Daly’s tables. That throws away its most useful part: the variogram, which is a description of the data.
• Its variance per cell is a real advantage, but not a complete one. For a given variogram it depends only on station geometry, not on how different the neighbouring values are. Thiemig et al. (2022) chose the Yamamoto (2000) estimator for this reason. Report the kriging variance together with cross-validation errors, not instead of them.

5. Read the variogram

The variogram is not just a step in the kriging machinery; it is information about the field. Its three parameters have a direct reading:

• Nugget: variability at scales below the station spacing, plus measurement error. A large nugget says the network does not resolve the field; denser sampling would be needed.
• Sill: the total variance of the field. A sill that keeps growing with distance signals a trend that should be modelled.
• Range: the distance over which values stay correlated. If it is longer than the typical distance to the stations used, the network resolves the dominant structure that day. If it is shorter, kriging returns towards the mean, and its variance says so.

Fitting a variogram every day is possible, and in our case it is done automatically with particle swarm optimization over five models. Two cautions:

• automatic daily fitting can be unstable; Thiemig et al. (2022) excluded kriging from EMO-5 partly for this reason. Keep fit diagnostics and flag bad days;
• a variogram fitted over a whole domain assumes the same structure everywhere. In large and varied domains, regional variograms may be needed.

6. Evaluate the errors honestly

Leave-one-out cross-validation (LOO) is the most useful error estimate we have: remove a station, estimate it from the others, repeat for all stations. Used carefully, it is sound. Its limits are well known (Daly, 2006):

• It sees only places with stations. Terrain without stations, often the highest, produces no error estimate at all.
• It is optimistic with clustered stations, because a twin station is always there to predict the one removed.
• Do not use it to tune and then to evaluate. If bandwidths or smoothing parameters are chosen by minimizing LOO error, the same LOO error is optimistic, and the selection favours smooth fields. Removing stations that the method predicts poorly lowers the reported error further.
• Compare methods only on the same network. Station-by-station comparison, with the same neighbours, is a valid paired test. MAE values from different networks, periods or wet-day thresholds are not comparable: in Daly’s example, removing the hard high-elevation stations lowered the MAE from 49.7 to 17.2 mm.

Complements that make an error assessment convincing:

1. Stratify LOO by station support and elevation. In the Po dataset, below 500 m, precipitation RMSE rises from 5.4 mm where the seventh-nearest station is within 10 km to 13.9 mm where it is 50–75 km away.
2. Check that the kriging variance is calibrated: standardized LOO errors, (z − ẑ)/σ, should have mean near 0 and variance near 1.
3. Withhold a block of stations, for example all those above 1,500 m, to test extrapolation (Daly’s Table V).
4. Thin the network in a data-rich period, deleting half the stations at random many times, to measure sensitivity to missing data.
5. Read residuals against predicted values, not observed ones. Kriging always flattens peaks, so residuals against observations always tilt and prove nothing.
6. Validate at the time step you distribute. Monthly or annual statistics let errors of opposite sign cancel; a daily product needs daily validation.
7. Remember what no gauge comparison can show. Snow undercatch is shared by all gauges, so the reference carries the same bias as the estimate. Only independent data, such as discharge, snow or radar, can reveal it, and each is confounded with something else.

7. Report what the grid contains

The last guideline is about communication. Daly’s advice to dataset developers is to state strengths and limits frankly. In practice:

• say “provided on a 1 km grid”, not “1 km resolution”;
• release daily station counts and support maps (distance to the nearest stations);
• release the per-cell uncertainty and the station-wise cross-validation errors, as separate products;
• release the daily variogram parameters and lapse rates;
• release station metadata and, where licences allow, the quality-controlled series;
• state the known biases: undercatch, high-elevation underestimation, valleys where the field is smoother than reality.

A user who knows, for every cell and every day, how much local information the value contains can decide whether it is good enough. That is more useful than any single number for “effective resolution”.

A short checklist

• Daily active stations and distance maps computed
• Station versus DEM elevation distribution plotted
• Time of observation verified
• Elevation used daily for temperature, climatologically for precipitation
• Variogram parameters read and saved, fit diagnostics kept
• LOO not used for tuning, or tuning declared
• LOO stratified by support and elevation; blocked and thinning tests done
• Validation at the distributed time step
• Uncertainty, errors, variograms and metadata released with the fields

References

• Bancheri, M., Serafin, F., Bottazzi, M., Abera, W., Formetta, G., Rigon, R. (2018). The design, deployment, and testing of kriging models in GEOframe with SIK-0.9.8. Geosci. Model Dev., 11, 2189–2207. doi:10.5194/gmd-11-2189-2018
• Crespi, A., Brunetti, M., Lentini, G., Maugeri, M. (2018). 1961–1990 high-resolution monthly precipitation climatologies for Italy. Int. J. Climatol., 38, 878–895. doi:10.1002/joc.5217
• Crespi, A., Matiu, M., Bertoldi, G., Petitta, M., Zebisch, M. (2021). A high-resolution gridded dataset of daily temperature and precipitation records (1980–2018) for Trentino-South Tyrol. Earth Syst. Sci. Data, 13, 2801–2818. doi:10.5194/essd-13-2801-2021
• Daly, C. (2006). Guidelines for assessing the suitability of spatial climate data sets. Int. J. Climatol., 26, 707–721. doi:10.1002/joc.1322
• Daly, C., et al. (2008). Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous United States. Int. J. Climatol., 28, 2031–2064. doi:10.1002/joc.1688
• Isotta, F. A., et al. (2014). The climate of daily precipitation in the Alps: development and analysis of a high-resolution grid dataset from pan-Alpine rain-gauge data. Int. J. Climatol., 34, 1657–1675. doi:10.1002/joc.3794
• Kirshbaum, D. J., Adler, B., Kalthoff, N., Barthlott, C., Serafin, S. (2018). Moist orographic convection: physical mechanisms and links to surface-exchange processes. Atmosphere, 9, 80. doi:10.3390/atmos9030080
• Manara, V., et al. (2026). A new daily high-resolution gridded precipitation dataset for the Italian Greater Alpine Region (1951–2023). J. Hydrol.: Reg. Stud., 67, 103776. doi:10.1016/j.ejrh.2026.103776
• Pavan, V., et al. (2019). High resolution climate precipitation analysis for north-central Italy, 1961–2015. Clim. Dyn., 52, 3435–3453. doi:10.1007/s00382-018-4337-6
• Pavan, V., et al. (2026). A new operational dataset of gridded minimum and maximum temperature over north-central Italy 1991 to present. Climate Services, 43, 100689. doi:10.1016/j.cliser.2026.100689
• Salehi, H., et al. (2026). A 30-year 1-km daily precipitation and air temperature dataset for the Po River District (Italy). Earth Syst. Sci. Data Discuss. doi:10.5194/essd-2026-467
• Thiemig, V., et al. (2022). EMO-5: a high-resolution multi-variable gridded meteorological dataset for Europe. Earth Syst. Sci. Data, 14, 3249–3272. doi:10.5194/essd-14-3249-2022
• Yamamoto, J. K. (2000). An alternative measure of the reliability of ordinary kriging estimates. Math. Geol., 32, 489–509. doi:10.1023/A:1007577916868

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)