Skip to main content

Numerical considerations

We want to solve the following initial value problem: \[ \begin{cases} \frac{\text{d}\mathbf{u}}{\text{d}t} = \mathbf{f}(\mathbf{u},t) \quad t_0 < t < t_\text{end} \\ \mathbf{u}(t_0) = \mathbf{u}_0 \end{cases}, \tag{1}\]

where \(\mathbf{f}\) denotes water_balance! and \(\mathbf{u_0} = \mathbf{0}\) the initial cumulative flows (and the PID integrals which also start out at \(0\)).

In general \(\mathbf{f}\) is a non-linear function in \(\mathbf{u}\). These non-linearities are introduced by e.g.:

The problem Equation 1 can be solved by various numerical time-integration methods. To do this the time interval \([t_0,t_\text{end}]\) is discretized into a finite number of time points \(t_0 < t_1 < \ldots < t_N = t_\text{end}\) for which approximate solutions \(\mathbf{w}_n \approx \mathbf{u}(t_n)\) are computed. In general we do not assume a fixed timestep (the interval between successive points in time). Rather, the solver attempts to make as large a step as possible while keeping error tolerances within requirements. The solver settings section details the available configuration options.

1 Example numerical methods

This section discusses two relatively simple numerical methods, mainly to demonstrate the difference between explicit and implicit methods.

1.1 Euler forward

The simplest numerical method is Euler forward: \[ \mathbf{w}_{n+1} = \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{w}_n, t_n). \tag{2}\]

Here \(\mathbf{w}_{n+1}\) is given as a simple explicit function of \(\mathbf{w}_n\).

1.2 Euler backward

Euler backward is formulated as follows: \[ \mathbf{w}_{n+1} = \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{w}_{n+1},t_{n+1}). \tag{3}\]

Note that this is an implicit equation for \(\mathbf{w}_{n+1}\), which is non-linear because of the non-linearity of \(\mathbf{f}\).

Generally one of the following iterative methods is used for finding solutions to non-linear equations like this:

  • Picard iteration for fixed points. This method aims to approximate \(\mathbf{w}_{n+1}\) as a fixed point of the function \[ \mathbf{g}(\mathbf{x}) = \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{x},t_{n+1}) \] by iterating \(\mathbf{g}\) on an initial guess of \(\mathbf{w}_{n+1}\);
  • Newton iterations: approximate \(\mathbf{w}_{n+1}\) as a root of the function \[ \mathbf{h}(\mathbf{x}) = \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{x},t_{n+1}) - \mathbf{x}, \] by iteratively finding the root of its linearized form:

\[\begin{align} \mathbf{0} =& \mathbf{h}(\mathbf{w}_{n+1}^k) + \mathbf{J}(\mathbf{h})(\mathbf{w}_{n+1}^k)(\mathbf{w}_{n+1}^{k+1}-\mathbf{w}_{n+1}^k) \\ =& \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{w}_{n+1}^k,t_{n+1}) - \mathbf{w}_{n+1}^k \\ +&\left[(t_{n+1}-t_n)\mathbf{J}(\mathbf{f})(\mathbf{w}_{n+1}^k)-\mathbf{I}\right](\mathbf{w}_{n+1}^{k+1}-\mathbf{w}_{n+1}^k). \end{align}\] Note that this thus requires an evaluation of the Jacobian of \(\mathbf{f}\) and solving a linear system per iteration.

2 The advantage of implicit methods

The implicit method Equation 3 is a coupled system of equations for \(\mathbf{w}_{n+1}\), while the explicit Equation 2 is fully decoupled. This means in general that in a time integration step with an implicit method the basins communicate information with eachother, while in an explicit method this does not happen. A consequence of this is that local events (e.g. a pump turns on) propagate slowly through the model using an explicit method but quickly using an implicit method, making implicit methods more stable.

3 The Jacobian

The Newton iterations above require the Jacobian of \(\mathbf{f}\), an \(N \times N\) matrix where \(N\) is the number of states. Its entries are given by \[ J_{i,j} = \frac{\partial f_i}{\partial u_j}, \] i.e. \(J_{i,j}\) quantifies how \(f_i\), the time derivative of state \(i\), changes with respect to changes in state \(j\). Most of these entries are \(0\), because flows in distant parts of the model do not depend on each other.

