Forward Structure

In smash a forward/direct spatially distributed model is obtained by chaining differentiable hydrological-hydraulic operators via simulated fluxes:

• (optional) a descriptors-to-parameters mapping \(\phi\) either for parameters imposing spatial constrain and/or regional mapping between physical descriptor and model conceptual parameters, see mapping section.

• (optional) a snow operator \(\mathcal{M}_{snw}\) generating a melt flux \(m_{lt}\) which is then summed with the precipitation flux to feed the hydrological operator \(\mathcal{M}_{rr}\).

• A hydrological production operator \(\mathcal{M}_{rr}\) generating an elementary discharge \(q_t\) which feeds the routing operator.

• A routing operator \(\mathcal{M}_{hy}\) simulating propagation of discharge \(Q)\).

The operators chaining principle is presented in section forward and inverse problems statement (cf. Eq. 2 ) and the chaining fluxes are explicitated in the diagram below. The forward model obtained reads \(\mathcal{M}=\mathcal{M}_{hy}\left(\,.\,,\mathcal{M}_{rr}\left(\,.\,,\mathcal{M}_{snw}\left(.\right)\right)\right)\) .

This section describes the various operators available in smash with mathematical/numerical expression, input data \(\left[\boldsymbol{I},\boldsymbol{D}\right](x,t)\), tunable conceptual parameters \(\boldsymbol{\theta}(x,t)\) and simulated state and fluxes \(\boldsymbol{U}(x,t)=\left[Q,\boldsymbol{h},\boldsymbol{q}\right](x,t)\).

These operators are written below for a given pixel \(x\) of the 2D spatial domain \(\Omega\) and for a time \(t\) in the simulation window \(\left]0,T\right]\).

Diagram of input data, hydrological-hydraulic operators, simulated quantities of a forward model \(\mathcal{M}=\mathcal{M}_{hy}\left(\,.\,,\mathcal{M}_{rr}\left(\,.\,,\mathcal{M}_{snw}\left(.\right)\right)\right)\) (cf. Eq. 2); recall the composition principle is explained in section forward and inverse problems statement.

Snow operator \(\mathcal{M}_{snw}\)

zero (Zero Snow)

This snow operator simply means that there is no snow operator.

\[m_{lt}(x, t) = 0\]

with \(m_{lt}\) the melt flux.

ssn (Simple Snow)

This snow operator is a simple degree-day snow operator. It can be expressed as follows:

\[m_{lt}(x, t) = f\left(\left[S, T_e\right](x, t), k_{mlt}(x), h_s(x, t)\right)\]

with \(m_{lt}\) the melt flux, \(S\) the snow, \(T_e\) the temperature, \(k_{mlt}\) the melt coefficient and \(h_s\) the state of the snow reservoir.

Note

Linking with the forward problem equation Eq. 1

• Internal fluxes, \(\{m_{lt}\}\in\boldsymbol{q}\)

• Atmospheric forcings, \(\{S, T_e\}\in\boldsymbol{\mathcal{I}}\)

• Parameters, \(\{k_{mlt}\}\in\boldsymbol{\theta}\)

• States, \(\{h_s\}\in\boldsymbol{h}\)

The function \(f\) is resolved numerically as follows:

• Update the snow reservoir state \(h_s\) for \(t^* \in \left] t-1 , t\right[\)

\[h_s(x, t^*) = h_s(x, t-1) + S(x, t)\]

• Compute the melt flux \(m_{lt}\)

\begin{eqnarray} m_{lt}(x, t) = \begin{cases} 0 &\text{if} \; T_e(x, t) \leq 0 \\ \min\left(h_s(x, t^*), k_{mlt}(x)\times T_e(x, t)\right) &\text{otherwise} \end{cases} \end{eqnarray}

• Update the snow reservoir state \(h_s\)

\[h_s(x, t) = h_s(x, t^*) - m_{lt}(x, t)\]

Hydrological operator \(\mathcal{M}_{rr}\)

Hydrological processes can be described at pixel scale in smash with one of the availabe hydrological operators adapted from state-of-the-art lumped models.

gr4 (Génie Rural 4)

This hydrological operator is derived from the GR4 model [Perrin et al., 2003].

Hint

Helpful links about GR:

• Brief history of GR models

• Scientific papers

• GR models in a R package

Diagram of the gr4 like hydrological operator

It can be expressed as follows:

