Document ID: EPA-HQ-OW-2023-0222-0255
Agency: epa
Document Type: Supporting & Related Material
Title: 
Posted Date: 2023-12-22T05:00Z

# ----- Code Purpose and Authors -----
# Purpose: Perform WQI calculations for the Delaware River Dissolved Oxygen Water Quality Standards
# Authors: Alyssa Le and Kate Munson
# Date created: 03/31/23
# Date modified: 05/31/23. Isabelle Morin adjusted variables to point to the user-specific folders and read in the correct data
# Version 5: Added additional scenarios combining unadjusted DO and estimated BOD changes
# 6/12/23: Modified the code to update the WQS values
# Version 6: Added 2012 and 2018 calculations (6/27/23)
# Version 7: Corrected the sensitivity scenario #6 to include adjusted BOD data (7/12/23) 
# Later also added code to recover some BOD values below the detection limit to use 1/2 detection limit (7/13/23) 
# version 8: Changed median to 6.1 mg/L as per EPA email (7/17/23).
# version 9: Added code to also use NARS-derived scores (8/8/23)
# Clean version for posting to docket: deleted old outdated comments and tested that the code runs through (11/6/23)
# version 10: Updated proposed WQ standard as per email dated 12/4/23 (10th percentile changed from 5.3 to 5.4 mg/L)

rm(list=ls())

options(stringsAsFactors = FALSE)
pacman::p_load(readxl, dplyr, tidyr, zoo, stringr, reshape2, stringi, sf, tidyverse, openxlsx)

# VARIABLES TO ADJUST -------------------------------------------------------------

      data_path <-"C:/Users/39304/ICF/CWA Anniversary & Econ Support - TD3_Delaware River WQS/Water Quality Data/"
      WQI_path <- "C:/Users/39304/ICF/CWA Anniversary & Econ Support - TD3_Delaware River WQS/WQI/R_programs/"
      gis_path <- "C:/Users/39304/ICF/CWA Anniversary & Econ Support - TD3_Delaware River WQS/GIS/Data/"
      eco_path <- "C:/Users/39304/ICF/CWA Anniversary & Econ Support - TD3_Delaware River WQS/GIS/Data/us_eco_l3/"
      
      Scen_RunDate <- "23-12-04"
      Scen_Year <- "2019"       
      Scen_Flow <- "Design"     # Options are: "Actual" or "Design"
      
      Scen_Name <- paste0(Scen_Year,"_",Scen_Flow) 
      
      if (Scen_Flow == "Actual") {
          # Baseline Future with actual flows and full LTCP implementation
          # HADO Actual Flows 2019
        Scen_Baseline <- paste0("Baseline_Actual_LTCP_",Scen_Year,"_")  #"Baseline_Actual_LTCP_2019_"
        Scen_Policy1 <- paste0("HADO_Actual_",Scen_Year,"_")   #"HADO_Actual_2019_"
        } else {
          # Baseline Design (Max) Flows inc LTCPs 2019
          # HADO Design (Max) Flows inc LTCPs 2019
          Scen_Baseline <- paste0("LTCP_Design_",Scen_Year,"_") #"LTCP_Design_2019_"
          Scen_Policy1 <- paste0("HADO_Design_",Scen_Year,"_")  #"HADO_Design_2019_"
        }
      
      # Subset data to summer months (July 1 to October 31)
      start_date <- paste0(Scen_Year,"-07-01")  # for 2019, we get "2019-07-01"
      end_date <- paste0(Scen_Year,"-10-31")    # for 2019, we get "2019-10-31"
      
      # Median of lowest 10% of DO values
      Med_10 <- 5.4
      
      # Median of time series of DO values
      Med_50 <- 6.1
      
      # Average BOD change at 12 sites with ammonia reduction to 1.5 mg/L, based on Table 7-2 of the Kleinfelder Nitrogen Reduction Report
      # https://www.nj.gov/drbc/library/documents/NitrogenReductionCostEstimates_KleinfelderJan2021.pdf
      BOD_PercCh <- -0.48

# GENERAL INPUT DATA -------------------------------------------------------------

      # Determine zones of interest to WQI calculations
      
      # Zone-to-model cell crosswalk (based on WASP_90percent_2019.DO.kavg.Rdata file delivered to EPA from Jim Hagy 2/27/2023)
      
      zones <- read.csv(paste0(data_path,"ICF_Rdata_Outputs/ij.csv"), header = TRUE)
      zones$IJ <- with(zones, paste0(I,"_",J))
      
      # The “fish maintenance area” (FMA) includes zones 3, 4, and upper 5
      
      fma_zones <- as.data.frame(unique(zones$IJ[zones$Zone == "3"|zones$Zone == "4"|zones$Zone == "5-upper"]))
      names(fma_zones) <- "IJ"
      
      length(unique(zones$IJ)) #1701 
      length(unique(fma_zones$IJ)) #579
      
        # QA - confirm that we've correctly selected the cells corresponding to the three zones
        fma_zones_QA <- zones %>%
          filter(Zone %in% c("3", "4", "5-upper")) %>%
          distinct(IJ)
      
      # FMA boundary and ecoregion data for subindex specification
      
      FMA_boundary <- st_read(paste0(gis_path, "Fish Maintenance Area/delriv_fma.shp")) 
      
      Ecoregions <- st_read(paste0(eco_path, "us_eco_l3.shp")) 

# BASELINE DATA -------------------------------------------------------------------

      #  * Modeled data -------------
      # Baseline scenario, see top of the file for specifications. Depth-averaged DO, TN, and TP concentrations
      
      model_param <- c("DO","TN","TP")
      
      for(i in 1:length(model_param))
      {
        param <- model_param[i]
        
        temp_bsln_mod <- read.csv(paste0(data_path,"Water-Quality Data/",Scen_Baseline,param,"_kavg.csv"), header = TRUE)
        
        bsln_mod_FMA <- merge(temp_bsln_mod[temp_bsln_mod$IJ %in% fma_zones$IJ,],zones[,c("IJ","Zone")],all.x=TRUE)
        
        melt_bsln_mod_FMA <- melt(bsln_mod_FMA,id = c("IJ","Zone"), variable.name = "datetime",value.name = "kavg")
        
        melt_bsln_mod_FMA$Date <- with(melt_bsln_mod_FMA,as.Date(paste0(substr(datetime,8,11),"/",substr(datetime,2,3),"/",substr(datetime,5,6))), "%m/%d/%Y")
        
        # Average to get daily values for each model cell
        bsln_mod_FMA_daily <- melt_bsln_mod_FMA[melt_bsln_mod_FMA$Date <= end_date & melt_bsln_mod_FMA$Date >= start_date,] %>%
          select(c("IJ","Zone","kavg","Date")) %>%
          group_by(IJ,Zone,Date) %>%
          summarise_all(mean)
        
        bsln_mod_FMA_daily$Param <- param
        
        if (exists("bind_bsln_mod_FMA")) {
          bind_bsln_mod_FMA <- rbind(bind_bsln_mod_FMA, bsln_mod_FMA_daily)
        } else {
          bind_bsln_mod_FMA <- bsln_mod_FMA_daily
        }
        
      }
      
      names(bind_bsln_mod_FMA) <- c("IJ","Zone","Date","Value","Param")