There are several ways to compute the Jacobian:

  • Using finite difference methods; however these are relatively expensive to compute and introduce truncation errors.
  • Hardcoding analytical expressions for the partial derivatives; these have the advantages over finite difference methods that they are cheaper to compute and don’t introduce truncation errors. However, computing and coding these partial derivatives is complex and thus error prone, and requires maintenance for every adjustment of \(\mathbf{f}\).
  • Using automatic differentiation; this has the advantages of being fast to compute and not requiring much extra implementation other than implementing \(\mathbf{f}\) in a way that is compatible with an automatic differentiation package, in the case of Ribasim that is ForwardDiff.jl.

Ribasim uses automatic differentiation by default. With forward mode automatic differentiation, one evaluation of \(\mathbf{f}\) yields one column of the Jacobian, or more precisely the Jacobian applied to one chosen direction. A naive computation of the Jacobian thus costs \(N\) evaluations of \(\mathbf{f}\), which is prohibitively expensive for large models.

3.1 Sparse matrix coloring

If two columns of the Jacobian have no nonzero entries in the same row, they can be computed in a single evaluation of \(\mathbf{f}\) by perturbing both corresponding states at once: the contributions can be told apart afterwards based on the sparsity pattern. The columns of a sparse Jacobian can be partitioned into groups (colors) of such structurally orthogonal columns, and only one evaluation per color is needed. Ribasim detects the sparsity pattern automatically with SparseConnectivityTracer.jl and computes the coloring with SparseMatrixColorings.jl.

The number of colors depends on the structure of the Jacobian. For the Jacobian of \(\mathbf{f}\) this structure is unfavorable: a perturbation of a cumulative flow state changes the storage of the Basins it connects, and that affects all flows connected to these Basins. Hence all flow states connected to the same Basin end up in different colors, and the number of colors is at least the number of connections of the most connected Basin. Therefore Ribasim does not compute the Jacobian of \(\mathbf{f}\) directly, but the Jacobian of an intermediate function which can be colored with very few colors, as explained below.

3.2 Structure of the right hand side

We split the state vector into the cumulative flows \(\mathbf{u}^\text{flow}\) and the PID integrals \(\mathbf{u}^\text{pid}\). The Basin storages depend linearly on the cumulative flows (see the equations): \[ \mathbf{S} = \mathbf{S}_0(t) + M\mathbf{u}^\text{flow}, \] where \(M\) is the \(n_\text{basin} \times n_\text{flow}\) signed incidence matrix: \(M_{b,k} = -1\) if Basin \(b\) is upstream of flow \(k\), \(M_{b,k} = 1\) if Basin \(b\) is downstream of flow \(k\) and \(M_{b,k} = 0\) otherwise. Evaporation and infiltration are flows with their Basin upstream and nothing downstream. The flows over the links out of the Junction nodes in Junction / static are not states: they are weighted sums of the flows that enter the subgraph of Junction nodes, with the products of the fractions along the way as weights, minus the abstractions of the UserDemand nodes in the subgraph. The column of \(M\) of a flow into such a Junction therefore holds these weights for the Basins downstream of the subgraph instead of a \(1\) for a Basin directly downstream, and the same holds for the UserDemand abstractions with negative weights. See Junction below.

For every flow \(k\) we then take a copy of the storage of its upstream Basin \(S^\text{up}_k\) and its downstream Basin \(S^\text{down}_k\): \[ \mathbf{S}^\text{up} = P^\text{up}\mathbf{S}, \quad \mathbf{S}^\text{down} = P^\text{down}\mathbf{S}, \] where \(P^\text{up}\) and \(P^\text{down}\) are constant selection matrices (note that \(M = (P^\text{down} - P^\text{up})^T\)). The flows are computed from these storage copies, the PID integrals \(\mathbf{u}^\text{pid}\) and the compound variables \(\mathbf{c}\) of the ContinuousControl nodes: \[ \mathbf{f}^\text{flow} = \mathbf{q}(\mathbf{v}), \quad \mathbf{v} = \left(\mathbf{S}^\text{up}, \mathbf{S}^\text{down}, \mathbf{u}^\text{pid}, \mathbf{c}\right). \]

The Jacobian \(\partial\mathbf{q}/\partial\mathbf{v}\) is what Ribasim computes with automatic differentiation. It is not the Jacobian of \(\mathbf{f}\), but it is very cheap to compute:

  • \(\partial\mathbf{q}/\partial\mathbf{S}^\text{up}\) and \(\partial\mathbf{q}/\partial\mathbf{S}^\text{down}\) are diagonal, since flow \(k\) only depends on its own storage copies \(S^\text{up}_k\) and \(S^\text{down}_k\). The exception are flows controlled by a PidControl node with a derivative term, which depend on all flows of the listened Basin;
  • \(\partial\mathbf{q}/\partial\mathbf{u}^\text{pid}\) and \(\partial\mathbf{q}/\partial\mathbf{c}\) have at most one nonzero per row, for the flows that are PID controlled or continuously controlled respectively.

