Equations

Ribasim currently simulates the following “natural” water balance terms:

Precipitation
Evaporation
Infiltration
Drainage
Urban runoff
Upstream and downstream flow

Additionally, Ribasim simulates the following “allocated” water balance terms:

UserDemand
Flushing

Depending on the type of boundary conditions, Ribasim requires relation between storage volume and wetted area \(A\), and between the storage volume and the water level \(h\). These are (currently) represented by piecewise linear relationships.

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}], \tag{1}\]

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 3 PID controller.

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.

Note

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)}. \]

Note

The presence of division by the basin area means that areas of size zero are not allowed.

2 Natural water balance terms

2.1 The reduction factor

At several points in the equations below a reduction factor is used. This is a term that makes certain transitions more smooth, for instance when a pump stops providing water when its source basin dries up. The reduction factor is given by

\[\begin{align} \phi(x; p) = \begin{cases} 0 &\text{if}\quad x < 0 \\ -2 \left(\frac{x}{p}\right)^3 + 3\left(\frac{x}{p}\right)^2 &\text{if}\quad 0 \le x \le p \\ 1 &\text{if}\quad x > p \end{cases} \end{align}\]

Here \(p > 0\) is the threshold value which determines the interval \([0,p]\) of the smooth transition between \(0\) and \(1\), see the plot below.

Code

import numpy as np
import matplotlib.pyplot as plt

def f(x, p = 3):
    x_scaled = x / p
    phi = (-2 * x_scaled + 3) * x_scaled**2
    phi = np.where(x < 0, 0, phi)
    phi = np.where(x > p, 1, phi)

    return phi

fontsize = 15
p = 3
N = 100
x_min = -1
x_max = 4
x = np.linspace(x_min,x_max,N)
phi = f(x,p)

fig,ax = plt.subplots(dpi=80)
ax.plot(x,phi)

y_lim = ax.get_ylim()

ax.set_xticks([0,p], [0,"$p$"], fontsize=fontsize)
ax.set_yticks([0,1], [0,1], fontsize=fontsize)
ax.hlines([0,1],x_min,x_max, color = "k", ls = ":", zorder=-1)
ax.vlines([0,p], *y_lim, color = "k", ls = ":")
ax.set_xlim(x_min,x_max)
ax.set_xlabel("$x$", fontsize=fontsize)
ax.set_ylabel(r"$\phi(x;p)$", fontsize=fontsize)
ax.set_ylim(y_lim)

fig.tight_layout()
plt.show()

2.2 Precipitation

The precipitation term is given by

\[ Q_P = P \cdot A. \tag{2}\]

Here \(P = P(t)\) is the precipitation rate and \(A\) is the maximum area given in the Basin / profile table. Precipitation in the Basin area is assumed to be directly added to the Basin storage. The modeler needs to ensure all precipitation enters the model, and there is no overlap in the maximum profile areas, else extra water is created. If a part of the catchment is not in any Basin profile, the modeler has to verify that water source is not forgotten. It can for instance be converted to a flow rate and added to a Basin as a FlowBoundary.

2.3 Evaporation

The evaporation term is given by

\[ Q_E = E_\text{pot} \cdot A(u) \cdot \phi(d;0.1). \tag{3}\]

Here \(E_\text{pot} = E_\text{pot}(t)\) is the potential evaporation rate and \(A\) is the wetted area. \(\phi\) is the reduction factor which depends on the depth \(d\). It provides a smooth gradient as \(u \rightarrow 0\).

A straightforward formulation \(Q_E = \mathrm{max}(E_\text{pot} A(u), 0)\) is unsuitable, as \(\frac{\mathrm{d}Q_E}{\mathrm{d}u}(u=0)\) is then not well-defined.