# * Observed data -------------
      # Physical chemistry results downloaded from the USGS water quality portal for DBRC boat run sites: CBOD, FC, TSS, NH3
      # (Delivered from EPA to ICF 3/20/2023)
      
      # br = boat run
      # obs = observed
      
      br_zones <- read_excel(paste0(data_path, "Water-Quality Data/G7pt2_3D_LookUp_Tables.xlsx"),
                             sheet = "BoatRun", 
                             col_types  = c(rep("skip",2),"text",rep("skip",5), rep("text", 3), "skip"))
      
      names(br_zones) <- c("MonitoringLocationIdentifier","Zone","I","J")
      
      br_zones$IJ <- with(br_zones, paste0(I,"_",J))
      
      temp_br_obs <- read.csv(paste0(data_path,"Water-Quality Data/resultphyschem.csv"), header = TRUE)
      
      br_obs <-
        temp_br_obs[temp_br_obs$CharacteristicName %in% c(
          "Carbonaceous biochemical oxygen demand, standard conditions",
          "Fecal Coliform",
          "Total suspended solids",
          "Ammonia-nitrogen",
          "Ammonia"
        ),]
      
      br_obs$ResultMeasureValue <- as.numeric(br_obs$ResultMeasureValue)
      br_obs$Date <- as.Date(br_obs$ActivityStartDate, "%m/%d/%Y")
      
      for(i in 1:nrow(br_obs)){
        if((br_obs$ResultDetectionConditionText[i] == "Not Detected"))
          #if((br_obs$ResultDetectionConditionText[i] == "Not Detected")&&(br_obs$CharacteristicName[i] == "Carbonaceous biochemical oxygen demand, standard conditions"))
          {
          br_obs$ResultMeasureValue[i] = as.numeric(br_obs$DetectionQuantitationLimitMeasure.MeasureValue[i])/2
        }
      }
      
      # --- QA Start --- #
      
      # Consider filters to include/exclude
      
      br_obs_summer <- br_obs[(br_obs$Date <= end_date & br_obs$Date >= start_date) & !is.na(br_obs$ResultMeasureValue),]
      
      unique(br_obs_summer$ActivityTypeCode)

      
      # [1] "Sample-Routine"                         "Quality Control Sample-Field Replicate" "Quality Control Sample-Equipment Blank"
      # [4] "Quality Control Sample-Field Blank"  
      
      # View(br_obs_summer[!(br_obs_summer$ActivityTypeCode == "Sample-Routine"),])
      # View(br_obs_summer[!(br_obs_summer$MonitoringLocationIdentifier == "31DELRBC_WQX-091017"),])
      #@KM: Usually, the location has routine samples AND QC samples. We are OK to remove QCs. 
      
      #@KM: Decision - select "Sample-Routine" only
      
      unique(br_obs_summer$ActivityMediaName)
      # [1] "Water"
      unique(br_obs_summer$ActivityMediaSubdivisionName)
      # [1] "Surface Water"
      
      unique(br_obs_summer$ResultValueTypeName)
      # [1] "Actual"    "Estimated"
      
      nrow(br_obs_summer[br_obs_summer$ResultValueTypeName == "Actual",]) # 303
      nrow(br_obs_summer[br_obs_summer$ResultValueTypeName == "Estimated",]) # 33
      unique(br_obs_summer$CharacteristicName[br_obs_summer$ResultValueTypeName == "Estimated"]) # [1] "Ammonia-nitrogen"

      # --- QA End --- #
      
      # Filter data to routine samples taken during the summer months
      
      # Note: Both ResultValueTypeName == "Actual" and ResultValueTypeName == "Estimated" data retained for Ammonia-nitrogen
      # because estimated data are a small portion of the available data and certain boat run sites (e.g., 31DELRBC_WQX-091002)
      # have no "actual" Ammonia-nitrogen data (only have estimated). We need estimated Ammonia-nitrogen data to estimate BOD below.
      
      bsln_br_obs <-
        merge(br_obs[(br_obs$Date <= end_date &
                        br_obs$Date >= start_date) &
                       !is.na(br_obs$ResultMeasureValue) &
                       (br_obs$ActivityTypeCode %in% "Sample-Routine")
                     ,], 
              br_zones[, c("MonitoringLocationIdentifier", "IJ")], all.x =
                TRUE) %>%
        select(
          c(
            "MonitoringLocationIdentifier",
            "CharacteristicName",
            "ResultMeasure.MeasureUnitCode",
            "ResultMeasureValue",
            "Date",
            "IJ"
          )
        ) %>%
        group_by(
          MonitoringLocationIdentifier,
          CharacteristicName,
          ResultMeasure.MeasureUnitCode,
          Date,
          IJ
        ) %>%
        summarise_all(mean, na.rm = TRUE)
      
      
      temp_bsln_obs_FMA <- merge(bsln_br_obs[bsln_br_obs$IJ %in% fma_zones$IJ,],zones[,c("IJ","Zone")],all.x=TRUE)
      
      # Calculate BOD based on nitrogenous BOD and carbonaceous BOD relationships found in the literature
      
      BOD_bsln_obs_FMA <- merge(temp_bsln_obs_FMA[temp_bsln_obs_FMA$CharacteristicName %in% "Ammonia-nitrogen",],
                                temp_bsln_obs_FMA[temp_bsln_obs_FMA$CharacteristicName %in% "Carbonaceous biochemical oxygen demand, standard conditions",],
                                by=c("IJ","MonitoringLocationIdentifier","ResultMeasure.MeasureUnitCode","Date","Zone"),all=TRUE)
      
      
      names(BOD_bsln_obs_FMA) <-
        c(
          "IJ",
          "MonitoringLocationIdentifier",
          "ResultMeasure.MeasureUnitCode",
          "Date",
          "Zone",
          "Ammonia",
          "Value_Ammonia",
          "CBOD",
          "Value_CBOD"
        )
      
      # BOD = NBOD + CBOD, where NBOD = 4.57*Ammonia
      
      BOD_bsln_obs_FMA$ResultMeasureValue <- with(BOD_bsln_obs_FMA,Value_Ammonia*4.57 + Value_CBOD)
      BOD_bsln_obs_FMA$Param <- "BOD"
      
      # Confirm at least 1 BOD value per monitoring location
      
      length(unique(BOD_bsln_obs_FMA$MonitoringLocationIdentifier)) # 8 
      length(unique(BOD_bsln_obs_FMA$MonitoringLocationIdentifier[!is.na(BOD_bsln_obs_FMA$ResultMeasureValue)])) # 8 
      
      temp_bsln_obs_FMA$Param <- with(temp_bsln_obs_FMA,
                                      ifelse(CharacteristicName == "Ammonia-nitrogen","NH3",
                                             ifelse(CharacteristicName == "Carbonaceous biochemical oxygen demand, standard conditions","CBOD",
                                                    ifelse(CharacteristicName == "Fecal Coliform","FC",
                                                           ifelse(CharacteristicName == "Total suspended solids","TSS","CHECK"
                                                           )))))
      
      bind_bsln_obs_FMA <- rbind(
        temp_bsln_obs_FMA[temp_bsln_obs_FMA$CharacteristicName %in% c("Fecal Coliform", "Total suspended solids"),
                          c("IJ","Zone","Date","ResultMeasureValue","Param")],
        BOD_bsln_obs_FMA[!is.na(BOD_bsln_obs_FMA$ResultMeasureValue),
                         c("IJ","Zone","Date","ResultMeasureValue","Param")]
      )
      
      names(bind_bsln_obs_FMA) <- c("IJ","Zone","Date","Value","Param")
      
      #  * Combined modeled and observed data -------------
      
      # Bind the modeled and observed baseline data together
      
      bind_bsln_FMA <- rbind(bind_bsln_mod_FMA,bind_bsln_obs_FMA)
      
      # --- QA Start --- #
      
      # Determine how many data points and stations per param. Determine range in dates of data.
      
      param_counts <- bind_bsln_FMA %>% 
        select(c("Param","IJ","Date")) %>% 
        group_by(Param, Zone) %>% 
        summarise(Stations = length(unique(IJ)),
                  DataPoints = n(),
                  MinDate = min(Date),
                  MaxDate = max(Date))
      
      write.csv(
        param_counts,
        paste0(
          WQI_path,"Outputs/QA/",Scen_Name,"DERiver_ObsParamCounts_",Scen_RunDate,".csv"
        ),
        row.names = F
      )
      
      # --- QA End --- #

