Routing in RORB – II

The last post discussed routing in the RORB model and discussed an implicit approach to numerically solving the routing equation.  Here I look at an alternative approach.

As discussed previously the key mathematical aspects of RORB are the continuity equation and a relationship between storage and discharge. The continuity equation can be written as a differential equation:

$I - Q = \frac{dS}{dt}$           (1)

Where $I$ is inflow, $Q$ is outflow, $S$ is storage $t$ is time.

The relationship between storage and discharge used in RORB is:

$S = kQ^m$           (2)

Looking at the differential on the right hand side of (1)

$\frac{dS}{dt} = \frac{dS}{dQ} \frac{dQ}{dt}$

From (2)

$\frac{dS}{dQ} = k m Q^{m - 1}$           (3)

Combining (3) and (1)

$\frac{dQ}{dt} = \frac{I - Q}{k m Q^{m - 1}}$           (4)

This is just a standard differential equation
$\frac{dQ}{dt} = f(Q, t, parameters)$

The parameters are two constants, $k$ and $m$ and $I$, the inflow hydrograph, which depends on $t$.

In this form, as a differential equation, it is possible to use any standard numerical approach to develop a solution.

There are many FORTRAN and C libraries for solving differential equations and these are available to R in the deSolve package.  There is an example in the code below.

This is an ‘initial value’ problem so an initial value must be provided for Q.  I’ve set this to a small positive value.  Its best not to use zero because the derivative (equation 4) will be infinite if $m$ is less than 1.

The figure below compares DE solution with that from the implicit solution in the previous post.  The answers aren’t exactly the same but are so close it’s not possible to see any differences on the graph.

Code is below, and as a gist.

library(deSolve)
# dQ/dt = (I - Q)/(kmQ^(m-1))

# Specify the inflow hydrograph, using a function as discussed at

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)
}
# set up the DE solver
parameters <- c(k = 20, m = 0.8)
state <- c(Q = 0.01) # Initial value
times <- seq(1, 101, 1) # Need to start the time sequence at 1 to allow a valid comparsion
# with the solution from the routing function

Routing_ode <- function(t, state, parameters){
with(as.list(c(state, parameters)),{
# rate of change
Inflow <- Ihydrograph(t)
dQ <- (Inflow - Q)/(k*m*Q^(m-1))

list(dQ)

})
}

out_de <- ode(y = state, times = times, func = Routing_ode, parms = parameters)

plot(out)

# Compare to routing function

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

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)
Q[1] <- 0.01

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

}

out_routing <- Routing(Ihydrograph, k = 20, m = 0.8, T = 100, plotit = FALSE, return.values = TRUE)
out_routing

x <- data.frame(out_de = out_de[ ,2], out_routing = out_routing)
graphics::matplot(x, type = 'l', lty = c(2,1), col = c('black', 'blue'), las = 1, ylab = 'Flow')
legend('topright', legend = c('DE','Routing Function'), lty = c(2,1), col = c('black', 'blue'), bty = 'n')