SIC^2: Simulation and Integration of Control for Canals

Exchange laws

Louis Poirel — 2012-06-21T14:35:03Z

An exchange law modifies one or more quality classes, based on parameters specific to the law, hydraulic variables, and possibly concentrations of some quality classes.

Thus, the exchange terms relative to law $i$ are calculated as follows:

$$\left(E_{k_1}^i, \dots, E_{k_n}^i\right)=L^i\left(C_{k_1}, \dots, C_{k_n}, p_1, \dots, V, h, \dots\right)$$

The exchange term of a class $k$ is the sum of its exchange term in every law where this class is present :

$$ E_k=E_{k}^{i_1}+E_{k}^{i_2}+\cdots+E_{k}^{i_n}$$

Engelund-Hansen (1967)

Louis Poirel — 2012-06-21T14:30:37Z

Engelund-Hansen formula

Transport capacity is calculated as follows :

$$C_{eq}=0.05\rho_S\frac{LU^2}{Q}\frac{(JR)^{\frac{3}{2}}}{\sqrt{g} (\rho_S/\rho-1)^2d}$$

where :

$ \rho_S$ is the sediment's density (kg/m³)
$ L $ is the stream's width (m)
$ U $ is the mean velocity (m/s)
$ Q $ is the water discharge (m³/s)
$ J $ is the slope (m/m)
$ R $ is the hydraulic radius (m)
$ g $ is the gravity (m/s²)
$ \rho $ is the density of water (kg/m³)
$ d $ is the sediment diameter

Specifications

Law's ID : 561
Number of acting classes : 4
Number of parameters : 6

Acting classes :

$ C_i $ : transported class modified by the law
$ C_j $ : fix class modified by the law
$ C_k $ : transported parameter class ($i=k$ for standard use)
$ T $ : water temperature

Parameters :

$ d $ : sediment diameter
$ \rho_S$ : sediment density
$ p $ : sediment porosity
$ \alpha$
$ i_{ech}$ : exchange formula : 1 = Han, 2 = Hazen
$\beta$

Sediment transport

Louis Poirel — 2012-06-21T14:08:26Z

Generalities

The calculation of sediment exchanges involves three steps:

Calculating an equilibrium concentration $C_{eq}$ according to a proposed law;
Calculating an adaptation time following of Han's formula$t_A=\beta\frac{Ru*}{VW}$ or Hazen's formula: $t_A=\beta\frac{R}{W}$ where $W$ is the fall velocity;
Calculation of the exchange term$E=\frac {\alpha C_{eq}-C }{t_A}$ ;
where $\alpha$ and $\beta$ are dimensionless parameters.

The fall rate is calculated according to Zanke's formula:

$$W=1.1*10\frac{\nu}{d}\left(\sqrt{1+\frac{0.01 g (\rho_S/\rho-1)d^3}{\nu^2}}-1\right)$$

This law is equal to Stokes' law for small diameters, and Newton's law for larger particles.

Algual detachment in response to flushing with reference shear stress equal to the first time step

Louis Poirel — 2012-06-21T13:49:15Z

Definition

The rate of accidental detachment of fix algae $ b_i $ towards drift algae $ a_i $ is calculated in each section as follows :

$$ \frac{\partial B_i}{\partial t}(x,t) = -S\frac{\partial A_i}{\partial t}(x,t) = -\frac{1}{\delta}\left ( \frac{\tau_{0}(x,t) - \tau_{0}(x,0)}{\tau_{0}(x,0)} - s_B \right )^\eta B_j(x,t) $$

if $ \frac{\tau_{0}(x) - \tau_{0}(x,0)}{\tau_{0}(x,0)} > s_B $

$ \frac{\partial A_i}{\partial t}(x,t) = \frac{\partial B_i}{\partial t}(x,t) = 0 $ atherwise.

Where :

$ A_i(x,t) $ : drift algae (kg/m³)
$ B_i(x,t) $ : fix algae (kg/m)
$ B_j(x,t) $ : fix algae ($i=j$ for standard applications).
$ S(x) $ : e cross sectional area (m²)
$ \tau(x,t) $ : shear stress (N m^-2)
$ \tau_{0}(x,0) $ : shear stress at $t=0$
$ s_B $ : sensitivity treshold
$ \delta $ : time constant (s)
$ \eta $ : adimensional exponent