These blocks can be colored with only a handful of colors, independent of the size of the model. The Jacobian of \(\mathbf{f}\) follows from the chain rule, where only the multiplications with the constant matrices \(P^\text{up}M\) and \(P^\text{down}M\) remain: \[ \frac{\partial \mathbf{f}^\text{flow}}{\partial \mathbf{u}^\text{flow}} = Q_S M, \quad Q_S = \frac{\partial\mathbf{q}}{\partial\mathbf{S}^\text{up}}P^\text{up} + \frac{\partial\mathbf{q}}{\partial\mathbf{S}^\text{down}}P^\text{down}, \] (ignoring control for now), where \(Q_S = \partial\mathbf{q}/\partial\mathbf{S}\) is the \(n_\text{flow} \times n_\text{basin}\) matrix with the dependence of the flows on the Basin storages. This matrix is never assembled; instead the structure is exploited directly in the linear solve.

3.3 The linear solve

Every Newton iteration of an implicit method requires solving a linear system of the form \[ \left(J - \gamma^{-1}I\right)\mathbf{x} = \mathbf{b}, \] where \(\gamma\) is proportional to the timestep (compare the Euler backward iteration above, where \(\gamma = t_{n+1} - t_n\)). This system has a dimension equal to the number of states, which is much larger than the number of Basins. Using \(J = Q_S M\) and defining \(\mathbf{s} = M\mathbf{x}\), the system reads \[ -\gamma^{-1}\mathbf{x} + Q_S\mathbf{s} = \mathbf{b} \quad \Rightarrow \quad \mathbf{x} = \gamma\left(Q_S\mathbf{s} - \mathbf{b}\right). \tag{4}\] Multiplying the left equation by \(M\) gives a system in terms of \(\mathbf{s}\) only: \[ \left(J_\text{inner} - \gamma^{-1}I_{n_\text{basin}}\right)\mathbf{s} = M\mathbf{b}, \quad J_\text{inner} = M Q_S. \tag{5}\] Hence the linear solve is done in two steps:

  1. Solve the \(n_\text{basin} \times n_\text{basin}\) system Equation 5 for the storage changes \(\mathbf{s}\). The sparsity pattern of \(J_\text{inner}\) is that of the adjacency matrix of the Basins, which is very sparse; by default this system is solved with a sparse LU factorization.
  2. Compute the solution in terms of the states \(\mathbf{x}\) with Equation 4, which is just a cheap matrix-vector product.

3.3.1 PID control

The PID integral state of a PidControl node \(j\) evolves with the error \(f^\text{pid}_j = e_j = h^\text{target}_j - h(S_{l_j})\), where \(l_j\) is the listened Basin, so \[ \frac{\partial f^\text{pid}_j}{\partial \mathbf{u}^\text{flow}} = -\frac{1}{A_{l_j}}M_{l_j,:}, \] where \(A_{l_j}\) is the area of the listened Basin. The listened Basin is directly up- or downstream of the controlled Pump or Outlet, so the dependence of the controlled flow on the listened Basin storage through the proportional term is already included in \(Q_S\). Let \(G = \partial\mathbf{q}/\partial\mathbf{u}^\text{pid}\) and let \(E\) be the matrix that selects the listened Basin per PidControl node and divides by its area. The PID rows of the linear system then give \[ \mathbf{x}^\text{pid} = -\gamma\left(\mathbf{b}^\text{pid} + E\mathbf{s}\right), \] and substituting this into the flow rows leads to \[ \mathbf{x}^\text{flow} = \gamma\left(Q_S\mathbf{s} + G\mathbf{x}^\text{pid} - \mathbf{b}^\text{flow}\right), \] with the inner system \[ \left(M Q_S - \gamma M G E - \gamma^{-1}I_{n_\text{basin}}\right)\mathbf{s} = M\mathbf{b}^\text{flow} + \gamma M G \mathbf{b}^\text{pid}. \] So the PID integral states are eliminated in the same way as the flow states, and the dimension of the inner system remains the number of Basins.

