##Depth to Flowrate (cfs) conversion
depth<-outfall3_clean_20212.20222$Depth
aco2dep<-acos((30-depth)/30)*2
sinangle<-sin(aco2dep)

# Depth is Lower than R
A<-30*30*0.5*(aco2dep-sinangle)
P<-30*aco2dep
R<-A/P
# Depth is higher than R
A2<-(3.14*30*30)-A
P2<-(2*3.14*30)-P
R2<-A2/P2

# Manning Equation
# Q = 1.49/n * A *R^(2/3) *S^0.5 #n=roughness coef; S=Slope 
# n= 0.015

## excluede rows with NA depth value of Outfall 5
depthdt5<-outfall5_clean_20212.20222[!is.na(outfall5_clean_20212.20222$Depth),]
depthlow5<-depthdt5[(depthdt5$Depth)<30,]
depthhigh5<-depthdt5[(depthdt5$Depth)>=30,]
depthl5<-depthlow5$Depth
depthh5<-depthhigh5$Depth

## when lower than R 
aco2dep5<-acos((30-depthl5)/30)*2
sinangle5<-sin(aco2dep5)
A5<-30*30*0.5*(aco2dep5-sinangle5)
P5<-30*aco2dep5
R5<-A5/P5
Ql5<- 106.4286*A5*(R5^(2/3))*(0.027^0.5)
Qlcft5<-Ql5/1728/0.14
## the unit would be cubic feet/sec

## when higher than R
acohdep5<-acos((30-depthh5)/30)*2
sinangleh5<-sin(acohdep5)
A5<-30*30*0.5*(acohdep5-sinangleh5)
P5<-30*acohdep5
A25<-(3.14*30*30)-A5
P25<-(2*3.14*30)-P5
R25<-A25/P25
Qh5<- 106.4286*A25*(R25^(2/3))*(0.027^0.5) ## unit cubic inch
Qhcfs5<-Qh5/1728/0.14 ## unit cubic feet/sec


## combine QL and Qh together
flowrate5<-append(Qlcft5,Qhcfs5)
## create dataframe for correlation test- figure 6
## combine corresponding qualitor 
turl<-depthlow5$Turbidity
turh<-depthhigh5$Turbidity
turb5<-append(turl,turh)

frtu5<-data.frame(flowrate5,turb5)
ggscatter(frtu5,x="flowrate5",y="turb5",
              add="reg.line",conf.int=TRUE,
              cor.coef=TRUE, cor.method="spearman",
              pch = 1, col =  "#00AFBB"          )

condl<-depthlow5$Conductivity
condh<-depthhigh5$Conductivity
condb5<-append(condl,condh)

frcon5<-data.frame(flowrate5,condb5)
ggscatter(frcon5,x="flowrate5",y="condb5",
          add="reg.line",conf.int=TRUE,
          cor.coef=TRUE, cor.method="spearman",
          pch = 1, col =  "#00AFBB"          )

DOl<-depthlow5$RDO
DOh<-depthhigh5$RDO
DO5<-append(DOl,DOh)

frdo5<-data.frame(flowrate5,DO5)
ggscatter(frdo5,x="flowrate5",y="DO5",
          add="reg.line",conf.int=TRUE,
          cor.coef=TRUE, cor.method="spearman",
          pch = 1, col =  "#00AFBB"          )

PHl<-depthlow5$pH
PHh<-depthhigh5$pH
PH5<-append(PHl,PHh)
frph5<-data.frame(flowrate5,PH5)
ggscatter(frph5,x="flowrate5",y="PH5",
          add="reg.line",conf.int=TRUE,
          cor.coef=TRUE, cor.method="spearman",
          pch = 1, col =  "#00AFBB"          )

### Outfall 3
## excluede rows with NA depth value of Outfall 3
depthdt3<-outfall3_clean_20212.20222[!is.na(outfall3_clean_20212.20222$Depth),]
depthlow3<-depthdt3[(depthdt3$Depth)<30,]
depthhigh3<-depthdt3[(depthdt3$Depth)>=30,]
depthl3<-depthlow3$Depth
depthh3<-depthhigh3$Depth