# SCENARIO DATA -------------------------------------------------------------------

      # * Modeled data -------------
      # HADO scenario, see top of the file for specifications
      
      for(i in 1:length(model_param))
      {
        param <- model_param[i]
        
        temp_scen_mod <- read.csv(paste0(data_path,"Water-Quality Data/",Scen_Policy1,param,"_kavg.csv"), header = TRUE)
        
        scen_mod_FMA <- merge(temp_scen_mod[temp_scen_mod$IJ %in% fma_zones$IJ,],zones[,c("IJ","Zone")],all.x=TRUE)
        
        melt_scen_mod_FMA <- melt(scen_mod_FMA,id = c("IJ","Zone"), variable.name = "datetime",value.name = "kavg")
        
        melt_scen_mod_FMA$Date <- with(melt_scen_mod_FMA,as.Date(paste0(substr(datetime,8,11),"/",substr(datetime,2,3),"/",substr(datetime,5,6))), "%m/%d/%Y")
        
        # Average to get daily values for each model cell
        
        scen_mod_FMA_daily <- melt_scen_mod_FMA[melt_scen_mod_FMA$Date <= end_date & melt_scen_mod_FMA$Date >= start_date,] %>%
          select(c("IJ","Zone","kavg","Date")) %>%
          group_by(IJ,Zone,Date) %>%
          summarise_all(mean)
        
        scen_mod_FMA_daily$Param <- param
        
        if (exists("bind_scen_mod_FMA")) {
          bind_scen_mod_FMA <- rbind(bind_scen_mod_FMA, scen_mod_FMA_daily)
        } else {
          bind_scen_mod_FMA <- scen_mod_FMA_daily
        }
        
      }
      
      names(bind_scen_mod_FMA) <- c("IJ","Zone","Date","Value","Param")
      
      # --- QA Start --- #
      
      # Compare HADO model TN and TP to calibration model TN and TP
      
      Cal_HADO_TN_TP_DO <- merge(bind_bsln_mod_FMA[bind_bsln_mod_FMA$Param %in% c("TN","TP","DO"),],
                              bind_scen_mod_FMA[bind_scen_mod_FMA$Param %in% c("TN","TP","DO"),],
                              by=c("IJ","Zone","Date","Param"),all.x=TRUE)
      
      names(Cal_HADO_TN_TP_DO) <- c("IJ","Zone","Date","Param","Value_Calib","Value_HADO")
      
      Cal_HADO_TN_TP_DO$HADO_Minus_Calib <- with(Cal_HADO_TN_TP_DO,Value_HADO-Value_Calib)
      
      min(Cal_HADO_TN_TP_DO$HADO_Minus_Calib[Cal_HADO_TN_TP_DO$Param == "TP"]) # 2019 design: -9.899167e-05 2019 actual: -6.128542e-05
      min(Cal_HADO_TN_TP_DO$HADO_Minus_Calib[Cal_HADO_TN_TP_DO$Param == "TN"]) # 2019 design: -0.0003863333 2019 actual: -0.000873875
      min(Cal_HADO_TN_TP_DO$HADO_Minus_Calib[Cal_HADO_TN_TP_DO$Param == "DO"]) # 2019 design: 0.1011168     2019 actual: 0.07712611
      
      max(Cal_HADO_TN_TP_DO$HADO_Minus_Calib[Cal_HADO_TN_TP_DO$Param == "TP"]) # 2019 design: 1.080556e-06  2019 actual: 9.958333e-07
      max(Cal_HADO_TN_TP_DO$HADO_Minus_Calib[Cal_HADO_TN_TP_DO$Param == "TN"]) # 2019 design: 0.07086509    2019 actual: 0.04614868
      max(Cal_HADO_TN_TP_DO$HADO_Minus_Calib[Cal_HADO_TN_TP_DO$Param == "DO"]) # 2019 design: 2.714549      2019 actual: 1.964259
      
        # QA - what do the changes in TN and TP look like on average?
        Cal_HADO_TN_TP_QA <- Cal_HADO_TN_TP_DO %>%
          filter(Param %in% c("TN", "TP", "DO")) %>%
          group_by(Param) %>%
          summarise(Avg_Change = mean(HADO_Minus_Calib, na.rm=T))
      
      # How often are DO, TN, and TP increasing and decreasing?
      nrow(Cal_HADO_TN_TP_DO[Cal_HADO_TN_TP_DO$Param == "DO" & Cal_HADO_TN_TP_DO$HADO_Minus_Calib <0,]) # 0
      nrow(Cal_HADO_TN_TP_DO[Cal_HADO_TN_TP_DO$Param == "DO" & Cal_HADO_TN_TP_DO$HADO_Minus_Calib >0,]) # 71217 - all cells show DO increasing (improvement)
      nrow(Cal_HADO_TN_TP_DO[Cal_HADO_TN_TP_DO$Param == "TN" & Cal_HADO_TN_TP_DO$HADO_Minus_Calib <0,]) # 23
      nrow(Cal_HADO_TN_TP_DO[Cal_HADO_TN_TP_DO$Param == "TN" & Cal_HADO_TN_TP_DO$HADO_Minus_Calib >0,]) # 71194 - mostly increases, degradation
      nrow(Cal_HADO_TN_TP_DO[Cal_HADO_TN_TP_DO$Param == "TP" & Cal_HADO_TN_TP_DO$HADO_Minus_Calib <0,]) # 71054
      nrow(Cal_HADO_TN_TP_DO[Cal_HADO_TN_TP_DO$Param == "TP" & Cal_HADO_TN_TP_DO$HADO_Minus_Calib >0,]) # 162 - mostly decreases (improvement)
      
      # --- QA End --- #