3.3.2 Continuous control

The compound variable \(c_m\) of a ContinuousControl node depends on Basin storages and on (listened) flows. Their derivatives \(\partial\mathbf{c}/\partial\mathbf{S}\) and \(\partial\mathbf{c}/\partial\mathbf{q}\) are also computed with automatic differentiation. For a continuously controlled flow this adds the term \[ \frac{\partial\mathbf{q}}{\partial\mathbf{c}}\left(\frac{\partial\mathbf{c}}{\partial\mathbf{S}} + \frac{\partial\mathbf{c}}{\partial\mathbf{q}}Q_S\right) \] to \(Q_S\), and thus \(M\) times this term to \(J_\text{inner}\). Here only the direct dependence of listened flows on their own Basin storages is taken into account, since chaining control (listening to a flow that is itself controlled by PidControl or ContinuousControl) is not supported.

3.3.3 Junction

Since the fractions of the Junction links are constant, the flow over every link out of a Junction in Junction / static is a constant weighted sum of the flows that enter the subgraph and of the inflows of the UserDemand nodes in it, and so is the contribution of these flows to the Basins downstream of the subgraph. These weights are in \(M\), so the links in the subgraphs need no states of their own, and a FlowBoundary that flows into such a Junction is integrated exactly like any other FlowBoundary. The only flows that need a computation are the inflows of the UserDemand nodes that abstract from a Junction, the minimum of the total inflow of that node and the demand, and their return flows. They are computed in the last step of the right-hand side as a function of the other flows \(\mathbf{q}^\text{np}\) (np: not computed in a Junction subgraph), which are computed as described above: \[ \mathbf{q} = K\mathbf{q}^\text{np}, \quad K = \frac{\partial\mathbf{q}}{\partial\mathbf{q}^\text{np}}. \] \(K\) is the identity except for the rows of these UserDemand flows, which hold the dependence on the flows entering the subgraph. It is computed with automatic differentiation like the other Jacobians, and only if the model has Junction nodes in Junction / static. For the derivative term of PID control, which uses the storage rate of the listened Basin, these UserDemand flows are computed once more before the PID controlled flows, from the flows known at that point.

Since \(K\) is applied after \(\mathbf{q}^\text{np}\), the flow rows of the linear system have \(Q_S\) and \(G\) replaced by \(KQ_S\) and \(KG\), and the solution of the flow states becomes \[ \mathbf{x}^\text{flow} = \gamma\left(K\left(Q_S\mathbf{s} + G\mathbf{x}^\text{pid}\right) - \mathbf{b}^\text{flow}\right). \] The inner system follows as before by multiplying by \(M\), which holds the routing weights: \[ \left(M K Q_S - \gamma M K G E - \gamma^{-1}I_{n_\text{basin}}\right)\mathbf{s} = M\mathbf{b}^\text{flow} + \gamma M K G \mathbf{b}^\text{pid}. \] For instance, a Pump upstream of a Basin \(b_1\) which flows into a Junction that divides its flow over Basins \(b_2\) and \(b_3\) with link fractions \(f\) and \(1-f\) has the column \(-e_{b_1} + f e_{b_2} + (1-f) e_{b_3}\) in \(M\). If a UserDemand abstracts from the Junction, its inflow has the column \(-f e_{b_2} - (1-f) e_{b_3}\), and the row of \(K\) of that inflow holds its dependence on the Pump flow. Note that \(\mathbf{b}^\text{flow}\) is not multiplied by \(K\), since the right-hand side of the linear system is given in terms of the states. The dimension of the inner system remains the number of Basins.

3.3.4 Exactness

The linear solve above is equivalent to solving the full linear system with the Jacobian of \(\mathbf{f}\), as long as control is not chained. If it is, the Jacobian used in the linear solve is approximate. An approximate Jacobian only affects the convergence speed of the Newton iterations, not the accuracy of the converged solution, since the latter is determined by \(\mathbf{f}\) itself.

The reduced linear solve can be disabled with the solver setting reduced_implicit_solve = false, see the solver settings. In that case the (sparse) Jacobian of \(\mathbf{f}\) is computed directly with automatic differentiation and sparse matrix coloring, and the full linear system is solved. This is mainly useful for verification, as it is considerably slower for larger models.

4 Convergence of the Newton iterations