## when lower than R 
aco2dep3<-acos((30-depthl3)/30)*2
sinangle3<-sin(aco2dep3)
A3<-30*30*0.5*(aco2dep3-sinangle3)
P3<-30*aco2dep3
R3<-A3/P3
Ql3<- 106.4286*A3*(R3^(2/3))*(0.017^0.5)
Qlcft3<-Ql3/1728/0.1
## the unit would be cubic feet/sec

## when higher than R
acohdep3<-acos((30-depthh3)/30)*2
sinangleh3<-sin(acohdep3)
A3<-30*30*0.5*(acohdep3-sinangleh3)
P3<-30*acohdep3
A23<-(3.14*30*30)-A3
P23<-(2*3.14*30)-P3
R23<-A23/P23
Qh3<- 106.4286*A23*(R23^(2/3))*(0.017^0.5) ## unit cubic inch
Qhcfs3<-Qh3/1728/0.1 ## unit cubic feet/sec


## combine QL and Qh together
flowrate3<-append(Qlcft3,Qhcfs3)
## create dataframe for correlation test 
## combine corresponding qualitor 
turl3<-depthlow3$Turbidity
turh3<-depthhigh3$Turbidity
turb3<-append(turl3,turh3)

frtu3<-data.frame(flowrate3,turb3)
ggscatter(frtu3,x="flowrate3",y="turb3",
          add="reg.line",conf.int=TRUE,
          cor.coef=TRUE, cor.method="spearman",
          pch = 1, col =  "#00AFBB"          )

condl3<-depthlow3$Conductivity
condh3<-depthhigh3$Conductivity
condb3<-append(condl3,condh3)

frcon3<-data.frame(flowrate3,condb3)
ggscatter(frcon3,x="flowrate3",y="condb3",
          add="reg.line",conf.int=TRUE,
          xlab = "Flowrate (cfs)",
          ylab = "Conductivity (uS/cm)",
          cor.coef=TRUE, cor.method="spearman",
          pch = 1, col =  "#00AFBB"          )

DOl3<-depthlow3$RDO
DOh3<-depthhigh3$RDO
DO3<-append(DOl3,DOh3)

frdo3<-data.frame(flowrate3,DO3)
ggscatter(frdo3,x="flowrate3",y="DO3",
          add="reg.line",conf.int=TRUE,
          cor.coef=TRUE, cor.method="spearman",
          pch = 1, col =  "#00AFBB"          )

PHl3<-depthlow3$pH
PHh3<-depthhigh3$pH
PH3<-append(PHl3,PHh3)
frph3<-data.frame(flowrate3,PH3)
ggscatter(frph3,x="flowrate3",y="PH3",
          add="reg.line",conf.int=TRUE,
          cor.coef=TRUE, cor.method="spearman",
          pch = 1, col =  "#00AFBB"          )

### Outfall 19
## excluede rows with NA depth value of Outfall 19
depthdt19<-outfall19_clean_20212.20222[!is.na(outfall19_clean_20212.20222$Depth),]
depthlow19<-depthdt19[(depthdt19$Depth)<20,]
depthhigh19<-depthdt19[(depthdt19$Depth)>=20,]
depthl19<-depthlow19$Depth
depthh19<-depthhigh19$Depth

## when lower than R 
aco2dep19<-acos((20-depthl19)/20)*2
sinangle19<-sin(aco2dep19)
A19<-20*20*0.5*(aco2dep19-sinangle19)
P19<-20*aco2dep19
R19<-A19/P19
Ql19<- 106.4286*A19*(R19^(2/3))*(0.031^0.5)
Qlcft19<-Ql19/1728/0.021
## the unit would be cubic feet/sec

## when higher than R
acohdep19<-acos((20-depthh19)/20)*2
sinangleh19<-sin(acohdep19)
A19<-20*20*0.5*(acohdep19-sinangleh19)
P19<-20*acohdep19
A219<-(3.14*20*20)-A19
P219<-(2*3.14*20)-P19
R219<-A219/P219
Qh19<- 106.4286*A219*(R219^(2/3))*(0.031^0.5) ## unit cubic inch
Qhcfs19<-Qh19/1728/0.021 ## unit cubic feet and normalized by drainage area