#  * Apply the DO standard to the HADO scenario DO data -------------

      HADO_DO <- bind_scen_mod_FMA[bind_scen_mod_FMA$Param == "DO",]
      
      # Write out unmodified HADO values
      HADO_DO_AsIs <- HADO_DO[,c("IJ","Zone","Date","Value","Param")]
      
      HADO_DO$Quantile = rank(HADO_DO$Value)/nrow(HADO_DO)
      
      quantile(HADO_DO$Value, probs = c(0.1, .5))
      
      summary_HADO_DO_AsIs <- HADO_DO %>% 
        group_by(Param) %>%
        summarise(n=n(),pct10=quantile(HADO_DO$Value,0.1),pct50=median(HADO_DO$Value))
      
      #               10%      50% 
      #               5.133340 6.224399   
      # 2019 design  5.194894 6.174386 
      # 2019 actual  5.269151 6.307272 
      
      # Calculate median over the entire time series
      median(HADO_DO$Value)
      # [1] 2019 design: 6.174386    2019 actual: 6.307272
      
      # # Confirm that median of lowest 10% of DO values is < Med_10
      HADO_DO_10 <- min(HADO_DO$Value[HADO_DO$Quantile>=0.1])
      # 
      # median(HADO_DO$Value[HADO_DO$Value < HADO_DO_10])
      # # [1] 2019 actual: 5.11333  2019 design: 5.035158
      # 
      # median(HADO_DO$Value[HADO_DO$Value < HADO_DO_10]) < Med_10
      # # [1] TRUE   
      
      HADO_DO_90 <- HADO_DO[HADO_DO$Value >= HADO_DO_10,]
      HADO_DO_90$Value_Rev10 <- with(HADO_DO_90,
                                     ifelse(
                                       Value <= Med_10, Med_10, Value
                                     ))
      
      # Calculate what DO would be if applying the water quality standard (WQS)
      
      temp_HADO_DO_WQS <- merge(HADO_DO,HADO_DO_90,all.x=TRUE)
      temp_HADO_DO_WQS$Value_Rev10 <- with(temp_HADO_DO_WQS,
                                        ifelse(
                                          is.na(Value_Rev10), Value, Value_Rev10
                                        ))
      temp_HADO_DO_WQS$Value_Diff10 <- temp_HADO_DO_WQS$Value_Rev10 -temp_HADO_DO_WQS$Value
      
      # Recalculate median with revision in Step 1
      median(temp_HADO_DO_WQS$Value_Rev10)
      # 2019 design: 6.174386
      # 2019 actual: 6.307272
      
      # Confirm that median of time series is < Med_50
      median(temp_HADO_DO_WQS$Value_Rev10) < Med_50
      # [1] 2019 actual: FALSE  2019 design: TRUE
      
      if (median(temp_HADO_DO_WQS$Value_Rev10) < Med_50) {
      
          temp_HADO_DO_WQS$Quantile = rank(temp_HADO_DO_WQS$Value_Rev10)/nrow(temp_HADO_DO_WQS)
          
          # IM: Appropriate thresholds vary and should be as small as possible so that there are no values within the range.
          if (!is_empty(temp_HADO_DO_WQS$Quantile[temp_HADO_DO_WQS$Value_Rev10 > (Med_50-0.00001) & temp_HADO_DO_WQS$Value_Rev10 < (Med_50+0.00001)])){
            Med_Quantile <- temp_HADO_DO_WQS$Quantile[temp_HADO_DO_WQS$Value_Rev10 > (Med_50-0.00001) & temp_HADO_DO_WQS$Value_Rev10 < (Med_50+0.00001)]
          } else {
              if (!is_empty(temp_HADO_DO_WQS$Quantile[temp_HADO_DO_WQS$Value_Rev10 > (Med_50-0.00002) & temp_HADO_DO_WQS$Value_Rev10 < (Med_50+0.00002)])){
              Med_Quantile <- temp_HADO_DO_WQS$Quantile[temp_HADO_DO_WQS$Value_Rev10 > (Med_50-0.00002) & temp_HADO_DO_WQS$Value_Rev10 < (Med_50+0.00002)]
              } else {
                if (!is_empty(temp_HADO_DO_WQS$Quantile[temp_HADO_DO_WQS$Value_Rev10 > (Med_50-0.00003) & temp_HADO_DO_WQS$Value_Rev10 < (Med_50+0.00003)])){
                  Med_Quantile <- temp_HADO_DO_WQS$Quantile[temp_HADO_DO_WQS$Value_Rev10 > (Med_50-0.00003) & temp_HADO_DO_WQS$Value_Rev10 < (Med_50+0.00003)]
                } else {
                  if (!is_empty(temp_HADO_DO_WQS$Quantile[temp_HADO_DO_WQS$Value_Rev10 > (Med_50-0.00004) & temp_HADO_DO_WQS$Value_Rev10 < (Med_50+0.00004)])){
                    Med_Quantile <- temp_HADO_DO_WQS$Quantile[temp_HADO_DO_WQS$Value_Rev10 > (Med_50-0.00004) & temp_HADO_DO_WQS$Value_Rev10 < (Med_50+0.00004)]
                  } else {
                    if (!is_empty(temp_HADO_DO_WQS$Quantile[temp_HADO_DO_WQS$Value_Rev10 > (Med_50-0.00005) & temp_HADO_DO_WQS$Value_Rev10 < (Med_50+0.00005)])){
                      Med_Quantile <- temp_HADO_DO_WQS$Quantile[temp_HADO_DO_WQS$Value_Rev10 > (Med_50-0.00005) & temp_HADO_DO_WQS$Value_Rev10 < (Med_50+0.00005)]
                    } else {
                      Med_Quantile <- temp_HADO_DO_WQS$Quantile[temp_HADO_DO_WQS$Value_Rev10 > (Med_50-0.00006) & temp_HADO_DO_WQS$Value_Rev10 < (Med_50+0.00006)]
                    }
                  }}}}
          temp_HADO_DO_WQS$Value_Rev50 <- with(temp_HADO_DO_WQS, 
                                               ifelse(
                                                 Quantile < Med_Quantile & Quantile >= 0.50, Med_50,  Value_Rev10
                                               ))
      } else {
          temp_HADO_DO_WQS$Value_Rev50 <-temp_HADO_DO_WQS$Value_Rev10
          }
        
      
      # Calculate total adjustment made relative to original value and summarize
      temp_HADO_DO_WQS$Value_Diff50 <- temp_HADO_DO_WQS$Value_Rev50-temp_HADO_DO_WQS$Value
      temp_HADO_DO_WQS$AdjFlag10 <- ifelse(!temp_HADO_DO_WQS$Value_Diff10==0,1,0)
      temp_HADO_DO_WQS$AdjFlag50 <- ifelse(!temp_HADO_DO_WQS$Value_Diff50==0,1,0)
      
      summary_HADO_DO_WQS_Adj <- temp_HADO_DO_WQS %>% 
        summarise(n=n(),sum10=sum(AdjFlag10),sum50=sum(AdjFlag50))
      
      summary_HADO_DO_WQS_Adj_cell <- temp_HADO_DO_WQS %>% 
        group_by(IJ) %>%
        summarise(n=n(),sum10=sum(AdjFlag10),sum50=sum(AdjFlag50))
      
      length(summary_HADO_DO_WQS_Adj_cell$IJ) # 579
      length(summary_HADO_DO_WQS_Adj_cell$IJ[summary_HADO_DO_WQS_Adj_cell$sum10>0]) # 303
      length(summary_HADO_DO_WQS_Adj_cell$IJ[summary_HADO_DO_WQS_Adj_cell$sum50>0]) # 303
      
      summary_HADO_DO_WQS <- temp_HADO_DO_WQS %>% 
        group_by(AdjFlag50) %>%
        summarise(n=n(),PctAdj = n()/nrow(temp_HADO_DO_WQS),mean=mean(Value_Diff50),min=min(Value_Diff50),max=max(Value_Diff50))
      
      temp_HADO_DO_WQS$AdjFlag10 <- NULL
      temp_HADO_DO_WQS$AdjFlag50 <- NULL
      temp_HADO_DO_WQS$Value_Diff10 <- NULL
      temp_HADO_DO_WQS$Value_Diff50 <- NULL
      
      # Write QA data sets and adjustment summary
      
      write.csv(
        summary_HADO_DO_WQS,
        paste0(
          WQI_path,"Outputs/QA/",Scen_Name,"DERiver_DOWQS_AdjSummary_",Scen_RunDate,".csv"
        ),
        row.names = F
      )
      
      write.csv(
        summary_HADO_DO_AsIs,
        paste0(
          WQI_path,"Outputs/QA/",Scen_Name,"DERiver_DOWQS_SummaryNoAdj_",Scen_RunDate,".csv"
        ),
        row.names = F
      )
      
      write.csv(
      summary_HADO_DO_WQS_Adj,
      paste0(
        WQI_path,"Outputs/QA/",Scen_Name,"DERiver_DOWQS_AdjSummaryN_",Scen_RunDate,".csv"
      ),
      row.names = F
      )
      
      # Confirm that median of revised time series is >= Med_50
      
      median(temp_HADO_DO_WQS$Value_Rev50) 
      # [1] 2019 actual = 6.307272 2019 design = 6.174386
      
      median(temp_HADO_DO_WQS$Value_Rev50) >= Med_50
      # [1] 2019 actual = TRUE
      # [1] 2019 design = TRUE
      
      # Select columns and rename
      HADO_DO_WQS <- temp_HADO_DO_WQS[,c("IJ","Zone","Date","Value_Rev50","Param")]
      names(HADO_DO_WQS) <- c("IJ","Zone","Date","Value","Param")
      
      # TEST 1: Change all cells by % change of average of what would need to change in the percentile range that would need to change to meet the standard
      
      #@KM: Testing out application of 10th percentile med; THEN 50th percentile med
      
        temp_HADO_DO_WQS_Q10 <- HADO_DO_90[HADO_DO_90$Value_Rev10 == Med_10,]
        
        temp_HADO_DO_WQS_Q10$T1_Ch_10 <- with(temp_HADO_DO_WQS_Q10, (Value_Rev10-Value)/Value)
        
        T1_Mean_Ch_10 <- mean(temp_HADO_DO_WQS_Q10$T1_Ch_10[temp_HADO_DO_WQS_Q10$T1_Ch_10>0])
        
        temp_HADO_DO_WQS$T1_Ch <- with(temp_HADO_DO_WQS, (Value_Rev50-Value)/Value)
        
        if (!is.null(nrow(temp_HADO_DO_WQS$T1_Ch[temp_HADO_DO_WQS$T1_Ch>0]))) 
          {
          T1_Mean_Ch <- mean(temp_HADO_DO_WQS$T1_Ch[temp_HADO_DO_WQS$T1_Ch>0])
          } else {
          T1_Mean_Ch<-0
          T1_Mean_Ch_10<-0
          }
        
        length(unique(temp_HADO_DO_WQS_Q10$IJ)) #411 cells
        summarize_changingdays_bycell <-
          temp_HADO_DO_WQS_Q10 %>%
          group_by(IJ) %>%
          summarise(n = n())
        
        write.csv(
          summarize_changingdays_bycell,
          paste0(
            WQI_path,"Outputs/QA/",Scen_Name,"_changingdays_bycell_",Scen_RunDate,".csv"
          ),
          row.names = F
        )
      
      # TEST 2: Change all cells by average of what would need to change in the percentile range that would need to change to meet the standard
        temp_HADO_DO_WQS$T2_Ch <- with(temp_HADO_DO_WQS, Value_Rev50-Value)
        
        temp_HADO_DO_WQS_Q10$T2_Ch_10 <- with(temp_HADO_DO_WQS_Q10, Value_Rev10-Value)
        
        T2_Mean_Ch_10 <- mean(temp_HADO_DO_WQS_Q10$T2_Ch_10[temp_HADO_DO_WQS_Q10$T2_Ch_10>0])
        
        if (!is.null(nrow(temp_HADO_DO_WQS$T2_Ch[temp_HADO_DO_WQS$T1_Ch>0]))) 
          {
          T2_Mean_Ch <- mean(temp_HADO_DO_WQS$T2_Ch[temp_HADO_DO_WQS$T1_Ch>0])
        } else {
          T2_Mean_Ch<-0
          T2_Mean_Ch_10<-0
        }
      
      # TEST 3: Change all of the cells by % change diff between median value and median value standard
      
      #@KM: Testing out application of 10th percentile med; THEN 50th percentile med
      
        T3_Med_Ch_10 <- (median(temp_HADO_DO_WQS_Q10$Value_Rev10) - median(temp_HADO_DO_WQS_Q10$Value)) / median(temp_HADO_DO_WQS_Q10$Value)
        
        T3_Med_Ch <- (median(HADO_DO_WQS$Value) - median(HADO_DO$Value)) / median(HADO_DO$Value)
      
      # TEST 4: Change all of the cells by diff between median value and median value standard
        T4_Med_Ch_10 <- median(temp_HADO_DO_WQS_Q10$Value_Rev10) - median(temp_HADO_DO_WQS_Q10$Value)
        
        T4_Med_Ch <- median(HADO_DO_WQS$Value) - median(HADO_DO$Value)
      
      # TEST 5: Change all of the cells by % change diff between average changes in HADO scenario compared to baseline
        temp_HADO_DO_WQS <- merge(temp_HADO_DO_WQS,bind_bsln_mod_FMA[bind_bsln_mod_FMA$Param == "DO",],by=c("IJ","Zone","Date","Param"))
        
        names(temp_HADO_DO_WQS) <- c("IJ","Zone","Date","Param","Value","Quantile","Value_Rev10","Value_Rev50","T1_Ch","T2_Ch","Value_Bsln")
        
        temp_HADO_DO_WQS$T5_Ch <- with(temp_HADO_DO_WQS, (Value-Value_Bsln)/Value_Bsln)
        
        T5_Mean_Ch <- mean(temp_HADO_DO_WQS$T5_Ch[temp_HADO_DO_WQS$T5_Ch>0])
      
      # TEST 6: Change all of the cells by diff between average changes in HADO scenario compared to baseline
      
        temp_HADO_DO_WQS$T6_Ch <- with(temp_HADO_DO_WQS, Value-Value_Bsln)
        
        T6_Mean_Ch <- mean(temp_HADO_DO_WQS$T6_Ch[temp_HADO_DO_WQS$T6_Ch>0])
      

        if (nrow(temp_HADO_DO_WQS_Q10)==0) 
        {
          T1_Mean_Ch_10<-0
          T2_Mean_Ch_10<-0
          T3_Med_Ch_10<-0
          T4_Med_Ch_10<-0
        }
        
      # Method adds both the avg change associated with the 10% portion of the standard and the 50% portion.
      
      temp_HADO_DO_WQS$T1 <-temp_HADO_DO_WQS$Value+T1_Mean_Ch*temp_HADO_DO_WQS$Value + T1_Mean_Ch_10*temp_HADO_DO_WQS$Value
                                  
      temp_HADO_DO_WQS$T2 <-temp_HADO_DO_WQS$Value+T2_Mean_Ch+T2_Mean_Ch_10
      
      temp_HADO_DO_WQS$T3 <-temp_HADO_DO_WQS$Value+T3_Med_Ch*temp_HADO_DO_WQS$Value+T3_Med_Ch_10*temp_HADO_DO_WQS$Value
      
      temp_HADO_DO_WQS$T4 <-temp_HADO_DO_WQS$Value+T4_Med_Ch+T4_Med_Ch_10
      
      temp_HADO_DO_WQS$T5 <-temp_HADO_DO_WQS$Value+T5_Mean_Ch*temp_HADO_DO_WQS$Value
      
      temp_HADO_DO_WQS$T6 <-temp_HADO_DO_WQS$Value+T6_Mean_Ch
      
      # --- QA Start: Check to see if meeting the standard --- #
      ###
      median(temp_HADO_DO_WQS$T1)
      median(temp_HADO_DO_WQS$T1) >= Med_50
      # [1] TRUE 
      
      temp_HADO_DO_WQS_T1 <- quantile(temp_HADO_DO_WQS$T1, probs = c(0.1))
      
      min(temp_HADO_DO_WQS$T1[temp_HADO_DO_WQS$Quantile >= 0.1]) # If < Med_10, then not meeting standard
      min(temp_HADO_DO_WQS$T1[temp_HADO_DO_WQS$Quantile >= 0.1]) >= Med_10
      # [1] FALSE
      
      ###
      median(temp_HADO_DO_WQS$T2)
      median(temp_HADO_DO_WQS$T2) >= Med_50
      # [1] TRUE
      
      temp_HADO_DO_WQS_T2 <- quantile(temp_HADO_DO_WQS$T2, probs = c(0.1))
      
      min(temp_HADO_DO_WQS$T2[temp_HADO_DO_WQS$Quantile >= 0.1]) # If < Med_10, then not meeting standard
      min(temp_HADO_DO_WQS$T2[temp_HADO_DO_WQS$Quantile >= 0.1]) >= Med_10
      # [1] FALSE
      
      ###
      median(temp_HADO_DO_WQS$T3)
      median(temp_HADO_DO_WQS$T3) >= Med_50
      # [1] TRUE
      
      temp_HADO_DO_WQS_T3 <- quantile(temp_HADO_DO_WQS$T3, probs = c(0.1))
      
      median(temp_HADO_DO_WQS$T3[temp_HADO_DO_WQS$T3 < temp_HADO_DO_WQS_T3])
      min(temp_HADO_DO_WQS$T3[temp_HADO_DO_WQS$Quantile >= 0.1]) # If < Med_10, then not meeting standard
      min(temp_HADO_DO_WQS$T3[temp_HADO_DO_WQS$Quantile >= 0.1]) >= Med_10
      # [1] 2019 actual = FALSE  2019 design = FALSE
      
      ###
      median(temp_HADO_DO_WQS$T4)
      median(temp_HADO_DO_WQS$T4) >= Med_50
      # [1] TRUE
      
      temp_HADO_DO_WQS_T4 <- quantile(temp_HADO_DO_WQS$T4, probs = c(0.1))
      
      median(temp_HADO_DO_WQS$T4[temp_HADO_DO_WQS$T4 < temp_HADO_DO_WQS_T4])
      min(temp_HADO_DO_WQS$T4[temp_HADO_DO_WQS$Quantile >= 0.1]) # If < Med_10, then not meeting standard
      min(temp_HADO_DO_WQS$T4[temp_HADO_DO_WQS$Quantile >= 0.1]) >= Med_10
      # [1] 2019 actual = FALSE  2019 design = FALSE
      
      # --- QA End --- #
      
      HADO_DO_WQS_Tests <- list()
      EPA_WQS_FMA_Tests <- list()

      for(i in 1:6) {

      HADO_DO_WQS_Tests[[i]] <- temp_HADO_DO_WQS[,c("IJ","Zone","Date",paste0("T",i),"Param")]
      names(HADO_DO_WQS_Tests[[i]]) <- c("IJ","Zone","Date","Value","Param")

      EPA_WQS_FMA_Tests[[i]] <- rbind(bind_bsln_mod_FMA[!(bind_bsln_mod_FMA$Param %in% "DO"),],HADO_DO_WQS_Tests[[i]],bind_bsln_obs_FMA)

      }

