In WDNR-Water-Use/CSLSchem: CSLS Water Chemistry Analysis

library(CSLSchem)
library(CSLSevap)
library(NISTunits)
library(dplyr)
library(lubridate)
library(ggplot2)
library(zoo)
library(reshape2)
library(extrafont)

text_size   <- 12

# start_date <- as_datetime("2018-08-21")
# end_date   <- as_datetime("2019-11-01")
start_date <- as_datetime("2018-05-22")
end_date   <- as_datetime("2019-11-15")
dt         <- "day"
parameter  <- "PHOSPHORUS TOTAL"

res_time    <- 200 #days
C_ice       <- 0
C_evap      <- 0
vs          <- 0  # m/d
vs_tadj     <- 0.2
kso         <- 0.005 #0.07 #g/m2-d
kso_tadj    <- 0 #1

lake_weather <- process_weather(start_date, end_date, dt)
lake_levels  <- process_levels(start_date, end_date, dt)
lake_temp    <- process_lake_temp(start_date, end_date, dt)
lake_chem    <- process_chem(parameter,  start_date, end_date, dt)

lake_ice <- NULL
dates    <-  seq(start_date, end_date, by = "1 day")
for (i in 1:length(dates)) {
  this_ice <- data.frame(date = dates[i],
                         I_mm = calculate_ice_thickness(dates[i]))
  lake_ice <- rbind(lake_ice, this_ice)
}
lake_ice$C_ice <- C_ice

lake_daily <- lake_weather %>%
              full_join(lake_levels, by = c("lake", "date")) %>%
              full_join(lake_temp, by = c("lake", "date")) %>%
              full_join(lake_ice, by = "date")  %>%
              full_join(lake_chem, by = c("lake", "date"))

# if (parameter %in% c("d18O", "d2H")) {
#   dt <- "month"
#   lake_weather <- process_weather(start_date, end_date, dt)
#   lake_temp    <- process_lake_temp(start_date, end_date, dt)
#   lake_chem    <- process_chem(parameter,  start_date, end_date, dt)
#   lake_Cevap   <- lake_weather %>%
#                   full_join(lake_temp, by = c("lake", "date")) %>%
#                   full_join(lake_chem, by = c("lake", "date"))
#   lake_Cevap   <- lake_Cevap %>%
#                   mutate(C_evap = calculate_Cevap(atmp = .data$atmp_C, 
#                                                   ltmp = .data$ltmp_surf_C, 
#                                                   RH = .data$RH_pct, 
#                                                   Cpcpn = .data$C_pcpn, 
#                                                   Clake = .data$C_lake, 
#                                                   parameter = parameter),
#                          date = .data$date + days(14)) %>%
#                   select(.data$lake, .data$date, .data$C_evap)
#   lake_Cevap   <- interpolate_values(lake_Cevap, 
#                                      group_vars = "lake", 
#                                      val_var = "C_evap", 
#                                      start_date, 
#                                      end_date)
#   lake_daily <- lake_daily %>% 
#                 full_join(lake_Cevap, by = c("lake", "date"))
# } else {
#   lake_daily$C_evap <- C_evap
# }

if (parameter %in% c("d18O", "d2H")) {
  dt <- "day"
  lake_weather <- process_weather(start_date, end_date, dt)
  lake_temp    <- process_lake_temp(start_date, end_date, dt)
  lake_chem    <- process_chem(parameter,  start_date, end_date, dt)
  lake_Cevap   <- lake_weather %>%
                  full_join(lake_temp, by = c("lake", "date")) %>%
                  full_join(lake_chem, by = c("lake", "date"))
  lake_Cevap   <- lake_Cevap %>%
                  mutate(C_evap = calculate_Cevap(atmp = .data$atmp_C, 
                                                  ltmp = .data$ltmp_surf_C, 
                                                  RH = .data$RH_pct, 
                                                  Cpcpn = .data$C_pcpn, 
                                                  Clake = .data$C_lake, 
                                                  parameter = parameter)) %>%
                  select(.data$lake, .data$date, .data$C_evap)
  lake_daily <- lake_daily %>% 
                full_join(lake_Cevap, by = c("lake", "date"))
} else {
  lake_daily$C_evap <- C_evap
}