## combine QL and Qh together
flowrate19<-append(Qlcft19,Qhcfs19)
## create dataframe for correlation test 
## combine corresponding qualitor 
turl19<-depthlow19$Turbidity
turh19<-depthhigh19$Turbidity
turb19<-append(turl19,turh19)

frtu19<-data.frame(flowrate19,turb19)
ggscatter(frtu19,x="flowrate19",y="turb19",
          add="reg.line",conf.int=TRUE,
          cor.coef=TRUE, cor.method="spearman",
          pch = 1, col =  "#00AFBB"          )

condl19<-depthlow19$Conductivity
condh19<-depthhigh19$Conductivity
condb19<-append(condl19,condh19)

frcon19<-data.frame(flowrate19,condb19)
ggscatter(frcon19,x="flowrate19",y="condb19",
          add="reg.line",conf.int=TRUE,
          cor.coef=TRUE, cor.method="spearman",
          pch = 1, col =  "#00AFBB"          )

DOl19<-depthlow19$RDO
DOh19<-depthhigh19$RDO
DO19<-append(DOl19,DOh19)

frdo19<-data.frame(flowrate19,DO19)
ggscatter(frdo19,x="flowrate19",y="DO19",
          add="reg.line",conf.int=TRUE,
          cor.coef=TRUE, cor.method="spearman",
          pch = 1, col =  "#00AFBB"          )

PHl19<-depthlow19$pH
PHh19<-depthhigh19$pH
PH19<-append(PHl19,PHh19)
frph19<-data.frame(flowrate19,PH19)
ggscatter(frph19,x="flowrate19",y="PH19",
          add="reg.line",conf.int=TRUE,
          cor.coef=TRUE, cor.method="spearman",
          pch = 1, col =  "#00AFBB"          )


try1<-2*acos((30-4.22)/30)
try2=sin(try1)
A381<-30*30*0.5*(try1-try2)
P381<-30*try1
R381<-A381/P381
Ql381<- 106.4286*A381*(R381^(2/3))*(0.017^0.5)
Qlcft381<-Ql381/1728
Qlcft381

try1<-2*acos((30-55.72)/30)
try2<-sin(try1)
A19<-30*30*0.5*(try1-try2)
P19<-30*try1
A219<-(3.14*30*30)-A19
P219<-(2*3.14*30)-P19
R219<-A219/P219
Qh19<- 106.4286*A219*(R219^(2/3))*(0.027^0.5) ## unit cubic inch
Qhcfs19<-Qh19/1728 ## unit cubic feet
Qhcfs19

try1<-acos((20-4.93)/20)*2
try2<-sin(try1)
A19<-20*20*0.5*(try1-try2)
P19<-20*try1
R19<-A19/P19
Ql19<- 106.4286*A19*(R19^(2/3))*(0.031^0.5)
Qlcft19<-Ql19/1728
Qlcft19


### 19 outfall Area 573992 sqft = 0.021 mi2
### 3 outfall area 2907053 sqft = 0.1 mi2
### 5 outfall area 3773190 sqft = 0.14 mi2


## Figure 5

## Figure 5, graph b
par(mfrow=c(1,3))
boxplot(peak.flow.rate, ylab="Peak Flow Rate (unit:cfs/mi2)", pch=23)
boxplot(discharge.rate, ylab="Rise Rate (unit:cfs/mi2/min)", pch=23)
boxplot(recession.rate, ylab="Recession Rate (unit:cfs/mi2min)", pch=23)

##linear mixed model exploration
library(dplyr)
library(ggplot2)
##correct the representation of categorical variables
behav<-mutate(norm.flow.rate,Outfall=as.factor(Outfall),
              BMP.Coverage=as.factor(BMP.Coverage),
              IMP=as.factor(IMP),
              Slope=as.factor(Slope))
head(behav)

##only event-level predictor(NAP)
mod1<-lm(formula = Norm.peak.flow.rate.cfs.mi2.~Precipitation_NAP,data=behav)
summary(mod1)
behav$fit1<-fitted(mod1)
gg<-ggplot(behav,aes(x=Precipitation_NAP, y=Norm.peak.flow.rate.cfs.mi2.))+
  geom_point()+
  geom_line(aes(y=fit1)
  )