# Bind the modeled and observed scenario data together considering six scenarios
      # ** HADO Unadjusted Scenario -------------
      HADO_AsIs_FMA <- rbind(bind_bsln_mod_FMA[!(bind_bsln_mod_FMA$Param %in% "DO"),],HADO_DO_AsIs,bind_bsln_obs_FMA)
      
      # ** EPA Standard Scenario -------------
      # EPA Standard for DO; calibration model TN and TP; observed BOD, TSS, FC
      
      EPA_WQS_FMA <- rbind(bind_bsln_mod_FMA[!(bind_bsln_mod_FMA$Param %in% "DO"),],HADO_DO_WQS,bind_bsln_obs_FMA)
      
      # ** Sensitivity Analysis Scenario 1 -------------
      
      # EPA Standard for DO; calibration model TN, TP, and observed BOD, TSS, FC reduced by percentage change in DO
      
      # Calculate average change in DO from WASP calibration to HADO scenario with applied EPA standard; 
      # apply to HADO scenario TN and TP and observed BOD, TSS, and FC data

# Average over time and space

      EPA_Standard <- HADO_DO_WQS
      
      names(EPA_Standard) <- c("IJ","Zone","Date","EPA_Standard_DO","Param") 
      
      temp_avg_ch_DO <- merge(EPA_Standard,bind_bsln_FMA[bind_bsln_FMA$Param=="DO",],all=TRUE)
      
      names(temp_avg_ch_DO) <- c("IJ","Zone","Date","Param","EPA_Standard_DO","Bsln_DO")
          
      avg_ch_DO <-
        temp_avg_ch_DO[, c("EPA_Standard_DO","Bsln_DO")] %>%
        mutate(perc_ch_DO = (EPA_Standard_DO-Bsln_DO)/Bsln_DO) %>%
        summarise(perc_ch_DO = mean(perc_ch_DO, na.rm=T))
      
      avg_ch_DO_byRegion <-
        temp_avg_ch_DO %>%
        group_by(Zone, Param) %>%
        summarise(n = n(),ch_avg = mean((EPA_Standard_DO/Bsln_DO)-1, na.rm=T))
      
      write.csv(
        avg_ch_DO_byRegion,
        paste0(
          WQI_path,"Outputs/QA/",Scen_Name,"DERiver_DOWQS_AvgDOChgByRegion_",Scen_RunDate,".csv"
        ),
        row.names = F
      )
      
      bind_bsln_obs_FMA_sens1 <- EPA_WQS_FMA[!EPA_WQS_FMA$Param == "DO",]
      bind_bsln_obs_FMA_sens1 <- merge(bind_bsln_obs_FMA_sens1,avg_ch_DO_byRegion[,c("Zone","ch_avg")],by = "Zone")
      bind_bsln_obs_FMA_sens1$Value <- bind_bsln_obs_FMA_sens1$Value*(1-(bind_bsln_obs_FMA_sens1$ch_avg))
      bind_bsln_obs_FMA_sens1$ch_avg <- NULL
      
      bind_sens1_FMA <- rbind(bind_bsln_obs_FMA_sens1,HADO_DO_WQS)
      
      # ** Sensitivity Analysis Scenario 2 -------------
      
      # EPA Standard for DO; calibration model TN, TP, and observed BOD reduced by percentage change in DO; TSS and FC remain constant
      
      # Calculate average change in DO from baseline to HADO scenario with applied EPA standard; 
      # apply to HADO scenario TN and TP and observed BOD data
      
      # Average over time and space
      
      bind_bsln_obs_FMA_sens2 <- EPA_WQS_FMA[!(EPA_WQS_FMA$Param == "DO"|EPA_WQS_FMA$Param == "TSS"|EPA_WQS_FMA$Param == "FC"),]
      bind_bsln_obs_FMA_sens2 <- merge(bind_bsln_obs_FMA_sens2,avg_ch_DO_byRegion[,c("Zone","ch_avg")],by = "Zone")
      bind_bsln_obs_FMA_sens2$Value <- bind_bsln_obs_FMA_sens2$Value*(1-(bind_bsln_obs_FMA_sens2$ch_avg))
      bind_bsln_obs_FMA_sens2$ch_avg <- NULL
      
      bind_sens2_FMA <- rbind(bind_bsln_obs_FMA_sens2,EPA_WQS_FMA[(EPA_WQS_FMA$Param == "TSS"|EPA_WQS_FMA$Param == "FC"),],HADO_DO_WQS)
      
      # ** Sensitivity Analysis Scenario 3 -------------
      
      # EPA Standard for DO; observed BOD reduced by percentage change in BOD at 12 effluent sites; TN, TP, TSS, FC remain constant
      
      bind_bsln_obs_FMA_sens3 <- EPA_WQS_FMA[EPA_WQS_FMA$Param == "BOD",]
      bind_bsln_obs_FMA_sens3$Value <- bind_bsln_obs_FMA_sens3$Value*(1+BOD_PercCh)
      
      bind_sens3_FMA <- rbind(bind_bsln_obs_FMA_sens3,EPA_WQS_FMA[!(EPA_WQS_FMA$Param == "BOD"),])
      
      # ** Sensitivity Analysis Scenario 4 -------------
      
      # EPA Standard for DO; observed BOD reduced by percentage change in DO; TN, TP, TSS and FC remain constant
      
      # Calculate average change in DO from baseline to HADO scenario with applied EPA standard; 
      # apply to observed BOD data
      
      # Average over time and space
      
      bind_bsln_obs_FMA_sens4 <- EPA_WQS_FMA[!(EPA_WQS_FMA$Param == "DO"|EPA_WQS_FMA$Param == "TN"|EPA_WQS_FMA$Param == "TP"|EPA_WQS_FMA$Param == "TSS"|EPA_WQS_FMA$Param == "FC"),]
      bind_bsln_obs_FMA_sens4 <- merge(bind_bsln_obs_FMA_sens4,avg_ch_DO_byRegion[,c("Zone","ch_avg")],by = "Zone")
      bind_bsln_obs_FMA_sens4$Value <- bind_bsln_obs_FMA_sens4$Value*(1-(bind_bsln_obs_FMA_sens4$ch_avg))
      bind_bsln_obs_FMA_sens4$ch_avg <- NULL
      
      bind_sens4_FMA <- rbind(bind_bsln_obs_FMA_sens4,EPA_WQS_FMA[(EPA_WQS_FMA$Param == "TN"|EPA_WQS_FMA$Param == "TP"|EPA_WQS_FMA$Param == "TSS"|EPA_WQS_FMA$Param == "FC"),],HADO_DO_WQS)
      
      # ** Sensitivity Analysis Scenario 5 -------------
      
      # Unadjusted HADO model for DO; observed BOD reduced by percentage change in BOD at 12 effluent sites; TN, TP, TSS and FC remain constant
      
      # Calculate average change in DO from baseline to HADO scenario "as is"; 
      # apply to observed BOD data
      
      # Average over time and space
      
      bind_bsln_obs_FMA_sens5 <- HADO_AsIs_FMA[!(HADO_AsIs_FMA$Param == "DO"|HADO_AsIs_FMA$Param == "TN"|HADO_AsIs_FMA$Param == "TP"|HADO_AsIs_FMA$Param == "TSS"|HADO_AsIs_FMA$Param == "FC"),]
      bind_bsln_obs_FMA_sens5$Value <- bind_bsln_obs_FMA_sens5$Value*(1+BOD_PercCh)
      
      bind_sens5_FMA <- rbind(bind_bsln_obs_FMA_sens5,HADO_AsIs_FMA[(HADO_AsIs_FMA$Param == "TN"|HADO_AsIs_FMA$Param == "TP"|HADO_AsIs_FMA$Param == "TSS"|HADO_AsIs_FMA$Param == "FC"),],HADO_DO_AsIs)
      
      # ** Sensitivity Analysis Scenario 6 -------------
      
      # Unadjusted HADO model for DO; observed BOD reduced by percentage change in DO; TN, TP, TSS and FC remain constant
      
      # Calculate average change in DO from baseline to HADO scenario "as is"; 
      # apply to observed BOD data
      
      # Average over time and space
      
      HADO_DO <- HADO_DO_AsIs
      names(HADO_DO) <- c("IJ","Zone","Date","HADO_DO","Param") 
      
      temp_avg_ch_DO <- merge(HADO_DO,bind_bsln_FMA[bind_bsln_FMA$Param=="DO",],all=TRUE)
      names(temp_avg_ch_DO) <- c("IJ","Zone","Date","Param","HADO_DO","Bsln_DO")
      
      avg_HADO_ch_DO <-
        temp_avg_ch_DO[, c("HADO_DO","Bsln_DO")] %>%
        mutate(perc_ch_DO = (HADO_DO-Bsln_DO)/Bsln_DO) %>%
        summarise(perc_ch_DO = mean(perc_ch_DO, na.rm=T))
      
      avg_ch_HADO_DO_byRegion <-
        temp_avg_ch_DO %>%
        group_by(Zone, Param) %>%
        summarise(n = n(),ch_avg = mean((HADO_DO/Bsln_DO)-1, na.rm=T))
      
      write.csv(
        avg_ch_HADO_DO_byRegion,
        paste0(
          WQI_path,"Outputs/QA/",Scen_Name,"DERiver_DOWQS_AvgHADODOChgByRegion_",Scen_RunDate,".csv"
        ),
        row.names = F
      )
      
      bind_bsln_obs_FMA_sens6 <- HADO_AsIs_FMA[!(HADO_AsIs_FMA$Param == "DO"|HADO_AsIs_FMA$Param == "TN"|HADO_AsIs_FMA$Param == "TP"|HADO_AsIs_FMA$Param == "TSS"|HADO_AsIs_FMA$Param == "FC"),]
      bind_bsln_obs_FMA_sens6 <- merge(bind_bsln_obs_FMA_sens6,avg_ch_HADO_DO_byRegion[,c("Zone","ch_avg")],by = "Zone")
      bind_bsln_obs_FMA_sens6$Value <- bind_bsln_obs_FMA_sens6$Value*(1-(bind_bsln_obs_FMA_sens6$ch_avg))
      bind_bsln_obs_FMA_sens6$ch_avg <- NULL
      
      bind_sens6_FMA <- rbind(bind_bsln_obs_FMA_sens6,HADO_AsIs_FMA[(HADO_AsIs_FMA$Param == "TN"|HADO_AsIs_FMA$Param == "TP"|HADO_AsIs_FMA$Param == "TSS"|HADO_AsIs_FMA$Param == "FC"),],HADO_DO_AsIs)

