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:

(1)

Where is inflow, is outflow, is storage is time.

The relationship between storage and discharge used in RORB is:

(2)

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

From (2)

(3)

Combining (3) and (1)

(4)

This is just a standard differential equation

The parameters are two constants, and and , the inflow hydrograph, which depends on .

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 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 # https://tonyladson.wordpress.com/2015/07/21/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) } # 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')