Specifications

Law's ID : 341
Number of acting classes : 3
Number of parameters : 3

Acting classes :

$ A_i $ : drift class modified by the law
$ B_i $ : fix class modified by the law
$ B_j $ : the parameter class of the law

Parameters :

$ s_B $ : sensitivity treshold
$ \delta $ : time constant
$ \eta $ : adimensional exponent

Algal growth (PhD Thesis, O. Fovet, 2010, p.101)

Louis Poirel — 2012-06-21T12:49:49Z

Definition

The biomass growth is calculated at time t, in each section x as follows :

$\frac{\partial B}{\partial t}(x,t) = \mu(x,t) B(x,t) F_{lim}(B'(x,t))$

where :

$F_{lim}(B'(x,t)) = \left ( 1 - \frac{B'(x,t)}{B_{Max}} \right )$
$ \mu(x,t) = \mu_{0} \theta^{T(t)-T_{0}} \frac{I(x,t)}{I_{opt}} e^{1- \frac{I(x,t)}{I_{opt}}} \textup{min} \left ( \frac{N_{i}(x,t)}{N_{i}(x,t) + K_{N_{I}}} \right )$
$ I(x,t) = I_{s}(x,t) e^{-k_{ext} h(x,t)} $
$ I_{s}(x,t) = (1-C_m(x))(1-a)R_N(t) $

Variables and parameters

$ B(x,t) $ : fix biomass modified by the law (kg m^-1)
$ B'(x,t) $ : fix biomass parameter in the law (kg m^-1)
$ B_{Max} $ : maximum value of fix biomass (kg m^-1)
$ \mu_{0} $ : reference growth rate (s ^-1)
$ \theta $ : growth coefficient
$ T(x,t) $ : water temperature (°C)
$ T_{0} $ : reference temperature (°C)
$ I(x,t) $ : light intensity at the section's bottom (W m^-2)
$ I_{s}(x,t) $ : solar light intensity (W m^-2)
$ C_{m}(x)$ : mask coefficient (cf. Temperature simulation)
$ a $ : albedo (cf. Temperature simulation)
$ R_{N}(t) $ solar radiation (W m^-2) (cf. Temperature simulation)
$ k_{ext} $ : extinction coefficient (due to turbidity)
$ h(x,t) $ : mean water level (m)
$ I_{opt}(x,t) $ : optimum light intensity (W m^-2)
$ N_{i}(x,t) $ : concentration of nutrient i (kg m^-3)
$ K_{N_{i}} $ : half-saturation constant of nutrient i (kg m^-3)

Limiting nutrient

This law can take into account up to 3 nutrients for algae growth. If less than three nutrients are used, $ K_{N_{i}} $ for unused nutrients should be set to 0.

Specifications

Law's ID : 301
Number of acting classes: 6
Number of meteo parameters : 3
Number of parameters : 9

Acting classes :

$ B(x,t) $ :the class modified by the law
$ B'(x,t) $ : the parameter class of the law
$ T(x,t) $ :Water temperature
$ N_{1}(x,t) $ : Nutrient 1
$ N_{2}(x,t) $ : Nutrient 2
$ N_{3}(x,t) $ : Nutrient 3

Meteo parameters :

$ C_{m}(x)$
$ a $
$ R_{N}(t) $

Parameters

$ B_{Max} $ : maximum value of fix biomass
$ \mu_{0} $ : reference growth rate
$ \theta $ : growth coefficient
$ T_{0} $ : reference temperature
$ k_{ext} $ : extinction coefficient
$ I_{opt} $ : optimum light intensity
$ K_{N_{1}} $ : half-saturation constant of nutrient 1
$ K_{N_{2}} $ : half-saturation constant of nutrient 2
$ K_{N_{3}} $ : half-saturation constant of nutrient 3

Exponential growth with fixed coefficient

Louis Poirel — 2012-06-21T06:19:57Z

Definition

This law changes the concentration of a class depending on an other class concentration (or eventually itself), and two fix coefficients $ k $ and $\alpha_k$ as shown here:

$ \frac{dC_i}{dt}=k C_j^{\alpha_k} $

Classic use

Numerous solutes have a first-order kinetics with an equation similar to :

$ \frac{dC_{Ni}}{dt}= - k_{Ni} C_{Ni} $

