Showing posts sorted by date for query Java. Sort by relevance Show all posts
Showing posts sorted by date for query Java. Sort by relevance Show all posts

Wednesday, October 7, 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

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.

Saturday, June 13, 2026

For future working: A note for a comprehensive object oriented implementation of hydrological dynamical systems

Here is an ambition worth stating plainly: a single engine that can solve any hydrological dynamical system, where building a new model means instantiating a few new classes rather than touching the engine at all. Not a model — a kit. The bricks are the storages and the fluxes; the engine assembles them, differentiates them, and solves them; and the same solvers are reused no matter what physics you bolt on. You add a new process by extending the kit, never by editing the machinery — open to extension, closed to modification — and you do it by writing to interfaces rather than to concrete classes, so the solver never needs to know, or care, what it is solving. The reward is a system that is at once generic, efficient, and expandable, with strikingly little code to maintain.
NOt forgetting that this is very preliminary material, food for thinking, probably bugged material, please find:

This is not a new creed. It is the generic-programming philosophy that Berti set out for scientific computing two decades ago — efficient, reusable components built on the right abstractions, so that code stops being rewritten for every variant — and it is exactly the spirit in which Niccolò and I built WHETGEO-1D (Tubini and Rigon, 2022), where a ClosureEquation interface and a factory let you swap van Genuchten for Brooks–Corey without disturbing a line of the solver. What I want to do here is take that same philosophy and lift it from the one-dimensional soil column up to the whole topology of a dynamical system, and the thing that makes it possible is an old friend.

The friend is the Extended Petri Net. A hydrological dynamical system is one: storages are places, fluxes are transitions, and a third set of objects, the controllers, complete the causal wiring. Follow that picture all the way down into the software, and the thing every modeller secretly dreads — assembling, by hand, the system of equations to hand to a solver — almost disappears.

The grammar, and why "universal" is not a boast

Before the software, the picture. An EPN is drawn with a small, fixed vocabulary, and once you have it the claim of universality stops being rhetorical. A place is a circle, a transition a square; a square carrying a black dot is a driven input, the rain falling into the net; a dashed square is the outlet to the world. A diamond is a splitter, where one quantity divides into branches whose fractions sum to one — precipitation parting into snow and rain. A dotted circle is a collector, a summing junction with no storage of its own. And a little histogram marks a flux computed by convolution, a unit hydrograph. That is very nearly the whole alphabet.

Why believe it is enough? Because Marialaura, Anna De Nardi and I drew all forty-six models of the MARRMoT collection (Knoben et al., 2019) in exactly this vocabulary — Collie, GR4J, TOPMODEL, HBV, Sacramento, VIC, the Tank cascades, the whole zoo — and each one turns out to be just a different choice of places, fluxes and wiring. If a single alphabet spells every word in the dictionary, it is the right alphabet. The classes that follow are nothing more than this grammar given types.

One class of storages, many kinds of flux

The first thing the Petri-net view tells you is an asymmetry. A place is universal: it is a state variable whose time variation we study, and one storage is, as an object, exactly like another. A bucket of soil water, a snowpack, a channel reach — all of them are just \(S_i\) with a balance to satisfy. So in an object-oriented design, storages are a single class.

Fluxes are the opposite. A flux can be a known time series (a forcing), or a constitutive law that declares a mathematical form and a dependence on the state. So a transition wants to be an interface with concrete implementations — an external flux here, a Darcy-Buckingham law there, a power-law discharge somewhere else. Each flux knows two things that matter for what follows: where it moves mass (its source and target places) and which state variables it reads. Those two are not the same, and keeping them apart is the whole trick.

Once you take that seriously, the "many kinds of flux" become a small, nameable family. There is the external flux (a forcing) and the constitutive flux (a law of the state); the controller, a quantity derived from the state that gates a flux without carrying any mass; the splitter and the collector of the grammar above; the stoichiometric flux, one rate metered into several currencies, which we will need the moment evapotranspiration appears; and the convolution flux, the one transition that carries a memory. They all implement the same interface. The only place the software must be careful is with the threshold laws that conceptual models love — the \(\text{if } S>S_{\max}\), the \(\max(\cdot,0)\) — which are not differentiable at the kink; those enter in regularised form, the same medicine the soil-water retention curve already takes.