print(gg+labs(y="peak flow rate(cfs.mi2)", x="Precipitation Size (inches in 24 hrs)"))
## event+BMP - peak flow rate
mod2<-lm(Norm.peak.flow.rate.cfs.mi2.~Precipitation_NAP+BMP.Coverage,data=behav)
summary(mod2)
behav$fit2<-fitted(mod2)
gg<-ggplot(behav,aes(x=Precipitation_NAP,y=Norm.peak.flow.rate.cfs.mi2.,color=BMP.Coverage))+
  geom_point()+
  geom_line(aes(y=fit2))
gg
gg+scale_x_continuous(name="Precipitation Size (unit: inches in 24 hrs)") +
  scale_y_continuous(name="Peak Flow Rate (unit:cfs/mi2)")

## event+imp - peak flow rate
mod3<-lm(Norm.peak.flow.rate.cfs.mi2.~Precipitation_NAP+IMP,data=behav)
summary(mod3)
behav$fit3<-fitted(mod3)
gg<-ggplot(behav,aes(x=Precipitation_NAP,y=Norm.peak.flow.rate.cfs.mi2.,color=IMP))+
  geom_point()+
  geom_line(aes(y=fit3))
gg
gg+scale_x_continuous(name="Precipitation Size (unit: inches in 24 hrs)") +
  scale_y_continuous(name="Peak Flow Rate (unit:cfs/mi2)")

## event+imp - Peak flow rise rate
mod4<-lm(Norm.Discharge.Rate.cfs.min.mi2.~Precipitation_NAP+IMP,data=behav)
summary(mod4)
behav$fit4<-fitted(mod4)
gg<-ggplot(behav,aes(x=Precipitation_NAP,y=Norm.Discharge.Rate.cfs.min.mi2.,color=IMP))+
  geom_point()+
  geom_line(aes(y=fit4))
gg
gg+scale_x_continuous(name="Precipitation Size (unit: inches in 24 hrs)") +
  scale_y_continuous(name="Peak Flow Rise Rate (unit:cfs/mi2/min)")

## event+bmp - peak flow rise rate
mod5<-lm(Norm.Discharge.Rate.cfs.min.mi2.~Precipitation_NAP+BMP.Coverage,data=behav)
summary(mod5)
behav$fit5<-fitted(mod5)
gg<-ggplot(behav,aes(x=Precipitation_NAP,y=Norm.Discharge.Rate.cfs.min.mi2.,color=BMP.Coverage))+
  geom_point()+
  geom_line(aes(y=fit5))
gg
gg+scale_x_continuous(name="Precipitation Size (unit: inches in 24 hrs)") +
  scale_y_continuous(name="Peak Flow Rise Rate (unit:cfs/mi2/min)")

## event+bmp - recession rate
mod6<-lm(Norm.Recession.Rate.cfs.min.mi2.~Precipitation_NAP+BMP.Coverage,data=behav)
summary(mod6)
behav$fit6<-fitted(mod6)
gg<-ggplot(behav,aes(x=Precipitation_NAP,y=Norm.Recession.Rate.cfs.min.mi2.,color=BMP.Coverage))+
  geom_point()+
  geom_line(aes(y=fit6))
gg
gg+scale_x_continuous(name="Precipitation Size (unit: inches in 24 hrs)") +
  scale_y_continuous(name="Peak Flow Recession Rate (unit:cfs/mi2/min)")
## event+imp - recession rate
mod7<-lm(Norm.Recession.Rate.cfs.min.mi2.~Precipitation_NAP+IMP,data=behav)
summary(mod7)
behav$fit7<-fitted(mod7)
gg<-ggplot(behav,aes(x=Precipitation_NAP,y=Norm.Recession.Rate.cfs.min.mi2.,color=IMP))+
  geom_point()+
  geom_line(aes(y=fit7))
gg
gg+scale_x_continuous(name="Precipitation Size (unit: inches in 24 hrs)") +
  scale_y_continuous(name="Peak Flow Recession Rate (unit:cfs/mi2/min)")
