The dynatop function in the dynamic TOPMODEL R package implements a set of numerical solutions to the governing equations of the dynamic TOPMODEL hillslope formulation. This returns the inflow to the channel reaches specified in the model. This vignette documents the formulation and the numerical solution used to route these channel inflows, along with other diffuse and point inputs to the model.

Notation and nomenclature

Table @ref(tab:ch-notation} outlines the notation used in this document.

Table 1: Outline of notation for describing the Channel HRU section of the hillslope HRU
Quantity type Symbol Description unit
Discharge from a reach at time \(t\) \(q\left(t\right)\) Instantaneous discharge m\(^3\)/s
Discharge from a reach at time \(t\) \(g\left(t\right)\) Instantaneous discharge m\(^3\)/s
Diffuse inflow \(d\left(t\right)\) Instantaneous inflow from diffuse inputs to the channel HRU m\(^3\)/s
Point inflow \(p\left(t\right)\) Instantaneous inflow from a point source to the head of the Channel HRU m\(^3\)/s
Velocity \(v_{ch}\) Effective velocity of water conveyance within the Channel HRU network m/s

Inputs to a Channel HRU

A Channel HRU is considered as representing a channel reach of length \(L\). Each Channel HRU has two types of inflow, labelled diffuse and point. Diffuse inflows are external inputs which occur at a constant rate over the length of the reach. These include precipitation, the inflow from Hillslope HRUs and time series specified in the diffuse_inputs table.

Point inflows are considered to occur at the head of the channel reach. These consist of the inflow from upstream channel HRUs and time series specified in the point_inputs table.

Time Delay Channel Routing

The time delay channel routing is selected by setting the channel_solver value in the model options vector to histogram. The method offers a basic representation of the movement of water in which a constant velocity to assigned to each channel reach.

For the diffuse inputs to a single channel reach of length \(L\) the contribution to the outflow at time \(t\) is given by \[ \frac{1}{L} \int\limits_{l=0}^{L} d\left(t - \frac{l}{v_{ch}}\right) dl \] where as for a point input to the same channel the contribution is \[ p\left(t - \frac{L}{v_{ch}}\right) \] such that \[ q\left(t\right) = \frac{1}{L} \int\limits_{l=0}^{L} d\left(t - \frac{l}{v_{ch}}\right) dl + p\left(t - \frac{L}{v_{ch}}\right) \]

Consider next a change of variables by taking \(\tau = l/v_{ch}\) and defining \(\tau_{L} = L/v_{ch}\) so that \[ q\left(t\right) = \frac{v_{ch}}{L} \int\limits_{\tau=0}^{\tau_{L}} d\left(t - \tau\right) d\tau + p\left(t - \tau_{L}\right) \]

Since the action on the point inputs is purely advective the flow at a gauged location can be written as the summation of the diffuse and external point inflows to the Channel HRUs (reaches) above the location. For the \(k\)th Channel HRU above the gauge let \(\tau_{0}^{\left[k\right]}\) be the advective time delay between water leaving the reach and arriving at the gauge. This is readily computed from the channel HRU connectivity, reach lengths and velocities. The gauged flow can thus be expressed as

\[ g\left(t\right) = \sum_{k} \frac{v_{ch}^{\left[k\right]}}{L^{\left[k\right]}} \int\limits_{\tau_{0}^{\left[k\right]}}^{\tau_{0}^{\left[k\right]} + \tau_{L}^{\left[k\right]}} d\left(t - \tau \right) d\tau + \hat{p}\left(t - \tau_{0}^{\left[k\right]} - \tau_{L}^{\left[k\right]}\right) \]

To solve for the gauge flow on a discrete time step \(\Delta t\) suppose that \(t = i\Delta t\) and that the flow index by the time subscript is valid over the preceding time step; for example \[ d\left(t\right) = d_{i\Delta t} \quad \left(i-1\right)\Delta t \leq t <i\Delta t \]

Then to evaluate \(g_{i\Delta t}\) while maintaining mass conservation in the formulation, the gauged flows must be integrated over the time period. Consider one of the contributing HRUs, then (dropping the superscript \(k\)) the contribution of the point inflow is \[ \frac{1}{\Delta t}\int\limits_{\left(i-1\right)\Delta t}^{i\Delta t} \hat{p}\left(t - \tau_{0} - \tau_{L}\right) dt = \frac{\left(r_{L}+1\right)\Delta t - \left(\tau_{0}+\tau_{L}\right)}{\Delta t} \hat{p}_{\left(i-r_L\right)\Delta t} + \frac{\tau_0 + \tau_L- r_L \Delta t}{\Delta t} \hat{p}_{\left(i-1-r_L\right)\Delta t} \] where \(r_L\) is the largest integer for which \(r_{L}\Delta t \leq \tau_0 + \tau_L\). The right hand side of this expression gives the time delay histogram for the point inflow.

