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}], \]
which is a system of coupled first order differential equations.
The model is given by a directed graph, consisting of a set of node IDs (vertices) \(V\) and edges \(E\), consisting of ordered pairs of node IDs. We denote the subset of the nodes given by the Basins \(B \subset V\), and the subset of nodes that prescribe flow \(N \subset V\).
The states \(\mathbf{u}\) of the model are given by cumulative flows since the start of the simulation as prescribed by the nodes \(N\): \[ u_n(t) = \int_{t_0}^t q_n\text{d}t' \quad \forall n \in N, \] as well as by the Basin forcings: \[ u_b^\text{forcing}(t) = \int_{t_0}^t q_b^\text{forcing}\text{d}t' \quad \forall b \in B. \]
Because of this definition, the initial conditions of all states are simple: \[ u_i(t_0) = 0 \quad \forall i. \]
From these cumulative flows, the storage in each Basin can be determined at each point in time: \[ S_b(t) = S_i(0) + S^\text{exact}(t) - u_b^\text{forcing}(t) + \sum_{n\;|\;(n,b)\in E} u(t) - \sum_{n\;|\;(b,n)\in E} u(t), \]
i.e. the storage is given by:
- the initial storage;
- plus the exactly integrated flows (more on that below);
- minus the cumulative outgoing forcings;
- plus the cumulative horizontal inflows;
- minus the cumulative horizontal outflows.
From these storages in combination with the Basin profiles the Basin levels \(h\) are computed. The relationship between the profile and the storage is given by \[ S_b = \int_{h_0}^h A_b(\ell)\text{d}\ell, \]
where \(A_b\) is the linear interpolation of the area as a function of the level. These levels are then inputs for determining the flows prescribed by the nodes \(N\). From this relation it also follows that
\[ \frac{\text{d}h}{\text{d}t} = \frac{1}{A_b}, \] and so areas of zero are not allowed in the Basin profiles.
1.1 The PID control integral state
There’s one other type of state, which is not a cumulative flow but a cumulative error. This is the error integral for PID control, further explained in PID equations.
1.2 Exactly integrating flows to minimize the number of states
The more states the problem has, the more time it takes to solve it. Therefore we want to minimize the number of states. Flows do not have to be states when they can be integrated over time exactly because they do not depend on the other states. This is true for FlowBoundary nodes, and Basin precipitation (which uses a fixed basin area) and drainage.
1.3 The Jacobian
The Jacobian is an \(N \times N\) matrix where \(N\) is the number of states in the simulation. It is computed as part of implicit time stepping methods. There are 2 different methods available for computing this matrix: finite difference 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}, \] i.e. \(J_{i,j}\) quantifies how \(f_j\). the time derivative of state \(j\), changes with respect to changes in state \(i\). Most of these entries are \(0\), because flows in distant parts of the model do not depend on each other.
1.4 The water balance error
The water balance error quantifies how well the water volume in the model is conserved for each Basin over an output save period, i.e. whether no water erroneously appears or disappears. It looks at the storage rate \[ \text{storage rate} = \frac{\Delta S_b}{\Delta t} \]
in a Basin over a time period \(\Delta t\) and compares that to the total inflows and outflows of that Basin over that period. More precisely, we first compute the total inflow and outflow, where:
- \(\text{total inflow}\): the precipitation, drainage and horizontal flows into the Basin;
- \(\text{total outflow}\): the evaporation, infiltration and horizontal flows out of the Basin.
Whether a flow is an inflow or an outflow depends on whether the flow contributes to or takes from the Basin storage, which means that this is independent of the edge direction. This is determined for each solver timestep individually.
Then from this we compute the errors:
\[ \begin{align} \text{balance error} =&& \text{storage rate} - (\text{total inflow} - \text{total outflow}) \\ \text{relative error}=&& \frac{\text{absolute error}}{0.5(\text{total inflow} + \text{total outflow})} \end{align} \] Hence the reference used for computing the relative error is the average of the total inflow and total outflow of the Basin (which are both non-negative).
The default tolerances are \(0.001 \text{ m}^3\) for the balance error and \(0.01\) for the relative error, which should not be exceeded for realistic models.
In extreme cases where the storage rate is many orders of magnitude smaller than the storage itself, these computations can have floating point truncation errors which can lead to large relative errors. This is however only when the storage is roughly \(\geq 10^{15}\) times bigger than the storage rate.
1.4.1 Example calculation
Say we have the following model:
and we want to calculate the water balance error for Basin 6. We have the following data:
- Time period length: \(10.0 \text{ s}\)
- Basin storage start: \(100.0 \text{ m}^3\)
- Basin storage end: \(50.0 \text{ m}^3\)
- UserDemand #11 inflow average: \(10.0 \text{ m}^3/\text{s}\)
- UserDemand #11 outflow average: \(5.0 \text{ m}^3/\text{s}\)
- Outlet #7 flow average: \(- 3.5 \text{ m}^3/\text{s}\)
- Outlet #11 flow average: \(4.0 \text{ m}^3/\text{s}\)
And so we get
\[ \begin{align} \text{storage rate} = && \frac{50.0 - 100.0}{10.0} &= & -6.0 \text{ m}^3/\text{s} \\ \text{total inflow} = && 5.0 + 3.5 &= & 8.5 \text{ m}^3/\text{s}\\ \text{total outflow} = && 10.0 + 4.0 &= & 14.0 \text{ m}^3/\text{s}\\ \text{balance error} = && -6.0 - (8.5 - 14.0) &= & -0.5 \text{ m}^3/\text{s}\\ \text{relative error} = && \frac{-0.5}{8.5 + 14.0} &\approx & -0.022 \end{align} \] Note that the balance error and relative error are negative, but we use their absolute value to compare to the respective tolerances.
1.5 Why this formulation
You might wonder why in the above explanation the states are given by the cumulative flows and not by the Basin storages, which is arguably conceptually simpler. The reason is that we do not just want to model the storages in the Basins over time, but we also want accurate output of each individual flow, e.g. to model the spread of pollutants.
When the states are given by the storages, generally the individual flows can not accurately be computed from that as a post processing step, because there are more flows than storages. Also, we can only compute flows at individual points in time explicitly, not over a whole interval. When the states are given by the cumulative flows however, the output of the problem solve gives these flows directly, and from those the storage over time can be computed accurately. Hence in short, the formulation above gives more information than a formulation with Basin storages as states.
1.6 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.
2 Performance
There are many things that can influence the calculations times, for instance:
- Solver tolerance: By default both the absolute and relative tolerance is
1e-7. - 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.