## event+slope - peak flow rate
mod8<-lm(Norm.peak.flow.rate.cfs.mi2.~Precipitation_NAP+Slope,data=behav)
summary(mod8)
behav$fit8<-fitted(mod8)
gg<-ggplot(behav,aes(x=Precipitation_NAP,y=Norm.peak.flow.rate.cfs.mi2.,color=Slope))+
  geom_point()+
  geom_line(aes(y=fit8))
gg
gg+scale_x_continuous(name="Precipitation Size (unit: inches in 24 hrs)") +
  scale_y_continuous(name="Peak Flow Rate (unit:cfs/mi2)")

## Figure 5, graph c d 
mod0<-lm(Norm.Recession.Rate.cfs.min.mi2.~Precipitation_NAP++Slope+IMP,data=norm.flow.rate)
summary(mod0)
behav$fit0<-fitted(mod0)
gg<-ggplot(behav,aes(x=Precipitation_NAP,y=Norm.Discharge.Rate.cfs.min.mi2.,color=Slope))+
  geom_point()+
  geom_line(aes(y=fit0))
gg
gg2<-ggPredict(mod0,se=TRUE,interactive=FALSE)
gg2+scale_x_continuous(name="Precipitation Size (unit: inches in 24 hrs)") +
  scale_y_continuous(name="Flow Rise Rate (unit:cfs/mi2/min)")

mod0<-lm(Norm.peak.flow.rate.cfs.mi2.~Precipitation_NAP++Slope+IMP,data=norm.flow.rate)
summary(mod0)
library(ggplot2)
library(ggeffects)
library(car)
# Create a partial regression plot
ggplot(norm.flow.rate, aes(x = Precipitation_NAP, y = Norm.peak.flow.rate.cfs.mi2.)) +
  ggpredict(mod0, terms = "Precipitation_NAP") +
  geom_line(color = "blue") +
  geom_point() +
  labs(x = "Precipitation_NAP", y = "Norm.peak.flow.rate.cfs.mi2.", title = "Partial Regression Plot") +
  theme_bw()
library(GGally)

# Create a scatterplot matrix
ggpairs(norm.flow.rate, columns = c("Norm.peak.flow.rate.cfs.mi2.", "Precipitation_NAP", "Slope", "IMP", "BMP.Coverage")) +
  theme_bw()

## EXTRA - tranformation of reponse variable (outside the paper scope)
##behavio<-mutate(behav,Discharge=log10(Discharge.Rate.cfs.min.))
##hist(behavio$Discharge)
##only event-level predictor(NAP)
##Dmod1<-lm(formula = Discharge~Precipitation_NAP,data=behavio)
##summary(Dmod1)
##behav$fit1<-fitted(mod1)
##ggplot(behav,aes(x=Precipitation_NAP, y=Discharge.Rate.cfs.min))+
  ##geom_point()+
  ##geom_line(aes(y=fit1))
## event+BMP
##Dmod2<-lm(Discharge~Precipitation_NAP*BMP.Coverage,data=behavio)
##summary(Dmod2)
##behavio$fit2<-fitted(Dmod2)
##gg2<-ggplot(behavio,aes(x=Precipitation_NAP,y=Discharge,color=BMP.Coverage))+
  ##geom_point()+
  ##geom_line(aes(y=fit2))
##gg2+scale_x_continuous(name="Precipitation Size (unit: inches in 24 hrs)") +
  ##scale_y_continuous(name="Log-Transformed Runoff Flow Rise Rate (unit:cfs/min)")

## EXTRA -tranformation of reponse variable
##behavior<-mutate(behavi,Recession=log10(Recession.Rate.in.min.))
##hist(behavior$Recession)
##only event-level predictor(NAP)
##Rmod1<-lm(formula = Recession~Precipitation_NAP,data=behavior)
##summary(Rmod1)
##behavi$fit1<-fitted(mod1)
##ggplot(behavi,aes(x=Precipitation_NAP, y=PDepth))+
  ##geom_point()+
  ##geom_line(aes(y=fit1))
## event+BMP
##Rmod2<-lm(Recession~Precipitation_NAP*BMP,data=behavior)
##summary(Rmod2)
##behavior$fit2<-fitted(Rmod2)
##ggplot(behavior,aes(x=Precipitation_NAP,y=Recession,color=BMP))+
  ##geom_point()+
  ##geom_line(aes(y=fit2))