As similar expression holds for each value of \(\tau\) along the HRU, such that the average inflow at that point over the time step is given by \[ \frac{1}{\Delta t}\int\limits_{\left(i-1\right)\Delta t}^{i\Delta t} d\left(t - \tau\right) dt = \frac{\left(r_{\tau}+1\right)\Delta t - \tau}{\Delta t} d_{\left(i-r_{\tau}\right)\Delta t} + \frac{\tau- r_{\tau} \Delta t}{\Delta t} d_{\left(i-1-r_{\tau}\right)\Delta t} \] where \(r_\tau\) is the largest integer for which \(r_{\tau}\Delta t \leq \tau\). Next consider the integral of the diffuse inflows over \(\tau\); that is \[ \int\limits_{\tau_{0}}^{\tau_{0} + \tau_{L}} \frac{\left(r_{\tau}+1\right)\Delta t - \tau}{\Delta t} d_{\left(i-r_{\tau}\right)\Delta t} + \frac{\tau- r_{\tau} \Delta t}{\Delta t} d_{\left(i-1-r_{\tau}\right)\Delta t} d\tau \]

The range of \(\tau\) can be divided into intervals which share the same value of \(r_{\tau}\). Defining the minimum value of \(r_\tau\) as \(r_{0} = r_{\tau_{0}}\) then if \(r_{L}>r_{0}\) the intervals are given by \[ \left(\tau_0,\left(r_{0}+1\right)\Delta t\right), \left(\left(r_{0}+1\right)\Delta t,\left(r_{0}+2\right)\Delta t\right), \ldots, \left(\left(r_{0}+n\right)\Delta t,\tau_{0}+\tau_{L}\right) \]

For \(j=1,\ldots,n-1\) the integral over an interval of \(\Delta t\) gives \[ \int\limits_{\left(r_{0}+j\right)\Delta t}^{\left(r_{0}+j+1\right)\Delta t} \frac{\left(r_{0}+j+1\right)\Delta t - \tau}{\Delta t} d_{\left(i-r_{0}-j\right)\Delta t} + \frac{\tau- \left(r_{0}+j\right) \Delta t}{\Delta t} d_{\left(i-r_{0}-j-1\right)\Delta t} d\tau \\ = \frac{\Delta t}{2} \left(\left(r_{0}+j+1\right)-\left(r_{0}+j\right)\right)^2 d_{\left(i-r-j\right)\Delta t} + \frac{\Delta t}{2} \left(\left(r_{0}+j+1\right)-\left(r_{0}+j\right)\right)^2 d_{\left(i-r-j-1\right)\Delta t} \\ = \frac{\Delta t}{2} d_{\left(i-r_{0}-j\right)\Delta t} + \frac{\Delta t}{2} d_{\left(i-r_{0}-j-1\right)\Delta t} \]

For the initial interval \[ \int\limits_{\tau_{0}}^{\left(r_{0}+1\right)\Delta t} \frac{\left(r_{0}+1\right)\Delta t - \tau}{\Delta t} d_{\left(i-r_{0}\right)\Delta t} + \frac{\tau- r_{0}\Delta t}{\Delta t} d_{\left(i-r_{0}-1\right)\Delta t} d\tau \\ = \frac{1}{2\Delta t} \left(\left(r_0+1\right)\Delta t - \tau_{0}\right)^2 d_{\left(i-r_{0}\right)\Delta t} + \left\{ \frac{1}{2\Delta t} \left(\left(r_0+1\right)\Delta t - \tau_{0}\right)^2 + \frac{ \left(\tau_{0} - r_{0}\Delta t\right)\left(\left(r_0+1\right)\Delta t - \tau_{0}\right)}{\Delta t} \right\} d_{\left(i-r_{0}-1\right)\Delta t} \]

and the final interval and noting that \(r_{0}+n = r_{L}\) \[ \int\limits_{r_{L}\Delta t}^{\tau_0+\tau_L} \frac{\left(r_{L}+1\right)\Delta t - \tau}{\Delta t} d_{\left(i-r_{L}\right)\Delta t} + \frac{\tau- r_{L} \Delta t}{\Delta t} d_{\left(i-r_{L}-1\right)\Delta t} d\tau \\ = \left\{ \frac{1}{2\Delta t} \left(\tau_0 + \tau_{L} - r_L\Delta t\right)^2 + \frac{ \left(\tau_0 + \tau_{L} - r_L\Delta t\right)\left(\left(r_L + 1 \right)\Delta t - \tau_{0} - \tau_L\right)}{\Delta t} \right\}d_{\left(i-r_{L}\right)\Delta t} + \frac{1}{2\Delta t} \left(\tau_0 + \tau_{L} - r_L\Delta t\right)^2 d_{\left(i-r_{L}-1\right)\Delta t} \]

