Solving the routing equation using a differential equation solver. See https://tonyladson.wordpress.com/2016/03/15/routing-in-rorb-ii/

Routing_DE.R

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