The HRU
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.
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. Appendices provide supporting information.
Notation
Table @ref(tab:hs-x-notation) 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⋆ → ⋆ 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⋆ and q⋆ are considered to be positive flows travelling downslope.
The superscripts + and − are used to denote variables for the outflow and inflow to the hillslope.
y|t indicates the value of the variable y at time t
| Quantity type | Symbol | Description | unit |
|---|---|---|---|
| Storage | ssf | Surface excess storage | m3 |
| srz | Root zone storage | m3 | |
| suz | Unsaturated zone storage | m3 | |
| ssz | Saturated zone storage deficit | m3 | |
| Vertical fluxes | rsf → rz | Flow from the surface excess store to the root zone | m3/s |
| rrz → uz | Flow from the root zone to the unsaturated zone | m3/s | |
| ruz → sz | Flow from unsaturated to saturated zone | m3/s | |
| Lateral fluxes | qsf | Lateral flow in the hillslope surface zone | m3/s |
| qsz | Lateral inflow to the hillslope surface zone | m3/s | |
| Vertical fluxes to the HRU | p | Precipitation rate | m3/s |
| ep | Potential Evapotranspiration rate | m3/s | |
| Other HRU properties | A | Plan area | m2 |
| w | Width of a hill slope cross section | m | |
| Δx | Effective length of the hillslope HRU (x = A/w$) | m | |
| Θsf | Further properties and parameters for the solution of the surface zone | :-: | |
| Θrz | Further properties and parameters for the solution of the root zone | :-: | |
| Θuz | Further properties and parameters for the solution of the saturated zone | :-: | |
| Θ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 + Δ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 rsf → 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 rsf → rz.
The outflow qsf+ is considered a function of qsf−, ssf, rsf → rz and parameters Θsf denoted by qsf+ = ℱ(ssf, qsf−, rsf → rz, Θsf) For the purposes the numerical scheme we limit the set to possible functions ℱ to those that are non-decreasing with respect to ssf and take the value 0 when ssf = 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 ℱ(0, qsf−, rsf → rz, Θsf) = 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
with the governing ODE
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 srz ≤ srzmax then rrz → uz ≤ 0. Negative values of rrz → uz may occur only when water is returned from the unsaturated zone due to saturation caused by lateral inflow to the saturated zone.
When srz = srzmax then In this case rrz → uz may be positive if so long as the unsaturated zone can receive the water. If rrz → uz is ‘throttled’ by the rate at which the unsaturated zone can receive water, then rsf → 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
The governing ODE is written as
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 ssz = suz the subsurface below the root zone can be considered saturated, as in there is no further available storage for water, but separated into two 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 ruz → sz occurs at the maximum possible rate and rrz → uz is limited to ensure that suz ≤ ssz.
Saturated Zone
For a HRU the alteration of the storage in the saturated zone is given in terms of the storage deficit ssz as: $$ \frac{ds_{sz}}{dt} = q^{+}_{sz} - r_{sf \rightarrow rz} - q^{-}_{sz} $$
The relationship between the outflow and storage deficit is given qsz+ = 𝒢(ssz, qsz−, Θsz) 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 𝒢(D, qsz−, Θsz) = 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 Δt. Hence on the first downward pass through the stores the maximum downward volumes of the vertical fluxes v̂* → * = Δtr̂* → * are determined. These then moderated by the solution of the saturated zone, with an upward pass giving the final states.
Downward Pass
Surface excess
The implicit solution of mass balance equation gives
ssf|t + Δt = ssf|t + Δt(qsf−|t + Δt − ℱ(ssf, qsf−, rsf → rz, Θsf)|t + Δt) − vsf → rz|t + Δt
Since the outflow is 0 when the ssf = 0 the maximum downward flux volume is v̂sf → rz|t + Δt = ssf|t + Δtqsf−|t + Δt
Root Zone
The implicit formulation for the root zone gives
Since flow to the unsaturated zone occur only to keep srz ≤ srzmax then
Unsaturated Zone
The implicit expression for the unsaturated zone is given in terms of suz|t + Δt and ssz|t + Δt as
Since suz|t + Δt ≤ ssz|t + Δt this places an additional condition of
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 ssz|t + Δt.
Saturated Zone
Building on the generic solution the implicit scheme for the saturated zone gives
ssz|t + Δt = ssz|t + Δt(𝒢(ssz, qsz−, Θsz)|t + Δt − qsz−|t + Δt) − vuz → sz|t + Δ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 ℋ(z) is a monotonic increasing function of z > 0 and ℋ(D) ≥ 0. Hence if ℋ(0) < 0 unique solution for ssz|Δt can be found. If ℋ(0) ≥ 0 then the incoming fluxes are enough to produce saturation in the subsurface and ssz|Δt = 0.
Maintaining a water (mass) balance and correct limits on the state values depends upon the accuracy of the solution of ℋ(z) = 0. Since each evaluation of ℋ(z) requires an evaluation of 𝒢 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 ℋ(z) = 0 allows scope for compromises based on computational time.
The upward pass presented below is mass conservative so long as ℋ(ssz|Δt) ≥ 0. We therefore seek z which satifies for some tolerance ϵ > 0 ϵ ≥ ℋ(z) ≥ 0 One approach to finding such a value is use a bracketing algorithm. At each iteration the lower and upper limits zl < zu satisify ℋ(zl) ≤ 0 ≤ ℋ(zu), thus bracket the solution. Iterating (using for example bisection) until ℋ(zl) ≤ ϵ give the approximation ssz|Δt = zu.
Upward Pass
Given an estimate of ssz|Δt such that ℋ(ssz|Δt) > 0 compute the vertical flux volumes in the following fashion
qsz+|t + Δt = min (qszmax, max (0, 2𝒢(ssz|t + Δt, Θsz) − qsz−|t + Δt))
vuz → sz = ssz|t + Δt(qsz+|t + Δt − qsz−|t + Δt) − ssz|t + Δt
suz|t + Δt = min (ssz|t + Δt, suz|t + v̂rz → uz − vuz → sz)
vrz → uz = suz|t + Δt − suz|t + vuz → sz
vsf → rz = min (v̂sf → rz, srzmax − s̄rz|0 − Δt(p − ep) + vrz → uz)
Solve for srz|t + Δt using the equation already given.
Compute the volume loss to evapotranspiration as vep|t + Δt = srz|t + vsf → rz − vrz → uz + Δtp − srz|t + Δ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 ≥ 0 with 𝒮(0) ≥ 0.
A bracketing solution algorithm gives at any iteration wl < ssf|t + Δt < wu. An initial value of wu can be computed as ssf|t + Δtqsf−|t + Δt − vsf → rz|t + Δt. Iterating until; for some tolerance ϵ > 0; wu − wl < ϵ take ssf|Δt = wl. In this case 𝒮(ssf|Δt) > 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
Recall that 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 rsf → 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).
The following subsections present various options for the surface zone solution based on the Muskingham approximation method.
The Muskingham method offers a parsimonious for capturing a range of surface dynamics. This 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 ssf in a reach or hillslope of length Δx can be expressed in terms of the representative flow qsf and velocity vsf as $$ s_{sf} = \frac{q_{sf} \Delta x}{v_{sf}} $$ The representative flow is related to the inflow qsf− and outflow qsf+ using the parameter η through qsf = ηqsf− + (1 − η)qsf+ which 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 substitution of $$ q_{sf}^{+} = \max\left(0, \frac{1}{1-\eta}\left(\frac{sv}{\Delta x} - \eta q^{-} \right) \right) $$ gives $$ \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.
Approximating Tank Models
Taking η = 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 = ηq− + (1 − η)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 η is limited to between 0 and $\frac{1}{2}$. Firstly η = 1/2 results from D = 0; that is the Kinematic wave equation. Taking η = 0 produces the maximum dispersion of $D = \frac{c\Delta x}{2}$, which is given by the linear tank.
Compound Channels
Compound channels are considered as channels with two regimes corresponding to different reltionships between the ssf, qsf and η depending upon whether ssf ≤ sc, for some storage threshold s1. This allows a basic representation of surface storage such as runoff attenuation features, or out of bank channel flow.
In solving this inflow is partitioned to ensure, where possible, ssf > sc
Available Formulations
The different formulation available in the model are given in Table @ref(tab:surface) with their parameters outlined in Table @ref(tab:surface-param)
| Type | Description | Parameters | ℱ(ssf, qsf−, rsf → rz, Θsf) |
|---|---|---|---|
| cnst | Linear Tank coupled below a constant parameter diffusive wave solutions | traf, sraf, csf, dsf | |
| kin | Linear tank with Kinematic Wave for remaining flow | traf, sraf, n, wsf, gsf | |
| comp | Two constant parameter diffusive wave solutions | vsf, 1, dsf, 1, s1, vsf, 2, dsf, 2 |
| Symbol | Description | unit |
|---|---|---|
| traf | Time constant of the runoff attenuation feature | s |
| csf | Free flowing surface water velocity | ms−1 |
| sraf | storage of the runoff attenuation feature | m3 |
| dsf | Free flowing surface water diffusion rate | m2s−1 |
| wsf | Width of surface channel | m |
| gsf | Gradient of surface channel | - |
| n | Manning n coefficent for roughness | sm−1/3 |
Saturated Zone representations
The function $q_{sz$ =
(s_{sz},_{sz})$ dscribing the representative lateral flow is related to
the HRU width through the transmissivity function. Table
@ref(tab:hs-transmissivity) gives 𝒢(ssz, Θsz)
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 @ref(tab:hs-transmissivity-param).
Since the subsurface has a kinematic approximation $$ q_{sz}^{+} = \max\left(0, 2\left( q_{sz} - \frac{1}{2} q_{sf}^{-}\right)\right) $$
| Name | G(s̄sz, Θsz) | Θsz | type value |
Notes | |
|---|---|---|---|---|---|
| Exponential | $T_{0}w\sin\left(\beta\right)\exp\left(-\frac{\cos\beta}{Am}s_{sz}\right)$ | T0, β, 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)$ | T0, β, m, D | bexp | Originally given in Beven & Freer 2001 | |
| Constant Celerity | $\frac{c_{sz}w}{A}\left(AD-s_{sz}\right)$ | csz, 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)$ | T0, β, m, m2, ω | dexp |
| Symbol | Description | unit |
|---|---|---|
| T0 | Transmissivity at saturation | m2/s |
| m,m2 | Exponential transmissivity constants | m−1 |
| β | Angle of hill slope | rad |
| csz | Saturated zone celerity | m/s |
| D | Storage deficit at which zero lateral flow occurs | m |
| ω | weighting parameter | - |
Appendix A - Monotonicity of ℋ(z)
By definition suz|t, Δt, A, Td, v̂rz → 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.