\[q_{t}(x, t) = f\left(\left[P, E\right](x, t), m_{lt}(x, t), \left[c_i, c_p, c_t, k_{exc}\right](x), \left[h_i, h_p, h_t\right](x, t)\right)\]

with \(q_{t}\) the elemental discharge, \(P\) the precipitation, \(E\) the potential evapotranspiration, \(m_{lt}\) the melt flux from the snow operator, \(c_i\) the maximum capacity of the interception reservoir, \(c_p\) the maximum capacity of the production reservoir, \(c_t\) the maximum capacity of the transfer reservoir, \(k_{exc}\) the exchange coefficient, \(h_i\) the state of the interception reservoir, \(h_p\) the state of the production reservoir and \(h_t\) the state of the transfer reservoir.

Note

Linking with the forward problem equation Eq. 1

• Internal fluxes, \(\{q_{t}, m_{lt}\}\in\boldsymbol{q}\)

• Atmospheric forcings, \(\{P, E\}\in\boldsymbol{\mathcal{I}}\)

• Parameters, \(\{c_i, c_p, c_t, k_{exc}\}\in\boldsymbol{\theta}\)

• States, \(\{h_i, h_p, h_t\}\in\boldsymbol{h}\)

The function \(f\) is resolved numerically as follows:

Interception

• Compute interception evaporation \(e_i\)

\[e_i(x, t) = \min(E(x, t), P(x, t) + m_{lt}(x, t) + h_i(x, t - 1)\times c_i(x))\]

• Compute the neutralized precipitation \(p_n\) and evaporation \(e_n\)

\begin{eqnarray} &p_n(x, t)& &=& &\max \left(0, \; P(x, t) + m_{lt}(x, t) - c_i(x) \times (1 - h_i(x, t - 1)) - e_i(x, t) \right)\\ &e_n(x, t)& &=& &E(x, t) - e_i(x, t) \end{eqnarray}

• Update the interception reservoir state \(h_i\)

\[h_i(x, t) = h_i(x, t - 1) + \frac{P(x, t) + m_{lt}(x, t) + e_i(x, t) - p_n(x, t)}{c_i(x)}\]

Production

• Compute the production infiltrating precipitation \(p_s\) and evaporation \(e_s\)

\begin{eqnarray} &p_s(x, t)& &=& &c_p(x) (1 - h_p(x, t - 1)^2) \frac{\tanh\left(\frac{p_n(x, t)}{c_p(x)}\right)}{1 + h_p(x, t - 1) \tanh\left(\frac{p_n(x, t)}{c_p(x)}\right)}\\ &e_s(x, t)& &=& &h_p(x, t - 1) c_p(x) (2 - h_p(x, t - 1)) \frac{\tanh\left(\frac{e_n(x, t)}{c_p(x)}\right)}{1 + (1 - h_p(x, t - 1)) \tanh\left(\frac{e_n(x, t)}{c_p(x)}\right)} \end{eqnarray}

• Update the production reservoir state \(h_p\)

\[h_p(x, t^*) = h_p(x, t - 1) + \frac{p_s(x, t) - e_s(x, t)}{c_p(x)}\]

• Compute the production runoff \(p_r\)

\begin{eqnarray} p_r(x, t) = \begin{cases} 0 &\text{if} \; p_n(x, t) \leq 0 \\ p_n(x, t) - (h_p(x, t^*) - h_p(x, t - 1))c_p(x) &\text{otherwise} \end{cases} \end{eqnarray}

• Compute the production percolation \(p_{erc}\)

\[p_{erc}(x, t) = h_p(x, t^*) c_p(x) \left(1 - \left(1 + \left(\frac{4}{9}h_p(x, t^*)\right)^4\right)^{-1/4}\right)\]

• Update the production reservoir state \(h_p\)

\[h_p(x, t) = h_p(x, t^*) - \frac{p_{erc}(x, t)}{c_p(x)}\]

Exchange

• Compute the exchange flux \(l_{exc}\)

\[l_{exc}(x, t) = k_{exc}(x) h_t(x, t - 1)^{7/2}\]

Transfer

• Split the production runoff \(p_r\) into two branches (transfer and direct), \(p_{rr}\) and \(p_{rd}\)

\begin{eqnarray} &p_{rr}(x, t)& &=& &0.9(p_r(x, t) + p_{erc}(x, t)) + l_{exc}(x, t)\\ &p_{rd}(x, t)& &=& &0.1(p_r(x, t) + p_{erc}(x, t)) \end{eqnarray}

• Update the transfer reservoir state \(h_t\)