where $ k_{Ni} $ is the reaction constant (which is the inverse of a time). For instance, bacteriological oxygen demand (DBO_5 $ $), has a constant around 0.3 days ^-1. This degradation models the process of water auto-purification in streams.

Specifications

Law's ID : 201
Number of acting classes : 2
Number of parameters : : 2

Acting classes :

$ C_i $ : the class modified by the law
$ C_j $ : the parameter class of the law

Parameters :

$ k $ : reaction constant
$ \alpha_k $ : reaction order

Evolution of the temperature

Louis Poirel — 2012-06-21T06:19:52Z

Definition

This law changes water temperature depending on meteorology and water temperature as explained in «Temperature simulation».

Specifications

Law's ID : 101
Number of acting classes : 1
Number of parameters : 0

Acting classes :

Water temperature

Temperature simulation

Louis Poirel — 2012-06-20T08:28:48Z

Transport and exchange equation

Water temperature in the network is a fonction of upstream temperature and thermal exchange with atmosphere and substratum.

The general equation is:

$\frac{\partial ST}{\partial t}+ \frac{\partial QT}{\partial x}= \frac{\partial}{\partial x}\left(DS\frac{\partial T}{\partial x}\right) + SE_T$

where $E_T$ is an exchange rate (°C/s) given by the heat exchange rate $P$ (W/m²). For a water colomn with a width of $\delta x$, energy balance is

$\rho S \delta x \ c_p E_T = P L_m \delta x$

where $\rho$ is the mass density of water, $c_p$ is the heat capacity of water, $L_m$ is the reach's width, $P$ is the heat power per unit area, yielding

$E_T = \frac{L_m P}{\rho c_p S}$

Numerical values are :

$\rho$ = 1000 kg/m³

$c_p$ = 4186 J/kg/°C

Calculation of the heat power term

The heat power per unit area is a balance of direct incoming solar and atmospheric radiation, minus the water surface radiation, and sensible and latent heat flow. The heat exchange between water and soil is difficult to evaluate since it depends on the soil's temperature and thermal parameters.

A mask coefficient $C_m$ is defined as a function of abscissa, which takes into account the potential coverage of the canal (passage in galleries in particular). Full coverage will cancel all exchange between water and atmosphere. Neglecting interactions with the substrate amounts to considering the environment as adiabatic or to assume that these exchanges are low with respect to the convection process of the water mass.

The heat power term is thus calculated :

$P = (1 - C_m)((1 - a) R_N + R_a - R_e - H_s - H_e)$

Solar radiation is a measured data that can be estimated using the solar radiation above atmosphere, the atmospheric attenuation and the time and location.Albedo $a$ is generally low for water, around 0.03.

Atmospheric radiation is given by Stefan Boltzmann formula, and takes into account water reverberation, clouds and air relative humidity ($w_a$, between 0 and 1) using Brutsaert formula (1982) [1]. We therefore need to calculate the equilibrium vapor pressure $e_s$ (Pa) and vapor pressure $e_e$ (en Pa) at air temperature $T_a$ (°C) with a humidity $w_a$:

$e_s = 101300 \exp \left (13.7 - \frac{5120}{273.15 + T_a}\right )$

$e_{a} = w_{a} e_{s}$

$c_a = 1,24 (1 - a)\left (\frac{e_a / 100}{273.15 + T_a}\right )^{1 / 7}$

$R_a = c_a \sigma (273.15 + T_a)^4$

where $\sigma = 5,67 10^{-8}$ W m^-2K^-4 is Stefan-Boltzmann constant.

Similarly, water radiation is :

$R_e = \epsilon \sigma (273.15 + T_e)^4$

where $T_e$ is the water temperature (°C) and $\epsilon$ is water emissivity ($\epsilon=0,97$).

Sensible heat flow is linked to thermal gradient between air and water and to the wind speed $U_V$.
We use De Bruin's formula (1978), in W/m$^2$/Pa, as cited in [2]:

$f_V =0.029+0.021 U_V$

and

$H_s = C_B f_V (T_e - T_a)$

where $C_B = 63$ Pa/° C (Bowen coefficient).

Then, latent heat flow is calculated as follows :

$H_e = f_V (e_s - e_a)$

[1] W. Brutsaert. Evaporation into the atmosphere. Kluwer Acad. Pubs., 1982.

