# Routing in RORB

I’ve been reviewing some of the background that explains the way RORB works, particularly Appendix A in the RORB manual and the 1974 paper by Mein, Laurenson and McMahon Simple Nonlinear Model for Flood Estimation.

The key mathematical aspects of RORB are the continuity equation and a relationship between storage and discharge.

For any storage, the continuity equation can be written as: “inflow in a time period, minus outflow in that time period equals change in storage”.

Or as used in RORB:

$\frac{I_t+I_{t+1}}{2}-\frac{Q_t+Q_{t+1}}{2} = \frac{S_{t+1}-S_t}{\Delta t}$        (1)

Where $I$ is inflow, $Q$ is outflow, $S$ is storage and $\Delta t$ is the time step, $t$ is one time period and $t+1$ is the next time period.

The relationship between storage and discharge used in RORB is:

$S = kQ^m$        (2)

This is justified in Appendix A of the RORB manual and in Mein et al. (1974).

Substituting equation (2) into (1):

$\frac{I_t+I_{t+1}}{2} + \frac{Q_t+Q_{t+1}}{2} = \frac{kQ_{t+1}^m-kQ_t^m}{\Delta t}$        (3)

In the general case, we are solving for the outflow $Q_{t+1}$ and know all the other parameters i.e. the inflow ($I_t$ and $I_{t+1}$), the outflow at the previous time step $Q_t$, along with $m$ and $k$ which are properties of a catchment and/or a waterway.

Equation (3) can be rearranged as follows

$Q_t + Q_{t+1} - I_t - I_{t+1} + \frac{2k}{\Delta t} \left( Q_{t+1}^m-Q_t^m \right) = 0$        (4)

I’ve written that equation (4) is equal to zero, but this will only be the case if we have the correct value of $Q_{t+1}$.

This opens the possibility of a numerical solution; guess a value of $Q_{t+1}$ and plug it into equation (4), if that doesn’t equal zero,  guess again. In applied mathematics, this is referred to as root finding i.e. finding where a function equals zero. There are a number of procedures to do this efficiently (see for example, Numerical Recipes).

Let’s start with the Newton-Raphson procedure which the method discussed in Mein et al. (1974)

Think of equation (4) as a function of $Q_{t+1}$ i.e.

$f \left( Q_{t+1}\right) = Q_t + Q_{t+1} - I_t - I_{t+1} + \frac{2k}{\Delta t} \left( Q_{t+1}^m-Q_t^m \right)$      (5)

We guess a value of $Q_{t+1}$ say $Q_{t+1.1}$ and would like to calculate a better guess that brings $f \left( Q_{t+1}\right)$ closer to zero. Consider the derivative of $f \left( Q_{t+1}\right)$ at our initial guess, this is just the slope of the tangent at $Q_{t+1.1}$ as shown in the figure below which hits the axes at a point closer to the required root so would be a good next guess.  Call this point$Q_{t+1.2}$ .

The slope is rise/run so in this case.

$f' \left( Q_{t+1}\right) = \frac{ f \left( Q_{t+1}\right)}{Q_{t+1.1} - Q_{t+1.2}}$        (6)

Now, rearrange (6) to give an expression for $Q_{t+1.2}$