The solver reuses the Jacobian and the factorization of \(W = J - \gamma^{-1}I\) over many timesteps, so the Newton iterations converge linearly, at a rate \(\theta\) that is estimated from successive increments. A solve is accepted when the error estimated from \(\theta\) and the last increment is small compared to the tolerance. OrdinaryDiffEq uses the criterion from Section IV.8 of Hairer and Wanner; the CVODE documentation describes the same idea. The first iteration of a solve uses the \(\theta\) of an earlier solve, and a solve that stops after one iteration does not measure it again, so a small \(\theta\) could otherwise be used for many timesteps after \(W\) has changed. Ribasim therefore resets \(\theta\) to 1 whenever \(W\) is rebuilt, like CVODE does. Without this reset, a stiff Basin can oscillate slightly every timestep, which keeps the timestep small for hours, as in the pump_pit test model.

5 Continuity considerations

The convergence of the Newton method can be proven given certain properties of \(\mathbf{f}\) around the initial guess and root to find. An important aspect is the smoothness of \(\mathbf{f}\). The basin profiles and \(Q(h)\) relations are given by interpolated data, and thus we have some control over the smoothness of these functions by the choice of interpolation method. This is discussed further below. the Manning resistance is not mentioned here since it is given by an analytical expression.

Control mechanisms can change parameters of \(\mathbf{f}\) discontinuously, leading to discontinuities of \(\mathbf{f}\). This however does not yield problems for the time integration methods in DifferentialEquations.jl, since the callback mechanisms used to change these parameters make the solver take these discontinuities into account.

5.1 Basin profiles

The basin profiles affect \(\mathbf{f}\) in many ways, anywhere where a basin level or area is required.

Note

This section needs to be updated and extended after once this issue is resolved.

5.2 Q(h) relations

TabulatedRatingCurve nodes contribute to \(\mathbf{f}\) with terms of the following form:

\[ Q(h(u)) \]

where the continuity of this term is given by the least continuous of \(Q\) and \(h\).

5.3 Empty basins

Reduction factors are introduced at several points in the definition of \(\mathbf{f}\) to smooth out otherwise discontinuous transitions (e.g. the flow rate of a pump going to zero when the source basin dries out). If flows are not too large with respect to basin storage, this will prevent basins from reaching 0. Rather, the basin gets a very small storage. The reduction factors help with performance, but are also an important tool to avoid getting negative storage in basins. Negative storage needs to be avoided since it is not a real solution, and would introduce water into the model that doesn’t exist. Another tool used to avoid negative storage is the isoutoutofdomain option, which Ribasim makes use of. This rejects timesteps that lead to negative storage, instead retrying with a smaller timestep.

6 Error control for cumulative flows

The states in \(\mathbf{u}\) are cumulative quantities that grow monotonically over time. With the default SciML error control the tolerance per state would be roughly \(\text{abstol} + |\mathbf{u}| * \text{reltol}\), see Basic Stepsize Control. That means the tolerance is based on the cumulative flow so far rather than on the flow over the current timestep, so it keeps loosening as the simulation progresses: after a thousand seconds at \(1\text{ m}^3/\text{s}\) an effective reltol of \(10^{-5}\) has degraded to \(10^{-2}\). It also makes the tolerance depend on arbitrary offsets: a large Basin with little throughflow would get a very loose tolerance purely because of its storage.

Ribasim therefore overrides calculate_residuals so that each state is scaled by the magnitude of its change over the timestep instead of its absolute magnitude:

\[ \text{residual}_i = \frac{\tilde{u}_i}{\text{abstol} + |u_{1,i} - u_{0,i}| * \text{reltol}} \]

This keeps the meaning of abstol and reltol constant over the whole simulation, and makes them independent of the size of the Basin storages. The same scaling is applied to the increment of the nonlinear solver, so that the nonlinear solve and the error estimate weigh the states in the same way.

Note that this changes the units in which the tolerances should be interpreted. abstol is an absolute tolerance on the volume passed by a state during a single timestep (in \(\text{m}^3\)), and reltol is relative to that same per-step volume. It acts as a floor for states that carry almost no flow during a timestep, where the reltol term vanishes; setting it too low makes such states dictate the timestep for the whole model. Since it is a volume per timestep and not a rate, abstol should not be scaled with dt.

The default of \(10^{-4}\ \text{m}^3\) is the largest value that is accurate across model scales, from a handful of Basins up to a national model. Because it is an absolute volume it does not scale with the model, so a large model can often raise it to \(10^{-2}\ \text{m}^3\) for a substantial speedup at no measurable loss of accuracy, while the same setting on a small model is roughly two orders of magnitude less accurate. If you raise it, verify against a run with a tighter tolerance.