By grouping the terms for the individual inputs together and multiply through by \(v_{ch}/L\) time delay histograms can be constructed for each HRU draining to the gauge. Table 2 shows the solution for \(r_{L}>r_0+1\) and Table 3 for the case \(r_{L}=r_0+1\).

Table 2: Outline of the time delay histogram for the diffuse input in the case \(r_{L}>r_0+1\)
Delay Inflow Weight
\(r=0,\ldots,r_0-1\) \(d_{\left(i-r\right)\Delta t}\) 0
\(r_0\) \(d_{\left(i-r_0\right)\Delta t}\) \(\frac{v_{ch}}{2L\Delta t} \left(\left(r_0+1\right)\Delta t - \tau_{0}\right)^2\)
\(r_0+1\) \(d_{\left(i-r_0-1\right)\Delta t}\) \(\frac{v_{ch}\Delta t}{2L} + \frac{v_{ch}}{2L\Delta t} \left(\left(r_0+1\right)\Delta t - \tau_{0}\right)^2 + \frac{ v_{ch}\left(\tau_{0} - r_{0}\Delta t\right)\left(\left(r_0+1\right)\Delta t - \tau_{0}\right)}{L\Delta t}\)
\(r=r_{0}+2,\ldots,r_{L}-1\) \(d_{\left(i-r\right)\Delta t}\) \(\frac{v_{ch}\Delta t}{L}\)
\(r_{L}\) \(d_{\left(i-r_{L}\right)\Delta t}\) \(\frac{v_{ch}\Delta t}{2L} + \frac{v_{ch}}{2L\Delta t} \left(\tau_0 + \tau_{L} - r_L\Delta t\right)^2 + \frac{ v_{ch}\left(\tau_0 + \tau_{L} - r_L\Delta t\right)\left(\left(r_L + 1 \right)\Delta t - \tau_{0} - \tau_L\right)}{L\Delta t}\)
\(r_{L}+1\) \(d_{\left(i-r_{L}-1\right)\Delta t}\) \(\frac{v_{ch}}{2L\Delta t} \left(\tau_0 + \tau_{L} - r_L\Delta t\right)^2\)
Table 3: Outline of the time delay histogram for the diffuse input in the case \(r_{L}=r_0+1\)
Delay Inflow Weight
\(r=0,\ldots,r_0-1\) \(d_{\left(i-r\right)\Delta t}\) 0
\(r_0\) \(d_{\left(i-r_0\right)\Delta t}\) \(\frac{v_{ch}}{2L\Delta t} \left(\left(r_0+1\right)\Delta t - \tau_{0}\right)^2\)
\(r_0+1 = r_{L}\) \(d_{\left(i-r_0-1\right)\Delta t}\) \[\frac{v_{ch}}{2L\Delta t}\left(\left(r_0+1\right)\Delta t - \tau_{0}\right)^2 + \frac{v_{ch}\left(\tau_{0} - r_{0}\Delta t\right)\left(\left(r_0+1\right)\Delta t- \tau_{0}\right)}{L\Delta t} +\\ \frac{v_{ch}}{L\Delta t}\left(\tau_0 + \tau_{L} - r_L\Delta t\right)^2 + \frac{ v_{ch}\left(\tau_0 + \tau_{L} - r_L\Delta t\right)\left(\left(r_L + 1 \right)\Delta t - \tau_{0} - \tau_L\right)}{L\Delta t}\]
\(r_{L}+1\) \(d_{\left(i-r_{L}-1\right)\Delta t}\) \(\frac{v_{ch}}{2L\Delta t} \left(\tau_0 + \tau_{L} - r_L\Delta t\right)^2\)

The final case to consider is \(r_{0} = r_{L}\). In this case there is a single integral to evaluate: \[ \int\limits_{\tau_{0}}^{\tau_{0}+\tau_{L}} \frac{\left(r_{0}+1\right)\Delta t - \tau}{\Delta t} d_{\left(i-r_{0}\right)\Delta t} + \frac{\tau- r_{0} \Delta t}{\Delta t} d_{\left(i-r_{0}-1\right)\Delta t} d\tau \\ = \frac{\tau_{L}}{\Delta t} \left(\left(r_{0}+1\right)\Delta t - \tau_{0} - \tau_{L}/2 \right) d_{\left(i-r_{0}\right)\Delta t} + \frac{\tau_{L}}{\Delta t} \left(\tau_{)} + \tau_{L}/2 - r_{0}\Delta t\right) d_{\left(i-r_{0}-1\right)\Delta t} \]

