This report is automatically generated with the R package knitr (version 1.39) .

---
title: '03_Wastewater Treatment Plants'
subtitle: 'MeHg Load Calculations'
author: "Mercury Program and Basin Planning Unit"
date: "3/1/2022"
output:
  html_document:
    code_folding: show
    toc: TRUE
    toc_float: TRUE
    toc_depth: 3
runtime: shiny
assets:
  css:
    - "http://fonts.googleapis.com/css?family=Raleway:300"
    - "http://fonts.googleapis.com/css?family=Oxygen"
---
---

<style>
body{
  font-family: 'Oxygen', sans-serif;
  font-size: 16px;
  line-height: 24px;
}

h1,h2,h3,h4 {
  font-family: 'Raleway', sans-serif;
}

.container { width: 1250px; }
h3 {
  background-color: #D4DAEC;
  text-indent: 50px; 
}
h4 {
  text-indent: 75px;
  margin-top: 35px;
  margin-bottom: 5px;
}
</style>

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE, fig.width=9.5)
```


```{r Libraries, echo=FALSE}
library(kableExtra) # better formatting of tables
library(shiny)

# Haf issue trying to set WD with Shiny in R project, reset working directory of rproj
wd <- rstudioapi::getActiveProject()
setwd(wd)

source("R Functions/functions_estimate NDDNQ values.R")
source("R Functions/functions_QA data.R")

# Read MeHg Data set
mehg <- readxl::read_xlsx("Reeval_Source_Analysis/Source Data/03a_Municipal WWTPs (NPDES)/Data/03.1_WWTP MeHg Data Prep_Clean_2023-02-13.xlsx", guess_max = 300000, sheet = 1) %>% 
  arrange(Discharger, Facility, SampleTypeCode, SampleDate) %>% 
  mutate(Year = year(SampleDate)) %>%
  filter(Discharger != 'Davis, City of' | SampleTypeCode != 'EFF-001') # Remove Discharges not in Delta

mgd <- readxl::read_xlsx("Reeval_Source_Analysis/Source Data/03a_Municipal WWTPs (NPDES)/Data/03.1_WWTP Flow Data Prep_Clean_2023-02-13.xlsx", guess_max = 300000, sheet = 1) %>% 
  arrange(Discharger, Facility, SampleTypeCode, SampleDate) %>% 
  mutate(Year = year(SampleDate))

StaffReport2010 <- readxl::read_xlsx("Reeval_Source_Analysis/Source Data/03a_Municipal WWTPs (NPDES)/Data/2010StaffReport-FacilityMeHgTotHgLoads.xlsx", guess_max = 300000, sheet = '2010 MeHg Loads ActivIn2021') %>% 
  select('Discharger','Facility','MeHg_Conc_2010','Flow_mgd_2010','MeHg_Load_2010')
```

## Check If WWTP MeHg & Flow Sampling Locations Match

### Years of MeHg Sampling
```{r}
mehg %>% 
  group_by(Discharger, Facility, SampleTypeCode) %>% 
  summarize(Years = paste(sort(unique(Year)), collapse=", "),
            .groups = 'drop') %>% 
  kbl(align='lcc', caption="MeHg All Years of Sampling") %>%
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)

# Get 5 most recent MeHg years to determine loading
mehg_5yrs <- NULL
Dischargers <- unique(mehg$GraphName)
for(discharger in Dischargers){
  temp1 <- mehg %>% 
    filter(GraphName == discharger) %>% 
    arrange(Year)
  years <- unique(temp1$Year)
    for(k in length(years):ifelse(length(years)>=5, length(years)-4, 1)){
      temp3 <- temp1 %>% 
        filter(Year == years[k])
      mehg_5yrs <- rbind(mehg_5yrs, temp3)
    }
}

mehg_yrsummary <- mehg_5yrs %>% 
  group_by(Discharger, Facility, SampleTypeCode, Subarea, Permit, Receiving_water) %>% 
  summarize(MeHg_Years = paste(sort(unique(Year)), collapse=", "),
            .groups = 'drop')
```


### Years of Flow Sampling Narrowed Down by MeHg Sampling Locations
```{r}
mgd_match <- mgd %>% 
  semi_join(., mehg_5yrs, by=c('Discharger', 'Facility', 'SampleTypeCode'))

