Runoff Routing based on Mein, R. G., E. M. Laurenson and T. A. McMahon (1974) see https://tonyladson.wordpress.com/2016/03/10/routing-in-rorb/

Routing.R

# Routing
# Based on discussion in 
# Mein, R. G., E. M. Laurenson and T. A. McMahon (1974). 
# "Simple nonlinear model for flood estimation." Journal of the Hydraulics Division, American Society of Civil Engineers 100(HY11): 1507-1518.


# Time stepping solution - ideas from  Section 23.2, p490 of
# Jones, O., Maillardet, R. and Robinson, A. (2014) Introduction to scientific programming and simulation using R.
# CRC Press



# Routing Function
# Inputs
#   Inflow hydrograph (as a function)
#   k - routing coefficient
#   m - routing parameter
#   T - Number of timesteps


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)

  # Continuing 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 # could keep track of storage if necessary
  
  }
  
  
  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)

}


# Specify an inflow hydrograph as a function
# See 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)
}

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