A non-smooth derivative results in extremely small timesteps and long computation time: ModelingToolkit identifies the singular behavior and adjusts its timestepping. In a physical interpretation, evaporation is switched on or off per individual droplet of water. In general, the effect of the reduction term is negligible, or not even necessary. As a surface water dries, its wetted area decreases and so does the evaporative flux. However, for (simplified) cases with constant wetted surface (a rectangular profile), evaporation only stops at \(u = 0\).

2.4 Infiltration and Drainage

Infiltration is provided as a lump sum for the basin. If Ribasim is coupled with MODFLOW 6, the infiltration is computed as the sum of all positive flows of the MODFLOW 6 boundary conditions in the basin:

\[ Q_\text{inf} = \sum_{i=1}^{n} \sum_{j=1}^{m} \max(Q_{\mathrm{mf6}_{i,j}}, 0.0) \tag{4}\]

Where \(i\) is the index of the boundary condition, \(j\) the MODFLOW 6 cell index, \(n\) the number of boundary conditions, and \(m\) the number of MODFLOW 6 cells in the basin. \(Q_{\mathrm{mf6}_{i,j}}\) is the flow computed by MODFLOW 6 for cell \(j\) for boundary condition \(i\).

Drainage is a lump sump for the basin, and consists of the sum of the absolute value of all negative flows of the MODFLOW 6 boundary conditions in the basin.

\[ Q_\text{drn} = \sum_{i=1}^{n} \sum_{j=1}^{m} \left| \min(Q_{\mathrm{mf6}_{i,j}}, 0.0) \right| \tag{5}\]

The interaction with MODFLOW 6 boundary conditions is explained in greater detail in the the iMOD Coupler docs.

2.5 Upstream and downstream flow

Ribasim’s basins can be connected to each other, and each basin expects an explicit connection. These connections are currently available for inter-basin flows:

Pump
TabulatedRatingCurve
LinearResistance
ManningResistance

The flow direction of the basin is not pre-determined: flow directions may freely reverse, provided the connection allows it. Currently, a LinearResistance allows bidirectional flow, but the

Additionally, three additional “connections” area available for the “outmost” basins (external nodes) in a network.

Terminal
LevelBoundary
FlowBoundary

2.5.1 Pump

The behaviour of pumps is very straight forward if these nodes are not PID controlled. Their flow is given by a fixed flow rate \(q\), multiplied by a reduction factor: \[ Q_\text{pump} = \phi(u; 10.0)q \]

Here \(u\) is the storage of the upstream basin. The reduction factor \(\phi\) makes sure that the flow of the pump goes smootly to \(0\) as the upstream basin dries out.

2.5.2 Outlet

The outlet is very similar to the pump, but it has a few extra reduction factors for physical constraints: \[ Q_\text{outlet} = \phi(u_a; 10.0)\phi(\Delta h; 0.1) \phi(h_a-h_\text{min};0.1)q. \] The subscript \(a\) denotes the upstream node and \(b\) the downstream node. The first reduction factor is equivalent to the one for the pump. The second one makes sure that the outlet flow goes to zero as the head difference \(\Delta h = h_a - h_b\) goes to zero. The last one makes sure that the outlet only produces flow when the upstream level is above the minimum chrest level \(h_\text{min}\).

Not all node types upstream or downstream of the outlet have a defined level. If this is the case, and therefore the reduction factor cannot be computed, it is defined to be \(1.0\).

2.5.3 TabulatedRatingCurve

The Tabulated Rating Curve is a tabulation of a basin’s discharge behavior. It describes a piecewise linear relationship between the basin’s level and its discharge. It can be understood as an empirical description of a basin’s properties. This can include an outlet, but also the lumped hydraulic behavior of the upstream channels.

Note

Currently, the discharge relies only on the basin’s level; it could also use the volume of both connected basins to simulate backwater effects, submersion of outlets, or even reversal of flows for high precipitation events.

2.5.4 LinearResistance

A LinearResistance connects two basins together. The flow between the two basins is determined by a linear relationship, up to an optional maximum flow rate:

