The Hillslope 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 hillslope HRU. The hillslope HRU is a representation of
an area of hill slope with, for example, similar topographic, soil and upslope area
characteristics.
The hillslope HRU is considered to be a area of hill slope 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.
In the following section the governing equations of the hillslope 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 a wide range of transmissivity profiles,
including those present in the code and presented in the third sections of the
vignette.
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 fluxes \(r_{\star \rightarrow \star}\) are considered to be positive when
travelling in the direction of the arrow, but can, in some cases be negative.
All fluxes \(l_{\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\)
Spatial averages are denoted with a \(\bar{}\) so for example integrating
over the length of the hillslope we find
\[\begin{equation}
\int w y dx = A\bar{y}
\end{equation}\]
- Average values over a time step are denoted with a \(\hat{}\). If the
temporal average is of a spatial average the \(\bar{}\) is dropped. For
example the temporal average of the spatial average \(\bar{y}\) over a time step is \(\hat{y}\).
Table 1: Outline of notation for describing a cross
section of the hillslope HRU
| Storage |
\(s_{sf}\) |
Surface excess storage per unit width |
m |
|
\(s_{rz}\) |
Root zone storage per unit width |
m |
|
\(s_{uz}\) |
Unsaturated zone storage per unit width |
m |
|
\(s_{sz}\) |
Saturated zone storage deficit per unit width |
m |
| Vertical fluxes |
\(r_{sf \rightarrow rz}\) |
Flow from the surface excess store to the root zone per unit width |
m/s |
|
\(r_{rz \rightarrow uz}\) |
Flow from the root zone to the unsaturated zone per unit width |
m/s |
|
\(r_{uz \rightarrow sz}\) |
Flow from unsaturated to saturated zone per unit width |
m/s |
| Lateral fluxes |
\(l_{sf}\) |
Lateral flow in the surface zone per unit outflow width |
m\(^2\)/s |
|
\(l_{sz}\) |
Lateral flow in the saturated zone per unit outflow width |
m\(^2\)/s |
|
\(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/s |
|
\(e_p\) |
Potential Evapotranspiration rate |
m/s |
| Other HRU properties |
\(A\) |
Plan area |
m\(^3\) |
|
\(w\) |
Width of a hill slope cross section |
m |
|
\(\Delta x\) |
Effective length of the hillslope HRU (Area divided by outflow width) |
m |
|
\(k_{sf}\) |
Maximum flow rate from surface excess to root zone per unit width |
m/s |
|
\(c_{sf}\) |
Surface excess storage celerity |
m/s |
|
\(s_{raf}\) |
Surface runoff attenuation feature storage |
m |
|
\(t_{raf}\) |
Surface runoff attenuation feature time constant |
s |
|
\(s_{rzmax}\) |
Maximum root zone storage depth |
m |
|
\(T_d\) |
Unsaturated zone time constant per m of saturated storage deficit. |
s/m |
|
\(\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.
Surface zone
The storage in the surface zone satisfies \(s_{sf} \geq 0\).
Surface storage is increased by lateral downslope flow from upslope HRUs.
Water flows to the root zone at a constant rate \(k_{sf}\), 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 governing ODE for a cross section of the hillslope can be written as
\[\begin{equation}
\frac{dws_{sf}}{dt} = - \frac{dq_{sf}}{dx} - w r_{sf \rightarrow rz}
\end{equation}\]
where \(r_{sf \rightarrow rz} \leq k_{sf}\).
Integration over the hillslope length and substituting in the spatial averages gives
\[\begin{equation}
A\frac{d \bar{s}_{sf}}{dt} = q_{sf}^{-} - q_{sf}^{+} - A\bar{r}_{sf \rightarrow rz}
\end{equation}\]
The outflow comes from two sources; the runoff-attenuation feature
(RAF) and surface flow; giving
\[
q_{sf}^{+} = w^{+}l_{sf}^{+} + \frac{A}{t_{raf}}\min\left(s_{raf},\bar{s}_{sf}\right)
\]
Since \(\Delta x = A/w^{+}\) and \(l_{sf}^{-} = q_{sf}^{-}/w^{+}\) re arranging gives
\[\begin{equation}
\frac{d \bar{s}_{sf}}{dt} = \frac{1}{\Delta x}\left(l_{sf}^{-} -
l_{sf}^{+}\right) - \frac{1}{t_{raf}}
\min\left(s_{raf},\bar{s}_{sf}\right) - \bar{r}_{sf \rightarrow rz}
\end{equation}\]
This equation represents the evolution in time of \(\bar{s}_{sf}\) expressed in
terms of the unit width fluxes, the runoff attenuation feature outflow and an effective length \(A/w^+\)
To complete the formulation assume that
\(l_{sf}^{+} = c_{sf}\max\left(0,\bar{s}_{sf} - s_{raf}\right)\).
which gives
\[\begin{equation}
\frac{d \bar{s}_{sf}}{dt} = \frac{1}{\Delta x}\left(l_{sf}^{-} -
c_{sf}\max\left(0,\bar{s}_{sf}-s_{raf}\right)\right)
- \frac{1}{t_{raf}} \min\left(s_{raf},\bar{s}_{sf}\right)
- \bar{r}_{sf \rightarrow rz}
\end{equation}\]
Root zone
The root zone gains water from precipitation and the surface zone. It loses water through
evaporation and to the unsaturated zone. Since
all the vertical fluxes are spatially uniform the evolution can be evaluated
in terms of the spatial averages. The root zone storage satisfies
\[\begin{equation} 0 \leq \bar{s}_{rz} \leq s_{rzmax} \end{equation}\]
with the governing ODE
\[\begin{equation}
\frac{d\bar{s}_{rz}}{dt} = p - \frac{e_p}{s_{rzmax}} \bar{s}_{rz} +
\bar{r}_{sf\rightarrow rz} - \bar{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 \(\bar{s}_{rz} \leq s_{rzmax}\) then \(\bar{r}_{rz \rightarrow uz} \leq 0\). Negative
values of \(\bar{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 \(\bar{s}_{rz} = s_{rzmax}\) then
\[\begin{equation}
p - e_p + \bar{r}_{sf\rightarrow rz} - \bar{r}_{rz \rightarrow uz} \leq 0
\end{equation}\]
In this case \(\bar{r}_{rz \rightarrow uz}\) may be positive if
\[\begin{equation}
p - e_p + \bar{r}_{sf\rightarrow rz} > 0
\end{equation}\]
so long as the unsaturated zone can receive the water. If \(\bar{r}_{rz \rightarrow uz}\) is ‘throttled’ by the rate at which the unsaturated zone can
receive water, then \(\bar{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 \bar{s}_{uz} \leq \bar{s}_{sz} \end{equation}\]
The governing ODE is written as
\[\begin{equation}
\frac{d\bar{s}_{uz}}{dt} = \bar{r}_{rz \rightarrow uz} - \bar{r}_{uz \rightarrow sz}
\end{equation}\]
If water is able to pass freely to the saturated zone, then it flows at the rate
\(\frac{\bar{s}_{uz}}{T_d \bar{s}_{sz}}\). If \(\bar{s}_{sz}=\bar{s}_{uz}=0\) this is interpreted to mean
that the flow rate is \(\frac{1}{T_d}\).
In this situation where \(\bar{s}_{sz}=\bar{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 \(\bar{r}_{uz \rightarrow sz}\) is
constrained by the ability of the saturated zone to receive water. If
this is the case \(\bar{r}_{uz \rightarrow sz}\) occurs at the maximum
possible rate and \(\bar{r}_{rz \rightarrow uz}\) is limited to ensure that
\(\bar{s}_{uz} \leq \bar{s}_{sz}\).
Saturated Zone
For a HRU the alteration of the storage in the saturated zone is
given by the kinematic equation
\[\begin{equation}
\frac{dw\left(D-s_{sz}\right)}{dt} = -\frac{dq_{sz}}{dx} + w r_{uz \rightarrow sz}
\end{equation}\]
Integration over the length of the hillslope and substitution of the
spatial averaged values gives
\[\begin{equation}
-A\frac{d\bar{s}_{sz}}{dt} = q_{sz}^{-} - q_{sz}^{+} + A\bar{r}_{uz \rightarrow sz}
\end{equation}\]
Rearranging this gives
\[\begin{equation}
\frac{d\bar{s}_{sz}}{dt} = \frac{1}{\Delta x}\left(l_{sz}^{+} - l_{sz}^{-}\right) - \bar{r}_{uz \rightarrow sz}
\end{equation}\]
The Kinematic approximation then depends upon the specification of a relationship
between \(\bar{s}_{sz}\) and \(l_{sz}^{+}\). This usually termed the transmissivity
profile however the use of the average, rather than cross sectional saturated
storage deficit means this term is used only loosely. More generally the
relationship is taken
to be a one-to-one, continuously differentiable function
\(g: \mathcal{R^{+}}\rightarrow \mathcal{R}^{+}\) which returns the lateral flow on
a unit width such that \(l_{sz}^{+} = g\left(\bar{s}_{sz},\Theta_{sz}\right) \geq 0\) which satisfies
\[\begin{equation}
-\frac{dl_{sz}^{+}}{d\bar{s}_{sz}} = -\frac{d}{d\bar{s}_{sz}}g\left(\bar{s}_{sz},\Theta_{sz}\right) \geq 0
\end{equation}\]
and
\[\begin{equation}
g\left(\bar{s}_{sz},\Theta_{sz}\right) \rightarrow 0 \quad \mathrm{as} \quad
\bar{s}_{sz}\rightarrow \infty
\end{equation}\]
The first condition ensures a non negative celerity (wave speed) which is in
keeping with the conceptualisation of the model. The first condition also implies
that \(l_{sz}^{+}\) decreases with increasing \(\bar{s}_{sz}\). This combined with
the second condition ensures that the model represents a state where no
lateral inflow is generated.
Using the relationship between \(l_{sz}^{+}\) and \(\bar{s}_{sz}\) gives
\[\begin{equation}
\frac{d\bar{s}_{sz}}{dt} = \frac{1}{\Delta x}\left(g\left(\bar{s}_{sz},\Theta_{sz}\right) - l_{sz}^{-}\right) - \bar{r}_{uz \rightarrow sz}
\end{equation}\]
Approximating Equations
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\).
Recalling that the temporal and spatial average vertical fluxes are denoted with
a \(\hat{}\) the implicit approximations are stated in the following sections using the unknown vertical fluxes
\(\hat{r}_{\star \rightarrow \star}\) which are then determined by maximising the
downward flux.
Surface excess
The implicit scheme gives
\[\begin{equation}
\left. \bar{s}_{sf} \right\rvert_{\Delta t}
=
\left. \bar{s}_{sf} \right\rvert_{0}
+ \frac{\Delta t}{\Delta x} \left(
\hat{l}_{sf}^{-} - c_{sf} \max\left(0, \left. \bar{s}_{sf}
\right\rvert_{\Delta t} - s_{raf}\right)
\right)
- \frac{\Delta t}{t_{raf}} \min\left(s_{raf},\left. \bar{s}_{sf}
\right\rvert_{\Delta t}\right)
- \Delta t \hat{r}_{sf \rightarrow rz}
\end{equation}\]
To ensure that the surface storage does not become negative
\[\begin{equation}
\hat{r}_{sf \rightarrow rz} \leq \min\left( k_{sf},
\left. \frac{1}{\Delta t}\left(\bar{s}_{sf} \right\rvert_{0} + \frac{\Delta
t}{\Delta x} \hat{l}_{sf}^{-}
\right)\right)
\end{equation}\]
The surface outflow over the time step is evaluated as
\[\begin{equation}
\hat{q}_{sf}^{+} = w^{+}c_{sf} \max\left(0, \left. \bar{s}_{sf}
\right\rvert_{\Delta t} -s_{raf}\right) +
\frac{A}{t_{raf}} \min\left(s_{raf},\left. \bar{s}_{sf}
\right\rvert_{\Delta t}\right)
\end{equation}\]
The solution for \(\left. \bar{s}_{sf} \right\rvert_{\Delta t}\) is conditional upon which side of \(s_{raf}\)
it falls.
If
\[
s_{raf} \lt
\left. \bar{s}_{sf} \right\rvert_{0}
+ \frac{\Delta t}{\Delta x} \hat{l}_{sf}^{-}
- \frac{\Delta t}{t_{raf}} s_{raf}
- \Delta t \hat{r}_{sf \rightarrow rz}
\]
then \(\left. \bar{s}_{sf} \right\rvert_{\Delta t} \gt s_{raf}\) and is given by
\[\begin{equation}
\left. \bar{s}_{sf} \right\rvert_{\Delta t}
= \left(1 + \frac{c_{sf}\Delta t}{\Delta x}\right)^{-1}
\left( \left. \bar{s}_{sf} \right\rvert_{0}
+ \frac{\Delta t}{\Delta x} \hat{l}_{sf}^{-}
+ \frac{c_{sf}\Delta t}{\Delta x} s_{raf}
- \frac{\Delta t}{t_{raf}} s_{raf}
- \Delta t \hat{r}_{sf \rightarrow rz}
\right)
\end{equation}\]
Otherwise
\[\begin{equation}
\left. \bar{s}_{sf} \right\rvert_{\Delta t}
= \left(1 + \frac{\Delta t}{t_{raf}}\right)^{-1}
\left( \left. \bar{s}_{sf} \right\rvert_{0}
+ \frac{\Delta t}{\Delta x} \hat{l}_{sf}^{-}
- \Delta t \hat{r}_{sf \rightarrow rz}
\right)
\end{equation}\]
To ensure that the surface storage does not become negative
\[\begin{equation}
\hat{r}_{sf \rightarrow rz} \leq \min\left( k_{sf},
\left. \frac{1}{\Delta t}\left(\bar{s}_{sf} \right\rvert_{0} + \frac{\Delta
t}{\Delta x} \hat{l}_{sf}^{-}
\right)\right)
\end{equation}\]
The surface outflow over the time step is evaluated as
\[\begin{equation}
\hat{q}_{sf}^{+} = w^{+}c_{sf} \left. \bar{s}_{sf} \right\rvert_{\Delta t}
\end{equation}\]
Root Zone
The implicit formulation for the root zone gives
\[\begin{equation}
\left. \bar{s}_{rz} \right.\rvert_{\Delta t} =
\left( 1 + \frac{e_{p}\Delta t}{s_{rzmax}} \right)^{-1}
\left(
\left. \bar{s}_{rz} \right.\rvert_{0}
+ \Delta t \left(
p + \hat{r}_{sf \rightarrow rz} - \hat{r}_{rz \rightarrow uz} \right)
\right)
\end{equation}\]
Since flow to the unsaturated zone occur only to keep
\(\bar{s}_{rz} \leq s_{rzmax}\) then
\[\begin{equation}
\hat{r}_{rz \rightarrow uz} \leq \max\left(0,
\frac{1}{\Delta t}\left(
\left. \bar{s}_{rz} \right.\rvert_{0}
+ \Delta t \left(
p + \hat{r}_{sf \rightarrow rz} - e_{p} \right)
- s_{rzmax}
\right)\right)
\end{equation}\]
Unsaturated Zone
The expression for the unsaturated zone is given in terms of
\(\left. \bar{s}_{uz} \right.\rvert_{\Delta t}\) and
\(\left. \bar{s}_{sz} \right.\rvert_{\Delta t}\) as
\[\begin{equation}
\left. \bar{s}_{uz} \right.\rvert_{\Delta t} =
\frac{T_{d}\left. \bar{s}_{sz} \right.\rvert_{\Delta t}}{
T_{d}\left. \bar{s}_{sz} \right.\rvert_{\Delta t} + \Delta t}
\left( \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \hat{r}_{rz
\rightarrow uz}\right)
\end{equation}\]
Since \(\left. \bar{s}_{uz} \right.\rvert_{\Delta t} \leq \left. \bar{s}_{sz} \right.\rvert_{\Delta t}\) this places an additional condition of
\[\begin{equation}
\hat{r}_{rz \rightarrow uz} \leq \frac{1}{\Delta t} \left(
\left. \bar{s}_{sz} \right.\rvert_{\Delta t} + \frac{\Delta t}{T_d} - \left. \bar{s}_{uz} \right.\rvert_{0}\right)
\end{equation}\]
This condition also implies that \(\hat{r}_{uz \rightarrow sz} \leq \frac{1}{T_d}\)
Saturated Zone
Initial evaluation and substitution of the flux from the unsaturated zone gives
\[\begin{equation}
\left. \bar{s}_{sz}\right.\rvert_{\Delta t} =
\left. \bar{s}_{sz}\right.\rvert_{\Delta 0} +
\frac{\Delta t}{\Delta x} \left(g\left(\left.\bar{s}_{sz}\right\rvert_{\Delta t},\Theta_{sz}\right) -
l_{sz}^{-}\right) - \Delta t \frac{
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \hat{r}_{rz \rightarrow uz} }
{ T_{d}\left. \bar{s}_{sz} \right.\rvert_{\Delta t} + \Delta t}
\end{equation}\]
Limiting the inflow by reaching the saturation level gives
\[\begin{equation}
\hat{r}_{rz \rightarrow uz} \leq \frac{1}{\Delta t}\left(
\left. \bar{s}_{sz}\right.\rvert_{\Delta 0} +
\frac{\Delta t}{\Delta x} \left(g\left(\left.\bar{s}_{sz}\right\rvert_{\Delta t},\Theta_{sz}\right) -
l_{sz}^{-}\right) - \left. \bar{s}_{uz} \right.\rvert_{0} \right)
\end{equation}\]
Selecting \(\bar{s}_{sz}\) and \(\hat{r}_{rz \rightarrow uz}\) to satisfy the above equation describes a non-linear
problem with two unknowns a multiple solutions. A further condition, that the
downward flux is maximised, is imposed to ensure a unique solution.
Numerical Algorithm
An implementation of the approximating equations which is consistent with
maximising the downward flux produces a fully implicit numeric scheme.
The maximum downward flux from the surface is given by
\[
\check{r}_{sf \rightarrow rz} = \min\left( k_{sf},
\left. \frac{1}{\Delta t}\left(\bar{s}_{sf} \right\rvert_{0} + \frac{\Delta
t}{\Delta x} \hat{l}_{sf}^{-}
\right)\right)
\]
which results in a maximum downward flux from the root zone of
\[
\check{r}_{rz \rightarrow uz} = \max\left(0, \frac{1}{\Delta t}\left(
\left. \bar{s}_{rz} \right.\rvert_{0} + \Delta t \left( p + \check{r}_{sf
\rightarrow rz} - e_{p} \right) - s_{rzmax} \right)\right)
\]
Consider the pseudo-function
\[
f\left(z\right) =
z -
\left. \bar{s}_{sz}\right.\rvert_{0} -
\frac{\Delta t}{\Delta x} \left(g\left(z,\Theta_{sz}\right) -
l_{sz}^{-}\right) + \Delta t \min\left(\frac{1}{T_d},\frac{
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} }
{ T_{d} z + \Delta t}\right)
\]
adapted from the approximating equation for the saturated zone storage with the minima
ensuring the limit on the flux from the unsaturated zone. If \(f(0) \geq 0\) the
downward flux is enough to produce saturation in the subsurface and
\(\left. \bar{s}_{sz}\right.\rvert_{\Delta t} = 0\). If \(f(0) < 0\) then
\(\left. \bar{s}_{sz}\right.\rvert_{\Delta t}\) is given by a solution of
\(f(z) = 0\). In Appendix A it is
shown that \(f\left(z\right)\) is a monotonic increasing function of
\(z > 0\); hence if \(f\left(0\right) \leq 0\) and \(f\left(D\right) \geq 0\) a
unique solution for \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t}\) can be found.
A direct implementation of the above results takes the form of Algorithm 1.
Algorithm 1: Direct Implementation
- Compute \(\check{r}_{sf \rightarrow rz}\) and \(\check{r}_{rz \rightarrow uz}\)
- Test for saturation. If \(f\left(0\right) \geq 0\) then go to Step 3 else go to Step 4
- Saturated. Set \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t} = 0\) and
\[\hat{r}_{rz \rightarrow uz} = \frac{1}{\Delta t}\left(
\left. \bar{s}_{sz}\right.\rvert_{\Delta 0} +
\frac{\Delta t}{\Delta x} \left(g\left(0,\Theta_{sz}\right) -
l_{sz}^{-}\right) - \left. \bar{s}_{uz} \right.\rvert_{0} \right)
\]
Go to Step 5.
- Not saturated. Find \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t}\) such that
\(f\left(\left. \bar{s}_{sz}\right.\rvert_{\Delta t}\right) = 0\) and take
\[
\hat{r}_{rz \rightarrow uz} = \min\left(\check{r}_{rz \rightarrow uz},\frac{1}{\Delta t} \left(
\left. \bar{s}_{sz} \right.\rvert_{\Delta t} + \frac{\Delta t}{T_d} -
\left. \bar{s}_{uz} \right.\rvert_{0}\right)\right)\]
Go to Step 5.
- Compute
\[
\hat{r}_{sf \rightarrow rz} = \min\left(\check{r}_{sf \rightarrow
rz}, \frac{1}{\Delta t}\left( s_{rzmax} -
\left. \bar{s}_{rz} \right.\rvert_{0} - \Delta t \left( p - e_p - \hat{r}_{rz
\rightarrow uz} \right) \right)\right)\]
- Solve for \(\left. \bar{s}_{uz}\right.\rvert_{\Delta t}\),
\(\left. \bar{s}_{rz}\right.\rvert_{\Delta t}\) and
\(\left. \bar{s}_{sf}\right.\rvert_{\Delta t}\) using the computed vertical
fluxes.
The challenge with Algorithm 1 is that maintaining a water (mass) balance and
correct limits on the state values
depends upon the accuracy of the solution of \(f\left(z\right)=0\). Since
finding \(z\) requires multiple evaluations of the transmissivity profile it can
rapidly become the most expensive part of the code.
To limit the computational cost suppose that \(f\left(z\right)=0\) is solved
approximately to give \(\tilde{z}\)
such that for some numerical tolerance \(\epsilon\) we have \(\left\lvert f\left(\tilde{z}\right) \right\rvert \leq \epsilon\). Using \(\tilde{z}\) to
evaluate the outward lateral flux from the saturated zone gives the
pseudo-function
\[
h\left(z\right) =
z -
\left. \bar{s}_{sz}\right.\rvert_{0} -
\frac{\Delta t}{\Delta x} \left(g\left(\tilde{z},\Theta_{sz}\right) -
l_{sz}^{-}\right) + \Delta t \min\left(\frac{1}{T_d},\frac{
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} }
{ T_{d} z + \Delta t}\right)
\]
Since \(h\left(z\right)\) is a monotonic increasing function of
\(z > 0\) the solution for \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t}\) is
similar to that using \(f\left(z\right)\). Taking \(z \in \left[0,D\right]\) the
solution for \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t}\) is given by
- If \(h(0) \geq 0\) the inward flux is enough to produce saturation in the subsurface and
\(\left. \bar{s}_{sz}\right.\rvert_{\Delta t} = 0\).
- If \(h(z_{max}) < 0\) (possible only due to the outflow flux approximation) then
\(\left. \bar{s}_{sz}\right.\rvert_{\Delta t} = D\) and the lateral outflow is
adjusted such that
\[
\frac{\Delta t}{\Delta x} l_{sz}^{+} =
D -
\left. \bar{s}_{sz}\right.\rvert_{0} +
\frac{\Delta t}{\Delta x} l_{sz}^{-} +
\Delta t \min\left(\frac{1}{T_d},\frac{
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} }
{ T_{d} D + \Delta t}\right)
\]
- If \(h(0) < 0\) and If \(h(D) \geq 0\) then with \(z_{C} = \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t\left( \check{r}_{rz \rightarrow uz} - \frac{1}{T_{d}}\right)\)
- If \(h\left( z_{C} \right) \geq 0\) then the inflow from the unsaturated
zone is \(1/T_{d}\) and
\[
\left. \bar{s}_{sz}\right.\rvert_{\Delta t} =
\left. \bar{s}_{sz}\right.\rvert_{0} +
\frac{\Delta t}{\Delta x} \left(g\left(\tilde{z},\Theta_{sz}\right) -
l_{sz}^{-}\right) - \frac{\Delta t}{T_d}
\]
- else (see Appendix C)
\[
z = \frac{1}{2}\left(-b + \sqrt{ b^2 - 4c }\right) \\
= \frac{1}{2} \left( z_{C} - h\left(z_{C}\right) + \sqrt{
\left(h\left(z_{C}\right)-z_{C}\right)^2 - 4 \frac{\Delta t}{T_{d}}
h\left(z_{C}\right) } \right)
\]
The implementation of this approach is presented in Algorithm 2. it can be seem
that the exchange for decreasing the accuracy of the solution of the zero
finding problem is a potential to evaluate a power and a square root.
Algorithm 2: Stable Implementation
- Compute \(\check{r}_{sf \rightarrow rz}\) and \(\check{r}_{rz \rightarrow uz}\)
- Compute the lateral outflow
- If \(f\left(0\right) \geq 0\) then \(l_{sz}^{+} = l_{szmax}\)
- else find \(\tilde{z}\) such that \(\left\lvert f\left(\tilde{z}\right) \right\rvert \leq \epsilon\). Set \(l_{sz}^{+} = g\left(\tilde{z},\Theta\right)\)
- Solve for \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t}\)
- If \(h(0) \geq 0\) then
\(\left. \bar{s}_{sz}\right.\rvert_{\Delta t} = 0\).
- If \(h(z_{max}) < 0\) then \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t} = D\) and
\[
\frac{\Delta t}{\Delta x} l_{sz}^{+} =
D -
\left. \bar{s}_{sz}\right.\rvert_{0} +
\frac{\Delta t}{\Delta x} l_{sz}^{-} +
\Delta t \min\left(\frac{1}{T_d},\frac{
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} }
{ T_{d} D + \Delta t}\right)
\]
- else with \(z_{C} = \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t\left( \check{r}_{rz \rightarrow uz} - \frac{1}{T_{d}}\right)\)
- If \(h\left( z_{C} \right) \geq 0\) then
\[
\left. \bar{s}_{sz}\right.\rvert_{\Delta t} =
\left. \bar{s}_{sz}\right.\rvert_{0} +
\frac{\Delta t}{\Delta x} \left(g\left(\tilde{z},\Theta_{sz}\right) -
l_{sz}^{-}\right) - \frac{\Delta t}{T_d}
\]
- else
\[
\left. \bar{s}_{sz}\right.\rvert_{\Delta t} =
= \frac{1}{2} \left( z_{C} - h\left(z_{C}\right) + \sqrt{
\left(h\left(z_{C}\right)-z_{C}\right)^2 - 4 \frac{\Delta t}{T_{d}}
h\left(z_{C}\right) } \right)
\]
- Solve for \(\hat{r}_{rz \rightarrow uz}\)
- If \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t} = 0\) then
\[\hat{r}_{rz \rightarrow uz} = \frac{1}{\Delta t}\left(
\left. \bar{s}_{sz}\right.\rvert_{\Delta 0} +
\frac{\Delta t}{\Delta x} \left(g\left(0,\Theta_{sz}\right) -
l_{sz}^{-}\right) - \left. \bar{s}_{uz} \right.\rvert_{0} \right)
\]
- else
\[
\hat{r}_{rz \rightarrow uz} = \min\left(\check{r}_{rz \rightarrow uz},\frac{1}{\Delta t} \left(
\left. \bar{s}_{sz} \right.\rvert_{\Delta t} + \frac{\Delta t}{T_d} -
\left. \bar{s}_{uz} \right.\rvert_{0}\right)\right)\]
- Compute
\[
\hat{r}_{sf \rightarrow rz} = \min\left(\check{r}_{sf \rightarrow
rz}, \frac{1}{\Delta t}\left( s_{rzmax} -
\left. \bar{s}_{rz} \right.\rvert_{0} - \Delta t \left( p - e_p - \hat{r}_{rz
\rightarrow uz} \right) \right)\right)\]
- Solve for \(\left. \bar{s}_{uz}\right.\rvert_{\Delta t}\),
\(\left. \bar{s}_{rz}\right.\rvert_{\Delta t}\) and
\(\left. \bar{s}_{sf}\right.\rvert_{\Delta t}\) using the computed vertical
fluxes.
Transmissivity Profiles
The solution and numerical scheme have been presented with a general
transmissivity function \(g\left(\bar{s}_{sz},\Theta_{sz}\right)\). Table
2 present the transmissivity profiles present as options within the dynatop
package, the corresponding value to use for the transmissivity_profile value
in the model options vector and the additional parameters or properties
required. The additional parameter and properties are defined in Table 3.
Table 2: Transmissivity profiles available with the
dynatop package.
| Exponential |
\(T_{0}\sin\left(\beta\right)\exp\left(-\frac{\cos\beta}{m}s_{sz}\right)\) |
\(\frac{\cos\beta}{m}g\left(\bar{s}_{sz},\Theta_{sz}\right)\) |
\(T_0,\beta,m\) |
"exp" |
Originally given in Beven & Freer 2001 |
| Bounded Exponential |
\(T_{0}\sin\left(\beta\right)\left(\exp\left(-\frac{\cos\beta}{m}s_{sz}\right) - \exp\left(-\frac{\cos\beta}{m}D\right)\right)\) |
\(\frac{cos \beta}{m}T_{0}\sin\left(\beta\right)\exp\left(-\frac{\cos\beta}{m}s_{sz}\right)\) |
\(T_0,\beta,m,D\) |
"bexp" |
Originally given in Beven & Freer 2001 |
| Constant Celerity |
\(c_{sz}\left(D-s_{sz}\right)\) |
\(c_{sz}\) |
\(c_{sz},D\) |
"cnst" |
|
| Double Exponential |
\(T_{0}\sin\left(\beta\right)\left( \omega\exp\left(-\frac{\cos\beta}{m}s_{sz}\right) + \left(1-\omega\right)\exp\left(\frac{\cos\beta}{m_2}s_{sz}\right)\right)\) |
\(\omega\frac{cos \beta}{m}T_{0}\sin\left(\beta\right)\exp\left(-\frac{\cos\beta}{m}s_{sz}\right) + \left(1-\omega \right)\frac{cos \beta}{m_2}T_{0}\sin\left(\beta\right)\exp\left(-\frac{\cos\beta}{m_2}s_{sz}\right)\) |
\(T_0,\beta,m,m_2,\omega\) |
“dexp” |
|
Table 3: Additional parameters used in the transmissivity profiles.
| \(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\) |
Depth at which zero lateral flow occurs |
m |
| \(\omega\) |
weighting parameter |
- |
Appendix A - Monotonicity of \(f\left(z\right)\)
By definition \(\left. \bar{s}_{uz} \right.\rvert_{0}\), \(\Delta t\), \(\Delta X\),
\(T_d\), \(\check{r}_{rz \rightarrow uz}\) and \(z\) are all greater or equal to 0. Rewrite the expression for \(f\left(z\right)\) as
\[
f\left(z\right) = \left\{ \begin{array}{cl}
z -
\left. \bar{s}_{sz}\right.\rvert_{\Delta 0} -
\frac{\Delta t}{\Delta x} \left(g\left(z,\Theta_{sz}\right) -
l_{sz}^{-}\right) + \frac{\Delta t}{T_d} &
\frac{
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} }
{ T_{d} z + \Delta t} \geq \frac{1}{T_d} \\
z -
\left. \bar{s}_{sz}\right.\rvert_{\Delta 0} -
\frac{\Delta t}{\Delta x} \left(g\left(z,\Theta_{sz}\right) -
l_{sz}^{-}\right) + \Delta t \frac{
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} }
{ T_{d} z + \Delta t} &
\frac{
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} }
{ T_{d} z + \Delta t} < \frac{1}{T_d}
\end{array}\right.
\]
Differentiating with respect to \(z\) gives
\[
\frac{d}{dz} f\left(z\right) = \left\{ \begin{array}{cl}
1 -
\frac{\Delta t}{\Delta x} \frac{d}{dz}g\left(z,\Theta_{sz}\right) &
\frac{
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} }
{ T_{d} z + \Delta t} \geq \frac{1}{T_d} \\
1 -
\frac{\Delta t}{\Delta x} \frac{d}{dz}g\left(z,\Theta_{sz}\right)
- T_d \Delta t \frac{
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} }
{ \left(T_{d} z + \Delta t\right)^2} &
\frac{
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} }
{ T_{d} z + \Delta t} < \frac{1}{T_d}
\end{array}\right.
\]
The definition of the function \(g\) states that
\[
-\frac{d}{dz}g\left(z,\Theta_{sz}\right) \geq 0
\]
Using this relationship and the expression for the range of \(z\) gives
\[
\frac{d}{dz} f\left(z\right) \geq \left\{ \begin{array}{cl}
1 &
\frac{
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} }
{ T_{d} z + \Delta t} \geq \frac{1}{T_d} \\
1 - \frac{\Delta t}{T_{d} z + \Delta t} &
\frac{
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} }
{ T_{d} z + \Delta t} < \frac{1}{T_d}
\end{array}\right.
\]
From this it can be seen that \(\frac{d}{dz} f\left(z\right) \geq 0\) with
equality possible when \(z=0\).
Appendix B - Comparison to the Infinitesimal slab derivation
Surface Zone
A similar expression for the evolution can be derived by considering the
infinitesimal slab at the downslope end of the HRU. Taking a discrete
approximation for the change of flow with length the governing ODE can be written as
\[\begin{equation}
\frac{dw^{+}s_{sf}^{+}}{dt} = \frac{q_{sf}^{-} - q_{sf}^{+}}{\Delta x} - w^{+}r_{sf \rightarrow rz}^{+}
\end{equation}\]
Dividing through by \(w^{+}\) gives
\[\begin{equation}
\frac{ds_{sf}^{+}}{dt} = \frac{q_{sf}^{-}/w^{+} - l_{sf}^{+}}{\Delta x} - r_{sf \rightarrow rz}^{+}
\end{equation}\]
Note the similarities between this and the finite volume solution. In fact if
\(r_{sf \rightarrow rz}\) is taken to be spatially uniform and
\(\Delta x = A/w^{+}\) then this solution for \(s_{sf}^{+}\) will evolve in a similar fashion to
\(\bar{s}_{sf}\)
Root Zone
The same formulation follows for the infinitesimal slab at the bottom of the
hillslope with the averaged terms (using \(\bar{}\)) replaced with those valid
at the bottom of the hillslope (using \(^{+}\)).
Unsaturated Zone
The same formulation follows for the infinitesimal slab at the bottom of the
hillslope with the averaged terms (using \(\bar{}\)) replaced with those valid
at the bottom of the hillslope (using \(^{+}\)).
Saturated Zone
Evaluating for the foot of the hillslope and taking a discrete approximation
of the spatial gradient gives
\[\begin{equation}
\frac{ds_{sz}^{+}}{dt} = \frac{l_{sz}^{+} - q_{sz}^{-}/w^{+}}{\Delta x} - r_{uz \rightarrow sz}^{+}
\end{equation}\]
As with the surface zone we see the similarity in form, indicating that,
subject to the evaluation of the temporal integral, the two formulations
should return the same results when \(\Delta x = A/w^{+}\)
Appendix C Solution to \(h\left(z\right)\) when \(h(0) < 0\) and \(h(D) \geq 0\)
In this case we see that \(0 < z < D\). Taking \(z_{C} = \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t\left( \check{r}_{rz \rightarrow uz} - \frac{1}{T_{d}}\right)\) the function \(h\) can
be rewritten as
\[
h\left(z\right) = \left\{ \begin{array}{cl}
z -
\left. \bar{s}_{sz}\right.\rvert_{\Delta 0} -
\frac{\Delta t}{\Delta x} \left(l_{sz}^{+} -
l_{sz}^{-}\right) + \frac{\Delta t}{T_d}
& z_{C} \geq z \\
z -
\left. \bar{s}_{sz}\right.\rvert_{0} -
\frac{\Delta t}{\Delta x} \left(l_{sz}^{+} -
l_{sz}^{-}\right) + \Delta t \frac{
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} }
{ T_{d} z + \Delta t}
& z_{C} < z
\end{array} \right.
\]
Using the monotonicity of \(h\left(z\right)\) we solve each part
separately. Firstly if \(h\left( z_{C} \right) \geq 0\) then the inflow from the unsaturated
zone is \(1/T_{d}\) and
\[
\left. \bar{s}_{sz}\right.\rvert_{\Delta t} =
\left. \bar{s}_{sz}\right.\rvert_{0} +
\frac{\Delta t}{\Delta x} \left(g\left(\tilde{z},\Theta_{sz}\right) -
l_{sz}^{-}\right) - \frac{\Delta t}{T_d}
\]
if \(h\left( z_{C} \right) < 0\) the value of \(z\) satisfies
\[
z -
\left. \bar{s}_{sz}\right.\rvert_{0} -
\frac{\Delta t}{\Delta x} \left(l_{sz}^{+} -
l_{sz}^{-}\right) + \Delta t \frac{
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} }
{ T_{d} z + \Delta t} =
0
\]
Rearrangement to a quadratic gives
\[
0 = T_{d} z^2 - T_{d} z \left(\left. \bar{s}_{sz}\right.\rvert_{0} +
\frac{\Delta t}{\Delta x} \left(l_{sz}^{+} -
l_{sz}^{-}\right) \right) +
\Delta t z -
\Delta t \left(\left. \bar{s}_{sz}\right.\rvert_{0} +
\frac{\Delta t}{\Delta x} \left(l_{sz}^{+} -
l_{sz}^{-}\right) \right)
+
\Delta t \left(
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz}
\right)
\]
\[
0 = z^2 - z \left(\left. \bar{s}_{sz}\right.\rvert_{0} +
\frac{\Delta t}{\Delta x} \left(l_{sz}^{+} -
l_{sz}^{-}\right) - \frac{\Delta t}{T_{d}}\right) +
\frac{\Delta t}{T_{d}} \left(
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz
\rightarrow uz} -
\left. \bar{s}_{sz}\right.\rvert_{0} -
\frac{\Delta t}{\Delta x} \left(l_{sz}^{+} -
l_{sz}^{-}\right) \right)
\]
The commonly quoted solution for the quadratic
\[
0 = z^2 + b z +c
\]
is
\[
z = \frac{-b \pm\sqrt{ b^2 - 4c }}{2}
\]
Comparing terms we see that
\[
-b = \left. \bar{s}_{sz}\right.\rvert_{0} +
\frac{\Delta t}{\Delta x} \left(l_{sz}^{+} -
l_{sz}^{-}\right) - \frac{\Delta t}{T_{d}}
\]
and
\[
c =
\frac{\Delta t}{T_{d}} \left(
\left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz
\rightarrow uz} -
\left. \bar{s}_{sz}\right.\rvert_{0} -
\frac{\Delta t}{\Delta x} \left(l_{sz}^{+} -
l_{sz}^{-}\right) \right)
\]
Since \(-b = z_{C} - h\left(z_{C}\right)\) gives \(-b > 0\).
Since \(z > z_{C}\) the subsurface (\(s_{uz}\) and
\(s_{sz}\)) does not fill so
\[
\left. \bar{s}_{sz}\right.\rvert_{0} - \left. \bar{s}_{uz}
\right.\rvert_{0} \gt
\Delta t \check{r}_{rz \rightarrow uz} -
\frac{\Delta t}{\Delta x} \left(l_{sz}^{+} -
l_{sz}^{-}\right)
\]
and \(c < 0\)
Combining these conditions results in
\(\sqrt{b^2 - 4c} > \left\lvert b \right\rvert\). Since \(z>0\) only one root of
the quadratic is valid (as always positive) so
\[
z = \frac{1}{2}\left(-b + \sqrt{ b^2 - 4c }\right)
\]