\[h_t(x, t^*) = \max\left(0, h_t(x, t - 1) + \frac{p_{rr}(x, t)}{c_t(x)}\right)\]

• Compute the transfer branch elemental discharge \(q_r\)

\begin{eqnarray} q_r(x, t) = h_t(x, t^*)c_t(x) - \left(\left(h_t(x, t^*)c_t(x)\right)^{-4} + c_t(x)^{-4}\right)^{-1/4} \end{eqnarray}

• Update the transfer reservoir state \(h_t\)

\[h_t(x, t) = h_t(x, t^*) - \frac{q_r(x, t)}{c_t(x)}\]

• Compute the direct branch elemental discharge \(q_d\)

\[q_d(x, t) = \max(0, p_{rd}(x, t) + l_{exc}(x, t))\]

• Compute the elemental discharge \(q_t\)

\[q_t(x, t) = q_r(x, t) + q_d(x, t)\]

gr5 (Génie Rural 5)

This hydrological operator is derived from the GR5 model [Le Moine, 2008]. It consists in a gr4 like model stucture (see diagram above) with a modified exchange flux with two parameters to account for seasonal variaitons.

Hint

Helpful links about GR:

• Brief history of GR models

• Scientific papers

• GR models in a R package

Diagram of the gr5 like hydrological operator

It can be expressed as follows:

\[q_{t}(x, t) = f\left(\left[P, E\right](x, t), m_{lt}(x, t), \left[c_i, c_p, c_t, k_{exc}, a_{exc}\right](x), \left[h_i, h_p, h_t\right](x, t)\right)\]

with \(q_{t}\) the elemental discharge, \(P\) the precipitation, \(E\) the potential evapotranspiration, \(m_{lt}\) the melt flux from the snow operator, \(c_i\) the maximum capacity of the interception reservoir, \(c_p\) the maximum capacity of the production reservoir, \(c_t\) the maximum capacity of the transfer reservoir, \(k_{exc}\) the exchange coefficient, \(a_{exc}\) the exchange threshold, \(h_i\) the state of the interception reservoir, \(h_p\) the state of the production reservoir and \(h_t\) the state of the transfer reservoir.

Note

Linking with the forward problem equation Eq. 1

• Internal fluxes, \(\{q_{t}, m_{lt}\}\in\boldsymbol{q}\)

• Atmospheric forcings, \(\{P, E\}\in\boldsymbol{\mathcal{I}}\)

• Parameters, \(\{c_i, c_p, c_t, k_{exc}, a_{exc}\}\in\boldsymbol{\theta}\)

• States, \(\{h_i, h_p, h_t\}\in\boldsymbol{h}\)

The function \(f\) is resolved numerically as follows:

Interception

Same as gr4 interception, see GR4 Interception

Production

Same as gr4 production, see GR4 Production

Exchange

• Compute the exchange flux \(l_{exc}\)

\[l_{exc}(x, t) = k_{exc}(x) \left(h_t(x, t - 1) - a_{exc}(x)\right)\]

Transfer

Same as gr4 transfer, see GR4 Transfer

grd (Génie Rural Distribué)

This hydrological operator is derived from the GR models and is a simplified strucutre used in [Jay-Allemand et al., 2020].

Diagram of the grd hydrological operator, a simplified GR like

It can be expressed as follows:

\[q_{t}(x, t) = f\left(\left[P, E\right](x, t), m_{lt}(x, t), \left[c_p, c_t\right](x), \left[h_p, h_t\right](x, t)\right)\]

with \(q_{t}\) the elemental discharge, \(P\) the precipitation, \(E\) the potential evapotranspiration, \(m_{lt}\) the melt flux from the snow operator, \(c_p\) the maximum capacity of the production reservoir, \(c_t\) the maximum capacity of the transfer reservoir, \(h_p\) the state of the production reservoir and \(h_t\) the state of the transfer reservoir.

Note

Linking with the forward problem equation Eq. 1

• Internal fluxes, \(\{q_{t}, m_{lt}\}\in\boldsymbol{q}\)

• Atmospheric forcings, \(\{P, E\}\in\boldsymbol{\mathcal{I}}\)

• Parameters, \(\{c_p, c_t\}\in\boldsymbol{\theta}\)

• States, \(\{h_p, h_t\}\in\boldsymbol{h}\)

The function \(f\) is resolved numerically as follows:

Interception

• Compute the interception evaporation \(e_i\)

\[e_i(x, t) = \min(E(x, t), P(x, t) + m_{lt}(x, t))\]