\[ Q = \mathrm{clamp}(\frac{h_a - h_b}{R}, -Q_{\max}, Q_{\max}) \tag{6}\]

Here \(h_a\) is the water level in the first basin and \(h_b\) is the water level in the second basin. \(R\) is the resistance of the link, and \(Q_{\max}\) is the maximum flow rate. A LinearResistance makes no assumptions about the direction of the flow: water flows from high to low.

2.5.5 Terminal

This only allows outflow from a basin into a terminal node.

2.5.6 LevelBoundary

This can be connected to a basin via a LinearResistance. This boundary node will then exchange water with the basin based on the difference in water level between the two.

2.5.7 FlowBoundary

This can be connected directly to a basin and prescribes the flow to or from that basin. We require that the edge connecting the flow boundary to the basin should point towards the basin, so that positive flow corresponds to water being added to the model.

2.5.8 Manning connection

Ribasim is capable of simulating steady flow between basins through a reach described by a trapezoidal profile and a Manning roughness coefficient.

We describe the discharge from basin \(a\) to basin \(b\) solely as a function of the water levels in \(a\) and \(b\).

\[ Q = f(h_a, h_b) \]

where:

The subscripts \(a,b\) denote basins
\(h\) is the hydraulic head, or water level

The energy equation for open channel flow is:

\[ H = h + \frac{v^2}{2g} \]

Where

\(H\) is total head
\(v\) is average water velocity
\(g\) is gravitational acceleration

The discharge \(Q\) is defined as:

\[ Q = Av \]

where \(A\) is cross-sectional area.

We use conservation of energy to relate the total head at \(a\) to \(b\), with \(H_a > H_b\) as follows:

\[ H_a = H_b + h_{\text{loss}} \]

Or:

\[ h_a + \frac{v_a^2}{2g} = h_b + \frac{v_b^2}{2g} + h_{\text{loss}} \]

Where \(v\) is the average water velocity. \(h_{\text{loss}}\) is a combination of friction and contraction/expansion losses:

\[ h_{\text{loss}} = S_f L + \frac{C}{2g} \left(v_b^2 - v_a^2\right) \]

Where:

\(L\) is the reach length
\(S_f\) is the representative friction slope
\(C\) is the expansion or contraction coefficient, \(0 \le C \le1\)

We assume velocity differences in a connection are negligible (\(v_a = v_b\)):

\[ h_a = h_b + S_f L \]

Friction losses are computed with the Gauckler-Manning formula:

\[ Q = \frac{A}{n} R_h^\frac{2}{3} \sqrt{S_f} \]

Where:

\(A\) is the representative area.
\(R_h\) is the representative wetted radius.
\(S_f\) is the representative friction slope.
\(n\) is Manning’s roughness coefficient.

We can rewrite to express \(S_f\) in terms of Q:

\[ S_f = Q^2 \frac{n^2}{A^2 R_h^{4/3}} \]

No water is added or removed in a connection:

\[ Q_a = Q_b = Q \]

Substituting:

\[ h_a = h_b + Q^2 \frac{n^2}{A^2 R_h^{4/3}} L \]

We can then express \(Q\) as a function of head difference \(\Delta h\):

\[ Q = \textrm{sign}(\Delta h) \frac{A}{n} R_h^{2/3}\sqrt{\frac{|\Delta h|}{L} } \]

The \(\textrm{sign}(\Delta h)\) term causes the direction of the flow to reverse if the head in basin \(b\) is larger than in basin \(a\).

This expression however leads to problems in simulation since the derivative of \(Q\) with respect to \(\Delta h\) tends to \(\pm \infty\) as \(\Delta h\) tends to 0. Therefore we use the slightly modified expression

\[ Q = \textrm{sign}(\Delta h) \frac{A}{n} R_h^{2/3}\sqrt{\frac{\Delta h}{L} s(\Delta h)} \]

