Calculate filling rule curve based on mSPA

library("pacman")
library(hydroTSM)
library(reshape2)
pacman::p_load(
  pacman,
  dplyr,
  GGally,
  ggplot2,
  ggthemes,
  ggvis,
  httr,
  lubridate,
  plotly,
  rio,
  rmarkdown,
  shiny,
  stringr,
  tidyr
)


data <- import("D:/INRAE-montpllier/Code/Test1/inputdata_b.txt")  #[Download here ](https://github.com/infohquan/link/blob/master/inputdata_b.txt)

time <-
  data.frame(seq(as.Date("1961/1/1"), as.Date("2007/12/31"), "month"))
month <-
  data.frame(format(
    as.Date(
      time$seq.as.Date..1961.1.1....as.Date..2007.12.31.....month..,
      format = "%Y-%m-%d"
    ),
    "%m"
  ))
year <-
  data.frame(format(
    as.Date(
      time$seq.as.Date..1961.1.1....as.Date..2007.12.31.....month..,
      format = "%Y-%m-%d"
    ),
    "%Y"
  ))
n <- NROW(na.omit(data))
Inflow <- data[, 1]
dem <- data[, 2]
prec <- data[, 3]
eva  <- data[, 4]
netev <- (data[, 4] - data[, 3]) / 1000
rule_curve <- c(34, 50, 63, 73, 77, 79, 80, 63, 46, 29, 12, 21)
DS <- rep(8.5, 12)
capacity <-  rep(82.5, 12)

# Calculate storage without net Evaporation Loss
Kt1 <- 0
Zt1 <- numeric(length = n)
Zt  <- numeric(length = n)
for (i in 1:n) {
  Kt1[i + 1] <- max(0, (Kt1[i] + dem[i] - Inflow[i]))
}
S <- max(Kt1)
for (i in 2:n + 1) {
  Zt1[i - 1] <- S - Kt1[i]
}
for (i in 2:n + 1) {
  Zt[i - 1] <- S - Kt1[i - 1]
}

# Estimate net evaporation using Piecewise Linear with 3 segments
b <- c(0.1105, 0.0514, 0.0043)
a <- c(0.138, 1.006 , 3.501)
s <- c(14.69, 52.91)

st_ini <- 82.5 - 8.5 # Assume storage is full
t1 <- 1

for (t in 1:n) {
  if (Zt1[t] >= s[2]) {
    ev <- netev * ((0.5 * b[3] * (Zt1 + Zt)) + a[3])
  } else if (Zt1[t] >= s[1]) {
    ev <- netev * ((0.5 * b[2] * (Zt1 + Zt)) + a[2])
  } else
    (Zt1[t] < s[1])
  ev <- netev * ((0.5 * b[1] * (Zt1 + Zt)) + a[1])
}

# Calculate storage with the incorporation of net evaporation
for (i in 1:n) {
  Kt1[i + 1] <- max(0, (Kt1[i] + dem[i] + ev[i] - Inflow[i]))
}
S <- max(Kt1)
for (i in 2:n + 1) {
  Zt1[i - 1] <- S - Kt1[i]
}
for (i in 2:n + 1) {
  Zt[i - 1] <- S - Kt1[i - 1]
}

stora <-
  data.frame(year, month, Zt) #  [ Storage = reservoir capacity - Deficit ]
names(stora)[1] <- "Year"
names(stora)[2] <- "Month"
ST <- reshape(
  data = stora,
  idvar = "Year",
  v.names = "Zt",
  timevar = "Month",
  direction = "wide"
)

storage = colMeans(ST[2:13])
St <- ST[2:13]
URC <- numeric(length = length(St[1]))
LRC <- numeric(length = length(St[1]))
SD <- numeric(length = length(St[1]))

for (ii in 1:12) {
  URC[ii] <- max(St[, ii])
  LRC[ii] <- min(St[, ii])
  SD[ii] <- sd(St[, ii])
}

URC_UPPER <- pmin(capacity, (URC + SD + DS))
URC_LOWER <- (URC - SD) + DS
LRC_UPPER <- (LRC + SD) + DS
LRC_LOWER <- pmax(DS, (LRC - SD) + DS)

Urc <-  (URC_UPPER + URC_LOWER) / 2
Lrc <-  (LRC_UPPER + LRC_LOWER) / 2

result <- data.frame(capacity, Urc, Lrc, storage, rule_curve, DS)
names(result)[1] <- "Capacity"
names(result)[4] <- "Estimated storage"
names(result)[5] <- "Filling curve"
names(result)[6] <- "Dead Storage"
result$Month <- c(1:12)
result$Month = factor(
  result$Month,
  labels = c(
    "Jan",
    "Feb",
    "Mar",
    "Apr",
    "May",
    "Jun",
    "Jul",
    "Aug",
    "Sep",
    "Oct",
    "Nov",
    "Dec"
  )
)
result.gather <-
  gather(result, key = "Rule Curve", value = IndexValue, -Month)

ggplot(result.gather,
       aes(Month, IndexValue, group = `Rule Curve`, color = `Rule Curve`)) +
  ylim(0, 85) +
  geom_line(aes(y = IndexValue), lwd = 1.2) +
  labs(title = "Rule Curve (1961-2007)",
       y = "Storage (MCM)",
       x = "") + theme_few(base_size = 20) +
  scale_color_brewer(palette = "Set2") +
  theme(plot.title = element_text(hjust = 0.5))



RESULT

`![alt text]`(/uploads/2485213e712e07d0f8c5966131ad9cee/Rplot.png)