Conceptual model
One of the underlying principles of dynamic TOPMODEL is that the landscape can be broken up into hydrologically similar regions, or Hydrological Response Units (HRUs), where all the area within a HRU behaves in a hydrologically similar fashion. Further discussion and the connection of HRUs is outlined elsewhere.
This document outlines the conceptual structure and computational implementation of the HRU. The HRU is a representation of an area of the catchment with, for example, similar topographic, soil and upslope area characteristics.
The HRU is considered to be a area of catchment with outflow occurring across a specified width of it boundary. It is formed of four zones representing the surface water, which passes water between HRUs and drains to the root zone. The root zone characterises the interactions between evapotranspiration and precipitation and when full spills to the unsaturated zone. This drains to the saturated zone which also interacts with other HRUs. The behaviour of the saturated zone is modelled using a kinematic wave approximation. The zones and variables used below are shown in the schematic diagram.
Figure 1: Schematic of the Hill slope HRU
In the following section the governing equations of the HRU are given in a finite volume form. Approximating equations for the solution of the governing equations are then derived with associated implicit and semi-implicit numerical schemes. These are valid for the wide range of surface zone and saturated zone transmissivity profiles presented in the vignette.
In this section a finite volume formulation of the Dynamic TOPMODEL equations is derived for the evolution of the states over the time step \(t\) to \(t + \Delta t\). The lateral flux rates are considered at the start and end of the time step. The vertical flux rates are considered as constant throughout the time step, allowing for a semi-analytical solution of some storage zones.
Two appendices provide supporting information justifying the numerical schemes and an alternative derivation of the governing equations based on the infinitesimal vertical slab used in the original derivation of Dynamic TOPMODEL.
Notation
Table 1 outlines the notation used for describing an infinitesimal slab across the hillslope HRU.
The following conventions are used:
All properties such as slope angles and time constants are considered uniform along the hillslope.
Precipitation and Potential Evapotranspiration occur at spatially uniform rates.
All vertical fluxes \(r_{\star \rightarrow \star}\) are considered to spatially uniform and to be positive when travelling in the direction of the arrow, but can, in some cases be negative.
All lateral fluxes \(q_{\star}\) and \(q_{\star}\) are considered to be positive flows travelling downslope.
The superscripts \(^+\) and \(^-\) are used to denote variables for the outflow and inflow to the hillslope.
\(\left.y\right\rvert_t\) indicates the value of the variable \(y\) at time \(t\)
In the derivation storage values at a given cross section of the HRU are expressed as a depth and denoted with a \(\bar{}\). So integrating over the length of the hillslope we find
\[\begin{equation} \int w \bar{y} dx = y \end{equation}\]
| Quantity type | Symbol | Description | unit |
|---|---|---|---|
| Storage | \(s_{sf}\) | Surface excess storage | m\(^3\) |
| \(s_{rz}\) | Root zone storage | m\(^3\) | |
| \(s_{uz}\) | Unsaturated zone storage | m\(^3\) | |
| \(s_{sz}\) | Saturated zone storage deficit | m\(^3\) | |
| Vertical fluxes | \(r_{sf \rightarrow rz}\) | Flow from the surface excess store to the root zone | m\(^3\)/s |
| \(r_{rz \rightarrow uz}\) | Flow from the root zone to the unsaturated zone | m\(^3\)/s | |
| \(r_{uz \rightarrow sz}\) | Flow from unsaturated to saturated zone | m\(^3\)/s | |
| Lateral fluxes | \(q_{sf}\) | Lateral flow in the hillslope surface zone | m\(^3\)/s |
| \(q_{sz}\) | Lateral inflow to the hillslope surface zone | m\(^3\)/s | |
| Vertical fluxes to the HRU | \(p\) | Precipitation rate | m\(^3\)/s |
| \(e_p\) | Potential Evapotranspiration rate | m\(^3\)/s | |
| Other HRU properties | \(A\) | Plan area | m\(^2\) |
| \(w\) | Width of a hill slope cross section | m | |
| \(\Delta x\) | Effective length of the hillslope HRU (x = A/w$) | m | |
| \(\Theta_{sf}\) | Further properties and parameters for the solution of the surface zone | :-: | |
| \(\Theta_{rz}\) | Further properties and parameters for the solution of the root zone | :-: | |
| \(\Theta_{uz}\) | Further properties and parameters for the solution of the saturated zone | :-: | |
| \(\Theta_{sz}\) | Further properties and parameters for the solution of the saturated zone | :-: |
Finite Volume Formulation
In this section a finite volume formulation of the Dynamic TOPMODEL equations is derived for the evolution of the states over the time step \(t\) to \(t + \Delta t\). The lateral flux rates are considered at the start and end of the time step. The vertical flux rates are considered as constant throughout the time step, allowing for a semi-analytical solution of some storage zones.
Surface zone
The storage in the surface zone satisfies the mass balance equation
\[ \frac{ds_{sf}}{dt} = q^{-}_{sf} - r_{sf \rightarrow rz} - q^{+}_{sf} \]
where surface storage is increased by lateral downslope flow from upslope HRUs. Water flows to the root zone at a constant rate \(r_{sf \rightarrow rz}\), unless limited by the available storage at the surface or the ability of the root zone to receive water (for example if the saturation storage deficit is 0 and the root zone storage is full). In cases where lateral flow in the saturated zone produce saturation storage deficits water may be returned from the root zone to the surface giving negative values of \(r_{sf \rightarrow rz}\).
The outflow \(q^{+}_{sf}\) is considered a function of \(q^{-}_{sf}\), \(s_{sf}\), \(r_{sf \rightarrow rz}\) and parameters \(\Theta_{sf}\) denoted by \[ q^{+}_{sf} = \mathcal{F}\left(s_{sf},q^{-}_{sf},r_{sf \rightarrow rz},\Theta_{sf}\right) \] For the purposes the numerical scheme we limit the set to possible functions \(\mathcal{F}\) to those that are non-decreasing with respect to \(s_{sf}\) and take the value \(0\) when \(s_{sf}=0\). More formally \[ \frac{\partial}{\partial s_{sf}} \mathcal{F}\left(s_{sf},q^{-}_{sf},r_{sf \rightarrow rz},\Theta_{sf}\right) \geq 0 \] and \[ \mathcal{F}\left(0,q^{-}_{sf},r_{sf \rightarrow rz},\Theta_{sf}\right) = 0 \]
Root zone
The root zone gains water from precipitation and the surface zone. It loses water through evaporation and to the unsaturated zone. All vertical fluxes are considered to be spatially uniform. The root zone storage satisfies
\[\begin{equation} 0 \leq s_{rz} \leq s_{rzmax} \end{equation}\]
with the governing ODE \[\begin{equation} \frac{ds_{rz}}{dt} = p - \frac{e_p}{s_{rzmax}} s_{rz} + r_{sf\rightarrow rz} - r_{rz \rightarrow uz} \end{equation}\]
Fluxes from the surface and to the unsaturated zone are controlled by the level of root zone storage along with the state of the unsaturated and saturated zones.
For \(s_{rz} \leq s_{rzmax}\) then \(r_{rz \rightarrow uz} \leq 0\). Negative values of \(r_{rz \rightarrow uz}\) may occur only when water is returned from the unsaturated zone due to saturation caused by lateral inflow to the saturated zone.
When \(s_{rz} = s_{rzmax}\) then \[\begin{equation} p - e_p + r_{sf\rightarrow rz} - r_{rz \rightarrow uz} \leq 0 \end{equation}\] In this case \(r_{rz \rightarrow uz}\) may be positive if \[\begin{equation} p - e_p + r_{sf\rightarrow rz} > 0 \end{equation}\] so long as the unsaturated zone can receive the water. If \(r_{rz \rightarrow uz}\) is ‘throttled’ by the rate at which the unsaturated zone can receive water, then \(r_{sf\rightarrow rz}\) is adjusted (potentially becoming negative) to ensure the equality is met.
Unsaturated Zone
The unsaturated zone acts as a non-linear tank subject to the constraint \[\begin{equation} 0 \leq s_{uz} \leq s_{sz} \end{equation}\]
The governing ODE is written as \[\begin{equation} \frac{ds_{uz}}{dt} = r_{rz \rightarrow uz} - \frac{A s_{uz}}{T_d s_{sz}} \end{equation}\]
If water is able to pass freely to the saturated zone, then it flows at the rate \(\frac{A s_{uz}}{T_d s_{sz}}\). In the situation where \(s_{sz}=s_{uz}\) the subsurface below the root zone can be considered saturated, as in there is no further available storage for water, but separated into parts: an upper part with vertical flow and a lower part with lateral flux.
It is possible that the downward flux is constrained by the ability of the saturated zone to receive water. If this is the case \(r_{uz \rightarrow sz}\) occurs at the maximum possible rate and \(r_{rz \rightarrow uz}\) is limited to ensure that \(s_{uz} \leq s_{sz}\).
Saturated Zone
For a HRU the alteration of the storage in the saturated zone is given in terms of the storage deficit \(s_sz\) as: \[ \frac{ds_{sz}}{dt} = q^{+}_{sz} - r_{sf \rightarrow rz} - q^{-}_{sz} \]
The relationship between the outflow and storage deficit is given \[ q^{+}_{sz} = \mathcal{G}\left(s_{sz},q^{-}_{sz},\Theta_{sz}\right) \] It is assummed that flow increases with storage so \[ - \frac{d}{ds_{sz}} \mathcal{G}\left(s_{sz},q^{-}_{sz},\Theta_{sz}\right) \geq 0 \] while at the maximum storage deficit \(\mathcal{G}\left(D,q^{-}_{sz},\Theta_{sz}\right) = 0\)
Numerical Solution
No analytical solution yet exists for simultaneous integration of the system of ODEs outlined above. In the following an implicit scheme, where fluxes between stores are considered constant over the time step and gradients are evaluated at the final state, is presented.
The basis of the solution is that a gravity driven system will maximise the downward flow of water within each timestep of size \(\Delta t\). Hence on the first downward pass through the stores the maximum downward volumes of the vertical fluxes \(\hat{v}_{* \rightarrow *} = \Delta t \hat{r}_{* \rightarrow *}\) are determined. These then moderated by the solution of the saturated zone, with an upward pass gives the final states.
Downward Pass
Surface excess
The implicit solution of mass balance equation gives
\[ \left. s_{sf} \right\rvert_{t+\Delta t} = \left. s_{sf} \right\rvert_{t} + \Delta t \left( \left. q_{sf}^{-} \right\rvert_{t+\Delta t} - \mathcal{F}\left(s_{sf},q^{-}_{sf},r_{sf \rightarrow rz},\Theta_{sf}\right) \right) - \left. v_{sf \rightarrow rz} \right\rvert_{t+\Delta t} \]
Since the outflow is 0 when the $s_sf =0 $ the maximum downward flux volume is \[ \left. \hat{v}_{sf \rightarrow rz} \right\rvert_{t+\Delta t} = \left. s_{sf} \right\rvert_{t} + \Delta t \left. q_{sf}^{-} \right\rvert_{t+\Delta t} \]
Root Zone
The implicit formulation for the root zone gives \[\begin{equation} \left. s_{rz} \right.\rvert_{t + \Delta t} = \left( 1 + \frac{e_{p}\Delta t}{s_{rzmax}} \right)^{-1} \left( \left. s_{rz} \right.\rvert_{0} + \Delta t p + v_{sf \rightarrow rz} - v_{rz \rightarrow uz} \right) \end{equation}\]
Since flow to the unsaturated zone occur only to keep \(s_{rz} \leq s_{rzmax}\) then \[\begin{equation} \hat{v}_{rz \rightarrow uz} = \max\left(0, \left. s_{rz} \right.\rvert_{0} + \Delta t \left( p -e_{p}\right) + \hat{v}_{sf \rightarrow rz} - s_{rzmax} \right) \end{equation}\]
Unsaturated Zone
The implicit expression for the unsaturated zone is given in terms of \(\left. s_{uz} \right.\rvert_{t+\Delta t}\) and \(\left. s_{sz} \right.\rvert_{t+\Delta t}\) as
\[\begin{equation} \left. s_{uz} \right.\rvert_{t+\Delta t} = \frac{T_{d}\left. s_{sz} \right.\rvert_{t+\Delta t}}{ T_{d}\left. s_{sz} \right.\rvert_{t+\Delta t} + A \Delta t} \left( \left. s_{uz} \right.\rvert_{t} + v_{rz \rightarrow uz}\right) \end{equation}\]
Since \(\left. s_{uz} \right.\rvert_{t+\Delta t} \leq \left. s_{sz} \right.\rvert_{t+\Delta t}\) this places an additional condition of \[\begin{equation} \left. s_{uz} \right.\rvert_{t} + v_{rz \rightarrow uz} \leq \left. s_{sz} \right.\rvert_{t+\Delta t} + \frac{A \Delta t}{T_d} \end{equation}\]
By mass balance this gives the maximum downward flux volume as \[ \hat{v}_{uz \rightarrow sz} = A\Delta t \min\left( \frac{\left( \left. s_{uz} \right.\rvert_{t} + \hat{v}_{rz \rightarrow uz}\right)} {T_{d}\left. s_{sz} \right.\rvert_{t+\Delta t} + A \Delta t}, \frac{1}{T_{d}} \right) \]
This is a function of \(\left. s_{sz} \right.\rvert_{t+\Delta t}\).
Saturated Zone
Building on the generic solution the implicit scheme for the saturated zone gives
\[ \left. s_{sz} \right\rvert_{t+\Delta t} = \left. s_{sz} \right\rvert_{t} + \Delta t \left( \mathcal{G}\left(s_{sz},q^{-}_{sz},\Theta_{sz}\right) - \left. q_{sz}^{-} \right\rvert_{t+\Delta t} \right) - \left. v_{uz \rightarrow sz} \right\rvert_{t+\Delta t} \]
Consider the pseudo-function \[ \mathcal{H}\left(z\right) = z - \left. s_{sz}\right.\rvert_{t} + \Delta t \left. q_{sz}^{-} \right\rvert_{t+\Delta t} + A\Delta t \min\left( \frac{\left( \left. s_{uz} \right.\rvert_{t} + \hat{v}_{rz \rightarrow uz}\right)} {T_{d}z + A \Delta t}, \frac{1}{T_{d}} \right) - \\ \Delta t \min\left(q_{sz}^{max},\max\left( 0, 2 \mathcal{G}\left(z,\Theta_{sz}\right) - \left. q_{sz}^{-} \right\rvert_{t+\Delta t} \right)\right) \]
adapted from the approximating equation for the saturated zone storage. In Appendix A it is shown that \(\mathcal{H}\left(z\right)\) is a monotonic increasing function of \(z > 0\) and \(\mathcal{H}\left(D\right) \geq 0\). Hence if \(\mathcal{H}\left(0\right) \lt 0\) unique solution for \(\left. s_{sz}\right.\rvert_{\Delta t}\) can be found. If \(\mathcal{H}\left(0\right) \geq 0\) then the incoming fluxes are enough to produce saturation in the subsurface and \(\left. s_{sz}\right.\rvert_{\Delta t} = 0\).
Maintaining a water (mass) balance and correct limits on the state values depends upon the accuracy of the solution of \(\mathcal{H}\left(z\right)=0\). Since each evaluation of \(\mathcal{H}\left(z\right)\) requires an evaluation of \(\mathcal{G}\) this rapidly become the most expensive part of the code. Use of a numerical scheme which allows mass balance to maintained regardless of the accuracy of the solution of \(\mathcal{H}\left(z\right)=0\) allows scope for compromises based on computational time.
The upward pass presented below is mass conservative so long as \(\mathcal{H}\left(\left. s_{sz}\right.\rvert_{\Delta t}\right)\geq0\). We therefore seek \(z\) which satifies for some tolerance \(\epsilon>0\) \[ \epsilon \geq $\mathcal{H}\left(z\right) \geq 0 \] One approach to finding such a value is use a bracketing algorithm. At each iteration the lower and upper limits \(z^{l} < z^{u}\) satisify \(\mathcal{H}\left(z^{l}\right)\leq 0 \leq \mathcal{H}\left(z^{u}\right)\), thus bracket the solution. Iterating (using for example bisection) until \(\mathcal{H}\left(z^{l}\right) \leq \epsilon\) give the approximation \(\left. s_{sz}\right.\rvert_{\Delta t} = z^{u}\).
Upward Pass
Given an estimate of \(\left. s_{sz}\right.\rvert_{\Delta t}\) such that \(\mathcal{H}\left(\left. s_{sz}\right.\rvert_{\Delta t}\right)>0\) compute the vertical flux volumes in the following fashion
- \(\left. q_{sz}^{+}\right.\rvert_{t+\Delta t} = \min\left(q_{sz}^{max},\max\left( 0, 2 \mathcal{G}\left(\left. s_{sz}\right.\rvert_{t+\Delta t},\Theta_{sz}\right) - \left. q_{sz}^{-} \right\rvert_{t+\Delta t} \right)\right)\)
- \(v_{uz \rightarrow sz} = \left. s_{sz}\right.\rvert_{t} + \Delta t \left(\left. q_{sz}^{+}\right.\rvert_{t+\Delta t} - \left. q_{sz}^{-} \right\rvert_{t+\Delta t}\right) - \left. s_{sz}\right.\rvert_{t+\Delta t}\)
- \(\left. s_{uz}\right.\rvert_{t+\Delta t} = \min\left(\left. s_{sz}\right.\rvert_{t+\Delta t},\left. s_{uz}\right.\rvert_{t} + \hat{v}_{rz \rightarrow uz} - v_{uz \rightarrow sz}\right)\)
- \(v_{rz \rightarrow uz} = \left. s_{uz}\right.\rvert_{t+\Delta t} -\left. s_{uz}\right.\rvert_{t} + v_{uz \rightarrow sz}\)
- \(v_{sf \rightarrow rz} = \min\left(\hat{v}_{sf \rightarrow rz}, s_{rzmax} - \left. \bar{s}_{rz} \right.\rvert_{0} - \Delta t \left( p - e_p \right) + {v}_{rz \rightarrow uz} \right)\)
- Solve for \(\left. s_{rz}\right.\rvert_{t+\Delta t}\) using the equation already given.
- Compute the volume loss to evapotranspiration as \(\left. v_{e_{p}}\right.\rvert_{t+\Delta t} = \left. s_{rz}\right.\rvert_{t} + v_{sf \rightarrow rz} - v_{rz \rightarrow uz} + \Delta t p - \left. s_{rz}\right.\rvert_{t+\Delta t}\)
This leaves the solution of the surface zone
Surface Zone
Proceeding in a similar fashion to the saturated zone consider the pseudo-function \[ \mathcal{S}\left(w\right) = \left. s_{sf} \right\rvert_{t} + \Delta t \left( \left. q_{sf}^{-} \right\rvert_{t+\Delta t} - \mathcal{F}\left(w,q_{sf}^{-},\frac{\left. v_{sf \rightarrow rz} \right\rvert_{t+\Delta t}}{\Delta t},\Theta_{sf}\right) \right) - \left. v_{sf \rightarrow rz} \right\rvert_{t+\Delta t} - w \] which is monotonically decreasing for \(w \geq 0\) with \(\mathcal{S}\left(0\right)\geq 0\).
A bracketing solution algorithm gives at any iteration \(w^{l} < \left. s_{sf} \right\rvert_{t+\Delta t} < w^{u}\). An initial value of \(w^{u}\) can be computed as \(\left. s_{sf} \right\rvert_{t} + \Delta t\left. q_{sf}^{-} \right\rvert_{t+\Delta t} - \left. v_{sf \rightarrow rz} \right\rvert_{t+\Delta t}\). Iterating until; for some tolerance \(\epsilon>0\); \(w^{u} - w^{l} < \epsilon\) take \(\left. s_{sf}\right.\rvert_{\Delta t} = w^{l}\). In this case \(\mathcal{S}\left(\left. s_{sf}\right.\rvert_{\Delta t}\right) > 0\) hence mass conservation can be maintained by taking \[ \left. q_{sf}^{+} \right\rvert_{t+\Delta t} = \left. q_{sf}^{-} \right\rvert_{t+\Delta t} + \frac{1}{\Delta t}\left( \left. s_{sf} \right\rvert_{t} - \left. v_{sf \rightarrow rz} \right\rvert_{t+\Delta t} - \left. s_{sf}\right.\rvert_{t + \Delta t}\right) \]
Surface Zone Representations
The storage in the surface zone satisfies the mass balance equation
\[ \frac{ds_{sf}}{dt} = q^{-}_{sf} - r_{sf \rightarrow rz} - q^{+}_{sf} \]
where surface storage is increased by lateral downslope flow from upslope HRUs. Water flows to the root zone at a constant rate \(r_{sf \rightarrow rz}\), unless limited by the available storage at the surface or the ability of the root zone to receive water (for example if the saturation storage deficit is 0 and the root zone storage is full).
In cases where lateral flows in the saturated zone produce saturation storage deficits water may be returned from the root zone to the surface giving negative values of \(r_{sf \rightarrow rz}\).
The following subsections present the various options for the surface zone solution starting from the first, that of the constant parameter Muskingham method.
Constant Parameter Muskingham solution
The constant parameter Muskingham method offers a parsimonious for capturing a range of surface dynamics. This sub-section outlines the key equations and shows how they may be used to capture a range of conceptual models.
Muskingham routing can be derived in a number of ways. The surface storage \(s_{sf}\) in a reach or hillslope of length \(\Delta x\) can be expressed in terms of the representative flow \(q_{sf}\) and velocity \(v_{sf}\) as \[ s_{sf} = \frac{q_{sf} \Delta x}{v_{sf}} \] The representative flow is realted to the inflow \(q_{sf}^{-}\) and outflow \(q_{sf}^{+}\) using the parameter \(\eta\) through \[ q_{sf} = \eta q_{sf}^{-} + \left(1-\eta\right) q_{sf}^{+} \] gives \[ s_{sf} = \frac{\Delta x}{v_{sf}} \left(\eta q_{sf}^{-} + \left(1-\eta\right) q_{sf}^{+} \right) \] Taking vertical inflow \(r\) the mass balance can be written \[ \frac{ds}{dt} = q^{-} + r - q^{+} \] where by substitution \[ \frac{ds}{dt} = q^{-} + r - \max\left(0, \frac{1}{1-\eta}\left(\frac{sv}{\Delta x} - \eta q^{-} \right) \right) \] where the maximum ensures that outflow is positive.
The numerical scheme is based on the implicit solution of the ODE. Let the initial state be \(s_{t}\) and the input rates \(q^{-}_{t+\Delta t}\) and \(r_{t+\Delta t}\) be known so
\[ s_{t+\Delta t} = \max\left(0, s_{t} + \Delta t \left( q^{-}_{t+\Delta t} + r_{t+\Delta t} - \frac{1}{1-\eta_{t+\Delta t}}\max\left(0,\frac{s_{t+\Delta t}v}{\Delta x} - \eta_{t} q^{-}_{t+\Delta t} \right) \right) \right) \]
where the outer maximum ensures that states are not negative. The maximum possible values of \(s_{t+\Delta t}\) is \(s_{m}=\max\left(0, s_{t} + \Delta t \left( q^{-}_{t+\Delta t} + r_{t+\Delta t}\right)\right)\) as the maximum possible values of \(s_{t+\Delta t}\) and outflow starts at \(s_{c} = \frac{\eta q^{-}_{t+\Delta t}\Delta x}{v}\)
If \(s_{m} \leq s_{c}\) then outflow cannnot occur and \(s_{t+\Delta t} = \max\left(0, s_{t} + \Delta t \left( q^{-}_{t+\Delta t} + r_{t+\Delta t} \right)\right)\) where for mass balance the value of \(r\) is limites so that \(\int_{t}^{t+\Delta t} r dt = s_{t+\Delta t} - s_{t} - \Delta t q^{-}_{t+\Delta t}\).
If \(s_{m} > s_{c}\) then outflow can occur and their exists a value \(s_{c} < s_{t+\Delta t} < s_{m}\) such that \[ s_{t+\Delta t} = s_{t} + \Delta t \left( q^{-}_{t+\Delta t} + r_{t+\Delta t} - \frac{1}{1-\eta_{t+\Delta t}}\left(\frac{s_{t+\Delta t}v}{\Delta x} - \eta_{t} q^{-}_{t+\Delta t} \right) \right) \]
which is given by
\[ s_{sf_{t+\Delta t}} = \frac{ \left(1-\eta\right) \Delta x}{ \left(1-\eta\right) \Delta x + v_{sf}\Delta t} \left( s_{_sf_t} + \Delta t \left( \frac{1}{1-\eta} q^{-}_{sf_{t+\Delta t}} + r_{t+\Delta t} \right)\right) \]
Compound channel
In a compound channel the convayance descriptors change from \(\left\{ v_{1}, \eta{1} \right\}\) to \(\left\{ v_{2}, \eta{2} \right\}\) when \(s>s_{1}\). The implicit solution for the mass balance can then be written as
\[ s_{t+\Delta t} = \max\left(0, s_{t} + \Delta t \left( q^{-}_{t+\Delta t} + r_{t+\Delta t} - \frac{1}{1-\eta_{1}}\max\left(0,\frac{\min\left(s_{1},s_{t+\Delta t}\right)v}{\Delta x} - \eta_{1} q^{-}_{1} \right) - \frac{1}{1-\eta_{2}}\max\left(0,\frac{\min\left(0,s_{t+\Delta t}-s_{1}\right)v}{\Delta x} - \eta_{2} q^{-}_{2} \right) \right) \right) \]
the storage exced
= (0, s_{t} + t ( q^{-}{t+t} + r{t+t} ))$ where for mass balance the value of \(r\) is limites so that \(\int_{t}^{t+\Delta t} r dt = s_{t+\Delta t} - s_{t} - \Delta t q^{-}_{t+\Delta t}\).
\[ \frac{ds}{dt} = q^{-}_{t+\Delta t} + r_{t+\Delta t} - \max\left(0, \frac{1}{1-\eta}\left(\frac{sv}{\Delta x} - \eta q^{-}_{t+\Delta t} \right) \right) \]
The solutions fall into four cases:
- Case 1: \(s_{t} \leq \frac{\eta q^{-}_{t+\Delta t}\Delta x}{v}\) and \(q^{-}_{t+\Delta t} + r_{t+\Delta t} \leq 0\). In this case there is no outflow at the insitial state and the storage volume it not increasing, hence outflow will not occur later in the time step so \(\int_{t}^{t+\Delta t} q^{+} dt = 0\), \(s_{t+\Delta t} = max\left(0,s_{t} + \Delta t\left(q^{-}_{t+\Delta t} + r_{t+ \Delta t}\right)\right)\) and to ensure mass balance \(\int_{t}^{t+\Delta t} r dt = s_{t+\Delta t} - s_{t} - \Delta t q^{-}_{t+\Delta t}\) where the maximum ensures that outflow is positive.
given by \[ s_{sf_{t+\Delta t}} = \frac{ \left(1-\eta\right) \Delta x}{ \left(1-\eta\right) \Delta x + v_{sf}\Delta t} \left( s_{_sf_t} + \Delta t \left( \frac{1}{1-\eta} q^{-}_{sf_{t+\Delta t}} + r_{t+\Delta t} \right)\right) \] subject to the constraint \(q_{sf_{t+\Delta t}}^{+} \geq 0\)
Approximating Tank Models
Taking \(\eta=0\) gives \[ \frac{ds}{dt} = q^{-} + r - \frac{sv}{\Delta x} \] which is the equation of a tank model with possibly varying time constant \(T = \frac{\Delta x}{v}\)
Approximating Diffusion Wave Routing
Diffuse routing with lateral inflow can be expressed as the parabolic equation (see for example doi:10.1016/j.advwatres.2005.08.008)
\[ \frac{dq}{dt} + c\frac{dq}{dx} - D \frac{d^2 q}{dx^2} = cl - D \frac{dl}{dx} \] where \(l\) is lateral inflow per unit length. Simplify this by considering that the lateral inflow \(r\) uniformly distributed so that \(l=\frac{r}{\Delta x}\) and \(\frac{dl}{dx}=0\) to give
\[ \frac{dq}{dt} + c\frac{dq}{dx} - D \frac{d^2 q}{dx^2} = cl \]
Let us relate this to the Muskingham Solution using the relationship \[ q = \eta q^{-} + \left(1-\eta\right) q^{+} \]
Taking Taylor series expansions for \(q^{-}\) and \(q^{+}\) based on \(q\) gives
\[ q^{+} \approx q + \eta \Delta x \frac{dq}{dx} + \frac{1}{2}\eta^2\Delta x^2 \frac{d^2 q}{dx^2} \] and \[ q^{-} \approx q - \left(1-\eta\right) \Delta x \frac{dq}{dx} + \frac{1}{2}\left(1-\eta\right)^2 \Delta x^2 \frac{d^2 q}{dx^2} \]
Subtracting the expression for \(q^{-}\) from that for \(q^{+}\) gives
\[ q^{-} - q^{+} \approx -\Delta x \frac{dq}{dx} + \frac{1}{2}\left(1-2\eta\right) \Delta x^2 \frac{d^2 q}{dx^2} \]
From the mass balance condition \[ \frac{ds}{dt} \approx \Delta x l -\Delta x \frac{dq}{dx} + \frac{1}{2}\left(1-2\eta\right)\Delta x^2 \frac{d^2 q}{dx^2} \]
Using the expression for storage and relationship between flow and area \[ \frac{dq}{dt} = \frac{dq}{da} \frac{da}{dt} = \frac{c}{\Delta x} \frac{ds}{dt} \]
Then combining this with the approximation produces \[ \frac{dq}{dt} + c\frac{dq}{dx} - \frac{c}{2}\left(1-2\eta\right)\Delta x \frac{d^2 q}{dx^2} \approx c l \]
Comparison to the original diffuse routing equation shows that \[ \eta \approx \frac{1}{2} - \frac{D}{c \Delta x} \]
The value of \(\eta\) is limited to between 0 and \(\frac{1}{2}\). Firstly \(\eta = 1/2\) results from \(D=0\); that is the Kinematic wave equation. Taking \(\eta=0\) produces the maximum dispersion of \(D = \frac{c\Delta x}{2}\), which is given by the linear tank.
Numerical scheme
The numerical scheme is based on the implicit solution of
\[ \frac{ds}{dt} = q^{-} + r - \frac{1}{1-\eta}\left(\frac{sv}{\Delta x} - \eta q^{-} \right) \]
given by \[ s_{t+\Delta t} = s_{t} + \Delta t \left( q^{-}_{t+\Delta t} + r_{t+\Delta t} - \frac{1}{1-\eta_{t+\Delta t}}\max\left(0,q_{t+\Delta t} - \eta_{t} q^{-}_{t+\Delta t} \right) \right) \]
where the maximum ensures that \(q^{+}_{t+\Delta t}\) does not become negative. The implementation outlines is equivalent to the addition of dispersion to ensure non-negative outflows.
For known \(\eta_{t + \Delta t}\) it is clear that the left hand side of the expression is increasing and the right hand side decreasing so long as \(\frac{dq}{ds} = \frac{1}{\Delta x}\frac{dq}{da} \geq 0\). If \(\eta_{t+ \Delta t}\) is a function of \(s_{t+ \Delta t}\) then it is required that \[ \frac{1}{\left(1-\eta\right)^{2}}\left(\frac{d\eta}{ds}\left(q-q^{-}\right) + \left(1-\eta\right)\frac{dq}{ds}\right) \geq 0 \] A simpler, more general, semi-implicit solution is implemented in the code where \(\eta_{t + \Delta t}\) is taken to be a function of \(s_{t}\), hence a known value.
The different formulation available in the model are given in Table 2 with their parameters outlined in Table 3
| Type | Description | Parameters | \(\mathcal{F}\left(s_{sf},\Theta_{sf}\right)\) | \(\mathcal{F}_{\eta}\left(s_{sf},\Theta_{sf}\right)\) | \(q_{sf}^{+}\) |
|---|---|---|---|---|---|
| cnst | Linear Tank coupled below a constant parameter diffusive wave solution | \(t_{raf}, s_{raf},c_{sf}, d_{sf}\) | \[ \left\{ \begin{array}{cl} \frac{s_{sf}}{t_{raf}} & s_{sf} \leq s_{raf} \\ \frac{s_{raf}}{t_{raf}} + \frac{c_{sf}}{\Delta x}(s_{sf}-s_{raf}) & \mathrm{otherwise} \end{array} \right.\] | \(\frac{1}{2} - \frac{d_{sf}}{c_{sf}\Delta x}\) | \[ \left\{ \begin{array}{cl} \frac{s_{sf}}{t_{raf}} & s_{sf} \leq s_{raf} \\ \frac{s_{raf}}{t_{raf}} + \frac{ \frac{c_{sf}}{\Delta x}(s_{sf}-s_{raf}) - \eta_{sf}\max\left(0,q_{sf}^{-}- \frac{s_{raf}}{t_{raf}}\right)}{1-\eta_{sf}}& \mathrm{otherwise} \end{array} \right.\] |
| kin | Linear tank with Kinematic Wave for remaining flow | \(t_{raf}, s_{raf}, n, w_{sf}, g_{sf}\) | \[ \left\{ \begin{array}{cl} \frac{s_{sf}}{t_{raf}} & s_{sf} \leq s_{raf} \\ \frac{s_{raf}}{t_{raf}} + \frac{g_{sf}^{1/2}}{n w_{sf}^{2/3}} \left(\frac{s_{sf}-s_{raf}}{\Delta x}\right)^{5/3} & \mathrm{otherwise} \end{array} \right.\] | \(\frac{1}{2}\) | \[ \left\{ \begin{array}{cl} \frac{s_{sf}}{t_{raf}} & s_{sf} \leq s_{raf} \\ \frac{s_{raf}}{t_{raf}} + 2\frac{g_{sf}^{1/2}}{n w_{sf}^{2/3}} \left(\frac{s_{sf}-s_{raf}}{\Delta x}\right)^{5/3} - \max\left(0,q_{sf}^{-}- \frac{s_{raf}}{t_{raf}}\right) & \mathrm{otherwise} \end{array} \right.\] |
| Symbol | Description | unit |
|---|---|---|
| \(t_{raf}\) | Time constant of the runoff attenuation feature | s |
| \(c_{sf}\) | Free flowing surface water velocity | ms\(^{-1}\) |
| \(s_{raf}\) | storage of the runoff attenuation feature | m\(^3\) |
| \(d_{sf}\) | Free flowing surface water diffusion rate | m\(^2\)s\(^{-1}\) |
| \(w_{sf}\) | Width of surface channel | m |
| \(g_{sf}\) | Gradient of surface channel | - |
| \(n\) | Manning \(n\) coefficent for roughness | sm\(^{-1/3}\) |
Saturated Zone represetnations
The function \(\mathcal{G}\left(s_{sz},\Theta_{sz}\right)\) is related to the
HRU width through the transmissivity function. Table
4 gives \(\mathcal{G}\left(s_{sz},\Theta_{sz}\right)\) for various the transmissivity profiles present as options within the dynatop
package, the corresponding value to use for the type value
in the saturated zone specification and the additional parameters
required. The additional parameters are defined in Table 5.
| Name | \(G\left(\bar{s}_{sz},\Theta_{sz}\right)\) | \(\Theta_{sz}\) |
type value |
Notes | |
|---|---|---|---|---|---|
| Exponential | \(T_{0}w\sin\left(\beta\right)\exp\left(-\frac{\cos\beta}{Am}s_{sz}\right)\) | \(T_0,\beta,m\) | exp | Originally given in Beven & Freer 2001 | |
| Bounded Exponential | \(T_{0}w\sin\left(\beta\right)\left(\exp\left(-\frac{\cos\beta}{Am}s_{sz}\right) - \exp\left(-\frac{\cos\beta}{m}D\right)\right)\) | \(T_0,\beta,m,D\) | bexp | Originally given in Beven & Freer 2001 | |
| Constant Celerity | \(\frac{c_{sz}w}{A}\left(AD-s_{sz}\right)\) | \(c_{sz},D\) | cnst | ||
| Double Exponential | \(T_{0}w\sin\left(\beta\right)\left( \omega\exp\left(-\frac{\cos\beta}{Am}s_{sz}\right) + \left(1-\omega\right)\exp\left(\frac{\cos\beta}{Am_2}s_{sz}\right)\right)\) | \(T_0,\beta,m,m_2,\omega\) | dexp |
| Symbol | Description | unit |
|---|---|---|
| \(T_0\) | Transmissivity at saturation | m\(^2\)/s |
| \(m\),\(m_2\) | Exponential transmissivity constants | m\(^{-1}\) |
| \(\beta\) | Angle of hill slope | rad |
| \(c_{sz}\) | Saturated zone celerity | m/s |
| \(D\) | Storage deficit at which zero lateral flow occurs | m |
| \(\omega\) | weighting parameter | - |
Appendix A - Monotonicity of \(\mathcal{H}\left(z\right)\)
By definition \(\left. s_{uz} \right.\rvert_{t}\), \(\Delta t\), A, \(T_d\), \(\hat{v}_{rz \rightarrow uz}\) and \(z\) are all greater or equal to 0. To show that the gradient of
\[ \mathcal{H}\left(z\right) = z - \left. s_{sz}\right.\rvert_{t} + \Delta t \left. q_{sz}^{-} \right\rvert_{t+\Delta t} + A\Delta t \min\left( \frac{\left( \left. s_{uz} \right.\rvert_{t} + \hat{v}_{rz \rightarrow uz}\right)} {T_{d}z + A \Delta t}, \frac{1}{T_{d}} \right) - \\ \Delta t \min\left(q_{sz}^{max},\max\left( 0, 2 \mathcal{G}\left(z,\Theta_{sz}\right) - \left. q_{sz}^{-} \right\rvert_{t+\Delta t} \right)\right) \]
is strictly positive consider first the gradient of the last term which takes the value
\[ \frac{d}{dz} \min\left(q_{sz}^{max},\max\left( 0, 2 \mathcal{G}\left(z,\Theta_{sz}\right) - \left. q_{sz}^{-} \right\rvert_{t+\Delta t} \right)\right) = \left\{ \begin{array}{cl} 0 & 2 \mathcal{G}\left(z,\Theta_{sz}\right) < \left. q_{sz}^{-} \right\rvert_{t+\Delta t}\\ 0 & 2 \mathcal{G}\left(z,\Theta_{sz}\right) - \left. q_{sz}^{-} \right\rvert_{t+\Delta t}> q_{sz}^{max} \\ 2 \frac{d}{dz} \mathcal{G}\left(z,\Theta_{sz}\right) & \mathrm{otherwise} \end{array} \right. \]
By definition \(-\frac{d}{dz} \mathcal{G}\left(z,\Theta_{sz}\right) \geq 0\) hence
\[ \frac{d}{dz} \mathcal{H}\left(z\right) \geq 1 + A\Delta t \frac{d}{dz} \min\left( \frac{\left( \left. s_{uz} \right.\rvert_{t} + \hat{v}_{rz \rightarrow uz}\right)} {T_{d}z + A \Delta t}, \frac{1}{T_{d}} \right) \]
Separating the range of \(z\) into two sections gives
\[ \frac{d}{dz} \mathcal{H}\left(z\right) \geq \left\{ \begin{array}{cl} 1 & \frac{\left( \left. s_{uz} \right.\rvert_{t} + \hat{v}_{rz \rightarrow uz}\right)} {T_{d}z + A \Delta t} \geq \frac{1}{T_{d}} \\ 1 - T_{d} A\Delta t \frac{\left. s_{uz} \right.\rvert_{t} + \hat{v}_{rz \rightarrow uz}} {\left(T_{d}z + A \Delta t\right)^2} & frac{\left( \left. s_{uz} \right.\rvert_{t} + \hat{v}_{rz \rightarrow uz}\right)}{T_{d}z + A \Delta t} \leq \frac{1}{T_{d}} \end{array} \right. \]
Further using the range of \(z\) valid for the second case
\[ \frac{d}{dz} \mathcal{H}\left(z\right) \geq \left\{ \begin{array}{cl} 1 & \frac{\left( \left. s_{uz} \right.\rvert_{t} + \hat{v}_{rz \rightarrow uz}\right)} {T_{d}z + A \Delta t} \geq \frac{1}{T_{d}} \\ 1 - \frac{A\Delta t}{T_{d}z + A \Delta t} & \frac{\left( \left. s_{uz} \right.\rvert_{t} + \hat{v}_{rz \rightarrow uz}\right)}{T_{d}z + A \Delta t} \leq \frac{1}{T_{d}} \end{array} \right. \]
From this it can be seen that \(\frac{d}{dz} \mathcal{H}\left(z\right) \geq 0\) with equality possible when \(z=0\).