to smooth out this problem. Here \(s(x) = \frac{2}{\pi}\arctan{1000x}\) can be thought of as a smooth approximation of the sign function.

Note

The computation of \(S_f\) is not exact: we base it on a representative area and hydraulic radius, rather than integrating \(S_f\) along the length of a reach. Direct analytic solutions exist for e.g. parabolic profiles (Tolkmitt), but other profiles requires relatively complicated approaches (such as approximating the profile with a polynomial).

We use the average value of the cross-sectional area, the average value of the water depth, and the average value of the hydraulic radius to compute a friction slope. The size of the resulting error will depend on the water depth difference between the upstream and downstream basin.

The cross sectional area for a trapezoidal or rectangular profile:

\[ A = w d + \frac{\Delta y}{\Delta z} d^2 \]

Where

\(w\) is the width at \(d = 0\) (A triangular profile has \(w = 0\))
\(\frac{\Delta y}{\Delta z}\) is the slope of the profile expressed as the horizontal length for one unit in the vertical (A slope of 45 degrees has \(\frac{\Delta y}{\Delta z} = 1\); a rectangular profile 0).

Accordingly, the wetted perimeter is:

\[ B = w + 2 d \sqrt{\left(\frac{\Delta y}{\Delta z}\right)^2 + 1} \]

3 UserDemand allocation

UserDemands have an allocated flow rate \(F^p\) per priority \(p=1,2,\ldots, p_{\max}\), which is either determined by allocation optimization or simply equal to the demand at time \(t\); \(F^p = d^p(t)\). The actual abstraction rate of a UserDemand is given by \[ Q_\text{userdemand, in} = \phi(u, 10.0)\phi(h-l_{\min}, 0.1)\sum_{p=1}^{p_{\max}} \min\left(F^p, d^p(t)\right). \]

From left to right:

The first reduction factor lets the UserDemand abstraction go smoothly to \(0\) as the source Basin dries out;
The second reduction factor lets the UserDemand abstraction go smoothly to \(0\) as the source Basin level approaches the minimum source Basin level (from above);
The last term is the sum of the allocations over the priorities. If the current demand happens to be lower than the allocation at some priority, the demand is taken instead.

UserDemands also have a return factor \(0 \le r \le 1\), which determines the return flow (outflow) of the UserDemand:

\[ Q_\text{userdemand, out} = r \cdot Q_\text{userdemand, in}. \]

Note that this means that the user_demand has a consumption rate of \((1-r)Q_\text{userdemand, in}\).

4 PID controller

The PID controller continuously sets the flow rate of a pump or outlet to bring the level of a certain basin closer to its setpoint. If we denote the setpoint by \(\text{SP}(t)\) and the basin level by \(y(t)\), then the error is given by \[ e(t) = \text{SP}(t) - y(t). \tag{7}\]

The output of the PID controller for the flow rate of the pump or outlet is then given by \[ Q_\text{PID}(t) = K_p e(t) + K_i\int_{t_0}^t e(\tau)\text{d}\tau + K_d \frac{\text{d}e}{\text{d}t}, \tag{8}\]

for given constant parameters \(K_p,K_i,K_d\). The pump or outlet can have associated minimum and maximum flow rates \(Q_{\min}, Q_{\max}\), and so \[ Q_\text{pump/outlet} = \text{clip}(\Phi Q_\text{PID}; Q_{\min}, Q_{\max}). \]

Here \(u_\text{us}\) is the storage of the basin upstream of the pump or outlet, \(\Phi\) is the product of reduction factors associated with the pump or outlet and

\[\begin{align} \text{clip}(Q; Q_{\min}, Q_{\max}) = \begin{cases} Q_{\min} & \text{if} \quad Q < Q_{\min} \\ Q & \text{if} \quad Q_{\min} \leq Q \leq Q_{\max} \\ Q_{\max} & \text{if} \quad Q > Q_{\max} \end{cases}. \end{align}\]