[2] C. Jacovides, G. Papaioannou, and P. Kerkides. Micro and large-scale parameters evaluation of evaporation from a lake. Agricultural Water Management, 13(2-4):263-272, 1988

Transiant calculation for fix classes

Louis Poirel — 2012-06-20T06:52:59Z

The transient equation for fix classes is :

$\frac{\partial C_k}{\partial t}(x,t)=E_k(x,t)$

where :

$C_k$ is the linear mass density ;
$E_k$ is the derivative of $C_k$ with respect to time.

The time discretization for $ j + 1 = t + \Delta t $ of this equation is
made by a semi-implicit Euler method :

$C_k^{j+1}-C_k^{j}-\frac{\Delta t}{2}\left(E_k^{j}+E_k^{j+1}\right)=0$

Resolution is then performed in the same Newton's method iteration as drifting quality classes.

The formula for iteration i is:

$C_k^{j+1,i}=C_k^{j+1,i-1}-\frac{f\left(C_k^{j+1,i-1}\right)}{f^\prime\left(C_k^{j+1,i-1}\right)}=C_k^{j+1,i-1}-\frac{b_{i}}{a_{i}}$

where :

$b_{i}=C_k^{j+1,i-1}-C_k^{j}-\frac{\Delta t}{2}\left(E_k^{j}+E_k^{j+1,i}\right)$
$a_{i}=1-\frac{\Delta t}{2}E_k^{\prime j+1,i}$

$E_k^{j + 1,i}$ and $E_k^{\prime j + 1,i}$ are fonctions of $C_k^{j + 1,i-1}$ and drifting classes concentrations $C_{k}^{\prime j + 1,i-1}$.

Transient calculation

Louis Poirel — 2012-06-20T06:42:45Z

Calculation at a section

Discretization for Preissmann Scheme

The advection equation (without diffusion) for drift classes is :

$\frac{\partial \mathit{CS}}{\partial t}+\frac{\partial \mathit{CQ}}{\partial x}=E(x,t,C)S$

Preissmann scheme discretization equations are :

$\frac{\partial U}{\partial t}\approx (1-\psi )\frac{U_{j}^{n+1}-U_{j}^{n}}{\Delta t}+\psi \frac{U_{j+1}^{n+1}-U_{j+1}^{n}}{\Delta t}$

$\frac{\partial U}{\partial x}\approx (1-\theta )\frac{U_{j+1}^{n}-U_{j}^{n}}{\Delta x_{j+1/2}}+\theta \frac{U_{j+1}^{n+1}-U_{j}^{n+1}}{\Delta x_{j+1/2}}$
$U(x,t)=\theta \left(\psi U_{j+1}^{n+1}+(1-\psi )U_{j}^{n+1}\right)+(1-\theta )\left(\psi U_{j+1}^{n}+(1-\psi )U_{j}^{n}\right)$

We note : $\Delta (U)_{i}=U_{i}^{n+1}-U_{i}^{n}$

These equations for section j+1 at time k+1 yield :

$\frac{(1-\psi )\Delta (\mathit{CS})_{j}+\psi .\Delta (\mathit{CS})_{j+1}}{\Delta t}+\frac{(\mathit{CQ})_{j+1}^{k}-(\mathit{CQ})_{j}^{k}+\theta \left(\Delta (\mathit{CQ})_{j+1}-\Delta (\mathit{CQ})_{j}\right)}{\Delta x}$
$=\left[(1-\psi )(\mathit{ES})_{j}^{k}+\psi (\mathit{ES})_{j+1}^{k}\right]+\theta \left[(1-\psi )\Delta (\mathit{ES})_{j}+\psi \Delta (\mathit{ES})_{j+1}\right]$

The transport equation can be rewritten as $F(C_{j+1}^{k+1})=0$ since $C_{j}^{k+1}$ is given by the upstream condition.

$\delta x\left[(1-\psi )\Delta (\mathit{CS})_{j}+\psi .\Delta (\mathit{CS})_{j+1}\right]+\Delta t\left[(\mathit{CQ})_{j+1}^{k}-(\mathit{CQ})_{j}^{k}+\theta \left(\Delta (\mathit{CQ})_{j+1}-\Delta (\mathit{CQ})_{j}\right)\right]$
$-\delta x\Delta t\left[(1-\psi )(\mathit{ES})_{j}^{k}+\psi (\mathit{ES})_{j+1}^{k}\right]-\delta x\Delta t\theta \left[(1-\psi )\Delta (\mathit{ES})_{j}+\psi \Delta (\mathit{ES})_{j+1}\right]=0$