7 Substance transport

The mass balance of a substance in Basin \(i\) is \[ \frac{\text{d}}{\text{d}t} \left( S_i C_i \right) = \sum_j q_{ji} C_j - q_{\text{out},i} C_i + f_i. \] It is solved after each timestep of the water balance, using the flows of that timestep. Integrating it over a timestep from \(t_n\) to \(t_{n+1}\) gives

\[ S_i^{n+1} C_i^{n+1} - S_i^n C_i^n = \sum_j \int_{t_n}^{t_{n+1}} q_{ji} C_j \,\text{d}t - \int_{t_n}^{t_{n+1}} q_{\text{out},i} C_i \,\text{d}t + F_i, \]

where the superscripts denote the time and \(F_i\) is the mass from the other sources over the timestep.

How the integrals over the flows are approximated determines the numerical scheme. The water balance gives the volumes \(V_{ji}\) and \(V_{\text{out},i}\) that flowed over the timestep. Ribasim combines these with one concentration per Basin for the entire timestep, a weighted average of its concentrations at the start and the end of the timestep:

\[ \int_{t_n}^{t_{n+1}} q_{ji} C_j \,\text{d}t \approx V_{ji} \left( \theta_j C_j^{n+1} + (1 - \theta_j) C_j^n \right). \]

Using this for all flows and moving the unknown concentrations to the left hand side gives a linear system for the concentrations of all Basins at the end of the timestep:

\[ \left( S_i^{n+1} + \theta_i V_{\text{out},i} \right) C_i^{n+1} - \sum_j \theta_j V_{ji} C_j^{n+1} = \left( S_i^n - (1 - \theta_i) V_{\text{out},i} \right) C_i^n + \sum_j (1 - \theta_j) V_{ji} C_j^n + F_i. \tag{6}\]

The weights that Ribasim uses are determined per Basin and per timestep:

\[ \theta_i = \max\left( \frac{1}{2}, 1 - \frac{S_i^n}{V_{\text{out},i}} \right). \]

The weight is \(\theta_i = 1/2\) wherever possible, as this is the trapezoidal rule, which is second order accurate. The concentrations can however only be guaranteed to stay non-negative if the coefficient of \(C_i^n\) in Equation 6 is non-negative, which requires \(\theta_i \geq 1 - S_i^n / V_{\text{out},i}\). So where more water flows out of a Basin within one timestep than twice its storage, the weight increases towards 1, using the smallest value that keeps the concentrations non-negative.

With these weights the matrix of Equation 6 has non-positive off-diagonal entries and is strictly diagonally dominant by column, since a Basin keeps a small positive storage. The system therefore always has a unique and non-negative solution, regardless of the size of the timestep. The mass that leaves one Basin over a flow is the mass that enters the other, so the scheme conserves mass, and a substance can travel through any number of Basins within one timestep. All substances share the same sparse matrix, which is factorized once per timestep.

If all Basins would use the same weight, the scheme reduces to these well known methods:

\(\theta\) Method Delwaq scheme Properties
\(0\) Forward Euler (explicit upwind) 1 First order. Needs no linear system. Non-negative if \(V_{\text{out},i} \leq S_i^n\). A substance travels at most one Basin per timestep.
\(1/2\) Trapezoidal rule (Crank-Nicolson) Second order. Non-negative if \(V_{\text{out},i} \leq 2 S_i^n\).
\(1\) Backward Euler (implicit upwind) 15 First order, with more numerical dispersion. Non-negative for any timestep.

The Delwaq coupling uses scheme 15. Adjusting the weight per Basin is also the idea of the local-theta schemes 21 and 22 of Delwaq, which however use \(\theta_i = 0\) wherever that is stable, and limit the fluxes to avoid new extremes. With instantaneously mixed Basins \(\theta_i = 1/2\) is more accurate, since an explicit step in which the full storage of a Basin flows out moves all of its water to the next Basin without mixing.

The ResidenceTime tracer has the integral of the storage over the timestep as a source in \(F_i\), which also uses the trapezoidal rule. The UserDemand return flows are part of \(F_i\) as well, with the concentrations of the Basins they came from at the start of the timestep.

Only the volume of a flow over a timestep is used, in its net direction, so effectively flows are constant over a timestep. The concentrations and loads from the input are constant over a timestep as well, since the solver stops at each time they change, so a pulse in one of these is followed exactly.