Examples

In this section some example histograms are generated based upon the length, velocity and \(\tau_0\) of the Channel HRU along with the time step \(\Delta t\).

Point input

First a function for defining the histogram for a point inflow can be coded as follows:

## L is length of channel
## vch is channel velocity
## Dt is the time step
## tau0 is the time delay between the gauge and foot of the reach
fp <- function(L,vch,Dt,tau0){
  tauL <- L/vch
  rL <- floor((tau0+tauL)/Dt)
  irL <- rL+1 ## index in vector since R starts at 1 not zero
  b <- rep(0,irL+1)
  b[irL] <- ((rL+1)*Dt - tau0-tauL)/Dt
  b[irL+1] <- (tau0+tauL - rL*Dt)/Dt
  return(b)
}

For the first example all the input should occur with one time step delay (second element in output vector) since \(L/v_ch\) exactly matches the time step and \(\tau_{0}=0\)

fp(L=100,vch=0.5,Dt=200,tau0=0)
#> [1] 0 1 0

In this example \(\tau_{0}+L/v_ch = \frac{7}{2}\Delta t\); so the first non zero value of the vector should be the 4th element with the 4th and 5th elements being equal.

fp(L=100,vch=1,Dt=200,tau0=600)
#> [1] 0.0 0.0 0.0 0.5 0.5

This build upon the above example, altering the velocity such that the only non-zero elements are the 4th and 5th, but with the 4th being larger then the 5th.

fp(L=100,vch=1.5,Dt=200,tau0=600)
#> [1] 0.0000000 0.0000000 0.0000000 0.6666667 0.3333333

Diffuse inflow

The function for the diffuse inflow is more complex since the special cases have to be handled:

fd <- function(L,vch,Dt,tau0){
  tauL <- L/vch
  r0 <- floor(tau0/Dt)
  ir0 <- r0+1 ## index in vector since R starts at 1 not zero
  rL <- floor((tau0+tauL)/Dt)
  irL <- rL+1 ## index in vector since R starts at 1 not zero
  b <- rep(0,irL+1)
  
  if(rL>r0){
    b[ir0:(irL-1)] <- Dt # inital values valid unless over written
    b[ir0] <- ( ((r0+1)*Dt - tau0)^2 ) / (2*Dt)
    b[ir0+1] <- b[ir0] + ( ((r0+1)*Dt - tau0)*(tau0 - (r0*Dt)) ) / Dt
    b[irL+1] <- ( (tau0+tauL - (rL*Dt))^2 ) / (2*Dt)
    b[irL] <- b[irL] + b[irL+1] + ( ((tau0+tauL - (rL*Dt))*((rL+1)*Dt - tau0-tauL)) / Dt ) ## added to self since rL could equal r0+1
    if( rL > (r0+1) ){
      b[ir0+1] <- b[ir0+1] + Dt/2
      b[irL] <- b[irL] + Dt/2
    }
  }else{
    b[ir0] <- (tauL/Dt)*( (r0+1)*Dt - tau0 - (tauL/2) )
    b[ir0+1] <- (tauL/Dt)*(tau0 + (tauL/2) - r0*Dt )
  }
  
  return(b*vch/L)
}

In the first example there is no advective time delay and \(L/v_{ch} = \Delta t\). This results in time delay histogram in which the first two values equal 0.5. This is an example of the case \(r_{L}=r_{0}+1\)

b <- fd(L=200,vch=1,Dt=200,tau0=0)
sum(b)
#> [1] 1
barplot(b)

Extending the previous example the velocity now altered to that \(r_{L}=r_{0}\). The resulting non-negative histogram values are now the first two in the vector.

b <- fd(L=200,vch=1.5,Dt=200,tau0=0)
sum(b)
#> [1] 1
barplot(b)

In contrast decreasing the channel velocity produces a more delayed histogram. In this case \(r_{L}>r_{0}+1\).

b <- fd(L=200,vch=0.5,Dt=200,tau0=0)
sum(b)
#> [1] 1
barplot(b)

Finally we look at two cases where \(\tau_0\) is positive; being either a direct multiple of the time step

b <- fd(L=600,vch=1,Dt=200,tau0=800)
sum(b)
#> [1] 1
barplot(b)

or not

b <- fd(L=600,vch=1,Dt=200,tau0=850)
sum(b)
#> [1] 1
barplot(b)