A Manual for dynatop
v0.5.0.1
1 Introduction
This manual presents details on the methods and applications contained
within the dynatop software for hydrological modelling. It is broken into the following sections:
- Conceptual Model provides a conceptual overview of the
dynatopdesign in relation to the original Dynamic TOPMODEL. - Hydrological Response Unit describes in detail the equations and and numerical solution used for each Hydrological Response Unit (HRU)
- Implementation provides further details of the software structure and data formats.
A seperate worked example walks through building and running a dynatop model for an example catchment.
1.1 Funding
The development of dynatop and this manual have benefited from funding from
- The NERC Q-NFM project (NE/R004722/1) at Lancaster University
- JBA Consulting through their Labs innovation funding.
1.2 Thanks
Thanks is due to the following
- Peter Metcalfe who developed the original port of Dynamic TOPMODEL to R during his PhD sponsored by the JBA Trust.
- The original developers of Dynamic TOPMODEL, Keith Beven and Jim Freer, who contributed to the original R port of Dynamic TOPMODEL.
1.3 Feedback
Feedback on this manual is always welcome via the issues page.
2 Conceptual Model
This a simplified, opinionated introduction. The original motivation for Dynamic TOPMODEL can be found in Beven and Freer (2001)
2.1 Motivation
Dynamic TOPMODEL was originally conceived as “A new version of the rainfall-runoff model TOPMODEL is described in which the assumption of a quasi-steady state saturated zone configuration is replaced by a kinematic wave routing of subsurface flow…” allowing for “…the simulation of dynamically variable up-slope contributing areas.”
Revisiting the idea of TOPMODEL Beven and Freer (2001) formulate models by
- Classifying the catchment up into areas that respond in a hydrologically similar manner
- Determining the connectivity between each of the classes
- Assigning each class to an appropriate Hydrological Representative Unit (HRU), which solves the hydrological processes for a standardised unit
As in earlier TOPMODEL work the identification of areas with a hydrological similar response was based on the classification of on topographic indices (but see Peters, Freer, and Beven (2003)). The use of topographic indices, such as the logarithm of up-slope area divided by the local gradient, can result HRUs representing disjoint spatial areas which may only partially capture the spatial ordering of hillslopes within the catchment.
Motivated by the desire to consider interventions at specific catchment
locations earlier versions of the dynatop package (as used in for example Beven et al. (2022)) included the concept of
hillslope order into the classification. This ensured that for a given HRU
only those upslope can provide inputs. This version of dynatop takes this concept a stage further with, by default, a HRU
being a single cell in the underlying digital terrain
model. Classification of HRUs by catchment properties is maintained to provide
a coherent method for assigning parameters. In the following these aspects are outlined with reference to the implementation of Dynamic TOPMODEL in the dynatop package.
2.2 Classification
The implementation of dynatop can assign HRUs to classes based on categorical
spatial data provided by the user. The model structure maintains the
values of categorical spatial data so they can be used in assigning parameters to the
HRUs. The classification creates a class for each unique combination
of categorical variable.
Consider a small catchment represented by a raster of spatial cells with the following soil classes
Figure 2.1: Map of soil classes
Since land use could influence the evapotranspiration and infiltration properties a second classification might be desirable. For the small catchment the land use classes are shown below
Figure 2.2: Map of land use classes
The landuse for a cell containing a channel may reflect the conditions in the parts of the cell not within the channel. The unique combinations of these classes (show below) provides the grouping used for parameterisation of the model.
Figure 2.3: Map of combined classes TODO
2.3 Connectivity
In keeping with earlier spatial analysis TOPMODEL and Dynamic TOPMODEL applications connectivity (flow directions) is assumed to be determined by the surface gradient with any down-slope region receiving flow bases on the method of Quinn et al. (1991). This gives a high degree of connectivity between the spatial cells. If a river intercepts the HRU it is assumed that both the surface and subsurface flow direction change to follow the river course. The following figure shows the connectivity for the small catchment, with the river cells highlighted in blue.
Figure 2.4: TODO - all river cells Map of combined topographic and land use classes along with connectivity
Each cell is represented by a single HRU, the inflows to which are the sum of all the inflow from the upslope cells. To represent the spatial ordering in the catchment the cells can be banded using the connectivity - all cells in a given band receive flows only from those in a higher band. For the small catchment this looks like the following
Figure 2.5: TODO include all channel cells Map of the band classification
Introducing the band into the classification has an advantage. Solving the HRUs starting with those in the highest band ensures that inflows at the current time step are available for use in the numeric solution. This allow the use of computation techniques that are more robust to the choice of solution time step.
The numeric scheme implemented in the
dynatoppackage makes use of an ordering stored in the model object so that the HRU inflows for the current time step are available. Although this ordering could be manually altered it is strongly recommended that the default values are used.
3 Hydrological Response Units
This chapter outlines the structure and computational implementation of the Hydrological Response Unit (HRU). We conceptualise the HRU as a cell, that is area where flux is vertical, and edges over which lateral fluxes pass.
The cell is formed of four zones representing the surface water, which passes over edges between cells 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 has edges to other cells. The zones and variables used below are shown in the schematic diagram.
In the following sections the governing equations of the cell are given in a finite volume form and the equations for transport over the edges are derived. Approximating equations for the solution of the governing equations are then derived with an associated semi-implicit numerical schemes.
3.1 Notation
The table below outlines the notation used for describing a 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.
| Quantity type | Symbol | Description | unit |
|---|---|---|---|
| Storage | \(s_{sf}\) | Surface excess storage | \(\rm{m}^3\) |
| \(s_{rz}\) | Root zone storage | \(\rm{m}^3\) | |
| \(s_{uz}\) | Unsaturated zone storage | \(\rm{m}^3\) | |
| \(s_{sz}\) | Saturated zone storage deficit | \(\rm{m}^3\) | |
| Vertical fluxes | \(r_{sf \rightarrow rz}\) | Flow from the surface excess store to the root zone | \(\rm{m}^3/\rm{s}\) |
| \(r_{rz \rightarrow uz}\) | Flow from the root zone to the unsaturated zone | \(\rm{m}^3\rm{s}\) | |
| \(r_{uz \rightarrow sz}\) | Flow from unsaturated to saturated zone | \(\rm{m}^3\rm{s}\) | |
| Lateral fluxes | \(q_{sf}\) | Lateral flow in the hillslope surface zone | \(\rm{m}^3/\rm{s}\) |
| \(q_{sz}\) | Lateral inflow to the hillslope surface zone | \(\rm{m}^3\rm{s}\) | |
| Vertical fluxes to the HRU | \(p\) | Precipitation rate | \(\rm{m}^3\rm{s}\) |
| \(e_p\) | Potential Evapotranspiration rate | \(\rm{m}^3\rm{s}\) | |
| Other HRU properties | \(A\) | Plan area | \(\rm{m}^2\) |
| \(w\) | Contour width of outflow | \(\rm{m}\) | |
| \(\beta\) | Surface slope angle | \(\rm{rad}\) | |
| Parameters | \(a_0\),\(a_1\),\(a_2\),\(b_0\),\(b_1\),\(b_2\),\(s_0\),\(s_1\) | Parameters for the solution of the surface zone | \(-\) |
| \(s_{rz_{max}}\) | Maximum rootzone depth | \(\rm{m}\) | |
| \(\lambda\),\(m\),\(K_0\),\(D\) | Parameters for the solution of the saturated zone | \(-\),\(\rm{m}\),\(\rm{m}/\rm{s}\),\(\rm{m}\) |
3.2 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 to occur over the cell edges and are considered as constant over 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.
3.2.1 Surface zone
Take flux over the surface to be positive when leaving the cell. The storage in the surface zone satisfies the mass balance equation
\[ \frac{ds_{sf}}{dt} = q_{sf}^{-} - q_{sf}^{+} - r_{sf \rightarrow rz} \]
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}\).
Surface outflow is given by
\[ q_{sf}^{+} = Q_{sf}\left(s_{sf}\right) \\ = a_{0}\min\left(s_{sf},s_{0}\right)^{b_{0}} + a_{1}\max\left(0,\min\left(s_{sf},s_{1}\right)-s_{0}\right)^{b_{1}} + a_{2}\max\left(0,s_{sf}-s_{1}\right)^{b_{2}} \]
This three layer structure can be made to match a variety of dynamics as outlined in [[Surface Zone Representations |surface_zone_representations]].
3.2.2 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
\[ 0 \leq s_{rz} \leq s_{rzmax} \]
with the governing ODE
\[ \frac{ds_{rz}}{dt} = p - \frac{e_p}{s_{rzmax}} s_{rz} + r_{sf\rightarrow rz} - r_{rz \rightarrow uz} \]
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
\[ p - e_p + r_{sf\rightarrow rz} - r_{rz \rightarrow uz} \leq 0 \]
In this case \(r_{rz \rightarrow uz}\) may be positive if
\[ p - e_p + r_{sf\rightarrow rz} > 0 \]
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.
3.2.3 Unsaturated & Saturated Zone
The unsaturated and saturated zones act in a coupled fashion subject to the constraint
\[ 0 \leq s_{uz} \leq s_{sz} \]
The continuity equations can be expressed as
\[ \frac{ds_{uz}}{dt} = r_{rz \rightarrow uz} - r_{uz \rightarrow sz} \]
\[ \frac{ds_{sz}}{dt} = q_{sz}^{+} - r_{uz \rightarrow uz} - q_{sz}^{-} \]
The outflow from the saturated zone is given in terms of a transmissivity \(T\left(s_{sz}\right)\) by
\[ q_{sz}^{+} = T\left(s_{sz}\right) w \tan\left(\beta\right) \]
In many TOPMODEL and dynamic TOPMODEL forumations the flux between the unsaturated ans saturated zones is given by
\[ r_{uz \rightarrow sz} = \frac{A s_{uz}}{T_d s_{sz}} \]
While it is acknowledged that this formaulation has no physical basis (Beven, The Primer, Box?? 20**) it has generally been interpreted as indicating that the unsaturated zone drains as a linear tank with time constant \(\frac{T_d s_{sz}}{A}\).
Here we develop an equally parsimonious alternative informed by the previsou analysis of (Beven 1984, Ion, 1989 paper, URBS model paper). Let \(K_0\) be the saturated conductivity of the surface layer. On the basis that saturated conductivity decreases with depth let, for for some storage defict \(z\) [m]
\[ K\left(z\right) = K_0 \exp\left(-\frac{z}{m}\right) \]
Assuming there is no preferential flow pathway the downslope transmissivity for a storage deficit \(z\) is given by
\[ T\left( z \right) = \tan\left(\beta\right) \int\limits_{z}^{\infty} K\left(x\right) dx \\ = \tan\left(\beta\right) K_0 \left[-m \exp\left(-\frac{x}{m}\right)\right]_{z}^{\infty} \\ = \tan\left(\beta\right) m K_0 \exp\left(-\frac{z}{m}\right) \]
Pragmatically we allow for preferential flow through the scaling factor \(\lambda\) which, noting that downslope flow is not perpendicular to the vertical boundary gives
\[ q_{sz}^{+} = w \lambda \cos\left(\beta\right) \tan\left(\beta\right) m K_0 \exp\left(-\frac{z}{zm}\right) \\ q_{sz}^{+} = w \lambda \sin\left(\beta\right) m K_0 \exp\left(-\frac{z}{m}\right) \\ = q_{0} \exp\left(-\frac{z}{m}\right) \]
Then using an approximation for z
\[ q_{sz}^{+} = q_{0} \exp\left(-\frac{s_sz}{Am}\right) \]
For vertical flow it is assumed that the flow interfaces saturate from above so the inflow to the unsaturated zone is \(\min\left(r_{rz \rightarrow uz},AK_{0}\right)\). However if the unsaturated zone is not full (i.e. not \(s_{uz} = s_{sz}\)) then the interface to the saturated zone may not be saturated from above. To approximate this the flux is taken to be
\[ r_{uz \rightarrow sz} = A K_0 \exp\left(-\frac{s_{sz}}{Am}\right) \frac{s_{uz}}{s_{sz}} \]
The ratio \(\frac{s_{uz}}{s_{sz}}\) approximates the gradient of the vertical hydraulic head.
The continuity equations can be rewritten as
\[ \frac{ds_{uz}}{dt} = \min\left(r_{rz \rightarrow uz},AK_{0}\right) - A K_0 \exp\left(-\frac{s_{sz}}{Am}\right) \frac{s_{uz}}{s_{sz}} \]
\[ \frac{ds_{sz}}{dt} = w\tan\left(\beta\right) m K_0 \exp\left(-\frac{s_{sz}}{Am}\right) - A K_0 \exp\left(-\frac{s_{sz}}{Am}\right) \frac{s_{uz}}{s_{sz}} - q_{sz}^{-} \]
or
\[ \frac{ds_{sz}}{dt} = \left( w\tan\left(\beta\right) m - A\frac{s_{uz}}{s_{sz}} \right) K_0 \exp\left(-\frac{s_{sz}}{Am}\right) - q_{sz}^{-} \]
3.3 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. A \(\tilde{}\) is used to note the values evaluated at the end the tiemstep.
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 giving the final states.
3.3.1 Downward Pass
3.3.1.1 Surface excess
The implicit solution of mass balance equation gives
\[ \tilde{s}_{sf} = s_{sf} + \Delta t \left( \tilde{q}_{sf}^{-} - Q_{sf}\left( \tilde{s}_{sf} \right) \right) - \tilde{v}_{sf \rightarrow rz} \]
Since downward flux is prioritised the maximum downward flux volume is \[ \hat{v}_{sf \rightarrow rz} = s_{sf} + \Delta t \tilde{q}_{sf}^{-} \]
3.3.1.2 Root Zone
The implicit formulation for the root zone gives \[\begin{equation} \tilde{s}_{rz} = \left( 1 + \frac{e_{p}\Delta t}{s_{rzmax}} \right)^{-1} \left( s_{rz} + \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}\) and at a maximum rate \(AK_{0}\) then \[\begin{equation} \hat{v}_{rz \rightarrow uz} = \min\left(A K_{0} \Delta t, \max\left(0, s_{rz} + \Delta t \left( p -e_{p}\right) + \hat{v}_{sf \rightarrow rz} - s_{rzmax} \right)\right) \end{equation}\]
3.3.1.3 Unsaturated & Saturated Zones
The unsaturated ans saturated zones are solved semi-implicitly by taking an implicit estimate the transmissivity of the interface between the unsaturated and saturated zones:
\[ \phi = AK_{0}\exp\left(-\frac{s_{sz}}{Am}\right) \]
Continuity for the unsaturated zone then gives \[ \tilde{s}_{uz} = \min\left( \tilde{s}_{sz}, s_{uz} + \hat{v}_{rz \rightarrow uz} - \Delta t \phi \frac{ \tilde{s}_{uz}} { \tilde{s}_{sz} } \right) \]
which can be rearranged to give \[ \frac{ \tilde{s}_{uz} }{ \tilde{s}_{sz} } = \min\left( 1, \frac{ s_{uz} + \hat{v}_{rz \rightarrow uz} }{\tilde{s}_{sz} + \phi \Delta t} \right) \]
and for the saturated zone
\[ \tilde{s}_{sz}= \max\left(0, s_{sz} + \Delta t q_{0} \exp\left(-\frac{\tilde{s}_{sz}}{Am}\right) - \Delta t \phi \min\left( 1, \frac{ s_{uz} + \hat{v}_{rz \rightarrow uz} }{\tilde{s}_{sz} + \phi \Delta t} \right) - q_{sz}^{-} \Delta t \right) \]
Consider the pseudo-function \[ \mathcal{H}\left(z\right) = s_{sz} + \Delta t q_{0} \exp\left(-\frac{z}{Am}\right) - \Delta t \phi \min\left( 1, \frac{ s_{uz} + \hat{v}_{rz \rightarrow uz} }{z + \phi \Delta t} \right) - q_{sz}^{-} \Delta t - z \]
adapted from the approximating equation for continuity in the saturated zone storage. Differentiating with respect to \(z\) gives
\[ \frac{\partial}{\partial z}\mathcal{H}\left(z\right) = \left\{ \begin{array}{cc} - \frac{\Delta t q_{0}}{Am}\exp\left(-\frac{z}{Am}\right) + \frac{\Delta t \phi}{z + \phi \Delta t} \frac{ \left. s_{uz} \right.\rvert_{t} + \hat{v}_{rz \rightarrow uz} } {z + \phi \Delta t} - 1 & z \geq s_{uz} + \hat{v}_{rz \rightarrow uz} - \phi \Delta t \\ -\frac{\Delta t q_{0}}{Am}\exp\left(-\frac{z}{Am}\right) - 1 & \rm{otherwise}\\ \end{array} \right. \]
from which in can be seen that \(\mathcal{H}\left(z\right)\) is a monotonic decreasing function of \(z > 0\) and \(\mathcal{H}\left(z\right) \rightarrow -\infty\) as \(z \rightarrow \infty\).
Hence if \(\mathcal{H}\left(0\right) \gt 0\) a unique solution \(\mathcal{H}\left(z\right) = 0\) can be found to provide a value for \(\left. s_{sz}\right.\rvert_{\Delta t}\). If \(\mathcal{H}\left(0\right) \leq 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\). The 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(\tilde{s}_{sz}\right)\leq 0\). We therefore seek \(z\) which satifies for some tolerance \(\epsilon>0\) \[ -\epsilon \leq \mathcal{H}\left(z\right) \leq 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)\geq 0 \geq \mathcal{H}\left(z^{u}\right)\), thus bracket the solution. Iterating (using for example bisection) until \(\mathcal{H}\left(z^{u}\right) \geq -\epsilon\) and give the approximation \(\tilde{s}_{sz} = z^{u}\).
3.3.2 Upward Pass
Given an estimate of \(\tilde{s}_{sz}\) such that \(\mathcal{H}\left(\tilde{s}_{sz}\right)>0\) compute the vertical flux volumes in the following fashion
- Compute \(\psi = \min\left( 1, \frac{ s_{uz} + \hat{v}_{rz \rightarrow uz} }{\tilde{s}_{sz} + \phi \Delta t} \right)\)
- \(v_{uz \rightarrow sz} = \Delta t \phi \psi\)
- \(q_{sf}^{+} = \frac{1}{\Delta t} \left( \tilde{s}_{sz} - s_{sz} + \Delta t q_{sf}^{-} + v_{uz \rightarrow sz} \right)\)
- \(\tilde{s}_{uz} = \min\left( \tilde{s}_{sz}, s_{uz} + \hat{v}_{rz \rightarrow uz} - v_{uz \rightarrow sz}\right)\)
- \(v_{rz \rightarrow uz} = \tilde{s}_{uz} + v_{uz \rightarrow sz} - s_{uz}\)
- Compute \(\omega = 1 + \frac{e_{p}\Delta t}{s_{rzmax}}\)
- Solve \(\tilde{s}_{rz} = \min\left(s_{rzmax}, \frac{1}{\omega}\left(s_{rz} + \Delta t p + \hat{v}_{sf \rightarrow rz} - v_{rz \rightarrow uz} \right)\right)\)
- \(v_{sf \rightarrow rz} = \tilde{s}_{rz} \omega + v_{rz \rightarrow uz} - \Delta t p - s_{rz}\)
- Compute the flux loss to evapotranspiration as \(e_a = \frac{1}{\Delta t} \left( s_{rz} + v_{sf \rightarrow rz} - v_{rz \rightarrow uz} + \Delta t p - \tilde{s}_{rz}\right)\)
This leaves the solution of the surface zone
3.3.2.1 Surface Zone
Proceeding in a similar fashion to the saturated zone consider the pseudo-function \[ \mathcal{S}\left(w\right) = s_{sf} + \Delta t \left( q_{sf}^{-} - Q_{sf}\left( w \right) \right) - v_{sf \rightarrow rz} - w \] which is monotonically decreasing for \(w \geq 0\).
If \(\mathcal{S}\left(0\right) \leq 0\) then \(\tilde{s}_{sf} = 0\). Otherwise a bracketing solution algorithm that gives at any iteration \(w^{l} < \tilde{s}_{sf} < w^{u}\) is used. An initial value of \(w^{u}\) can be computed as \(s_{sf} + \Delta t q_{sf}^{-} - v_{sf \rightarrow rz}\). Iterating until; for some tolerance \(\epsilon>0\); \(S\left(w^{l}\right) < \epsilon\) take \(\tilde{s}_{sf} = w^{l}\). In this case \(\mathcal{S}\left(\tilde{s}_{sf}\right) > 0\) hence mass conservation can be maintained by taking \[ q_{sf}^{+} = q_{sf}^{-} + \frac{1}{\Delta t}\left( s_{sf} - v_{sf \rightarrow rz} - \tilde{s}_{sf}\right) \]
3.4 Linking HRUs
Fluxes between the HRUs are represented as occurring at two levels, the surface and the saturated zones. For the \(i\)th HRU the outflow in \(m^3/s\) at time \(t+\Delta t\) is \(\left.q_{sf}^{\left[i\right]}\right\rvert_{t+\Delta t}\) at the surface and \(\left.q_{sz}^{\left[i\right]}\right\rvert_{t+\Delta t}\) in the saturated zone.
The fraction of outflow going from the \(i\) th to the \(j\) th HRU is \(f_{i,j}\). The redistribution is conservative as does not vary in time so \[ \sum\limits_{j} f_{i,j} = 1 \] Since there is no exchange between the surface and saturated zones during redistribution the inflows at time \(t+\Delta t\) to the \(k\)th HRU can be computed as \[ \sum\limits_{i>k} f_{i,k}\left.q_{sf}^{\left[i\right]}\right\rvert_{t+\Delta t} \] and \[ \sum\limits_{i>k} f_{i,k}\left.q_{sz}^{\left[i\right]}\right\rvert_{t+\Delta t} \]
4 Implementation
The dynatop package represents an implementation of the concepts behind a dynamic TOPMODEL.
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.
While a Hydrological Response Unit (HRU) has typically been thought of as an area of hill slope with, for example, similar topographic, soil and upslope area characteristics the idea may be generalised. In doing this a catchments is conceptualised as a collection of HRUs which may represent of different types of landscape (e.g. hill slope, channel, lake) which exchange fluxes at specified levels (e.g. saturated zone, surface). Currently a single type of HRU is present in the model, which may be parameterised to represent a number of different physical features.
The aim of this document is to outline
- the computational sequence and linking the HRUs
- conventions for the input data
- other miscellaneous information that may help developers
The following schematic outlines the stages resulting a working Dynamic
TOPMODEL for a catchment. It highlights where the dynatopGIS and dynatop
packages fit into this process.
Figure 4.1: Schematic of a Dynamic TOPMODEL construction
As we will see in the example the dynatopGIS package can help
- Compute standard variables for classification (such as the topographic index)
- Compute the spatial ordering of the HRUs
- Build classifications of the catchment
- Construct models suitable for
dynatop
The dynatop package allows simulation and visualisation of Dynamic TOPMODELs as well as
providing some helper for processing time series input data.
Although in this training presumes that model for dynatop are generated
using dynatopGIS this need not be the case so long as the inputs to
dynatop are of the correct format.
4.1 Implementation notes
The dynatopGIS and dynatop packages implement structured, object orientated, data
flows. The packages are written using the object orientated framework
provided by the R6 package. This means that some aspects of working with the
objects may appear idiosyncratic for some R users. In this training these problems are largely obscured, except for the
call structure. However, before adapting the code, or doing more complex analysis
users should read about R6 class objects (e.g. in the R6 package vignettes
or in the Advanced R book). One particular gotcha is when copying an object. Using
my_new_object <- my_objectcreates a pointer, that is altering my_new_object also alters
my_object. To create a new independent copy of my_object use
my_new_object <- my_object$clone()4.2 Input Series
It is expected that the precipitation and potential evapotranspiration inputs series are given in \(m\) accrued over the proceeding time step. So if the data has a time step \(\Delta t\) the value given at \(t+\Delta t\) is accrued in the interval between time \(t\) and time \(t+\Delta t\).
4.3 Relating the vignette notation and code variables
To aid the readability of the code the variables are labelled consistently with regards to the vignettes. For example:
- Single subscript: \(l_{sz}\) becomes
l_sz, \(s_{uz}\) becomess_uzetc. - Directional variables: \(q_{rz \rightarrow uz}\) becomes
q_rz_uz, etc. - Greek letters are spelt: \(\lambda\) becomes
lambda
However the notation used in the vignettes for intermediate values of fluxes or states (e.g. ) is dropped from the code where it can be inferred from the computational sequence (or comments)
4.4 Numerical solutions
Currently the code uses a purely implicit scheme for the solution of the HRUs. This requires a solving two zero finding problems for each HRU at each (sub)time step. Currently a bisection algorithm is used with a user specified tolerance and maximum number of iterations. The tolerance is defined as the difference between the upper and lower limits of the interval containing the zero point. Notionally faster algorithms (such as the TOMS-748 algorithm) could be used.
4.4.1 The model object
The purpose of this vignette is to document the structure of a Dynamic
TOPMODEL model ‘object’; that is the structure of the model that can be generated in
the dynatopGIS package and is
used in calls to the dynatop function in the current package. Currently
the model is stored as an R variable of the list class.
This should be read alongside the “Getting Started” vignette which demonstrates the manipulation and use of a model object.
The model structure will be documented with reference to the model in the test catchment data supplied with the package. This can be loaded with
library(dynatop)
data("Swindale")
model <- Swindale$modelThe top level of the model contains the following elements:
names(model)We consider each of these separately.
4.5 map
The map element gives the path to the raster map of HRU locations which can
be used in the creation of a new dynatop object to enable plotting of
the state variables.
4.6 output_flux
The output_flux element contains a data frame which defines the series to be output
during a simulation. This is passed to the sim method of the dynatop
object. The columns of output_flux are given in Table
4.1.
A single output value for each series is given for each data time step. If multiple rows of the data frame contain the same name the values of the series are summed. If the simulation is run with substeps then the outputs are averaged over the sub steps. Note that strictly only “q_sf”, “q_sf_in”, “q_sz” and “q_sz_in” are fluxes, the remainder are volumes.
| Name | Class | Unit | Description |
|---|---|---|---|
| name | character | - | name given to the output series |
| id | integer | - | id of the HRU from which the output is to be taken |
| flux | character | - | Flux to be recorded. One of “precip”, “pet”, “aet”, “q_sf”, “q_sf_in”, “q_sz”, “q_sz_in”, “s_sf”, “s_rz”, “s_uz”, “s_sz”, “v_sf_rz”, “v_rz_uz” or “v_uz_sz”. |
4.7 hru
The hru contains a list, where each list element is the description of a
single HRU. Each HRU is described by a list containing the elements outlined
in Table 4.2 and described in more detail in the sections
below. In many cases the naming of the variables within the vectors relates to
the variables names in the HRU vignette as given in
implementation document.
| Name | Class | Description |
|---|---|---|
| id | integer | The unique ID of the HRU |
| states | numeric | Named numeric vector of states |
| properties | numeric | Named numeric vector of properties |
| sf | list | Description of the surface zone solution |
| rz | list | Description of the surface zone solution |
| uz | list | Description of the unsaturated zone solution |
| sz | list | Description of the saturated solution |
| sf_flow_direction | list | Description of the surface zone outflows to other HRUs |
| sz_flow_direction | list | Description of the saturated zone outflows to other HRUs |
| initialisation | numeric | Named numeric vector of initialisation parameters |
| precip | list | Description of the precipitation input to the HRU |
| pet | list | Description of the PET input to the HRU |
| class | list | Classification information from dynatopGIS |
4.7.1 states
A numeric vector of states, see Table 4.3, such as
model$hru[[1]]$states| Name | Class | Unit | Description |
|---|---|---|---|
| s_sf | numeric | m | Surface Zone storage volume per unit area |
| s_rz | numeric | m | Root Zone storage volume per unit area |
| s_uz | numeric | m | Unsaturated Zone storage volume per unit area |
| s_sz | numeric | m | Saturated Zone storage deficit per unit area |
4.7.2 properties
A numeric vector of properties, see Table 4.4 such as
model$hru[[1]]$properties| Name | Class | Unit | Description |
|---|---|---|---|
| area | numeric | m\(^2\) | Surface area of the HRU |
| width | numeric | m | Width of the HRU |
| Dx | numeric | m | Effective length of the HRU |
| gradient | numeric | - | Average gradient |
4.7.3 sf, rz, uz,& sz
The list variables describe the surface (sf), root zone (rz), unsaturated zone
(uz) ans saturated zone (sz) representations. Each list contains a character
vector type giving the type of representation and parameters a numeric
vector of names parameters. For example
model$hru[[1]]$szThe description of the types and parameters are given in the HRU vignette
4.7.4 sf_flow_direction & sz_flow_direction
Each of these lists contains two vectors as specified in Table 4.5. In the special case that no outflow goes to another HRU, such as at the outlet of the catchment, the description of the flow direction takes the following form:
model$hru[[1]]$sf_flow_direction| Name | Class | Unit | Description |
|---|---|---|---|
| id | integer | - | Vector of id values for the HRUs to which the flow is going. Must be lower then the current id |
| frc | numeric | 0-1 | Fraction of the flow which goes from the corresponding id. Must sum to 1 |
4.7.5 initialisation
The named numeric vector contains the values to use in the initialisation of the HRU as given in Table 4.6 and shown below.
model$hru[[1]]$initialisation| Name | Class | Unit | Description |
|---|---|---|---|
| s_rz0 | numeric | 0-1 | Initial root zone depth expressed as fraction of maximum depth (parameter) |
| r_uz_sz0 | numeric | m/s | Initial recharge to the saturated zone per unit area (parameter) |
4.7.6 precip & pet
These list contain the description of the precipitation (precip) and PET (pet) inputs used by the HRU. Each list takes the form of two vectors as given in Table 4.7 and shown below. From this the average input value is calculated.
model$hru[[1]]$precip| Name | Class | Unit | Description |
|---|---|---|---|
| name | character | - | Character vector containing the names of the input series to use. |
| frc | numeric | 0-1 | The fraction of the HRU area for which the input series is valid. Must sum to 1. |
4.7.7 class
Class is a list variable. while not used in the simulation code it can be used to store information that might be of use in, for example, determining which parameter values to apply to the HRU (see the “Getting Started” vignette for an example).
Currently if the model is constructed using
dynatopGIS class is populated
based on the classification and channel data used in the model build process,
but the format of this is subject to change.
4.7.7.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1