Reservoir level pool routing

level-pool-routing.R

lagpad <- function(x, k) {
    c(rep(NA, k), x)[1 : length(x)] 
}

level_pool_routing <- function(lt, qh, area, delta_t, initial_outflow, initial_storage, linear.fit){
  
  lt$ii <- apply(cbind(lagpad(lt$inflow, 1), lt$inflow), 1, sum)

  if(is.null(qh$storage)){
    qh$storage <- area * qh$elevation
  }
  qh$stq     <- ((2 * qh$storage) / (delta_t)) + qh$discharge
  
  lt$sjtminq   <- NA
  lt$sj1tplusq <- NA
  lt$outflow   <- NA
  lt[1, c("sj1tplusq")] <- c(NA)
  lt[1, c("sjtminq")] <- ((2 * initial_storage) / delta_t) - initial_outflow
  
  lt[1, "outflow"] <- initial_outflow
  
  if(linear.fit == TRUE){
    fit <- lm(discharge ~ stq, data = qh)
  }else{
    fit <- mgcv::gam(discharge ~ s(stq, k = 2), data = qh)
    
    # plot(qh$stq, qh$discharge)
    # lines(qh$stq, predict(fit))
  }
  
  for(i in seq_len(nrow(lt))[-1]){
    lt[i, "sj1tplusq"] <- lt[i-1, "sjtminq"] + lt[i, "ii"]
    lt[i, "outflow"]   <- predict(fit, data.frame(stq = lt[i, "sj1tplusq"]))
    lt[i, "sjtminq"]   <- lt[i, "sj1tplusq"] - (lt[i, "outflow"] * 2)
  }

  lt
}

lt <- data.frame(time = seq(0, 210, by = 10),
                 inflow = c(seq(0, 360, by = 60), seq(320, 0, by = -40), rep(0, 6)))

qh <- data.frame(elevation = seq(0, 10, by = 0.5),
                 discharge = c(0, 3, 8, 17, 30, 43, 60, 78, 97, 117, 137, 156, 173, 190,
                               205, 218, 231, 242, 253, 264, 275))

res_ppt <- level_pool_routing(lt, qh, area = 43560, delta_t = 600, initial_outflow = 0, initial_storage = 0, linear.fit = FALSE)

plot(res_ppt$time, res_ppt$inflow)
lines(res_ppt$time, res_ppt$outflow)