The two graphs

An EPN actually carries two graphs over the same nodes, and conflating them is, I think, why generic model engines so often turn into a tangle.

The first is the incidence graph — the Petri net proper. It records which flux moves water between which storages, and it is summarised by the incidence matrix \(\mathbf{C}\), with \(+1\) where a flux enters a storage and \(-1\) where it leaves. From it, the entire model collapses into one line:

\[ \frac{\mathrm{d}\mathbf{S}}{\mathrm{d}t} \;=\; \mathbf{C}\,\mathbf{Q}(\mathbf{S},t). \]

That is conservation, and nothing more. The second graph is the causal one, directed on the storages alone: there is an arc from \(S_j\) to \(S_i\) whenever the equation for \(S_i\) contains a flux that reads \(S_j\). This is the graph that decides how the equations are coupled.

And here is where the controllers finally make sense. A controller is an arc that lives in the second graph but not in the first — a flux that is gated or modulated by some storage without any mass passing through that storage. The read arcs and inhibitor arcs of Petri-net theory. They carry no water, so they are absent from \(\mathbf{C}\); but they carry causation, so they are present in the dependency graph and therefore in the Jacobian. That is exactly what we mean when we say controllers "complete the causal structure of the model".

Assembling the system the solver can eat

If you want an implicit step — and for cyclically coupled, stiff budgets you almost always do — backward Euler turns the master equation into a residual,

\[ \mathbf{F}(\mathbf{S}^{\,n+1}) \;:=\; \bigl(\mathbf{S}^{\,n+1}-\mathbf{S}^{\,n}\bigr) -\Delta t\,\mathbf{C}\,\mathbf{Q}(\mathbf{S}^{\,n+1},t^{\,n+1}) \;=\; \mathbf{0}, \]

and a Newton method needs its Jacobian,

\[ \mathbf{J} \;=\; \mathbf{I} \;-\; \Delta t\,\mathbf{C}\,\frac{\partial \mathbf{Q}}{\partial \mathbf{S}}. \]

This is the passage people worry about: how do you build this nonlinear system automatically for an arbitrary model? The answer is that the structure is already in the objects. The sparsity of \(\partial\mathbf{Q}/\partial\mathbf{S}\) is just the list of dependencies each flux declared. So the assembler becomes a single scatter loop over the fluxes — exactly like assembling element matrices in a finite-element code. You walk the fluxes once; each one drops its value into the balance of its two endpoints and its derivative into the matching rows of the Jacobian. You never build \(\mathbf{C}\) or the flux Jacobian as dense matrices; they assemble themselves, entry by entry.

