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.