$\delta x\left[(1-\psi )\left((\mathit{CS})_{j}^{k+1}-(\mathit{CS})_{j}^{k}\right)+\psi .\left((\mathit{CS})_{j+1}^{k+1}-(\mathit{CS})_{j+1}^{k}\right)\right]$
$+\Delta t\left[(\mathit{CQ})_{j+1}^{k}-(\mathit{CQ})_{j}^{k}+\theta \left(\left((\mathit{CQ})_{j+1}^{k+1}-(\mathit{CQ})_{j+1}^{k}\right)-\left((\mathit{CQ})_{j}^{k+1}-(\mathit{CQ})_{j}^{k}\right)\right)\right]$
$-\delta x\Delta t\left[(1-\psi )(\mathit{ES})_{j}^{k}+\psi (\mathit{ES})_{j+1}^{k}\right]$
$-\delta x\Delta t\theta \left[(1-\psi )\left((\mathit{ES})_{j}^{k+1}-(\mathit{ES})_{j}^{k}\right)+\psi \left((\mathit{ES})_{j+1}^{k+1}-(\mathit{ES})_{j+1}^{k}\right)\right]=0$

Resolution by conjugate gradient method

This equation $F(C_{j+1}^{k+1})=0$ is solved using a conjugate gradient method, requiring the computation of a Jacobian matrix.

Let $S(C_{j+1}^{k+1})= ^{t}F(C_{j+1}^{k+1})F(C_{j+1}^{k+1})$. The equation is solved by minimizing $S$, using the conjugate gradient method. One need to calculate the partial derivatives of exchange terms with respect to every quality class in order to calculate $F$'s Jacobian matrix $J$ which gives $S$'s gradient $\nabla S$.

Discretization at a node

Mass balance at a node

The mass balance equation at a node is integrated between $t$ et $t+\Delta t$:

$\int _{t}^{t+\Delta t}{\left(\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})\right)\mathit{dt}}-\int _{C_{\mathit{nd}}}^{C_{\mathit{nd}}+\Delta C_{\mathit{nd}}}{V_{\mathit{cas}}\mathit{dC}_{\mathit{nd}}}$
$+\int _{t}^{t+\Delta t}{\left(C_{\mathit{nd}}\left[\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right]\right)\mathit{dt}}-\int _{V_{\mathit{cas}}}^{V_{\mathit{cas}}+\Delta V_{\mathit{cas}}}{C_{\mathit{nd}}\mathit{dV}_{\mathit{cas}}}=0$

This integral is calculated using the trapezoidal rule formula :

$\frac{\Delta t}{2}\left(\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})\right)^{k}$
$+{\frac{\Delta t}{2}}\left(\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})\right)^{k+1}$
$-{\frac{V_{\mathit{cas}}^{k+1}-V_{\mathit{cas}}^{k}}{2}}\left[C_{\mathit{nd}}^{k}+C_{\mathit{nd}}^{k+1}\right]-\frac{C_{\mathit{nd}}^{k+1}-C_{\mathit{nd}}^{k}}{2}\left[V_{\mathit{Cas}}^{k}+V_{\mathit{Cas}}^{k+1}\right]$
$+{\frac{\Delta t}{2}}C_{\mathit{nd}}^{k}\left[\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right]^{k}$
$+{\frac{\Delta t}{2}}C_{\mathit{nd}}^{k+1}\left[\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right]^{k+1}=0$

The source term E is the only non-linear term in this equation.

Newton's method for a node with a pond

Let $f\left(C_{\mathit{nd}}^{i}\right)=0$ be the result of the integration at the ith iteration.

Newton's method gives :
$C_{\mathit{nd}}^{i}=C_{\mathit{nd}}^{i-1}-\frac{f\left(C_{\mathit{nd}}^{i-1}\right)}{f'\left(C_{\mathit{nd}}^{i-1}\right)}=C_{\mathit{nd}}^{i-1}-\frac{b^{i}}{a^{i}}$

$b^{i}=\frac{\Delta t}{2}\left[\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})+C_{\mathit{nd}}^{k}\left(\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right)\right]^{k}$
$+{\frac{\Delta t}{2}}\left[\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}}^{i-1})+C_{\mathit{nd}}^{i-1}\left(\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right)\right]^{k+1}$
$-V_{\mathit{cas}}^{k+1}C_{\mathit{nd}}^{k}-V_{\mathit{cas}}^{k}C_{\mathit{nd}}^{k+1,i-1}$