For the integral term we denote \[ I(t) = \int_{t_0}^t e(\tau)\text{d}\tau, \]

where \(t_0\) is the last time the PID controller was made active. \(I(t)\) is treated as a state of the system and thus it has its own equation in the system in Equation 1: \[ \frac{\text{d}I}{\text{d}t} = e(t). \]

Note

In the case of the controlled outlet, the upstream node can also be a level boundary. In this case we define \(\phi = 1\).

4.1 The derivative term

When \(K_d \ne 0\) this adds a level of complexity. We can see this by looking at the error derivative more closely: \[ \frac{\text{d}e}{\text{d}t} = \frac{\text{d}\text{SP}}{\text{d}t} - \frac{1}{A(u_\text{PID})}\frac{\text{d}u_\text{PID}}{\text{d}t}, \] where \(A(u_\text{PID})\) is the area of the controlled basin as a function of the storage of the controlled basin \(u_\text{PID}\). The complexity arises from the fact that \(Q_\text{PID}\) is a contribution to \(\frac{\text{d}u_\text{PID}}{\text{d}t} = f_\text{PID}\), which makes Equation 8 an implicit equation for \(Q_\text{PID}\). We define

\[ f_\text{PID} = \hat{f}_\text{PID} \pm Q_\text{pump/outlet}, \]

that is, \(\hat{f}_\text{PID}\) is the right hand side of the ODE for the controlled basin storage state without the contribution of the PID controlled pump. The plus sign holds for an outlet and the minus sign for a pump, dictated by the way the pump and outlet connectivity to the controlled basin is enforced.

Using this, solving Equation 8 for \(Q_\text{PID}\) yields \[ Q_\text{pump/outlet} = \text{clip}\left(\phi(u_\text{us})\frac{K_pe + K_iI + K_d \left(\frac{\text{d}\text{SP}}{\text{d}t}-\frac{\hat{f}_\text{PID}}{A(u_\text{PID})}\right)}{1\pm\phi(u_\text{us})\frac{K_d}{A(u_\text{PID})}};Q_{\min},Q_{\max}\right), \] where the clipping is again done last. Note that to compute this, \(\hat{f}_\text{PID}\) has to be known first, meaning that the PID controlled pump/outlet flow rate has to be computed after all other contributions to the PID controlled basin’s storage are known.

4.2 The sign of the parameters

Note by Equation 7 that the error is positive if the setpoint is larger than the basin level and negative if the setpoint is smaller than the basin level.

We enforce the convention that when a pump is controlled, its edge points away from the basin, and when an outlet is controlled, its edge points towards the basin, so that the main flow direction along these edges is positive. Therefore, positive flows of the pump and outlet have opposite effects on the basin, and thus the parameters \(K_p,K_i,K_d\) of the pump and outlet must have oppositive signs to achieve the same goal.

5 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.

6 Performance

Ribasim needs to be sufficiently fast to make application on the national scale, with ~10000 basins, practical. Therefore, whilst developing, we need to be mindful that our approach can scale to such a size. We currently simulate a set of 40 connected free-flowing basins for a period of 2 years in about 10 seconds.

For a real scaling test we would need to do a national simulation, but we expect these computation times not to be problematic, considering that simulating those same 40 basins in MODFLOW 6 takes minutes.

There are many things that can influence the calculations times, for instance:

Solver tolerance we currently use a very conservative tolerance of 1e-10, whereas with larger tolerances step sizes can easily be 3x longer, leading to a 3x speedup.
ODE solvers: The Rosenbrock23 method 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.

When scaling up the model, we will need spend some time to strike the right balance between error tolerance, schematization and forcing. The SciML software we use has many showcases of large scale applications, as well as documentation on how to achieve the optimal performance for your system.

References

SciML Development Team. 2022. “Event Handling and Callback Functions.” https://docs.sciml.ai/DiffEqDocs/latest/features/callback_functions/.