mgd_match %>% 
  group_by(Discharger, Facility, SampleTypeCode) %>% 
  summarize(Years = paste(sort(unique(Year)), collapse=", "),
            .groups = 'drop') %>% 
  kbl(align='lcc', caption="Flow All Years of Sampling") %>%
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)

# Get 5 most recent Flow years to determine loading
mgd_5yrs <- NULL
StationNames <- unique(mgd_match$GraphName)
for(station in StationNames){
  temp1 <- mgd_match %>% 
    filter(GraphName == station) %>% 
    arrange(Year)
  years <- unique(temp1$Year)
    for(k in length(years):ifelse(length(years)>=5, length(years)-4, 1)){
      temp3 <- temp1 %>% 
        filter(Year == years[k])
      mgd_5yrs <- rbind(mgd_5yrs, temp3)
    }
}
```


### MeHg Sampling & Corresponding Flow Sampling
```{r}
mehg_freq <- mehg_5yrs %>% 
  group_by(Discharger, Facility, SampleTypeCode, Subarea, Permit, Receiving_water) %>% 
  summarize(Yrs_MeHg = paste(sort(unique(Year)), collapse=", "),
            n_MeHg = n(),
            .groups = 'drop')

mgd_freq <- mgd_5yrs %>% 
  group_by(Discharger, Facility, SampleTypeCode) %>% 
  summarize(Yrs_Flow = paste(sort(unique(Year)), collapse=", "),
            n_Flow = n(),
            .groups = 'drop')

yr_freq <- mehg_freq %>%  
  full_join(., mgd_freq, by=c('Discharger', 'Facility', 'SampleTypeCode'))

# Yearly Sampling Frequencies for All Sites
yr_freq %>% 
  arrange(Discharger, Facility, SampleTypeCode) %>% 
  kbl(align='lcc', caption="Yearly Sampling Frequency") %>%
  kable_paper('hover', full_width=T) %>%
  kable_styling(fixed_thead = T)
```



## Calculate Loads

### Calculate MeHg & Flow data Median & Mean by Pooling 5yr Data
```{r}
setwd(wd)

# Save Raw Data used to Calculated Loads
mun_ind_load_data <- mehg_5yrs %>%
  select(c("SourceID", "StationName", "SampleDate", "Result", "Unit", "MDL", "RL", "ResultQualCode", "Censored", "StationCode")) %>%
  bind_rows(., mgd_5yrs %>% select(c("SourceID", "StationName", "SampleDate", "Result", "Unit", "StationCode")) %>% mutate(MDL=NA_real_, RL=NA_real_, ResultQualCode=NA, Censored = NA)) %>%
  arrange(StationName, desc(Unit), SampleDate)

writexl::write_xlsx(list("NPDES Mun. & Ind." = mun_ind_load_data),
                    path='Reeval_Source_Analysis/Source Data/03a_Municipal WWTPs (NPDES)/Output/Appdx A - Data Compilation_Mun & Ind Sources.xlsx')



# Format Aq MeHg Data to calculate Loads
mehg_results <- mehg_5yrs %>% 
  group_by(Discharger, Facility, Subarea, Permit, SampleTypeCode, FacilityType, Treatment, Receiving_water) %>% 
  summarize(MeHgNDs = nrow(cur_data() %>% filter(ResultQualCode == 'ND')), # this operation needs to go 1st for some reason, maybe an issue with cur_data()??
            Locations = unique(SampleTypeCode),
            MeHg_n = n(),
            MeHg_StartDate = sprintf('%04d-%02d', year(min(SampleDate)), month(min(SampleDate))),
            MeHg_EndDate   = sprintf('%04d-%02d', year(max(SampleDate)), month(max(SampleDate))),
            MeHg_Years = paste(sort(unique(Year)), collapse=", "),
            MeHg_min = min(Result),
            MeHg_max = max(Result),
            MeHg_range = paste(ifelse(min(Result)<median(MDL, na.rm=T), paste0('<',median(MDL, na.rm=T)), min(Result)), max(Result), sep=' - '),
            ND_min = min(MDL, na.rm=T), 
            ND_max = max(MDL, na.rm=T), 
            
            ND_median = median(MDL, na.rm=T),
            MeHg_median = median(Result, na.rm=T),
            final_MeHg_median = NA, # setup column for later use
            uses_ND_median    = NA, # setup column for later use
            
            ND_mean   = mean(MDL, na.rm=T),
            MeHg_mean = mean(Result, na.rm=T),
            final_MeHg_mean   = NA, # setup column for later use
            uses_ND_mean      = NA, # setup column for later use
            .groups = 'drop') %>%
  arrange(Discharger, Facility)