$Q_{t+1.2} = Q_{t+1.1} - \frac{f \left( Q_{t+1}\right)}{f' \left( Q_{t+1}\right)}$      (7)

Next we need the derivative of equation (5) with respect to $Q_{t+1}$

$f' = 1 + \frac{2km}{\Delta t}Q^{m-1}_{t+1}$       (8)

Using this approach to find Q_{t+1.2}, substitute (8) and (5) into (7) to give.

$Q_{t+1.2} = Q_{t+1.1} - \frac{\Delta t \left(Q_{t+1.1} + Q_t - I_{t+1} - I_t \right) +2k \left( Q^m_{t+1.1} - Q^m_t \right)}{\Delta t + 2kmQ^{m-1}_{t+1.1}}$      (9)

Equation (9) is the same as equation (8) in Mein et al. (1974).

The iterative solution proceeds as follows.

1. Select an initial value of $Q$ for the first timestep $Q _{t}$,
2. Guess $Q _{t+1.1}$ (Mein et al. use $(I_{t+1} + Q_t)/2$
3. Use equation (9) to calculate $Q_{t+1.2}$
4. Set $Q_{t+1.1} = Q_{t+1.2}$
5. Calculate a new $Q_{t+1.2}$
6. Repeat steps 3 to 5 until $Q_{t+1.2}$ has a similar value to $Q_{t+1.1}$ or, equivalently, until equation (4) does in fact equal zero (or is nearly equal to zero).
7. Go to the next timestep
8. Set $Q _{t}$ (current timestep) equal to  $Q_{t+1}$ (from the previous timestep)
9. Repeat steps 2 to 8 for all time steps

This sounds straightforward but there can be numerical issues.  The initial guess needs to be quite close to the answer or the next guess provided by equation (9) could be a negative flow.  Putting this back into equation (5) will cause an error because it requires raising a negative number to a fractional power. Another challenge is that the derivative (equation 8) is undefined when $Q^m_{t+1}$ is zero if $m < 1$.

When reading the RORB manual, section 2.2.5 there is a comment that: “The Regular Falsi convergence algorithm was used”.  There is a good description of this root finding approach in Wikipedia.  The upshot is that RORB does not use Newton-Raphson as was discussed in Mein et al. 1974.

In R, the uniroot function provides a reliable way of root finding without the numerical issues of the Newton-Raphson procedure. uniroot is based on a C function zeroin as documented here.

A simple routing function, in R, is below.

Routing <- function(Ihydrograph, k, m, T, dt = 1, plotit = TRUE, return.values = FALSE){

# Ihydrograph is a function that describes the inflow hydrograph
# k and m are the routing parameters
# dt is the time step

I <- rep(0, T+1) # Initialise input vector
Q <- rep(0, T+1) # Initialise output vector
S <- rep(0, T+1) # Initialise storage vector

I[1] <- Ihydrograph(1) # initialise inflow

# Continuity equation that we seek to set to zero to determine Q at the next time step

myF <- function(Qip1, Qi, Ii, Iip1, k, m, dt){

Qip1 + Qi - Iip1 - Ii + (2*k/dt)*(Qip1^m - Qi^m)

}

for(i in 1:100){
I[i+1] <- Ihydrograph(i+1)
Q[i+1] <- uniroot(myF, interval = c(0,10), Qi = Q[i], Ii = I[i], Iip1 = I[i+1], k, m, dt)\$root
S[i+1] <- k*Q[i+1]^m

}

if(plotit){
plot(I, type = 'l', col = 'blue', las = 1, ylab = 'flow', xlab = 'Time step') # Inflow hydrograph
lines(Q) # outflow hydrograph
legend('topright', legend = c('Inflow', 'Outflow'), lty = 1, col = c('blue', 'black'), bty = 'n' )
}

if(return.values) return(Q)

}


Let’s test it using a hydrograph as a function (see this blog for details).

What about the influence of k and m on routing?

Influence of the m parameter on routing

Influence of the k parameter on routing

RORB starts with a rainfall excess hyetograph (the rainfall that is left over after losses have been subtracted) and then routes this through a series of elements to produce flow at the catchment outlet.  Here I’ve modelled a single element but that gives an indication of how the program works.

Code for the routing function and to produce these figures (also available as a gist).

Code to produce these figures.
# Specify an inflow hydrograph as a function

my.Qmin <- 0
my.Qmax <- 1
my.beta <- 10
my.t.peak <- 20

Hydro <- function(tt, t.peak = 3, Qmin = 1, Qmax = 2, beta = 5) {
Qmin + (Qmax - Qmin)*( (tt/t.peak) * (exp(1 - tt/t.peak)))^beta

}
# Initialise the inflow
Ihydrograph <- function(tt){
Hydro(tt, t.peak = my.t.peak, Qmin = my.Qmin, Qmax = my.Qmax, beta = my.beta)
}

# Routing
Routing(Ihydrograph, k = 20, m = 1, T = 100, plotit = TRUE, return.values = FALSE)
# Step function hydrograph (or possibly based on measured inflows)

inflow <- c(0,0, 1,1,1,1,1, rep(0, times = 94))
Ihydrograph <- approxfun(1:101, inflow)
Routing(Ihydrograph, k = 20, m = 1, T = 100, plotit = TRUE, return.values = FALSE)

############################################
# Influence of m
# Initialise the inflow

my.Qmin = 0
my.Qmax = 1
my.beta = 10
my.t.peak = 20

Hydro <- function(tt, t.peak = 3, Qmin = 1, Qmax = 2, beta = 5) {
Qmin + (Qmax - Qmin)*( (tt/t.peak) * (exp(1 - tt/t.peak)))^beta

}
Ihydrograph <- function(tt){
Hydro(tt, t.peak = my.t.peak, Qmin = my.Qmin, Qmax = my.Qmax, beta = my.beta)
}

# Wrap Routing as a function of m

Routing_m <- function(m, Ihydrograph, k, T, dt = 1, plotit = TRUE, return.values = FALSE){
Routing(Ihydrograph, k, m, T, dt = 1, plotit = plotit, return.values = return.values )
}
m.seq <- seq(0.3, 1, 0.1)
x <- lapply(m.seq, FUN = Routing_m, Ihydrograph = Ihydrograph, k = 20, T = 100, dt = 1, plotit = FALSE, return.values = TRUE)
x1 <- do.call(cbind, x)

plot(Ihydrograph(1:100), type = 'l', col = 'blue', las = 1, ylab = 'flow', xlab = 'Time step') # Inflow hydrograph
matlines(x1)
legend('topright', legend = m.seq, col = 1:6, lty = 1:5, bty = 'n')

############################################
# Influence of k

# Wrap Routing as a function of k

Routing_k <- function(k, Ihydrograph, m, T, dt = 1, plotit = TRUE, return.values = FALSE){
Routing(Ihydrograph, k, m, T, dt = 1, plotit = plotit, return.values = return.values )
}
k.seq <- c(1, 5, 10, 15, 20, 50)
x <- lapply(k.seq, FUN = Routing_k, Ihydrograph = Ihydrograph, m = 0.8, T = 100, dt = 1, plotit = FALSE, return.values = TRUE)
x1 <- do.call(cbind, x)

plot(Ihydrograph(1:100), type = 'l', col = 'blue', las = 1, ylab = 'flow', xlab = 'Time step') # Inflow hydrograph
matlines(x1)
legend('topright', legend = k.seq, col = 1:6, lty = 1:5, bty = 'n' )



## One thought on “Routing in RORB”

1. Pingback: Routing in RORB – II | tonyladson