Equations
1 Formal model description
In this section we give a formal description of the problem that is solved by Ribasim. The problem is of the form
\[ \frac{\text{d}\mathbf{u}}{\text{d}t} = f(\mathbf{u},p(t),t),\quad t \in [t_0,t_\text{end}], \]
i.e. a system of coupled first order ordinary differential equations, with initial condition \(\mathbf{u}(t_0)= \mathbf{u}_0\) and time dependent input data denoted by \(p(t)\).
The model is given by a directed graph, consisting of a set of nodes (or vertices) \(V\) and edges \(E\). Let \(V\) be the set of node IDs and let \(E\) be the set of ordered tuples \((i,j)\) meaning that node \(i\) is connected to node \(j\).
We can split the set of nodes into two subsets \(V = B \cup N\), where \(B\) is the set of basins and \(N\) is the set of non-basins. The basins have an associated storage state and the non-basins dictate how water flows to or from basins.
\(\mathbf{u}(t)\) is given by all the states of the model, which are (currently) the storage of the basins and the integral terms of the PID controllers, the latter being explained in PID equations.
Given a single basin with node ID \(i \in B\), the equation that dictates the change of its storage over time is given by
\[ \frac{\text{d}u_i}{\text{d}t} = \sum_{(i',j') \in E | j' = i} Q_{i',j'} - \sum_{(i',j') \in E | i' = i} Q_{i',j'} + F_i(p,t). \]
Here \(Q_{i,j}\) is the flow along an edge, where the graph direction dictates positive flow. So the first term denotes flow towards the basin, the second one denotes flow away from the basin, and the third term denotes external forcing. \(F_i(p,t)\) is given by input data, and \(Q_{i' ,j'}\) is determined by the type of nodes that connect to that edge.
The various node and forcing types that the model can contain are explained in the section Natural water balance terms.
In general a model has more nodes than states, so in the Julia core there is a distinction between node indices and state indices. For simplicity these are treated as equal in the documentation when it comes to basins and their storage.
1.1 The Jacobian
The Jacobian is a \(n\times n\) matrix where \(n\) is the number of states in the simulation. The Jacobian is computed either using finite difference methods or automatic differentiation. For more details on the computation of the Jacobian and how it is used in the solvers see numerical considerations.
The entries of the Jacobian \(J\) are given by \[ J[i,j] = \frac{\partial f_j}{\partial u_i}, \]
hence \(J[i,j]\) quantifies how \(f_j\), the derivative of state \(j\) with respect to time, changes with a change in state \(i\). If a node creates dependendies between basin storages (or other states), then this yields contributions to the Jacobian. If \(j\) corresponds to a storage state, then
\[ J[i,j] = \sum_{(i',j') \in E | j' = i} \frac{\partial Q_{i',j'}}{\partial u_i} - \sum_{(i',j') \in E | i' = i} \frac{\partial Q_{i',j'}}{\partial u_i}, \]
Most of these terms are always \(0\), because a flow over an edge only depends on a small number of states. Therefore the matrix \(J\) is very sparse.
For many contributions to the Jacobian the derivative of the level \(l(u)\) of a basin with respect to its storage \(u\) is required. To get an expression for this, we first look at the storage as a function of the level:
\[ u(l) = \int_{l_0}^l A(\ell)d\ell. \]
From this we obtain \(u'(l) = A(l)\) and thus \[ \frac{\text{d}l}{\text{d}u} = \frac{1}{A(u)}. \]
The presence of division by the basin area means that areas of size zero are not allowed.
2 Numerical solution
Ribasim uses OrdinaryDiffEq.jl to provide a numerical solution to the water balance equations. Changes to forcings or parameters such as precipitation, but also the allocated water abstraction is managed through the use of callback functions (SciML Development Team 2022). In a coupled run, the exchanges with MODFLOW 6 are also managed via the use of a callback function. For more a more in-depth discussion of numerical computations see Numerical considerations.
3 Performance
There are many things that can influence the calculations times, for instance:
- Solver tolerance: By default we use absolute and relative tolerances of
1e-6and1e-5respectively. - ODE solvers: The
QNDFmethod we use is robust to oscillations and massive stiffness, however other solvers should be tried as well. - Forcing: Every time new forcing data is injected into the model, it needs to pause. Moreover, the larger the forcing fluxes are, the bigger the shock to the system, leading to smaller timesteps and thus longer simulation times.
Similarly to other models that solve using a system of equations, like MODFLOW 6, if one Basin takes longer to converge due to extreme forcing, or bad schematization, the system as a whole need to iterate longer. It is important to be mindful of this, as poor schematization choices can slow down the model.