mehg_results %>% 
  arrange(Discharger, Facility, SampleTypeCode) %>% 
  kbl(align='lcc', caption="MeHg Data used to Calculate Loads") %>%
  kable_paper('hover', full_width=T) %>%
  kable_styling(fixed_thead = T)


# Format MGD Data to calculate Loads
mgd_results1 <- mgd_5yrs %>% 
  group_by(Discharger, Facility, SampleTypeCode) %>% 
  summarize(Flow_n = n(),
            Flow_StartDate = sprintf('%04d-%02d', year(min(SampleDate)), month(min(SampleDate))),
            Flow_EndDate   = sprintf('%04d-%02d', year(max(SampleDate)), month(max(SampleDate))),
            Flow_Years = paste(sort(unique(Year)), collapse=", "),
            Flow_min = min(Result),
            Flow_max = max(Result),
            .groups = 'drop')

# Calculate Mean & Median MGD based on Annual Flow Volume
mgd_results3 <- mgd_5yrs %>% 
  group_by(Discharger, Facility, SampleTypeCode, Year) %>% 
  summarize(Flow_AnnualSum = sum(Result),
            .groups = 'drop') %>% 
  group_by(Discharger, Facility, SampleTypeCode) %>% 
  summarize(Flow_median = ifelse(median(Flow_AnnualSum) == 0, mean(Flow_AnnualSum)/365, median(Flow_AnnualSum)/365),
            Flow_mean = mean(Flow_AnnualSum)/365,
            .groups = 'drop')

mgd_results <- mgd_results1 %>% 
  left_join(mgd_results3, by=c('Discharger', 'Facility', 'SampleTypeCode'))

mgd_results %>% 
  arrange(Discharger, Facility, SampleTypeCode) %>% 
  kbl(align='lcc', caption="MGD Data used to Calculate Loads") %>%
  kable_paper('hover', full_width=T) %>%
  kable_styling(fixed_thead = T)

# Merge MeHg & Flow Data
combine_mehg_mgd <- mehg_results %>% 
  left_join(., mgd_results, by=c('Discharger', 'Facility', 'SampleTypeCode'))

combine_mehg_mgd %>% 
  arrange(Discharger, Facility, SampleTypeCode) %>% 
  kbl(align='lcc', caption="Combined MeHg & MGD Data used to Calculate Loads") %>%
  kable_paper('hover', full_width=T) %>%
  kable_styling(fixed_thead = T)
```

### Add Lathrop R5-2022-0004 Permitted Flow of 2.5mgd
```{r}
lathrop_permit_mgd <- tibble(Discharger  = "Lathrop, City of",
                             Flow_Years  = "2022 Permitted Flow",
                             Flow_median = 2.5, # Permitted mgd from R5-2022-0004
                             Flow_mean = 2.5)   # Permitted mgd from R5-2022-0004
                              

combine_mehg_mgd <- combine_mehg_mgd %>% 
  rows_update(lathrop_permit_mgd, by = "Discharger")
```

### See if MeHg and MGD data from Sac Combined Eff-002 & Eff-006 should be pooled or not
```{r}
# Sac combined is the only NPDES facility with multiple discharge locations

sac_combined_mehg1 <- mehg_5yrs %>% 
  filter(Discharger == "Sacramento, City of") %>% 
  select(Discharger, Facility, SampleTypeCode, Result, Year)

# Combine SampleTypeCodes
sac_combined_mehg <- sac_combined_mehg1 %>% 
  mutate(SampleTypeCode = "EFF-002 & 006 Pooled") %>% 
  bind_rows(sac_combined_mehg1)