# Summarize results by region and parameter 

      avg_bsln_byRegion <-
        bind_bsln_FMA %>%
          group_by(Zone, Param) %>%
          summarise(n = n(),Bsln_avg = mean(Value, na.rm=T),Bsln_min = min(Value, na.rm=T),Bsln_max = max(Value, na.rm=T))
      
      avg_EPAWQ_byRegion <-
        EPA_WQS_FMA %>%
        group_by(Zone, Param) %>%
        summarise(n = n(),EPAWQ_avg = mean(Value, na.rm=T))

# SUBINDEX and WQI CALCULATIONS ------------------------------------------

      # * FMA Ecoregion Overlay -------------
      
      # Connect FMA with Level III Omernik ecoregions
      
      #*# Check to make sure datasets are in appropriate projections to calculate areas
      
      proj <- st_crs(Ecoregions)
      
      FMA_boundary_Albers <- st_transform(FMA_boundary, crs=proj$wkt)
      
      # Overlay ecoregion with FMA boundary
      FMA_eco <- st_intersection(FMA_boundary_Albers, Ecoregions) 
      # Note to user: You may receive the warning message 'attribute variables are assumed to be spatially constant throughout all geometries'.
      # This warning is only of concern if values are aggregates, such as population count. Because attribute variables in this analysis are
      # constant, you may ignore this warning.
      
      # Calculate area of intersection
      FMA_eco$Int_Area_m2 <- st_area(FMA_eco)
      
      FMA_eco <- FMA_eco %>% 
        st_drop_geometry
      
      FMA_eco_water <- FMA_eco[FMA_eco$TYPE=="water",]
      
      # FMA covers two ecoregions:
      unique(FMA_eco_water$US_L3NAME)
      # [1] "Middle Atlantic Coastal Plain" "Atlantic Coastal Pine Barrens"
      
      FMA_eco_water[,c("US_L3NAME","Int_Area_m2")]
      # US_L3NAME       Int_Area_m2
      # 1   Middle Atlantic Coastal Plain 90969173.50 [m^2]
      # 1.1 Atlantic Coastal Pine Barrens    13656.19 [m^2]
      
      # Note: because the Middle Atlantic Coastal Plain ecoregion makes up the majority of the FMA (>99%),
      # we rely on subindex (SI) curves from this region only

# * Subindex Curves -------------

      si.tn <- read.csv(paste0(WQI_path,"Inputs/TNSI.csv"), header = TRUE)
      si.tn <- si.tn %>% rename_with( ~ paste0("TN_", .x), !c("EcoID", "EcoName"))
      
      si.tp <- read.csv(paste0(WQI_path,"Inputs/TPSI.csv"), header = TRUE)
      si.tp <- si.tp %>% rename_with( ~ paste0("TP_", .x), !c("EcoID", "EcoName"))
      
      si.ssc <- read.csv(paste0(WQI_path,"Inputs/SSCSI.csv"), header = TRUE)
      si.ssc <- si.ssc %>% rename_with( ~ paste0("SSC_", .x), !c("EcoID", "EcoName"))
      
      temp_si <- Reduce(function(x,y) merge(x=x, y=y, by=c("EcoID","EcoName"), all.x=T),
                        list(si.tn, si.tp, si.ssc))
      
      # Subset to ecoregions of relevance to the FMA #*#
      
      si <- temp_si[temp_si$EcoName %in% "Middle Atlantic Coastal Plain",]
      # Note that this ecoregion corresponds to Coastal Plains when using the NARS aggregated ecoregions