df <- lake_daily %>%
      group_by(.data$lake) %>%
      mutate(P_m3 = .data$P_mm*.data$area_m2/1000,
             E_m3 = .data$E_mm*.data$area_m2/1000,
             I_m3 = .data$I_mm*.data$area_m2/1000,
             GW_m3 = res_time/mean(.data$vol_m3)) %>%
      ungroup() %>%
      select(.data$lake, .data$date, .data$day, .data$atmp_C, 
             .data$RH_pct, .data$irr_factor, .data$ltmp_bot_C, 
             .data$ltmp_surf_C, .data$area_m2, .data$vol_m3, 
             .data$dV_m3, .data$P_mm, .data$E_mm, .data$I_mm, 
             .data$P_m3, .data$E_m3, .data$I_m3, .data$GW_m3, 
             .data$C_lake, .data$C_pcpn, .data$C_GWin, .data$C_evap, 
             .data$C_ice)

df$GWin_m3_calc <- (df$P_m3*(df$C_pcpn - df$C_lake) +
                      df$E_m3*(df$C_lake - df$C_evap) -
                      df$vol_m3*(df$C_lake - lag(df$C_lake)))/
                    (df$C_lake - df$C_GWin)

lake_model$first_order <- NA
lake_model$second_order <- NA
lake_model$C_lake_calc <- NA
for (i in 1:nrow(lake_model)) {
  if (i == 1) {
    lake_model$C_lake_calc[i] <- lake_model$C_lake[i]
  } else {
    lake_model$first_order[i] <- (vs*(1+vs_tadj)^(lake_model$ltmp_C[i]-20))*
                                 lake_model$C_lake_calc[i-1]*lake_model$area_m2[i]
    lake_model$second_order[i] <- (kso*(1+kso_tadj)^(lake_model$ltmp_C[i]-20))*
                                 lake_model$C_lake_calc[i-1]*lake_model$area_m2[i]
    lake_model$C_lake_calc[i] <- lake_model$C_lake_calc[i-1] + (1/lake_model$vol_m3[i])*
                                 (lake_model$GWin_m3[i]*(lake_model$C_GWin[i] - 
                                                          lake_model$C_lake_calc[i-1]) +
                                    lake_model$P_m3[i]*(lake_model$C_pcpn[i] - 
                                                          lake_model$C_lake_calc[i-1]) +
                                    lake_model$E_m3[i]*(lake_model$C_lake_calc[i-1] -                                                                   C_evap) +
                                    lake_model$Ifrm_m3[i]*(lake_model$C_lake_calc[i-1] -                                                                   C_ice) + 
                                    lake_model$Imelt_m3[i]*(C_ice - 
                                                            lake_model$C_lake_calc[i-1]) -
                                    lake_model$first_order[i] - 
                                    lake_model$second_order[i])

  }
}

ggplot(lake_model) + 
  geom_line(aes(x = date, 
                y = C_lake_calc,
                color = "Calculated")) +
  geom_line(aes(x = date, 
                y = C_lake,
                color = "Mean Measured")) +
  scale_x_datetime(date_breaks = "6 months",
                   date_minor_breaks = "2 months",
                   date_labels = "%b '%y") +
  scale_color_manual(name = "",
                     breaks = c("Calculated", "Mean Measured"),
                     values = c("black", "grey70")) +
  labs(x = "",
       y = "Lake Concentration (mg/L)", 
       title = "Pleasant Lake") +
  theme_bw() +
  theme(text = element_text(family = "Segoe UI Semilight",
                                        size = text_size),
        legend.position = "top",
        legend.margin = margin(0,0,0,0))