Notation and nomenclature

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

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\).

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} \]

## 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)
}

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

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

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

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)
}

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

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

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

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

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