sac_combined_mehg %>% 
  ggplot(aes(x=SampleTypeCode, y=Result, fill=SampleTypeCode)) +
  geom_boxplot() +
  scale_fill_viridis_d(option="E", alpha=0.6) +
  geom_jitter(aes(color=Year), size=1.5, alpha=0.9) +
  stat_summary(fun=mean, geom="point", shape=21, size=3, color="red", fill="white") +
  theme_light() +
  theme(
    legend.position="none",
    plot.title = element_text(size=11)
    ) +
  ggtitle(unique(sac_combined_mehg$Discharger)) +
  xlab("Methylmercury")

# Randomization Test of MeHg Concentrations
mehg_ptest <- randomizationTest(sac_combined_mehg1, Result, SampleTypeCode, N=5000)
mehg_ptest$p_value_diff_medians
mehg_ptest$median_diff_dist_plot
cat("MeHg concentrations are not significantly different bewtween Eff-002 & Eff-006")


# MGD Flow
sac_combined_mgd1 <- mgd_5yrs %>% 
  filter(Discharger == "Sacramento, City of") %>% 
  select(Discharger, Facility, SampleTypeCode, Result, Year) %>% 
  mutate(Type = "MGD")

# Combine SampleTypeCodes
sac_combined_mgd <- sac_combined_mgd1 %>% 
  mutate(SampleTypeCode = "EFF-002 & 006 Pooled") %>% 
  bind_rows(sac_combined_mgd1)

sac_combined_mgd %>% 
  ggplot(aes(x=SampleTypeCode, y=Result, fill=SampleTypeCode)) +
  geom_boxplot() +
  scale_fill_viridis_d(option="E", alpha=0.6) +
  geom_jitter(aes(color=Year), size=1.5, alpha=0.9) +
  stat_summary(fun=mean, geom="point", shape=21, size=3, color="red", fill="white") +
  theme_light() +
  theme(
    legend.position="none",
    plot.title = element_text(size=11)
    ) +
  ggtitle(unique(sac_combined_mgd$Discharger)) +
  xlab("Million Gallons per Day")

# Randomization Test of Daily Flow
mgd_ptest <- randomizationTest(sac_combined_mgd1, Result, SampleTypeCode, N=5000)
mgd_ptest$p_value_diff_medians
mgd_ptest$median_diff_dist_plot
cat("Daily flow volumes are not significantly different between Eff-002 & Eff-006")


# Box Plot of Annual Flow Volumes by Location Compared to Median of EFF-002 & Eff-006 Median Annual Flow Volume
sac_combined_annualFlow <- mgd_5yrs %>% 
  filter(Discharger == "Sacramento, City of") %>%
  group_by(Discharger, Facility, SampleTypeCode, Year) %>% 
  summarize(Flow_AnnualSum = sum(Result),
            .groups = 'drop')

sac_combined_annualFlow_all <- mgd_5yrs %>% 
  filter(Discharger == "Sacramento, City of") %>%
  group_by(Discharger, Facility, SampleTypeCode, Year) %>% 
  summarize(Flow_AnnualSum = sum(Result),
            .groups = 'drop') %>% 
  group_by(Discharger, Facility, SampleTypeCode) %>% 
  summarize(Flow_AnnualSum = ifelse(median(Flow_AnnualSum) == 0, mean(Flow_AnnualSum), median(Flow_AnnualSum)),
            .groups = 'drop') %>% 
  mutate(SampleTypeCode = "EFF-002 & 006 Median",
         Year = "Median")

sac_combined_annualFlow <- sac_combined_annualFlow %>% 
  mutate(Year = as.character(Year)) %>% 
  bind_rows(sac_combined_annualFlow_all)

sac_combined_annualFlow %>% 
  ggplot(aes(x=SampleTypeCode, y=Flow_AnnualSum, fill=SampleTypeCode)) +
  geom_boxplot() +
  scale_fill_viridis_d(option="E", alpha=0.6) +
  geom_jitter(aes(color=Year, shape=Year), size=1.5, alpha=0.9) +
  stat_summary(fun=mean, geom="point", shape=21, size=3, color="red", fill="white") +
  theme_light() +
  theme(
    legend.position="none",
    plot.title = element_text(size=11)
    ) +
  ggtitle(paste(unique(sac_combined_annualFlow$Discharger), "- Median of EFF-002 & Eff-006 Flow Medians")) +
  xlab("Million Gallons per Year")