If $b^{i}=b_{1}+b_{2}^{i}$ where $b_1$ is constant over iterations et $b_2$ we get :

$b_{1}=\frac{\Delta t}{2}\left[\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})+C_{\mathit{nd}}^{k}\left(\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right)\right]^{k}$
$+{\frac{\Delta t}{2}}\left[\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}\right]^{k+1}-V_{\mathit{cas}}^{k+1}C_{\mathit{nd}}^{k}$

We define
$b_{11}^{k}=\left[\sum _{s_{b}Q_{b}>0}C_{b}s_{b}Q_{b}+\sum _{Q_{p}>0}C_{p}Q_{p}\right]^{k}$
and
$a_{1}^{k}=\left[\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}\right]^{k}$

yielding : $b_{1}=\frac{\Delta t}{2}\left(\left[b_{11}+V_{\mathit{Cas}}E(t,C_{\mathit{nd}})\right]^{k}+b_{11}^{k+1}\right)+C_{\mathit{nd}}^{k}\left(\frac{\Delta t}{2}\left[a_{1}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right]^{k}-V_{\mathit{cas}}^{k+1}\right)$

and $b_{2}^{i}=\frac{\Delta t}{2}\left[V_{\mathit{Cas}}E(t,C_{\mathit{nd}}^{i-1})+C_{\mathit{nd}}^{i-1}\left(a_{1}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right)\right]^{k+1}-V_{\mathit{cas}}^{k}C_{\mathit{nd}}^{k+1,i-1}$

The derivative with respect to $C_{nd}$ is :
$a^{i}=\frac{\Delta t}{2}\left[V_{\mathit{Cas}}E'(t,C_{\mathit{nd}}^{i-1})+\left(a_{1}-k_{\mathit{inf}}S_{\mathit{cas}}v_{\mathit{inf}}\right)\right]^{k+1}-V_{\mathit{cas}}^{k}$

As $C_{nd}$ is a vector, this equation is true for every quality class $C_{nd}(i)$.

The different equations may be interdependent due to the exchange rate E.

As a consequence, matrix a is non-diagonal due to E'. If this term is negligible compared to the others, the Newton's method can still solve the system.

Calculation at a node without a pond

A node without a pond is a nodes where the volume (or area) of the pod is zero at time k+1 and k.

The equation is then linear and the exact result is found within one Newton iteration. If distribution coefficients are heterogeneous, they are adjusted to ensure mass balance at the node.

The formula for $C_{nd}$ is :
$C_{\mathit{nd}}^{k+1}=C_{\mathit{nd}}^{k}-\frac{f\left(C_{\mathit{nd}}^{k}\right)}{f'\left(C_{\mathit{nd}}^{k}\right)}=C_{\mathit{nd}}^{k}-\frac{b}{a}$

Equations for a node with a pond are :
$b=\frac{\Delta t}{2}\left[b_{11}^{k}+b_{11}^{k+1}+C_{\mathit{nd}}^{k}\left(a_{1}^{k}+a_{1}^{k+1}\right)\right]$
et $a=\frac{\Delta t}{2}a_{1}^{k+1}$

And :
$C_{\mathit{nd}}^{k+1}=C_{\mathit{nd}}^{k}-\frac{b_{11}^{k}+b_{11}^{k+1}+C_{\mathit{nd}}^{k}\left(a_{1}^{k}+a_{1}^{k+1}\right)}{a_{1}^{k+1}}$

Given the adjustable coefficient :
$a_{1}^{k}=\left[\sum _{s_{b}Q_{b}<0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0}k_{p}Q_{p}\right]^{k}$ , we get :

$a_{1}^{k}=\left[\sum _{s_{b}Q_{b}<0,b_{i}=0}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0,b_{i}=0}k_{p}Q_{p}+k_{a}\left(\sum _{s_{b}Q_{b}<0,b_{i}=1}k_{b}s_{b}Q_{b}+\sum _{Q_{p}<0,b_{i}=1}k_{p}Q_{p}\right)\right]^{k}$