The derivatives I would get from automatic differentiation. If each flux is written over a differentiable scalar (in Java, Hipparchus' DerivativeStructure), one forward evaluation returns both the flux value and its partials. Concretely, a linear reservoir is just

DerivativeStructure evaluate(AdContext ctx, double t) {
    return ctx.value(from).divide(tau);   // Q = S_from / tau
}

and the engine differentiates it for you. The consequence is the thing I actually care about: you write the physics once, per flux, and never again touch the solver. Adding a new constitutive law costs you its forward formula and nothing else.

Loops, sequences, and where the parallelism hides

There is one more gift in the causal graph. Run Tarjan's algorithm on it and you get its strongly connected components. A component with more than one storage is a loop — a knot of mutually dependent equations that has to be solved simultaneously. Every other component is a single equation that can be solved in turn. The graph of components is a DAG, and its topological order is simply the order in which to solve. In the algebra of the Jacobian this is a block-triangular structure; it is the same decomposition that equation-based modelling languages call BLT. And components sitting at the same level of that DAG do not depend on one another, so they can go to different threads. The parallelism was never something we had to impose — it was sitting in the dependency graph all along, waiting to be read off.

The solver, finally, gets to be ignorant. It sees only a residual, a Jacobian, and a starting guess. That is enough to let an off-the-shelf Newton–Raphson handle the easy blocks and a Nested Newton (Casulli–Zanolli) handle the monotone Richards-type ones — without either of them ever knowing what they are solving.

Is it efficient? Yes — but structurally

The natural worry is that all this generality must cost something at runtime. It does not, or rather, the cost lands exactly where it should. The only expensive operation — the implicit Newton solve — happens only inside the loops; the rest of the network, which is most of it, is cheap explicit updates. And the whole structural analysis — finding the loops, the solution order, the sparsity of the Jacobian — depends only on the graph, which never changes in time. So you do it once, at the start, and pay nothing more for it on the millions of timesteps that follow. The efficiency is not a clever trick in the inner loop; it is a consequence of letting the structure tell you where the hard work actually is. The one honest caveat is that you must respect that structure: solve everything monolithically and you throw the gift away. Couple only what is genuinely stiff, and let the rest run loose and parallel.

Net3, upgraded

This is where Francesco Serafin's Net3 comes in, and the fit is almost too neat — Net3 grew out of the same EPN picture. Net3 takes a model as a directed acyclic graph and runs it in parallel, which is wonderful for a river network (a tree is a DAG) but awkward for coupled dynamics, because coupling makes loops, and a DAG cannot hold a loop. The repair is the one move we already have: take each loop, each strongly connected component, and crush it into a single super-node whose insides are the simultaneous solve. The graph of super-nodes is a DAG again, and Net3 schedules it exactly as before. The coupled solving lives inside the node; the parallel routing lives between nodes. Net3 keeps doing what it is good at, and inherits the one thing it could not do.

Three budgets, one net — and the trouble with ET

The original EPN paper already hints that we might want to solve more than water: the energy budget, the carbon budget, all at once. The beauty is that the equation is the same, \(\dot{\mathbf{S}}_b = \mathbf{C}_b\,\mathbf{Q}_b\); only the meaning of the parameters changes. Water places hold storages, energy places hold internal energy, carbon places hold pools. What ties them together are the controllers — quantities derived from one budget's state that reach into another. Temperature, born of the energy budget, governs evapotranspiration in the water budget; soil moisture, born of the water budget, governs the thermal conductivity in the energy budget and the stomata in the carbon budget. Each of these is, once again, an arc that lives in the causal graph but in nobody's incidence matrix.

And then there is the genuinely awkward case, the one that makes the whole thing interesting: a single quantity that belongs to two budgets at once. Evapotranspiration is the textbook offender. It is a loss of water, at rate \(E\); it is also a loss of energy, at rate \(\lambda E\), the latent heat carried away. How do you stop the two budgets from quarrelling about how much of it happened?

The clean answer is to stop thinking of separate budgets and to write one net over all the currencies at once, letting the incidence matrix carry not just \(\pm 1\) but real stoichiometric coefficients. Then ET is a single column with an entry of \(-1\) in the water rows and \(-\lambda\) in the energy rows. You evaluate it once — it has one rate — and drop that one rate into both budgets, scaled appropriately. The consequence is the thing I find quietly satisfying: the latent heat is exactly \(\lambda\) times the water loss at every step of the iteration, not just at the end. The two budgets cannot disagree, because there is only one number. Conservation across currencies stops being something you check and becomes something the structure guarantees.

If you have ever wondered what Penman–Monteith really is, this is it. Penman solves for the surface temperature and the evaporation rate together, between the energy and the mass budgets — which is precisely solving one of these little coupled super-nodes with a shared ET column and a shared temperature controller. The multi-currency net is just Penman–Monteith let off its leash: the same closure, now for any number of budgets, solved numerically rather than by hand. And when all the controllers descend from a single free energy, the cross-couplings ought to be symmetric — Onsager reciprocity — which gives a quiet, structural way to check that the coupled model is thermodynamically honest.

Routing, travel times, and the unit hydrograph

One flux refuses to be memoryless, and it rewards a closer look, because its kernel is not arbitrary — it is a travel-time distribution. The instantaneous unit hydrograph carries the effective rainfall to the outlet by convolution, and the cleanest way to hold it is not the IUH \(f\) itself but its integral, the S-function

\[ s_f(T) \;:=\; \int_0^T f(y)\,\mathrm{d}y, \]

which is, up to the catchment area, the cumulative distribution of travel times — the very residence-time object that runs through the age-ranked budget story (Rigon, Bancheri and Green, 2016). Write the discharge against travel time and a single rain record contributes differences of \(s_f\) over its own length; the records then simply superpose, because the routing is linear and time-invariant. That is the entire numerical recipe, worked out impulse by impulse in the illustrated guide (Rigon et al., 2022b; Rigon et al., 2016).

The interesting part is what the kernel's origin decides. A parametric hydrograph — an exponential, which is just a single linear reservoir; or a Nash cascade — reduces exactly to a little chain of reservoir places, so it folds back into the basic vocabulary and asks for no new type. A geomorphological hydrograph does not: the width function, read off a digital elevation model as the area of the catchment at each flow distance from the outlet and mapped to time by a velocity, gives \(s_f\) straight from the shape of the basin and is no finite chain of reservoirs. That is exactly why the convolution flux has to be a first-class citizen rather than sugar over a cascade. And the linearity is a genuine assumption: when the hydrograph changes shape with the storm — an event-specific GIUH — the kernel becomes a controller of the forcing, and the clean convolution gives way to the general, time-varying case.

There is a quiet bonus for anyone thinking of running the model in real time. Because each incoming rain record commits the discharge for the next several steps — water already fallen, travel times already fixed — the convolution hands you a short forecast for nothing, unmodifiable by rain that has not yet arrived. Feeding the engine live data is then only a matter of swapping the file behind a forcing for a streaming source; the engine never notices whether the rain fell last year or a minute ago.

Why I like this

What pleases me here is how much the architecture buys by simply refusing to mix two things up. Conservation lives in the incidence matrix. Causation lives in the dependency graph. The physics lives in the fluxes, written one at a time. And the assembler — the small piece of code that compiles all of it into a system of equations — is one loop. It is the kind of separation that, once you see it, makes you wonder why the equations ever felt like the hard part. They were never the hard part. The hard part was deciding what was a place, what was a flux, and what was only a controller.

And that, in the end, is the whole point of the kit. Every new model — a different catchment, a snow scheme, a coupled carbon budget — is a handful of new flux classes implementing the same interface, dropped into the same engine, solved by the same Newton. Nothing in the machinery changes; the machinery was closed to modification from the start. The code stays small not because we were clever line by line, but because we let the abstraction carry the weight. That is what generic programming promised for scientific computing, and it is satisfying to watch hydrology turn out to be such a natural place to collect on the promise.

References

Bancheri, M., Serafin, F., and Rigon, R. (2019). The Representation of Hydrological Dynamical Systems Using Extended Petri Nets (EPN). Water Resources Research, 55(11), 8895–8921. doi:10.1029/2019WR025099.

Berti, G. (2000). Generic Software Components for Scientific Computing. PhD thesis, BTU Cottbus. See also Berti, G. (2006), GrAL — the grid algorithms library, Future Generation Computer Systems, 22(1–2), 110–122.

Knoben, W. J. M., Freer, J. E., Fowler, K. J. A., Peel, M. C., and 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, 2463–2480. doi:10.5194/gmd-12-2463-2019.

Rigon, R., Bancheri, M., Formetta, G., and de Lavenne, A. (2016a). The geomorphological unit hydrograph from a historical-critical perspective. Earth Surface Processes and Landforms, 41(1), 27–37. doi:10.1002/esp.3855.

Rigon, R., Bancheri, M., and Green, T. R. (2016b). Age-ranked hydrological budgets and a travel time description of catchment hydrology. Hydrology and Earth System Sciences, 20(12), 4929–4947. doi:10.5194/hess-20-4929-2016.

Rigon, R., Formetta, G., Bancheri, M., Tubini, N., D'Amato, C., David, O., and Massari, C. (2022a). HESS Opinions: Participatory Digital eARth Twin Hydrology systems (DARTHs) for everyone — a blueprint for hydrologists. Hydrology and Earth System Sciences, 26, 4773–4800. doi:10.5194/hess-26-4773-2022.

Rigon, R., Franceschi, S., Formetta, G., Bancheri, M., and Tubini, N. (2022b). An illustrated guide to IUH/GIUH estimation. Authorea preprint. doi:10.22541/au.164192110.08629205/v1.

Serafin, F. (2019). Enabling Modeling Framework with Surrogate Modeling Capabilities and Complex Networks (the Net3 subsystem). PhD thesis, University of Trento. See also Serafin, F., David, O., Carlson, J. R., Green, T. R., and Rigon, R. (2021), Environmental Modelling & Software, 146, 105231.

Tubini, N. and Rigon, R. (2022). Implementing the Water, HEat and Transport model in GEOframe (WHETGEO-1D v.1.0): algorithms, informatics, design patterns, open science features, and 1D deployment. Geoscientific Model Development, 15, 75–104. doi:10.5194/gmd-15-75-2022.

Friday, February 6, 2026

The Hydrology 2026 Lab

 The lab component makes up nearly half of the course, following the motto:

"Learning by doing."

Throughout the lab, you will conduct at least three key numerical experiments:

  • Time Series Analysis – Exploring various data elaborations with various Jupyter Notebooks and a little of Python
  • Intensity-Duration-Frequency (IDF) Curves – Estimating rainfall intensity over different time scales with various Jupyter Notebooks and a little of Python, as well 
  • Infiltration Experiments – Investigating soil absorption dynamics using the WHETGEO system
  • Evaporation & Transpiration Experiments – Understanding water loss processes in different conditions using the GEOSPACE ssytem

Resources

🔹[Vimeo Showcase – General Lab Videos]

🔹[OSF Repository – Lab Materials]

🔹[Theory and Concepts here]

📌 Detailed videos and materials for each experiment are listed below.




📅 26 February 2026 – Installation and Introduction to working with Jupyter and Notebooks

📅 6 March 2026 Reading a file

📅 9 March 2026 How to read and plot data

📅 12 March 2026 -  San Martino Reprise


Interpolating the Gumbel distribution to annual precipitation maxima
Where to find the data
  • Finding the data for the Time Series Analysis (Vimeo2025)

📅 13 April 2026 

Below, you'll discover the resources for the second part of the lab, which encompass instructions for utilizing Object Modelling System version 3 (OMS3) models such as WHETGEO1D and GEOET. Within the OMS Projects, you'll find a directory named Jupyter_Notebook, housing sets of notebooks designed to guide you through handling both input and output data for these programs.

OMS and further Anaconda Installations
A general introduction to some OMS3 concepts

📅 16 April 2026 

📅 28 April 2026 

 Introduction to Infiltration with WHETGEO 1D. WHETGEO specifically focuses on simulating the coupled processes of water and heat transport, which are crucial for understanding: Soil moisture dynamics, Evapotranspiration processes, Energy balance at the land surface, Groundwater-surface water interactions, Snow and ice processes (where temperature is critical)

📅 21 May 2026 - Radiation estimation

is preliminary to any Evapotranspiration estimation since our models use it as input. The estimation of radiation, as it results from the theory, needs a few GIS operations to determine the sky view factor and the aspect of the terrain in the point of interest. These operation for the present exercise are skipped and pre-analyzed data are provided. 

📅 4 June 2026 - Estimation of Evapotranspiration


Please below find Notebooks and video related to the estimation of radiation. Theradiationin this class, will be functional to the estimation of evapotranspiration with the Priestley-Taylor and Penman-FAO models,
 Introduction to the estimation of evapotranspiration

Musical Coda





Tuesday, June 3, 2025

OMS Runner Library: Streamlining Hydrological Model Execution

 The OMS Runner Library v1.2.2 represents a significant advancement in hydrological modeling workflow automation, specifically designed to simplify the execution of OMS3 (Object Modeling System) simulations. For hydrologists and water resources engineers working with GEOframe and OMS3, this Python library addresses the seamless integration and execution of simulation models across different computing platforms. What follows assume a lot of knowlege that you can get by looking to some of our Winter Schools or some of our lab classes as  Physical Hydrology (in Italian) or  Biosphere Atmosphere and Climate Interactions. 

What is OMS3?

The Object Modeling System (OMS3) is a Java-based framework widely used in environmental and hydrological modeling. It provides a robust platform for developing, coupling, and executing complex simulation models. However, working with OMS3 often requires dealing with Java classpaths, configuration files, and platform-specific execution commands – tasks that can be time-consuming and error-prone, especially for researchers focused on scientific analysis rather than software engineering.

The Solution: Python Integration

The OMS Runner Library bridges this gap by providing a comprehensive Python interface for OMS3 operations. This is particularly valuable because Python has become the lingua franca of scientific computing, with most hydrologists already familiar with its ecosystem of tools like pandas, matplotlib, and Jupyter notebooks.

The library automatically handles the complexities of Java environment detection, ensuring that Java JDK 11 is properly configured across Windows, macOS, and Linux systems. This cross-platform compatibility is crucial for research teams working in diverse computing environments, from field laptops running Windows to high-performance computing clusters running Linux.

Please find:

Version 1.2.4

Version 1.2.2

Key Capabilities

One of the library's standout features is its intelligent simulation management. It can automatically discover simulation files within a project, maintain configuration databases, and execute models either individually or in sophisticated batch processing workflows. For hydrologists working with multiple scenarios – such as climate change impact assessments or calibration procedures – the parallel execution capabilities can reduce computational time.

The library supports various execution patterns: sequential processing for dependent simulations, parallel execution for independent model runs, and asynchronous background processing for long-running computations. This flexibility allows researchers to optimize their workflows based on available computational resources and modeling requirements.

Practical Applications

In practical hydrological applications, this translates to significant productivity gains. A researcher studying watershed responses to different precipitation scenarios can now set up dozens of model runs with just a few lines of Python code, monitor their progress through Jupyter notebooks, and automatically collect results for analysis. The library's integration with popular Python data analysis tools means results can be immediately processed, visualized, and shared.

Users can explore more about GEOframe's capabilities and latest developments at the GEOframe blog, where detailed tutorials and case studies demonstrate advanced hydrological modeling workflows.

The comprehensive logging and error handling features are particularly valuable in operational hydrology contexts, where model reliability and traceability are paramount. The library maintains detailed execution histories, facilitates debugging, and provides clear diagnostic information when issues arise.


Wednesday, May 14, 2025

Biosphere Atmosphere Climate Interactions lab 2025

Go to the Theory and Concepts page

This page contains the laboratory materials for the BACI2025 course, specifically covering the sections I oversee. The laboratory sessions build upon concepts from the hydrology class while extending into more advanced simulations, including applications of the GEOSPACE system.  GEOSPACE is based upon GEOframe is an open-source, component-based hydrological modeling framework developed primarily by researchers in Italy. It's designed to be modular and flexible, allowing researchers to combine different hydrological processes and models.

Installations
Generalities about the Object Modelling System

OMS3 (Object Modeling System 3) is a Java-based environmental modeling framework developed by the U.S. Geological Survey (USGS) and Colorado State University. It's designed to facilitate the development, integration, and deployment of environmental and agricultural models.

OMS Project that will be used
WHETGEO specifically focuses on simulating the coupled processes of water and heat transport, which are crucial for understanding: Soil moisture dynamics, Evapotranspiration processes, Energy balance at the land surface, Groundwater-surface water interactions, Snow and ice processes (where temperature is critical)
  • Exploring a WHETGEO project. Everything is explained in the WHETGEO1D_RichardsCoupled_Computational_grid notebook (Vimeo2025)
  •  Please also browse directly the files:
    • The grid file: ex00_grid.csv , Initial Conditions: ex00_ic.csv, Parameters: (for van Genuchten): Richards_VG.csv, Dictionary: dictionary.csv (Vimeo2025)
    • WHETGEO1D.sim (Vimeo2025)
  • Analysis of WHETGEO1D.sim (Vimeo2025)
  • To run WHETGEO use the Runner_WHETGEO.ipynb notebook (Vimeo 2025)
    • Please observe that the oms.py file needs to be in the same folder than the Notebook
    • During the course a more solid solution was provided here.
  • How to visualize and interpret the data through the Visualize_output.ipynb Notebook (Vimeo 2025)
  • Utilities:
    • To proper formatting (according to OMS3) rainfall files, peruse the TimeSeriesFormatter.ipynb Notebook (Vimeo2025)
    • How to create an empty TimeSeries using TimeSeriesCreator.ipynb (Vimeo2025) 
RADIATION ESTIMATION is preliminary to any Evapotranspiration estimation since our models use it as input. The estimation of radiation, as it results from the theory, needs a few GIS operations to determine the sky view factor and the aspect of the terrain in the point of interest. These operation for the present exercise are skipped and pre-analyzed data are provided. 
  • Please see the 05_NET.sim file in the project's simulation folder (Vimeo2025)
  • Inputs are in Input_Radiation.ipynb (Vimeo2025)
  • Outputs are in the Output_Radiation.ipynb (Vimeo2025)
  • The simulation runner is Runner_NetRadiation.ipynb (Vimeo2025)
GEOET  is a group of  components within GEOframe that focuse specifically on evapotranspiration (ET) modeling. At present GEOET includes three main estimator of ET: The Priestley-Taylor, the Penman-FAO and the Prospero one.  
  • Exploring a GEOET project
GEOSPACE, the Soil Plant Atmosphere Continuum Estimator in GEOframe (GEOSPACE) is presented in this paper, in particular the one-dimensional development, GEOSPACE-1D. GEOSPACE  (D'Amato et al., 2025) is a sophisticated modeling framework designed to simulate the complex interactions within the soil-plant-atmosphere continuum (SPAC) system.
  • Exploring a GEOSPACE project

REFERENCES

Tuesday, February 4, 2025

The Hydrology Class Lab 2025

 The lab component makes up nearly half of the course, following the motto:

"Learning by doing."

Throughout the lab, you will conduct at least three key numerical experiments:

  • Time Series Analysis – Exploring various data elaborations with various Jupyter Notebooks and a little of Python
  • Intensity-Duration-Frequency (IDF) Curves – Estimating rainfall intensity over different time scales with various Jupyter Notebooks and a little of Python, as well 
  • Infiltration Experiments – Investigating soil absorption dynamics using the WHETGEO system
  • Evaporation & Transpiration Experiments – Understanding water loss processes in different conditions using the GEOSPACE ssytem

Resources

🔹[Vimeo Showcase – General Lab Videos]

🔹[OSF Repository – Lab Materials]

🔹[Theory and Concepts here]

📌 Detailed videos and materials for each experiment are listed below.



 2025-03-03 Introduction to working with Jupyter and Notebooks

2025-03-10 How to read and plot data
2025-03-25 San Martino Reprise
2025-03-31 Gumbel derivation

Interpolating the Gumbel distribution to annual precipitation maxima
2025-04-08 Where to find the data
  • Finding the data for the Time Series Analysis (Vimeo2025)
2025-05-05 - Installations
  • Installing the GEOframe_vincennes environment
  • Installing the Java Development toolkit V11
2025-05-06  

Below, you'll discover the resources for the second part of the lab, which encompass instructions for utilizing Object Modelling System version 3 (OMS3) models such as WHETGEO1D and GEOET. Within the OMS Projects, you'll find a directory named Jupyter_Notebook, housing sets of notebooks designed to guide you through handling both input and output data for these programs.

A general introduction to some OMS3 concepts
2025-05-11
 Introduction to Infiltration with WHETGEO 1D. WHETGEO specifically focuses on simulating the coupled processes of water and heat transport, which are crucial for understanding: Soil moisture dynamics, Evapotranspiration processes, Energy balance at the land surface, Groundwater-surface water interactions, Snow and ice processes (where temperature is critical)
2025-06-03

    Radiation estimation is preliminary to any Evapotranspiration estimation since our models use it as input. The estimation of radiation, as it results from the theory, needs a few GIS operations to determine the sky view factor and the aspect of the terrain in the point of interest. These operation for the present exercise are skipped and pre-analyzed data are provided. 
  • Please see the 05_NET.sim file in the project's simulation folder (Vimeo2025)
  • Inputs are in Input_Radiation.ipynb (Vimeo2025)
  • Outputs are in the Output_Radiation.ipynb (Vimeo2025)
  • The simulation runner is Runner_NetRadiation.ipynb (Vimeo2025)
 

    Thursday, December 12, 2024

    Positions in snow modelling, Po River basin hydrology, soil-plant-atmosphere interactions and GEOframe system development @UniTrento

     Dear All,

    I am seeking motivated  master graduated interested in working in areas related to Snow modelling, Po River basin hydrology, Soil-plant-atmosphere interactions and GEOframe system development. Below are some exciting thesis opportunities, each with potential for continuation into a Ph.D. program. Post doc positions could be considered as well for appropriate persons.


    1. Snow Modelling (SUNSET PRIN Project)

    This topic focuses on snow dynamics modeling using GEOframe-NewAGE and GEOtop, within the SUNSET PRIN project (Details). Opportunities include fieldwork, guided by Prof. Stefano Ferraris (University of Turin), with Dr. John Mohd Wani as co-supervisor.

    2. Po River Basin Projects (ADBPo Collaboration)

    The Po River Basin thesis topics align with the long-term collaboration with the Basin Authority of River Po (Details). These projects could lead to Ph.D. opportunities and professional roles.

    Topics include:

    - Modeling Romagna Catchments for Drought and Flood Prevention
      Focused on hourly-scale modeling for water management.
    - Co-Supervisor: Ing. Gaia Roati (Po Basin Authority)
    - Includes periods at the Po Basin Authority in Parma.
    - Earth Observation for Po River Basin Calibration
      Systematic use of satellite data to validate and improve GEOframe-NewAGE models.
    - Part of: SpaceItUp! PNRR (Italian Space Agency) and an upcoming ESA project
    - Collaboration: Ing. Hossein Salehi, Fondazione Edmund Mach, and Prof. Manuela Girotto (UC Berkeley).

    3. Land-Surface Interactions with the use of the GEOSPACE system and its development (EPFL collaboration, ESA Projects, EU Projects)
    - Modeling Transpiration and Soil-Atmosphere Interactions
      Utilize GEOframe-NewAGE and GEOSPACE for basin-scale modeling.
    - Co-Supervisor: Dr. Concetta D’Amato (EPFL, Sion Campus)
    - Includes a potential study period in Sion.
    - Depending on the specific topic other collaboration should be envisioned

    Topics include:
    - Estimating the effects of evaporation and transpiration at Po scale, integrating GEOSPACE with Earth Observation.
    - Understanding the effects of soil evolution under the action of biota and under global warming.
    - New parameterizations of the atmosphere - plant interactions

    4. Informatics-Oriented GEOframe Development (SIM Project)
    - Integrating Large Language Models with GEOframe
      Explore the potential of Generative AI for improving user interaction and programming within GEOframe.
    - Transforming GEOframe into a DARTH
      Enhance the GEOframe infrastructure and Object Modelling System (OMS) codebase.
    - Co-Supervisors: Prof. Giuseppe Formetta and Dr. Olaf David (Colorado State University)
    - Includes a potential study period in Fort Collins, Colorado.


    Additional Information

    All theses involve using and extending GEOframe tools, requiring proficiency in Python and Java. Coding skills are especially critical for informatics-oriented topics.
    These projects provide an excellent foundation for doctoral studies and professional development in hydrology, environmental modeling, and computational science.

    Please feel free to share this information and contact me for further details.

    Best regards, 
    Riccardo Rigon
    riccardo<dot> rigon<@>unitn<dot>it

    The anticipated salary for pre-doctoral and Ph.D. students is €1,350 per month (net), plus an additional €3,000 annually for supplementary activities. For postdoctoral researchers, the net annual salary ranges from €24,000 to €30,000, depending on individual qualifications and experience. Additional income opportunities may also be available. The cost of living in the region is more affordable compared to many other European countries.


    IMPORTANT !!!!

    P.S. - In your response, please specify which of the above proposals you are interested in pursuing, along with a brief explanation of your motivation. Send your CV with your age and gender included. If you do not have prior experience with the GEOframe system,, we kindly request that you first enroll in our  GEOframe School (we are happy to waive your subscription fees). Please complete the enrollment and mention in your communication that you have done so.
    The school has already covered the installation process and some theoretical aspects, but all materials are available online for self-paced learning. The next session will take place in January, and individuals currently in Italy are encouraged to attend in person.