# Box Plot of Annual Flow Volumes by Location Compared to Pooling of EFF-002 & Eff-006 Flow Volumes
sac_combined_pooledannualFlow <- mgd_5yrs %>% 
  filter(Discharger == "Sacramento, City of") %>%
  group_by(Discharger, Facility, Year) %>% 
  summarize(Flow_AnnualSum = sum(Result),
            .groups = 'drop') %>% 
  mutate(SampleTypeCode = "EFF-002 & 006 Pooled",
         Year = as.character(Year))

sac_combined_annualFlow_pooled <- sac_combined_annualFlow %>% 
  filter(!grepl("&", SampleTypeCode)) %>% 
  bind_rows(sac_combined_pooledannualFlow)

sac_combined_annualFlow_pooled %>% 
  ggplot(aes(x=SampleTypeCode, y=Flow_AnnualSum, fill=SampleTypeCode)) +
  geom_boxplot() +
  scale_fill_viridis_d(option="E", alpha=0.6) +
  geom_jitter(aes(color=Year, shape=Year), size=1.5, alpha=0.9) +
  stat_summary(fun=mean, geom="point", shape=21, size=3, color="red", fill="white") +
  theme_light() +
  theme(
    legend.position="none",
    plot.title = element_text(size=11)
    ) +
  ggtitle(paste(unique(sac_combined_annualFlow$Discharger), "- Pooling of EFF-002 & Eff-006 Flow Data")) +
  xlab("Million Gallons per Year")

# Randomization Test of Annual Flow
annualflow_ptest <- randomizationTest(sac_combined_annualFlow %>% filter(!grepl("&", SampleTypeCode)), Flow_AnnualSum, SampleTypeCode, N=5000)
annualflow_ptest$p_value_diff_medians
annualflow_ptest$median_diff_dist_plot
cat("Annual flow volumes are significantly different between Eff-002 & Eff-006")
```
####Analysis Conclusion
Even though the median MeHg concentrations and MGD flow rates were not significantly different between Eff-002 & Eff-006, the annual flow volumes (which are used to calculate MeHg loads) are different enough (p < 0.05) that it would not be appropriate to pool the data of the two discharge locations. Therefore, the MeHg load will be calculated for each location and the two loads will be added together to give Sac City Combined one allocation rather than an allocation for each discharge.   
```