# * SI and WQI Calculations -------------

      # Loop through baseline and scenario dataframes to calculate SI and WQI values
      
      #*# Note to user: If adding/changing scenarios, must adjust SI_List and Scenarios
      
      SI_List <- list(bind_bsln_FMA,HADO_AsIs_FMA,EPA_WQS_FMA,bind_sens1_FMA,bind_sens2_FMA,bind_sens3_FMA,bind_sens4_FMA,bind_sens5_FMA,bind_sens6_FMA,
                      EPA_WQS_FMA_Tests[[1]],EPA_WQS_FMA_Tests[[2]],EPA_WQS_FMA_Tests[[3]],EPA_WQS_FMA_Tests[[4]],EPA_WQS_FMA_Tests[[5]],EPA_WQS_FMA_Tests[[6]]) 
      
      Scenarios <- c("Baseline","HADO_AsIs","EPA_Standard","Sensitivity_Scenario_1","Sensitivity_Scenario_2","Sensitivity_Scenario_3","Sensitivity_Scenario_4",
                     "Sensitivity_Scenario_5","Sensitivity_Scenario_6","Test_1","Test_2","Test_3","Test_4","Test_5","Test_6")  
      
      Scenario_Descr <- c("Baseline",
                          "HADO Model as is",
                          "EPA standard applied to HADO Model",
                          "EPA Standard for DO; calibration model TN, TP, and observed BOD, TSS, FC reduced by percentage change in DO by zone",
                          "EPA Standard for DO; calibration model TN, TP, and observed BOD reduced by percentage change in DO  by zone; TSS and FC remain constant",
                          "EPA Standard for DO; observed BOD reduced by percentage change in BOD at 12 effluent sites; TN, TP, TSS, FC remain constant",
                          "EPA Standard for DO; observed BOD reduced by percentage change in DO by zone; TN, TP, TSS, FC remain constant",
                          "HADO Model as is for DO; observed BOD reduced by percentage change in BOD at 12 effluent sites; TN, TP, TSS, FC remain constant",
                          "HADO Model as is for DO; observed BOD reduced by percentage change in DO by zone; TN, TP, TSS, FC remain constant",
                          "DO - Change all cells by % change of average of what would need to change in the percentile range that would need to change to meet the standard",
                          "DO - Change all cells by average of what would need to change in the percentile range that would need to change to meet the standard",
                          "DO - Change all of the cells by % change diff between median value and median value in the percentile range that would need to change to meet the standard",
                          "DO - Change all of the cells by diff between median value and median value in the percentile range that would need to change to meet the standard",
                          "DO - Change all of the cells by % change diff between average changes in HADO scenario compared to baseline",
                          "DO - Change all of the cells by diff between average changes in HADO scenario compared to baseline")
      
      
      for (i in seq_along(SI_List)) {
        Scen <- SI_List[[i]]
        
        # Cast Param as columns:
        
        cast_FMA <- dcast(Scen, IJ + Zone + Date ~ Param, value.var = "Value")
        
        SWQI <- merge(cast_FMA, si)
        
        # Dissolved oxygen sub-index value
        
        SWQI$si_do <- with(SWQI,
                           ifelse(
                             !is.na(DO) & DO <= 3.3,
                             10,
                             ifelse(
                               !is.na(DO) & DO > 3.3 & DO < 10.5,
                               (-80.29 + 31.88 * DO - 1.401 * (DO) ^
                                  2),
                               ifelse(!is.na(DO) &
                                        DO >= 10.5, 100,
                                      NA)
                             )
                           ))
        
        # Biological dissolved oxygen sub-index value
        
        SWQI$si_bod <- with(SWQI,
                            ifelse(!is.na(BOD) & BOD > 8, 10,
                                   ifelse(!is.na(BOD) & BOD <= 8,
                                          (
                                            100 * exp(BOD * -0.1993)
                                          ),
                                          NA)))
        
        # Total nitrogen sub-index value
        
        SWQI$si_tn <- with(SWQI,
                           ifelse(
                             !is.na(TN) & TN > TN_val10,
                             10,
                             ifelse(
                               !is.na(TN) & TN > TN_val100 & TN <= TN_val10,
                               (TN_coeff_a * exp(TN * TN_coeff_b)),
                               ifelse(!is.na(TN) &
                                        TN <= TN_val100, 100,
                                      NA)
                             )
                           ))
        
        # Total phosphorus sub-index value
        
        SWQI$si_tp <- with(SWQI,
                           ifelse(
                             !is.na(TP) & TP > TP_val10,
                             10,
                             ifelse(
                               !is.na(TP) & TP > TP_val100 & TP <= TP_val10,
                               (TP_coeff_a * exp(TP * TP_coeff_b)),
                               ifelse(!is.na(TP) &
                                        TP <= TP_val100, 100,
                                      NA)
                             )
                           ))
        
        SWQI$si_fc <- with(SWQI,
                           ifelse(
                             !is.na(FC) & FC > 1600,
                             10,
                             ifelse(
                               !is.na(FC) & FC <= 1600 & FC > 50,
                               (98 * exp((FC - 50) * -9.9178E-4)),
                               ifelse(!is.na(FC) &
                                        FC <= 50, 98,
                                      NA)
                             )
                           ))
        
        # TSS to SSC conversion is based on the Northeast SPARROW model output: https://www.usgs.gov/centers/maryland-delaware-d.c.-water-science-center/science/2012-sparrow-models-northeast-total
      
        SWQI$SSC <- with(SWQI, TSS * (1 / exp(-1.46)))
        SWQI$si_ssc <- with(SWQI,
                            ifelse(
                              !is.na(SSC) & SSC > SSC_val10,
                              10,
                              ifelse(
                                !is.na(SSC) & SSC > SSC_val100 & SSC <= SSC_val10,
                                (SSC_coeff_a * exp(SSC * SSC_coeff_b)),
                                ifelse(!is.na(SSC) &
                                         SSC <= SSC_val100, 100,
                                       NA)
                              )
                            ))
        
        # Total sub-index values, based on NARS curves for Coastal Plains ecoregions with median at 70 and 10th percentile at 100
        # These are the curves used in the final results instead of the previous SPARROW-based curves.
        
        NARS_TN_val10<-3.540210143
        NARS_TN_val100<-0.3214
        NARS_TN_coeff_a<-125.8492296
        NARS_TN_coeff_b<--0.715352876
          
        NARS_TP_val10<-0.456986278
        NARS_TP_val100<-0.0280582
        NARS_TP_coeff_a<-116.255818
        NARS_TP_coeff_b<--5.368231203
        
        NARS_TSS_val10<-59.90126612 #Note that curves are developed for TSS
        NARS_TSS_val100<-1.8
        NARS_TSS_coeff_a<-107.3940924
        NARS_TSS_coeff_b<--0.039630549

        SWQI$si_NARS_tn <- with(SWQI,
                           ifelse(
                             !is.na(TN) & TN > NARS_TN_val10,
                             10,
                             ifelse(
                               !is.na(TN) & TN > NARS_TN_val100 & TN <= NARS_TN_val10,
                               (NARS_TN_coeff_a * exp(TN * NARS_TN_coeff_b)),
                               ifelse(!is.na(TN) &
                                        TN <= NARS_TN_val100, 100,
                                      NA)
                             )
                           ))
        
        SWQI$si_NARS_tp <- with(SWQI,
                           ifelse(
                             !is.na(TP) & TP > NARS_TP_val10,
                             10,
                             ifelse(
                               !is.na(TP) & TP > NARS_TP_val100 & TP <= NARS_TP_val10,
                               (NARS_TP_coeff_a * exp(TP * NARS_TP_coeff_b)),
                               ifelse(!is.na(TP) &
                                        TP <= NARS_TP_val100, 100,
                                      NA)
                             )
                           ))
        
        SWQI$si_NARS_tss <- with(SWQI,
                                ifelse(
                                  !is.na(TSS) & TSS > NARS_TSS_val10,
                                  10,
                                  ifelse(
                                    !is.na(TSS) & TSS > NARS_TSS_val100 & TSS <= NARS_TSS_val10,
                                    (NARS_TSS_coeff_a * exp(TSS * NARS_TSS_coeff_b)),
                                    ifelse(!is.na(TSS) &
                                             TSS <= NARS_TSS_val100, 100,
                                           NA)
                                  )
                                ))
        
        # Bounds correction
        
        SWQI$si_do[which(!is.na(SWQI$si_do) & SWQI$si_do > 100)] <- 100
        SWQI$si_do[which(!is.na(SWQI$si_do) & SWQI$si_do < 10)] <- 10
        
        SWQI$si_bod[which(!is.na(SWQI$si_bod) & SWQI$si_bod > 100)] <- 100
        SWQI$si_bod[which(!is.na(SWQI$si_bod) & SWQI$si_bod < 10)] <- 10
        
        SWQI$si_tn[which(!is.na(SWQI$si_tn) & SWQI$si_tn > 100)] <- 100
        SWQI$si_tn[which(!is.na(SWQI$si_tn) & SWQI$si_tn < 10)] <- 10
        
        SWQI$si_tp[which(!is.na(SWQI$si_tp) & SWQI$si_tp > 100)] <- 100
        SWQI$si_tp[which(!is.na(SWQI$si_tp) & SWQI$si_tp < 10)] <- 10
        
        SWQI$si_ssc[which(!is.na(SWQI$si_ssc) & SWQI$si_ssc > 100)] <- 100
        SWQI$si_ssc[which(!is.na(SWQI$si_ssc) & SWQI$si_ssc < 10)] <- 10
        
        SWQI$si_fc[which(!is.na(SWQI$si_fc) & SWQI$si_fc > 98)] <- 98
        SWQI$si_fc[which(!is.na(SWQI$si_fc) & SWQI$si_fc < 10)] <- 10
        
        SWQI$si_NARS_tn[which(!is.na(SWQI$si_NARS_tn) & SWQI$si_NARS_tn > 100)] <- 100
        SWQI$si_NARS_tn[which(!is.na(SWQI$si_NARS_tn) & SWQI$si_NARS_tn < 10)] <- 10
        
        SWQI$si_NARS_tp[which(!is.na(SWQI$si_NARS_tp) & SWQI$si_NARS_tp > 100)] <- 100
        SWQI$si_NARS_tp[which(!is.na(SWQI$si_NARS_tp) & SWQI$si_NARS_tp < 10)] <- 10
        
        SWQI$si_NARS_tss[which(!is.na(SWQI$si_NARS_tss) & SWQI$si_NARS_tss > 100)] <- 100
        SWQI$si_NARS_tss[which(!is.na(SWQI$si_NARS_tss) & SWQI$si_NARS_tss < 10)] <- 10
        
          # --- QA Start --- #
        
        # Compare calculated DO SI values before averaging
        
        SWQI_do_comp <- SWQI[, c("IJ", "Zone", "Date", "si_do")]
        names(SWQI_do_comp) <- c("IJ", "Zone", "Date", paste0("si_do_",Scenarios[i]))
        
        # Bind together
        
        if (exists("bind_SWQI_do_comp")) {
          bind_SWQI_do_comp <- merge(bind_SWQI_do_comp, SWQI_do_comp)
        } else {
          bind_SWQI_do_comp <- SWQI_do_comp
        }
        
        # Compare calculated TN and TP SI values before averaging
        
        SWQI_tn_tp_comp <- SWQI[, c("IJ", "Zone", "Date", "si_tn", "si_tp")]
        names(SWQI_tn_tp_comp) <- c("IJ", "Zone", "Date", paste0("si_tn_",Scenarios[i]), paste0("si_tp_",Scenarios[i]))
        
        # Bind together
        
        if (exists("bind_SWQI_tn_tp_comp")) {
          bind_SWQI_tn_tp_comp <- merge(bind_SWQI_tn_tp_comp, SWQI_tn_tp_comp)
        } else {
          bind_SWQI_tn_tp_comp <- SWQI_tn_tp_comp
        }
        
        # --- QA End --- #
        
        SWQI_summer_FMA <-
          SWQI[, c("si_tn", "si_tp", "si_ssc", "si_do", "si_bod", "si_fc","si_NARS_tn", "si_NARS_tp", "si_NARS_tss")] %>%
          summarise_all(mean, na.rm = TRUE)
        
        DO_summer_FMA <-
          SWQI[, c("DO","si_do")] %>%
          summarise_all(mean, na.rm = TRUE)
        
        SWQI_summer_FMA$WQI <-
          with(SWQI_summer_FMA,
               exp(
                 .24 * log(si_do) + .22 * log(si_fc) + .15 * log(si_bod) +
                   .14 *
                   log(si_tn) + .14 * log(si_tp) + .11 * log(si_ssc)
               ))
        
        SWQI_summer_FMA$NARS_WQI <-
          with(SWQI_summer_FMA,
               exp(
                 .24 * log(si_do) + .22 * log(si_fc) + .15 * log(si_bod) +
                   .14 *
                   log(si_NARS_tn) + .14 * log(si_NARS_tp) + .11 * log(si_NARS_tss)
               ))
        
        SWQI_summer_FMA$Scenario <- Scenarios[i]
        DO_summer_FMA$Scenario <- Scenarios[i]
        
        SWQI_summer_FMA$Scenario_Description <- Scenario_Descr[i]
        DO_summer_FMA$Scenario_Description <- Scenario_Descr[i]
        
        # Bind together
        
        if (exists("bind_SWQI_summer_FMA")) {
          bind_SWQI_summer_FMA <- rbind(bind_SWQI_summer_FMA, SWQI_summer_FMA)
        } else {
          bind_SWQI_summer_FMA <- SWQI_summer_FMA
        }
        
        if (exists("bind_DO_summer_FMA")) {
          bind_DO_summer_FMA <- rbind(bind_DO_summer_FMA, DO_summer_FMA)
        } else {
          bind_DO_summer_FMA <- DO_summer_FMA
        }
        
      }

