The Kinetic Theory of water flow in porous media - Lausanne seminars

 Work in progress

Daniele Andreis' PhD.

Under construction

 

Gaia Roati's Ph.D.

Under construction 

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 open discussion 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.

A caution that applies to all products built from normals. IH-GAR, APGD and Crespi et al. (2021) obtain daily values by scaling station anomalies with elevation-aware normals. A normal is the number of wet days times their intensity, so this scaling can turn a higher frequency of rain at altitude into heavier storms. Monthly totals stay correct; daily wet-day frequency and intensity at altitude are worth checking before using these products as daily forcing.

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, and ordinary kriging of daily totals spreads small amounts across the wet/dry boundary. In the revision we are testing, on the gauges and within single events, whether elevation acts on the number of wet days or on storm intensity, and we will build precipitation from occurrence and intensity separately. 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 brings in the elevation dependence of precipitation where it acts, on wet-day frequency, on intensity or both, 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.

It also has a cost that is easy to miss when the goal is a daily field. A monthly normal is the number of wet days times their mean intensity. Scaling a station’s daily value by the ratio between the cell’s normal and the station’s normal multiplies each storm by the ratio of wet-day frequencies as well as of intensities. If a ridge is wetter than the valley because it rains there more often, the anomaly method gives the ridge the valley’s wet days with heavier storms. Monthly totals come out right by construction, so monthly validation cannot see the problem; wet-day frequency and intensity can. For monthly climatology this does not matter. For daily forcing of a hydrological model, where the intensity of each storm controls runoff and the split between rain and snow, it does. There is a further reason: the elevation gradient of a monthly normal is an average over the weather types of that month, each with its own gradient (Kirshbaum et al., 2018), so it describes no single storm.

The remedy is to ask the data where elevation acts. Group wet stations into events, and estimate the elevation effect on the probability of rain and on storm intensity from contrasts between stations in the same storm (a mixed model with one random effect per event does this). Then put elevation where it is found: if it acts on occurrence, interpolate occurrence with an elevation drift and intensity without, and combine them. To claim that intensity does not depend on elevation, use an equivalence test against a margin fixed in advance; a non-significant slope is not enough. One practical detail: the product of an occurrence probability and an intensity is an expected value, which smears the wet/dry boundary again. For a daily map, decide which cells are wet by thresholding the probability, or produce conditional simulations.

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
Occurrence–intensity (hurdle) kriging	Variograms of indicator and of log-intensity	In occurrence, intensity or both, as the data show	Kriging variance of each part	Daily precipitation

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. In kriging the spatial structure is estimated from the observations, not assumed through prescribed weights or tuned kernels; the estimate is the best linear unbiased one for that structure; and stations close to each other share their weight instead of doubling it, which no distance-weighting scheme does.
• 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, and score wet-day frequency and intensity separately. 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; for precipitation, placed in occurrence or intensity according to a within-event test
• Wet-day frequency and intensity validated separately, not only totals
• 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