• Compute the neutralized precipitation \(p_n\) and evaporation \(e_n\)

\begin{eqnarray} &p_n(x, t)& &=& &\max \left(0, \; P(x, t) + m_{lt}(x, t) - e_i(x, t) \right)\\ &e_n(x, t)& &=& &E(x, t) - e_i(x, t) \end{eqnarray}

Production

Same as gr4 production, see GR4 Production

Transfer

• Update the transfer reservoir state \(h_t\)

\[h_t(x, t^*) = \max\left(0, h_t(x, t - 1) + \frac{p_{r}(x, t)}{c_t(x)}\right)\]

• Compute the transfer branch elemental discharge \(q_r\)

\begin{eqnarray} q_r(x, t) = h_t(x, t^*)c_t(x) - \left(\left(h_t(x, t^*)c_t(x)\right)^{-4} + c_t(x)^{-4}\right)^{-1/4} \end{eqnarray}

• Update the transfer reservoir state \(h_t\)

\[h_t(x, t) = h_t(x, t^*) - \frac{q_r(x, t)}{c_t(x)}\]

• Compute the elemental discharge \(q_t\)

\[q_t(x, t) = q_r(x, t)\]

loieau (LoiEau)

This hydrological operator is derived from the GR model [Folton and Arnaud, 2020].

Hint

Helpful links about LoiEau:

• Database

Diagram of the loieau like hydrological operator

It can be expressed as follows:

\[q_{t}(x, t) = f\left(\left[P, E\right](x, t), m_{lt}(x, t), \left[c_a, c_c, k_b\right](x), \left[h_a, h_c\right](x, t)\right)\]

with \(q_{t}\) the elemental discharge, \(P\) the precipitation, \(E\) the potential evapotranspiration, \(m_{lt}\) the melt flux from the snow operator, \(c_a\) the maximum capacity of the production reservoir, \(c_c\) the maximum capacity of the transfer reservoir, \(k_b\) the transfer coefficient, \(h_a\) the state of the production reservoir and \(h_c\) the state of the transfer reservoir.

Note

Linking with the forward problem equation Eq. 1

• Internal fluxes, \(\{q_{t}, m_{lt}\}\in\boldsymbol{q}\)

• Atmospheric forcings, \(\{P, E\}\in\boldsymbol{\mathcal{I}}\)

• Parameters, \(\{c_a, c_c, k_b\}\in\boldsymbol{\theta}\)

• States, \(\{h_a, h_c\}\in\boldsymbol{h}\)

The function \(f\) is resolved numerically as follows:

Interception

Same as grd interception, see GRD Interception

Production

Same as gr4 production, see GR4 Production

Note

The parameter \(c_p\) is replaced by \(c_a\) and the state \(h_p\) by \(h_a\)

Transfer

• Split the production runoff \(p_r\) into two branches (transfer and direct), \(p_{rr}\) and \(p_{rd}\)

\begin{eqnarray} &p_{rr}(x, t)& &=& &0.9(p_r(x, t) + p_{erc}(x, t))\\ &p_{rd}(x, t)& &=& &0.1(p_r(x, t) + p_{erc}(x, t)) \end{eqnarray}

• Update the transfer reservoir state \(h_c\)

\[h_c(x, t^*) = \max\left(0, h_c(x, t - 1) + \frac{p_{rr}(x, t)}{c_c(x)}\right)\]

• Compute the transfer branch elemental discharge \(q_r\)

\begin{eqnarray} q_r(x, t) = h_c(x, t^*)c_c(x) - \left(\left(h_c(x, t^*)c_c(x)\right)^{-3} + c_c(x)^{-3}\right)^{-1/3} \end{eqnarray}

• Update the transfer reservoir state \(h_c\)

\[h_c(x, t) = h_c(x, t^*) - \frac{q_r(x, t)}{c_c(x)}\]

• Compute the direct branch elemental discharge \(q_d\)

\[q_d(x, t) = \max(0, p_{rd}(x, t))\]

• Compute the elemental discharge \(q_t\)

\[q_t(x, t) = k_b(x)\left(q_r(x, t) + q_d(x, t)\right)\]

vic3l (Variable Infiltration Curve 3 Layers)

This hydrological operator is derived from the VIC model [Liang et al., 1994].

Hint

Helpful links about VIC:

• Model overview

• References

• GitHub

Diagram of the vic3l like hydrological operator

It can be expressed as follows:

\[q_{t}(x, t) = f\left(\left[P, E\right](x, t), m_{lt}(x, t), \left[b, c_{usl}, c_{msl}, c_{bsl}, k_s, p_{bc}, d_{sm}, d_s, w_s\right](x), \left[h_{cl}, h_{usl}, h_{msl}, h_{bsl}\right](x, t)\right)\]

with \(q_{t}\) the elemental discharge, \(P\) the precipitation, \(E\) the potential evapotranspiration, \(m_{lt}\) the melt flux from the snow operator, \(b\) the variable infiltration curve parameter, \(c_{usl}\) the maximum capacity of the upper soil layer, \(c_{msl}\) the maximum capacity of the medium soil layer, \(c_{bsl}\) the maximum capacity of the bottom soil layer, \(k_s\) the saturated hydraulic conductivity, \(p_{bc}\) the Brooks and Corey exponent, \(d_{sm}\) the maximum velocity of baseflow, \(d_s\) the non-linear baseflow threshold maximum velocity, \(w_s\) the non-linear baseflow threshold soil moisture, \(h_{cl}\) the state of the canopy layer, \(h_{usl}\) the state of the upper soil layer, \(h_{msl}\) the state of the medium soil layer and \(h_{bsl}\) the state of the bottom soil layer.

Note

Linking with the forward problem equation Eq. 1

• Internal fluxes, \(\{q_{t}, m_{lt}\}\in\boldsymbol{q}\)

• Atmospheric forcings, \(\{P, E\}\in\boldsymbol{\mathcal{I}}\)

• Parameters, \(\{b, c_{usl}, c_{msl}, c_{bsl}, k_s, p_{bc}, d_{sm}, d_s, w_s\}\in\boldsymbol{\theta}\)

• States, \(\{h_{cl}, h_{usl}, h_{msl}, h_{bsl}\}\in\boldsymbol{h}\)

The function \(f\) is resolved numerically as follows:

Canopy layer interception

• Compute the canopy layer interception evaporation \(e_c\)

\[e_c(x, t) = \min(E(x, t)h_{cl}(x, t - 1)^{2/3}, P(x, t) + m_{lt}(x, t) + h_{cl}(x, t - 1))\]

• Compute the neutralized precipitation \(p_n\) and evaporation \(e_n\)

\begin{eqnarray} &p_n(x, t)& &=& &\max\left(0, P(x, t) + m_{lt}(x, t) - (1 - h_{cl}(x, t - 1)) - e_c(x, t)\right)\\ &e_n(x, t)& &=& &E(x, t) - e_c(x, t) \end{eqnarray}

• Update the canopy layer interception state \(h_{cl}\)

\[h_{cl}(x, t) = h_{cl}(x, t - 1) + P(x, t) - e_c(x, t) - p_n(x, t)\]

Upper soil layer evaporation

• Compute the maximum \(i_{m}\) and the corresponding soil saturation \(i_{0}\) infiltration

\begin{eqnarray} &i_{m}(x, t)& &=& &(1 + b(x))c_{usl}(x)\\ &i_{0}(x, t)& &=& &i_{m}(x, t)\left(1 - (1 - h_{usl}(x, t - 1))^{1/(1 - b(x))}\right) \end{eqnarray}

• Compute the upper soil layer evaporation \(e_s\)

\begin{eqnarray} e_s(x, t) = \begin{cases} e_n(x, t) &\text{if} \; i_{0}(x, t) \geq i_{m}(x, t) \\ \beta(x, t)e_n(x, t) &\text{otherwise} \end{cases} \end{eqnarray}

with \(\beta\), the beta function in the ARNO evaporation [Todini, 1996] (Appendix A)

• Update the upper soil layer reservoir state \(h_{usl}\)

\[h_{usl}(x, t) = h_{usl}(x, t - 1) - \frac{e_s(x, t)}{c_{usl}(x)}\]

Infiltration

• Compute the maximum capacity \(c_{umsl}\), the soil moisture \(w_{umsl}\) and the relative state \(h_{umsl}\) of the first two layers

\begin{eqnarray} &c_{umsl}(x)& &=& &c_{usl}(x) + c_{msl}(x)\\ &w_{umsl}(x, t - 1)& &=& &h_{usl}(x, t - 1)c_{usl}(x) + h_{msl}(x, t - 1)c_{msl}(x)\\ &h_{umsl}(x, t - 1)& &=& &\frac{w_{umsl}(x, t - 1)}{c_{umsl}(x)} \end{eqnarray}