# Calculate WQI changes

      bind_SWQI_summer_FMA$WQI_ch <-
        bind_SWQI_summer_FMA$WQI - bind_SWQI_summer_FMA$WQI[bind_SWQI_summer_FMA$Scenario == "Baseline"]
      
      bind_SWQI_summer_FMA$WQI_perc_ch <-
        (bind_SWQI_summer_FMA$WQI - bind_SWQI_summer_FMA$WQI[bind_SWQI_summer_FMA$Scenario == "Baseline"]) /
        bind_SWQI_summer_FMA$WQI[bind_SWQI_summer_FMA$Scenario == "Baseline"]
      
      bind_SWQI_summer_FMA$NARS_WQI_ch <-
        bind_SWQI_summer_FMA$NARS_WQI - bind_SWQI_summer_FMA$NARS_WQI[bind_SWQI_summer_FMA$Scenario == "Baseline"]
      
      bind_SWQI_summer_FMA$NARS_WQI_perc_ch <-
        (bind_SWQI_summer_FMA$NARS_WQI - bind_SWQI_summer_FMA$NARS_WQI[bind_SWQI_summer_FMA$Scenario == "Baseline"]) /
        bind_SWQI_summer_FMA$NARS_WQI[bind_SWQI_summer_FMA$Scenario == "Baseline"]
      
      bind_SWQI_summer_FMA <- bind_SWQI_summer_FMA %>% relocate(Scenario,Scenario_Description)

# Calculate DO changes

      bind_DO_summer_FMA$DO_ch <-
        bind_DO_summer_FMA$DO - bind_DO_summer_FMA$DO[bind_DO_summer_FMA$Scenario == "Baseline"]
      
      bind_DO_summer_FMA$DO_perc_ch <-
        (bind_DO_summer_FMA$DO - bind_DO_summer_FMA$DO[bind_DO_summer_FMA$Scenario == "Baseline"]) /
        bind_DO_summer_FMA$DO[bind_DO_summer_FMA$Scenario == "Baseline"]
      
      bind_DO_summer_FMA$si_do_ch <-
        bind_DO_summer_FMA$si_do - bind_DO_summer_FMA$si_do[bind_DO_summer_FMA$Scenario == "Baseline"]
      
      bind_DO_summer_FMA$si_do_perc_ch <-
        (bind_DO_summer_FMA$si_do - bind_DO_summer_FMA$si_do[bind_DO_summer_FMA$Scenario == "Baseline"]) /
        bind_DO_summer_FMA$si_do[bind_DO_summer_FMA$Scenario == "Baseline"]
      
      bind_DO_summer_FMA <- bind_DO_summer_FMA %>% relocate(Scenario,Scenario_Description)

# Write to files

      # QA files

      write.csv(
        bind_SWQI_do_comp,
        paste0(
          WQI_path,
          "Outputs/QA/",Scen_Name,"SI_DO_ScenarioComp_",Scen_RunDate,".csv"
        ),
        row.names = F
      )
      
      write.csv(
        bind_SWQI_tn_tp_comp,
        paste0(
          WQI_path,
          "Outputs/QA/",Scen_Name,"SI_TN_TP_ScenarioComp_",Scen_RunDate,".csv"
        ),
        row.names = F
      )
      
      dataset_names <-
        list(
          'Total Nitrogen' = Cal_HADO_TN_TP_DO[Cal_HADO_TN_TP_DO$Param == "TN",],
          'Total Phosphorus' = Cal_HADO_TN_TP_DO[Cal_HADO_TN_TP_DO$Param == "TP",],
          'Dissolved Oxygen' = Cal_HADO_TN_TP_DO[Cal_HADO_TN_TP_DO$Param == "DO",]
        )
      
      write.xlsx(
        dataset_names,
        file = paste0(WQI_path, "Outputs/QA/",Scen_Name,"_Comparison_",Scen_RunDate,".xlsx")
      )
      
      # Final output file
      
      write.csv(
        bind_SWQI_summer_FMA,
        paste0(
          WQI_path,"Outputs/",Scen_Name,"DERiver_DOWQS_WQI_PercChange_",Scen_RunDate,".csv"
        ),
        row.names = F
      )
      
      write.csv(
        bind_DO_summer_FMA,
        paste0(
          WQI_path,"Outputs/",Scen_Name,"DERiver_DOWQS_DO_PercChange_",Scen_RunDate,".csv"
        ),
        row.names = F
      )

      
      write.csv(
        avg_bsln_byRegion,
        paste0(
          WQI_path,"Outputs/QA/",Scen_Name,"DERiver_avgBlsn_",Scen_RunDate,".csv"
        ),
        row.names = F
      )
      
      
      write.csv(
        avg_EPAWQ_byRegion,
        paste0(
          WQI_path,"Outputs/QA/",Scen_Name,"DERiver_avgEPAWQ_",Scen_RunDate,".csv"
        ),
        row.names = F
      )