### Calculate MeHg Annual Median & Mean Load
```{r} 
annual_Loads_loc <- combine_mehg_mgd %>% 
  mutate(final_MeHg_median = if_else(round(MeHg_median, 3) <= ND_median, ND_median, round(MeHg_median, 3)),
         uses_ND_median    = if_else(round(MeHg_median, 3) <= ND_median,     T,        F),
         
         MedianLoad_grams = (final_MeHg_median * 1/1000000000) * (Flow_median * 1/0.2642 * 1000000 * 365),  # ng/L * L/d * g/1,000,000,000 ng * MGD * 1 L/0.2642 gal * 1,000,000 gal/1 Mgal * 365d/yr
         
         final_MeHg_mean   = if_else(round(MeHg_mean, 3) <= ND_mean,   ND_mean,   round(MeHg_mean, 3)),
         uses_ND_mean      = if_else(round(MeHg_mean, 3) <= ND_mean,     T,        F),
         
         MeanLoad_grams   = (final_MeHg_mean * 1/1000000000) * (Flow_mean * 1/0.2642 * 1000000 * 365)) %>%    # (ng/L * g/1,000,000,000 ng) * (MGD * 1 L/0.2642 gal * 1,000,000 gal/1 Mgal * 365d/yr)
  left_join(StaffReport2010, by=c('Discharger','Facility'))

annual_Loads_loc %>% 
  arrange(Discharger, Facility, SampleTypeCode) %>% 
  kbl(align='lcc', caption="Combined MeHg & MGD Data used to Calculate Loads") %>%
  kable_paper('hover', full_width=T) %>%
  kable_styling(fixed_thead = T)
  


annual_Loads <- annual_Loads_loc %>% 
  group_by(Discharger, Facility, Subarea, Permit, FacilityType, Treatment, Receiving_water) %>%
  summarize(Locations = paste(unique(Locations), collapse="; "),
            nLocations = n(),
            MeHg_n = sum(MeHg_n),
            MeHg_StartDate = min(MeHg_StartDate),
            MeHg_EndDate   = max(MeHg_EndDate),
            MeHg_Years = paste(sort(unique(MeHg_Years)), collapse="; "),
            MeHg_min = min(MeHg_min),
            MeHg_max = max(MeHg_max),
            MeHg_range = paste(ifelse(min(MeHg_min)<median(ND_median), paste0('<',median(ND_median)), min(MeHg_min)), max(MeHg_max), sep=' - '),
            MeHgNDs = sum(MeHgNDs), 
            ND_range = paste(min(ND_min), max(ND_max), sep=' - '),
            
            ND_median = median(ND_median),
            MeHg_median = median(MeHg_median),
            final_MeHg_median = if_else(nLocations > 1, NA_real_,
                                        if_else(round(median(final_MeHg_median), 3) <= ND_median, ND_median, round(median(final_MeHg_median), 3))),
            uses_ND_median    = if_else(round(median(final_MeHg_median), 3) <= ND_median,             T,        F),
            
            ND_mean   = mean(ND_mean),
            MeHg_mean = mean(MeHg_mean),
            final_MeHg_mean = if_else(nLocations > 1, NA_real_,
                                      if_else(round(mean(final_MeHg_mean), 3) <= ND_mean, ND_mean, round(mean(final_MeHg_mean), 3))),
            uses_ND_mean    = if_else(round(median(final_MeHg_median), 3) <= ND_median,     T,        F),
            
            
            Flow_n = sum(Flow_n),
            Flow_StartDate = min(Flow_StartDate),
            Flow_EndDate   = max(Flow_EndDate),
            Flow_Years = paste(sort(unique(Flow_Years)), collapse="; "),
            Flow_min = min(Flow_min),
            Flow_max = max(Flow_max),
            Flow_median = median(Flow_median),
            final_Flow_median = if_else(nLocations > 1, NA_real_, Flow_median),
            MedianLoad_grams = sum(MedianLoad_grams),
            
            Flow_mean = mean(Flow_mean),
            final_Flow_mean = if_else(nLocations > 1, NA_real_, Flow_mean),
            MeanLoad_grams  = sum(MeanLoad_grams),
            
            .groups = 'drop') %>%
  left_join(StaffReport2010, by=c('Discharger','Facility'))

annual_Loads %>% 
  arrange(Discharger, Facility) %>% 
  kbl(align='lcc', caption="Combined MeHg & MGD Data used to Calculate Loads") %>%
  kable_paper('hover', full_width=T) %>%
  kable_styling(fixed_thead = T)
```