• Compute the maximum \(i_{m}\) and the corresponding soil saturation \(i_{0}\) infiltration

\begin{eqnarray} &i_{m}(x, t)& &=& &(1 + b(x))c_{umsl}(x)\\ &i_{0}(x, t)& &=& &i_{m}(x, t)\left(1 - (1 - h_{umsl}(x, t - 1))^{1/(1 - b(x))}\right) \end{eqnarray}

• Compute the infiltration \(i\)

\begin{eqnarray} i(x, t) = \begin{cases} c_{umsl}(x) - w_{umsl}(x, t - 1) &\text{if} \; i_{0}(x, t) + p_n(x, t) > i_{m}(x, t) \\ c_{umsl}(x) - w_{umsl}(x, t - 1) - c_{umsl}(x)\left(1 - \frac{i_{0}(x, t) + p_n(x, t)}{i_m(x, t)}\right)^{b(x) + 1} &\text{otherwise} \end{cases} \end{eqnarray}

• Distribute the infiltration \(i\) between the first two layers, \(i_{usl}\) and \(i_{msl}\)

\begin{eqnarray} &i_{usl}(x, t)& &=& &\min((1 - h_{usl}(x, t - 1)c_{usl}(x), i(x, t))\\ &i_{msl}(x, t)& &=& &\min((1 - h_{msl}(x, t - 1)c_{msl}(x), i(x, t) - i_{usl}(x, t)) \end{eqnarray}

• Update the first two layers reservoir states, \(h_{usl}\) and \(h_{msl}\)

\begin{eqnarray} &h_{usl}(x, t)& &=& &h_{usl}(x, t - 1) + i_{usl}(x, t)\\ &h_{msl}(x, t)& &=& &h_{msl}(x, t - 1) + i_{msl}(x, t) \end{eqnarray}

• Compute the runoff \(q_r\)

\[q_r(x, t) = p_n(x, t) - (i_{usl}(x, t) + i_{msl}(x, t))\]

Drainage

• Compute the soil moisture in the first two layers, \(w_{usl}\) and \(w_{msl}\)

\begin{eqnarray} &w_{usl}(x, t - 1)& &=& &h_{usl}(x, t - 1)c_{usl}(x)\\ &w_{msl}(x, t - 1)& &=& &h_{msl}(x, t - 1)c_{msl}(x) \end{eqnarray}

• Compute the drainage flux \(d_{umsl}\) from the upper soil layer to medium soil layer

\[d_{umsl}(x, t^*) = k_s(x) * h_{usl}(x, t - 1)^{p_{bc}}\]

• Update the drainage flux \(d_{umsl}\) according to under and over soil layer saturation

\[d_{umsl}(x, t) = \min(d_{umsl}(x, t^*), \min(w_{usl}(x, t - 1), c_{msl}(x) - w_{msl}(x, t - 1)))\]

• Update the first two layers reservoir states, \(h_{usl}\) and \(h_{msl}\)

\begin{eqnarray} &h_{usl}(x, t)& &=& &h_{usl}(x, t - 1) - \frac{d_{umsl}(x, t)}{c_{usl}(x)}\\ &h_{msl}(x, t)& &=& &h_{msl}(x, t - 1) + \frac{d_{umsl}(x, t)}{c_{msl}(x)} \end{eqnarray}

Note

The same approach is performed for drainage in the medium and bottom layers. Hence the three first steps are skiped for readability and the update of the reservoir states is directly written.

• Update of the reservoirs states, \(h_{msl}\) and \(h_{bsl}\)

\begin{eqnarray} &h_{msl}(x, t)& &=& &h_{msl}(x, t - 1) - \frac{d_{mbsl}(x, t)}{c_{msl}(x)}\\ &h_{bsl}(x, t)& &=& &h_{bsl}(x, t - 1) + \frac{d_{mbsl}(x, t)}{c_{bsl}(x)} \end{eqnarray}

Baseflow

• Compute the baseflow \(q_b\)

\begin{eqnarray} q_b(x, t) = \begin{cases} \frac{d_{sm}(x)d_s(x)}{w_s(x)}h_{bsl}(x, t - 1) &\text{if} \; h_{bsl}(x, t - 1) > w_s(x) \\ \frac{d_{sm}(x)d_s(x)}{w_s(x)}h_{bsl}(x, t - 1) + d_{sm}(x)\left(1 - \frac{d_s(x)}{w_s(x)}\right)\left(\frac{h_{bsl}(x, t - 1) - w_s(x)}{1 - w_s(x)}\right)^2 &\text{otherwise} \end{cases} \end{eqnarray}

• Update the bottom soil layer reservoir state \(h_{bsl}\)

\[h_{bsl}(x, t) = h_{bsl}(x, t - 1) - \frac{q_b(x, t)}{c_{bsl}(x)}\]

Routing operator \(\mathcal{M}_{hy}\)

The following routing operators are grid based and adapted to perform on the same grid than the snow and production operators. They take as input a 8 direction (D8) drainage plan \(\mathcal{D}_{\Omega}\left(x\right)\) obtained by terrain elevation processing.

For all the following models, the 2D flow routing problem over the spatial domain \(\Omega\) reduces to a 1D problem by using the drainage plan \(\mathcal{D}_{\Omega}\left(x\right)\). The lattest, for a given cell \(x\in\Omega\) defines 1 to 7 upstream cells which surface discharge can inflow the current cell \(x\) - each cell has a unique downstream cell.

lag0 (Instantaneous Routing)

This routing operator is a simple aggregation of upstream discharge to downstream following the drainage plan. It can be expressed as follows:

\[Q(x, t) = f\left(Q(x', t), q_{t}(x, t)\right),\;\forall x'\in \Omega_x\]

with \(Q\) the surface discharge, \(q_t\) the elemental discharge and \(\Omega_x\) a 2D spatial domain that corresponds to all upstream cells flowing into cell \(x\), i.e. the whole upstream catchment. Note that \(\Omega_x\) is a subset of \(\Omega\), \(\Omega_x\subset\Omega\) and for the most upstream cells, \(\Omega_x=\emptyset\).

Note

Linking with the forward problem equation Eq. 1

• Surface discharge, \(Q\)

• Internal fluxes, \(\{q_{t}\}\in\boldsymbol{q}\)

The function \(f\) is resolved numerically as follows:

Upstream discharge

• Compute the upstream discharge \(q_{up}\)

\begin{eqnarray} q_{up}(x, t) = \begin{cases} 0 &\text{if} \; \Omega_x = \emptyset \\ \sum_{k\in\Omega_x} Q(k, t) &\text{otherwise} \end{cases} \end{eqnarray}

Surface discharge

• Compute the surface discharge \(Q\)

\[Q(x, t) = q_{up}(x, t) + \alpha(x) q_t(x, t)\]

with \(\alpha\) a conversion factor from \(mm.\Delta t^{-1}\) to \(m^3.s^{-1}\) for a single cell.

lr (Linear Reservoir)

This routing operator is using a linear reservoir to rout upstream discharge to downstream following the drainage plan. It can be expressed as follows:

\[Q(x, t) = f\left(Q(x', t), q_{t}(x, t), l_{lr}(x), h_{lr}(x, t)\right),\;\forall x'\in \Omega_x\]

with \(Q\) the surface discharge, \(q_t\) the elemental discharge, \(l_{lr}\) the routing lag time, \(h_{lr}\) the state of the routing reservoir and \(\Omega_x\) a 2D spatial domain that corresponds to all upstream cells flowing into cell \(x\). Note that \(\Omega_x\) is a subset of \(\Omega\), \(\Omega_x\subset\Omega\) and for the most upstream cells, \(\Omega_x=\emptyset\).

Note

Linking with the forward problem equation Eq. 1

• Surface discharge, \(Q\)

• Internal fluxes, \(\{q_{t}\}\in\boldsymbol{q}\)

• Parameters, \(\{l_{lr}\}\in\boldsymbol{\theta}\)

• States, \(\{h_{lr}\}\in\boldsymbol{h}\)

The function \(f\) is resolved numerically as follows:

Upstream discharge

Same as lag0 upstream discharge, see LAG0 Upstream Discharge

Surface discharge

• Update the routing reservoir state \(h_{lr}\)

\[h_{lr}(x, t^*) = h_{lr}(x, t) + \frac{1}{\beta(x)} q_{up}(x, t)\]

with \(\beta\) a conversion factor from \(mm.\Delta t^{-1}\) to \(m^3.s^{-1}\) for the whole upstream domain \(\Omega_x\).

• Compute the routed discharge \(q_{rt}\)

\[q_{rt}(x, t) = h_{lr}(x, t^*) \left(1 - \exp\left(\frac{-\Delta t}{60\times l_{lr}}\right)\right)\]

• Update the routing reservoir state \(h_{lr}\)

\[h_{lr}(x, t) = h_{lr}(x, t^*) - q_{rt}(x, t)\]

• Compute the surface discharge \(Q\)

\[Q(x, t) = \beta(x)q_{rt}(x, t) + \alpha(x)q_t(x, t)\]

with \(\alpha\) a conversion factor from from \(mm.\Delta t^{-1}\) to \(m^3.s^{-1}\) for a single cell.

kw (Kinematic Wave)

This routing operator is based on a conceptual 1D kinematic wave model that is numerically solved with a linearized implicit numerical scheme [Chow et al., 1998]. This is applicable given the drainage plan \(\mathcal{D}_{\Omega}\left(x\right)\) that enables reducing the routing problem to 1D.

The kinematic wave model is a simplification of 1D Saint-Venant hydraulic model. First the mass equation writes:

(1)\[\partial_{t}A+\partial_{x}Q =q\]

with \(\partial_{\square}\) denoting the partial derivation either in time or space, \(A\) the cross sectional flow area, \(Q\) the flow discharge and \(q\) the lateral inflows.

Assuming that the momentum equation reduces to

(2)\[S_0=S_f\]

with \(S_0\) the bottom slope and \(S_f\) the friction slope - i.e. a locally uniform flow with energy grade line parallel to the channel bottom. This momentum equation can be written as [Chow et al., 1998]:

(3)\[A=a_{kw} Q ^{b_{kw}}\]

with \(a_{kw}\) and \(b_{kw}\) two constants to be estimated - that can also be written using Manning friction law.

Injecting the momentum parameterization of Eq. 3 into mass equation Eq. 1 leads to the following one equation kinematic wave model [Chow et al., 1998]:

(4)\[\partial_{x}Q+a_{kw}b_{kw} Q^{b_{kw}-1}\partial_{t}Q=q\]

Hint

Helpful link about kinematic wave:

• Numerical Solution (page 294, section 9.6)

The solution of this equation can written as:

\[Q(x, t) = f\left(Q(x', t'), q_{t}(x, t'), \left[a_{kw}, b_{kw}\right](x)\right),\;\forall (x', t') \in \Omega_x\times[t-1, t]\]

with \(Q\) the surface discharge, \(q_t\) the elemental discharge, \(a_{kw}\) the alpha kinematic wave parameter, \(b_{kw}\) the beta kinematic wave parameter and \(\Omega_x\) a 2D spatial domain that corresponds to all upstream cells flowing into cell \(x\). Note that \(\Omega_x\) is a subset of \(\Omega\), \(\Omega_x\subset\Omega\) and for the most upstream cells, \(\Omega_x=\emptyset\).

Note

Linking with the forward problem equation Eq. 1

• Surface discharge, \(Q\)

• Internal fluxes, \(\{q_{t}\}\in\boldsymbol{q}\)

• Parameters, \(\{a_{kw}, b_{kw}\}\in\boldsymbol{\theta}\)

For the sake of clarity, the following variables are renamed for this section and the finite difference numerical scheme writting:

Renamed variables

Before	After
\(Q(x, t)\)	\(Q_i^j\)
\(Q(x, t - 1)\)	\(Q_{i}^{j-1}\)
\(q_t(x, t)\)	\(q_{i}^{j}\)
\(q_t(x, t - 1)\)	\(q_{i}^{j-1}\)

The function \(f\) is resolved numerically as follows:

Upstream discharge

Same as lag0 upstream discharge, see LAG0 Upstream Discharge

Note

\(q_{up}\) is denoted here \(Q_{i-1}^{j}\)

Surface discharge

• Compute the intermediate variables \(d_1\) and \(d_2\)

\begin{eqnarray} &d_1& &=& &\frac{\Delta t}{\Delta x}\\ &d_2& &=& &a_{kw} b_{kw} \left(\frac{\left(Q_i^{j-1} + Q_{i-1}^j\right)}{2}\right)^{b_{kw} - 1} \end{eqnarray}

• Compute the intermediate variables \(n_1\), \(n_2\) and \(n_3\)

\begin{eqnarray} &n_1& &=& &d_1 Q_{i-1}^j\\ &n_2& &=& &d_2 Q_{i}^{j-1}\\ &n_3& &=& &d_1 \frac{\left(q_i^{j-1} + q_{i}^{j}\right)}{2} \end{eqnarray}

• Compute the surface discharge \(Q_i^j\)

\[Q_i^j = Q(x, t) = \frac{n_1 + n_2 + n_3}{d_1 + d_2}\]