Impervious area calculations using the runoff coefficient model

Impervious_RoC.R

library(tidyverse)

rain_design <- 83.4
t_pattern <- c(4.1, 8.0, 11.0, 23.3, 15.3, 8.1, 7.0, 7.0, 5.1, 6.1, 3.1, 1.9)

if(!dplyr::near(sum(t_pattern), 100)) cat('The sum of the temporal pattern is not 100%')

t_step <-  0.5
t_step_count <- 1:length(t_pattern)
t_hour <- t_step_count * t_step
hyeto <- t_pattern * rain_design/100

hyeto
max(hyeto)
sum(hyeto)


hyeto_df <- data_frame(t_step_count, t_hour, rain = hyeto)
hyeto_df

# Runoff coefficent

# Function to reduce the rainfall using the runoff coefficient model
# Rainfall at each time step is multiplied by the RoC

Remove_RoC <- function(hyeto, RoC) {
  
  RoC * hyeto
}

Remove_RoC(hyeto, 0.9)

Remove_IL <- function(Hyeto, IL){
  
  # IL <- 1
  # Hyeto <- c(30, 7, 4, 3, 7, 20, 15, 3)
  # # 
  if(sum(Hyeto) < IL) {
    warning('Initial loss is not satisfied')
    Hyeto[] <- 0
    return(Hyeto)
  } 
  
  rain_cum <- cumsum(Hyeto)
  
  tzero <- min(which(rain_cum > IL)) - 1
  
  if(tzero == 0) {
    Hyeto[1] <- Hyeto[1] - IL
    return(Hyeto)
  }
  
  Hyeto[tzero + 1] <- Hyeto[tzero + 1] - (IL - rain_cum[tzero])
  Hyeto[1:tzero] <- 0
  Hyeto
  
}

Remove_CL <- function(Hyeto, CL, t_step) {
  
  CL_t_step = CL * t_step # Continuing loss per time step
  Hyeto <- Hyeto - CL_t_step
  Hyeto[Hyeto < 0] <- 0
  Hyeto
  
  
}


RoC <- 0.9
IL <-  0



# Example calculation

hyeto_df %>% 
  mutate(rain_IL = Remove_IL(rain, IL),
         rain_ILRoC = Remove_RoC(rain_IL, RoC)) %>% 
  select(rain_ILRoC) %>% 
  unlist(use.names = FALSE) # pull out the rainfall excess

# [1]  3.07746  6.00480  8.25660 17.48898 11.48418  6.07986  5.25420  5.25420  3.82806  4.57866  2.32686
# [12]  1.42614
# Rainfall in the 3rd time step is 19.4322 mm

(17.5 / 0.5) * 10 * 1/3.6

# Calculate and plot the rainfall excess hydrograph

RoC <- 0.9
IL <-  0
CL <-  0
sub_area <- 10

hyeto_df %>%
  mutate(rain_excess_RoC = Remove_RoC(Remove_IL(rain, IL), RoC),
         rain_excess_ILCL = Remove_CL(Remove_IL(rain, IL), CL, t_step)) %>% 
  add_row(t_step_count = 0,
          t_hour = 0,
          rain = 0,
          rain_excess_ILCL = 0,
          rain_excess_RoC = 0,
          .before = 1) %>% 
  add_row(t_step_count = 13,
          t_hour = 6.5,
          rain = 0,
          rain_excess_ILCL = 0,
          rain_excess_RoC = 0) %>% 
  mutate(flow_ILCL = (rain_excess_ILCL/t_step) * sub_area * 1/3.6,
         flow_RoC = (rain_excess_RoC/t_step) * sub_area * 1/3.6) %>% 
  ggplot(aes(x = t_hour)) + 
  geom_line(aes(y = flow_ILCL), alpha = 0.5, colour = 'black') +
  geom_point(aes(y = flow_ILCL), shape = 21, color = 'black', fill = 'grey', size = 3) +
  geom_line(aes(y = flow_RoC), alpha = 0.5, colour = 'blue') +
  geom_point(aes(y = flow_RoC), shape = 21, color = 'black', fill = 'blue', size = 3) +
  annotate(geom = 'rect', xmin = 3.1, xmax = 3.5, ymin = 100, ymax = 108, fill = 'grey', colour = 'black') +
  annotate(geom = 'text', x = 3.7, y = 104, hjust = 0, label = 'Initial loss / continuing loss' ) +
  annotate(geom = 'rect', xmin = 3.1, xmax = 3.5, ymin = 85, ymax = 93, fill = 'blue', colour = 'black') +
  annotate(geom = 'text', x = 3.7, y = 89, hjust = 0, label = 'Runoff coefficient' ) +
  theme_gray(base_size = 14) +
  scale_x_continuous(name = 'Time (hour)') +
  scale_y_continuous(name = 'Flow ('*m^-3~s^-1*')',
                     breaks = seq(0,110,10),
                     limits = c(0, 110))