### Calculate Summary Statistics & Add Rows to Bottom of Table
```{r}
# Calculate Summary Statistics for each FacilityType (Municipal, Groundwater, and All)
FacType <- c('All', 'Municipal', 'Groundwater')
# loads <- NULL
for(facility in FacType){
  if(facility == 'All'){
    #Calculate Summary Statistics based on Discharge location
    loads <- annual_Loads_loc

  }else{
    #Calculate Summary Statistics based on Facility
    loads <- annual_Loads %>%
      filter(grepl(facility, FacilityType))
  }
    
    colMin <- loads %>%
  add_row(Facility = "Minimum", summarize(., across(contains(c('_n','Date', 'min', 'max', 'median', 'mean', '2010')), ~min(.x, na.rm=T)))) %>%
  slice_tail(n=1)
    
    colMax <- loads %>%
  add_row(Facility = "Maximum", summarize(., across(contains(c('_n','Date', 'min', 'max', 'median', 'mean', '2010')), ~max(.x, na.rm=T)))) %>%
  slice_tail(n=1)
    
    colMedians <- loads %>%
  add_row(Facility = "Median of Medians", summarize(., across(contains(c('median', '2010')), ~median(.x, na.rm=T)))) %>%
  slice_tail(n=1)
    
    colMeans <- loads %>%
  add_row(Facility = "Mean of Means", summarize(., across(contains(c('mean', '2010')), ~mean(.x, na.rm=T)))) %>%
  slice_tail(n=1)
    
    colTotal <- loads %>%
  add_row(Facility = "Summation", summarize(., across(contains(c('median', 'mean', '2010')), ~sum(.x, na.rm=T)))) %>%
  slice_tail(n=1)
    
    assign(paste0(facility,"_result_summary"), 
       loads %>%
         bind_rows(., colMin, colMax, colMedians, colMeans, colTotal)  %>% 
         mutate(FlowMedian_to2010_PercDiff = (Flow_median - Flow_mgd_2010)/ Flow_mgd_2010 * 100,
                LoadMedian_to2010_PercDiff = (MedianLoad_grams - MeHg_Load_2010)/ MeHg_Load_2010 * 100,
                LoadMean_to2010_PercDiff   = (MeanLoad_grams - MeHg_Load_2010)/ MeHg_Load_2010 * 100)
    )
} 
```



## Export to excel
```{r}
writexl::write_xlsx(list(All=All_result_summary,
                         WWTP=Municipal_result_summary,
                         Ind=Groundwater_result_summary), paste0(wd, "/Reeval_Source_Analysis/Source Data/03a_Municipal WWTPs (NPDES)/Output/03a_WWTP MeHg Annual Loads_", today(), ".xlsx"))
```
## Error: <text>:20:1: unexpected '<'
## 19: 
## 20: <
##     ^

The R session information (including the OS info, R version and all packages used):

    sessionInfo()
## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] lubridate_1.8.0    plotly_4.10.0      readxl_1.4.0       actuar_3.2-2      
##  [5] NADA_1.6-1.1       forcats_0.5.1      stringr_1.4.0      dplyr_1.0.9       
##  [9] purrr_0.3.4        readr_2.1.2        tidyr_1.2.0        tibble_3.1.7      
## [13] ggplot2_3.3.6      tidyverse_1.3.1    fitdistrplus_1.1-8 survival_3.4-0    
## [17] MASS_7.3-58.1     
## 
## loaded via a namespace (and not attached):
##  [1] lattice_0.20-45    assertthat_0.2.1   digest_0.6.29      utf8_1.2.2        
##  [5] R6_2.5.1           cellranger_1.1.0   backports_1.4.1    reprex_2.0.1      
##  [9] evaluate_0.15      httr_1.4.3         highr_0.9          pillar_1.7.0      
## [13] rlang_1.0.2        lazyeval_0.2.2     rstudioapi_0.13    data.table_1.14.2 
## [17] Matrix_1.5-1       rmarkdown_2.14     labeling_0.4.2     splines_4.2.2     
## [21] htmlwidgets_1.5.4  munsell_0.5.0      broom_0.8.0        compiler_4.2.2    
## [25] modelr_0.1.8       xfun_0.31          pkgconfig_2.0.3    htmltools_0.5.2   
## [29] tidyselect_1.1.2   viridisLite_0.4.0  fansi_1.0.3        crayon_1.5.1      
## [33] tzdb_0.3.0         dbplyr_2.2.0       withr_2.5.0        grid_4.2.2        
## [37] jsonlite_1.8.0     gtable_0.3.0       lifecycle_1.0.1    DBI_1.1.2         
## [41] magrittr_2.0.3     scales_1.2.0       writexl_1.4.0      cli_3.3.0         
## [45] stringi_1.7.6      farver_2.1.0       fs_1.5.2           xml2_1.3.3        
## [49] ellipsis_0.3.2     generics_0.1.2     vctrs_0.4.1        expint_0.1-7      
## [53] RColorBrewer_1.1-3 tools_4.2.2        glue_1.6.2         hms_1.1.1         
## [57] yaml_2.3.5         fastmap_1.1.0      colorspace_2.0-3   rvest_1.0.2       
## [61] knitr_1.39         haven_2.5.0
    Sys.time()
## [1] "2023-12-27 10:07:06 PST"